Permutation Invariant Gaussian Matrix Models

Permutation invariant Gaussian matrix models were recently developed for applications in computational linguistics. A 5-parameter family of models was solved. In this paper, we use a representation theoretic approach to solve the general 13-parameter Gaussian model, which can be viewed as a zero-dimensional quantum field theory. We express the two linear and eleven quadratic terms in the action in terms of representation theoretic parameters. These parameters are coefficients of simple quadratic expressions in terms of appropriate linear combinations of the matrix variables transforming in specific irreducible representations of the symmetric group $S_D$ where $D$ is the size of the matrices. They allow the identification of constraints which ensure a convergent Gaussian measure and well-defined expectation values for polynomial functions of the random matrix at all orders. A graph-theoretic interpretation is known to allow the enumeration of permutation invariants of matrices at linear, quadratic and higher orders. We express the expectation values of all the quadratic graph-basis invariants and a selection of cubic and quartic invariants in terms of the representation theoretic parameters of the model.

In the context of distributional semantics [1,2], the meaning of words is represented by vectors which are constructed from the co-occurrences of a word of interest with a set of context words. In tensorial compositional distributional semantics [3,4,5,6,7], different types of words, depending on their grammatical role, are associated with vectors, matrices or higher rank tensors. In [8,9] we initiated a study of the statistics of these tensors in the framework of matrix/tensor models. We focused on matrices associated with adjectives or verbs, constructed by a linear regression method, from the vectors for nouns and for adjective-noun composites or verb-noun composites.
We developed a 5-parameter Gaussian model, The parameters J S , J 0 , a, b, Λ are coefficients of five linearly independent linear and quadratic functions of the D 2 random matrix variables M i,j which are permutation invariant, i.e. obey the equation for σ ∈ S D , the symmetric group of all permutations of D distinct objects. This S D invariance implements the notion that the meaning represented by the word-matrices is independent of the ordering of the D context words. General observables of the model are polynomials f (M) obeying the condition (1.2). At quadratic order there are 11 linearly independent polynomials, which are listed in Appendix B of [8]. A three dimensional subspace of quadratic invariants was used in the model above. The most general Gaussian matrix model compatible with S D symmetry considers all the eleven parameters and allows coefficients for each of them. What makes the 5-parameter model relatively easy to handle is that the diagonal variables M ii are each decoupled from each other and from the off- where EXP is the product of exponentials in (1.1). Representation theory of S D offers the techniques to solve the general permutation invariant Gaussian model. The D 2 matrix elements M ij transform as the tensor product V D ⊗ V D of two copies of the natural representation V D . We first decompose V D ⊗ V D into irreducible representations of the diagonal S D .
The trivial (one-dimensional) representation V 0 occurs with multiplicity 2. The (D − 1)dimensional irreducible representation (irrep) V H occurs with multiplicity 3. occurs with multiplicity 1. As a result of these multiplicities, the 11 parameters can be decomposed as 11 = 1 + 1 + 3 + 6 (1.5) 3 is the size of a symmetric 2 × 2 matrix. 6 is the size of a symmetric 3 × 3 matrix. More precisely the parameters form where R + is the set of real numbers greater or equal to zero, M + r is the space of positive semi-definite matrices of size r. Calculating the correlators of this Gaussian model amounts to inverting a symmetric 2 × 2 matrix, inverting a symmetric 3 × 3 matrix, and applying Wick contraction rules, as in quantum field theory, for calculating correlators. There is a graph basis for permutation invariant functions of M. This is explained in Appendix B of [8] which gives examples of graph basis invariants and representation theoretic counting formulae which make contact with the sequence A052171 -directed multi-graphs with loops on any number of nodes -of the Online Encyclopaedia of Integer Sequences (OEIS) [10].
In this paper we show how all the linear and quadratic moments of the graph-basis invariants are expressed in terms of the representation theoretic parameters of (1.6). We also show how some cubic and quartic graph basis invariants are expressed in terms of these parameters. These results are analytic expressions valid for all D.
The paper is organised as follows. Section 2 introduces the relevant facts from the representation theory of S D we need in a fairly self-contained way, which can be read with little prior familiarity of rep theory, but only knowledge of linear algebra. This is used to define the 13-parameter family of Gaussian models (equations (2.71) ,(2.72), (2.73)). Section 3 calculates the expectation values of linear and quadratic graph-basis invariants in the Gaussian model. Sections 4 and 5 describe calculations of expectation values of a selection of cubic and quartic graph-basis invariants in the model.

General permutation invariant Gaussian Matrix models
We solved a permutation invariant Gaussian Matrix model with 2 linear and 3 quadratic parameters [8], obtaining analytic expressions for low order moments of permutation in-variant polynomial functions of a matrix variable as a function of the 5 parameters (section 6 of [8]). The linear parameters are coefficients of linear permutation invariant functions of M and the quadratic parameters (denoted Λ, a, b) are coefficients of quadratic functions. We explained the existence of a 2 + 11 parameter family of models, based on the fact that there are 11 linearly independent quadratic permutation invariant functions of a matrix. The general 2 + 11-parameter family of models can be solved with the help of techniques from the representation theory of S D . We give a brief informal description of the key concepts we will use here. Further information can be found in [12,11,13,14], and we will give more precise references below. A representation of a finite group G is a pair (V, D V ) consisting of a vector space V and a homomorphism D V from G to the space of invertible linear operators acting on V.
Physicists often speak of a representation V of G, when the accompanying homomorphism is left implicit. The homomorphism associates to each g ∈ G a linear operator D V (g). Distinguished among the representations of G are the irreducible representations (irreps). It is known that any representation of G is isomorphic to a direct sum of irreducible representations. For further explanations of these statements see Lecture 1 of [12]. When a representation V is a direct sum of V 1 , V 2 , · · · , V k , we express this as This implies that the linear operators D V (g) corresponding to group elements g ∈ G can, after an appropriate choice of basis in V , be put in a block diagonal form where the blocks are D V 1 (g), D V 2 (g), · · · , D V k (g). The problem of finding this change of basis is called "reducing the representation V into a direct sum of irreducibles". Given two representations (V 1 , D V 1 ) and (V 2 , D V 2 ) of G, the tensor product space V 1 ⊗ V 2 is a representation of the product group G × G, which consists of pairs (g 1 , g 2 ) with g 1 , g 2 ∈ G. The product group G × G has a subgroup of pairs (g, g) which is called the diagonal subgroup of G, denoted Diag(G). The tensor product space V 1 ⊗ V 2 is also a representation of this diagonal subgroup (see for example Chapter 1 of [14]). The linear transformation which reduces V 1 ⊗ V 2 into a direct sum of irreducibles of Diag(G) is called the Clebsch-Gordan decomposition. The matrix elements of the transformation are called Clebsch-Gordan coefficients. More details on these can be found in Chapter 5 of [11]. These can be used to construct projection operators for the subspaces of the tensor product space corresponding to particular irreducible representations.
In section 2.1, we introduce the natural representation V D of S D . We note that the space of linear combinations of the matrix variables M ij is isomorphic as a vector space to V D ⊗V D . We recall the known fact that V D is isomorphic to a direct sum of two irreducible representations and give the explicit change of basis which demonstrates this isomorphism. The tensor product is thus isomorphic to a direct sum This leads to the definition (Equation (2.23)) of S D ×S D covariant variables S 00 , S 0H a , S H,0 a , S HH ab , which correspond to the four terms in the expansion (2.2).
In section 2.2, we describe the space of linear combinations of M ij as a representation of Diag(S D ): The multiplicity index α keeps track of the fact the same irrep appears multiple times in the decomposition into irreducibles of Span (M ij ). The isomorphism of representations of S D above implies the identity relating the dimensions This decomposition leads to the definition, in equations (2.51), (2.52) , (2.53), of variables S V i ;α transforming according to the decomposition (2.3). The next key observation is to think about the vector space of quadratic polynomials in indeterminates {x 1 , x 2 , · · · , x N } in a way which is amenable to the methods of representation theory. Consider a vector space V N spanned by x 1 , x 2 , · · · , x N . The quadratic polynomials are spanned by the set of monomials x i x j which contains N(N + 1)/2 elements. The vector space can be identified with the subspace of the tensor product V N ⊗V N which is invariant under the exchange of the two factors using the map In section (2.2) we apply this observation to the space of quadratic polynomials in the matrix variables M ij . They form a vector space which is isomorphic to Sym 2 (V D ⊗ V D ). Using the decomposition (2.3), we are able to find the S D invariants by using a general theorem about invariants in tensor products of irreducible representations. For two irreps V R , V S , the tensor product V R ⊗ V S contains the trivial representation of the diagonal S D only if R = S, i.e. V R is isomorphic to V S , and further it is also known that this invariant appears in the symmetric subspace Sym 2 (V R ) ⊂ (V R ⊗ V R ). For further information on this useful fact, the reader is referred to Chapter 5 of [11].
This culminates in section (2.3) in an elegant representation theoretic description of the quadratic invariants in the matrix variables, using the linear combinations S V i ;α . With this description in hand, we introduce a set of representation theoretic parameters for the 13-parameter Gaussian matrix models, see equations (2.71) and (2.72). In terms of these parameters, the linear and quadratic expectation values of S V i ;α are simple (see equations (2.75), (2.76), (2.77)). The computation of the correlators of low order polynomial invariant functions of the matrices then follows using Wick's theorem from quantum field theory (see for example Appendix A of [26]).

Matrix variables
on the basis vectors and extended by linearity. With this definition, ρ V D is a homomorphism from S D to linear operators acting on V D We introduce an inner product (. , .) where the e i form an orthonormal basis (e i , e j ) = δ ij . (2.8) We can form the following linear combinations . . .
since, for any σ, we have e σ −1 (1) + e σ −1 (2) + · · · + e σ −1 (D) = e 1 + e 2 + · · · + e D . (2.11) Thus the one-dimensional vector space spanned by E 0 is an S D invariant vector subspace of V D . We can call this vector space V 0 . The vector space spanned by E a , where 1 ≤ a ≤ (D − 1), which we call V H , is also an S D -invariant subspace We have a matrix D H (σ) with matrix elements D H ab (σ) such that These matrices are obtained by using the action on the e i and the change of basis coefficients. The vectors E A for 0 ≤ A ≤ D − 1 are orthonormal under the inner product (2.8) (2.14) All the above facts are summarised by saying that the natural representation V D of S D decomposes as an orthogonal direct sum of irreducible representations of S D as By reading off the coefficients in the expansion of the E 0 , E a in V H , we can define the coefficients using the inner product (2.8). They are . (2.17) The orthonormality means that The last equation implies that As we will see, this function F (i, j) will play an important role in calculations of correlators in the Gaussian model. It is the projector in V D for the subspace V H , obeying Now we will use these coefficients C A,i to build linear combinations of the matrix elements M i,j which have well-defined transformation properties under S D × S D . Define The a, b indices range over 1 · · · (D − 1). These variables are irreducible under S D × S D , Under the diagonal S D , the first three transform as V 0 , V H , V H while S HH ab form a reducible representation. Conversely, we can write these M variables in terms of the S variables, using the orthogonality properties of the C 0,i , C a,i , The next step is to consider quadratic products of these S-variables, and identify the products which are invariant. In order to do this we need to understand the transformation properties of the above S variables in terms of the diagonal action of S D . It is easy to see that S 00 is invariant. S 0H a and S H0 a both have a single a index running over {1, 2, · · · , (D − 1)}, and they transform in the same way as V H . The vector space spanned by S HH ab form a space of dimension (D − 1) 2 which is Permutations act on this as • The representation space V H ⊗ V H can be decomposed into irreducible representations (irreps) of the diagonal S D action as (2.27) In Young diagram notation for irreps of S D , listing the row lengths of the Young diagram, we have • The vector D−1 a=1 E a ⊗ E a is invariant under the diagonal action of σ on V H ⊗ V H . Using the fact that V H is a subspace of V D described by the coefficients C a,i defined in (2.16), the action of σ on V H is given by These can be verified to satisfy the homomorphism property We also have D H ab (σ −1 ) = D H ba (σ). Using these properties, we can show that a E a ⊗ E a is invariant under the diagonal action. The vector has unit norm, using the inner product on V D ⊗ V D obtained from (2.8), and defines a normalized vector in the V 0 subspace of the direct sum decompsotion of V H ⊗ V H given in (2.27). From this expression, we can read off the Clebsch-Gordan coefficients for the trivial representation (2.32) Using these we define S HH→V 0 as (2.34) The coefficients C H,H→H b, c ; a are some representation theoretic numbers (called Clebsch-Gordan coefficients ) which satisfy the orthonormality condition As shown in Appendix B, these Clebsch-Gordan coefficients are proportional to It is a useful fact that the Clebsch-Gordan coefficients for V H ⊗ V H → V H can be usefully written in terms of the C a,i describing V H as a subspace of the natural representation. This has recently played a role in the explicit description of a ring structure on primary fields of free scalar conformal field theory [15]. It would be interesting to explore the more general construction of explicit Clebsch-Gordan coefficients and projectors in the representation theory of S D in terms of the C a,i .
• Similarly for V 2 , V 3 we have corresponding vectors and Clebsch-Gordan coefficients where a ranges from 1 to Dim We have the orthogonality property Here the a, a 1 , a 2 runs over 1 to (D−1)(D−2) 2 .
• The projector for the subspace of (2.40) The quadratic invariant corresponding to V 2 is The quadratic invariant corresponding to V 3 is similar. We just have to calculate • The following is an important fact about invariants. Every irreducible representation of S D , let us denote it by V R has the property that contains the trivial irrep once. This invariant is formed by taking the sum over an orthonormal basis A e V A ⊗ e V A . The invariance is proved as follows In the first equality we have used the definition of the diagonal action of σ on the tensor product space.
• To summarize the matrix variables M ij can be linearly transformed to the following variables, organised according to representations of the diagonal S D Trivial rep: (2.50) • For convenience, we will also use simpler names where we introduced labels 1, 2 to distinguish two occurrences of the trivial irrep V 0 in the space spanned by the M ij . The variables S 0,0 , S H,H→V 0 were first introduced in (2.23) and (2.33) respectively. We will also use where we introduced labels 1, 2, 3 to distinguish the three occurrences of V H in the space spanned by M ij . The variables S 0H a , S H0 a were introduced earlier in (2.23). For the multiplicity-free cases, we introduce (2.53) The M ij variables can be written as linear combinations of the S variables. Rep-basis expansion of M ij is (2.54) In going from first to second line, we have used the fact that the transition from the natural representation to the trivial representation is given by simple constant coefficients In the third line, we have used the Clebsch-Gordan coefficients for For V = V 0 , which is one dimensional, we just have in accordance with (2.31). The index c ranges over a set of orthonormal basis vectors for the irrep V , i.e. extends over a range equal to the dimension of V , denoted DimV . It is now useful to collect together the terms corresponding to each (2.58) Using the notation of (2.51), (2.52), (2.53) , we write this as (2.59) • The discussion so far has included explicit bases for V H inside V D which are easy to write down. A key object in the above discussion is the projector F (i, j) defined in (2.21). For the irreps V 2 , V 3 which appear in V H ⊗V H , we will not need to write down explicit bases. Although Clebsch-Gordan coefficients for H, H → V 2 and H, H → V 3 appear in some of the above formulae, we will only need some of their orthogonality properties rather than their explicit forms. The projectors for V 2 , V 3 in V D ⊗ V D can be written in terms of the F (i, j), and it is these projectors which play a role in the correlators we will be calculating.

Representation theoretic description of quadratic invariants
With the above background of facts from representation theory at hand, we can give a useful description of quadratic invariants. Quadratic invariant functions of M ij form the invariant subspace of  (2.64) We introduced parameters (Λ H ) αβ forming a symmetric 3 × 3 matrix. When we define the general Gaussian measure, we will see that this matrix will be required to be a positive definite matrix. The quadratic invariants constructed from the V 2 , V 3 variables are When we define the general Gaussian measure, we will take the parameters

Definition of the Gaussian models
The measure dM for integration over the matrix variables M ij is taken to be the Euclidean measure on R D 2 parametrised by the D 2 variables Now the S B variables are obtained from M A by an orthogonal basis change, and symmetric group properties also allow the matrix to be chosen to be real. This implies that the matrix is orthogonal Hence det J has magnitude 1, and we have the claimed identity (2.67). The model is defined by integration. The partition function is where the action is a combination of linear and quadratic functions.
The expectation values of permutation invariant polynomials f (M) are defined by These expectation values can be computed using standard techniques from quantum field theory, specialised to matrix fields in zero space-time dimensions (See Appendix A for some explanations). Textbook discussions of these techniques are given, for example in [25], [26]. For linear functions, the non-vanishing expectation values are those of the invariant variables, which transform as V 0 under the S D action We have defined variables µ 1 , µ 2 for convenience. The variables transforming according The quadratic expectation values are The delta function means that these expectation values vanish unless the two irreps V i , V j are equal. While δ ab is the identity in the state space for each V i . The fact that the mixing matrix in the multiplicity indices α, β is the inverse of the coupling matrix Λ V is a special (zero-dimensional) case of a standard result in quantum field theory, where the propagator is the inverse of the operator defining the quadratic terms in the action. The decoupling between different irreps follows because of the factorised form of the measure dMe −S in (2.71).
The requirement of an S D invariant Gaussian measure has led us to define variables S V,α , transforming in irreducible representations of S D . The action is simple in terms of these variables. This is reflected in the fact that the above one and two-point functions are simple in terms of the parameters of the model.
When Λ V 2 > 0, Λ V 3 > 0 and Λ H , Λ V 0 are positive-definite real symmetric matrices (i.e real symmetric matrices with positive eigenvalues), then the partition function Z is well defined as well as the numerators in the definition of f (M) . We can relax these conditions, allowing Λ V 2 , Λ V 3 ≥ 0 and Λ H , Λ V 0 positive semi-definite, by appropriately restricting the f (M) we consider. For example, if Λ V 2 = 0, we consider functions f (M) which do not depend on S V 2 , which ensures that the ratios defining f (M) are welldefined.
Thus the complete set of constraints on the representation theoretic parameters are With these linear and quadratic expectation values of representation theoretic matrix variables S available, the expectation value of a general polynomial function of M ij can be expressed in terms of finite sums of products involving these linear and quadratic expectation values. This is an application of Wick's theorem in the context of QFT. We will explain this for the integrals at hand in Appendix A and describe the consequences of Wick's theorem explicitly for expectation values of functions up to quartic in the matrix variables. We will be particularly interested in the expectation values of polynomial functions of the M ij which are invariant under S D action and can be parametrised by graphs. While the mixing between between the S variables in the quadratic action is simple, there are non-trivial couplings between the D 2 variables M ij if we expand the action in terms of the M variables. This will lead to non-trivial expressions for the expectation values of the graph-basis polynomials.
These expectation values were computed for the 5-parameter Gaussian model in [8]. They were referred to as theoretical expectation values f (M) , which were compared with experimental expectation values f (M) EXP T . These experimental expectation values were calculated by considering a list of words labelled by an index A ranging from 1 to N, and their corresponding matrices M A , We will now proceed to explicitly apply Wick's theorem to calculate the expectation values of permutation invariant functions labelled by graphs for the case of quadratic functions (2-edge graphs), cubic (3-edge graphs) and quartic functions (4-edge graphs). We will leave the comparison of the results of this 13-parameter Gaussian model to linguistic data for the future.

3
Graph basis invariants in terms of rep theory parameters will play an important role in the following. Its meaning is that it is the projector for the hook representation in the natural representation. Deriving expressions for expectation values of permutation invariant polynomial functions of the matrix variable M amounts to doing appropriate sums of products of F factors, with the arguments of these F factors being related to each other according to the nature of the polynomial under consideration. In terms of the variables defined in Section 2, repeated here for convenience, All the terms can be expressed in terms of the F -function defined in (3.1) We will refer to the terms depending on Λ V 0 as V 0 -channel contributions to the 2-point functions, those on Λ H as V H -channel (or H-channel) contributions, those on Λ V 2 as V 2channel and those on Λ V 3 as V 3 -channel contributions. It will be convenient to denote these different channel contributions as In arriving at the expressions for the last two terms in (3.6), we used the fact that these terms in (3.5) can be expressed as . (3.10) The factor D D−2 is explained in Appendix C.

Calculation of i,j M ij M ij
Following (3.7) the expectation value M ij M ij can be written as a sum over V -channel contributions, where V ranges over the four irreps.

Contributions from
The projector has eigenvalue 1 on the subspace transforming in the irrep V 2 and zero elsewhere, hence the (Dim V 2 ). Similarly we can also write the trace in V H ⊗ V H and express this in terms of irreducible characters In the last line we have used the fact that the Kronecker coefficient for V H ⊗ V H → V 2 is 1.

Contribution from V 0 channel
This is (3.17)

Summing all channels
The disconnected piece is (3.20)

3.2
Calculation of i,j M ij M ji As in (3.7) the expectation value M ij M ij can be written as a sum over V -channel contributions, where V ranges over the four irreps.

3.2.1
Contribution from multiplicity 1 channels V 2 , V 3 The τ is the swop which acts on the two factors of V H . We have used the fact that V 2 appears in the symmetric part of We use the fact that V 3 is the antisymmetric part of V H ⊗ V H .

Contribution from V 0 channel
The first term from (3.5) is The second term is The third term is which vanishes using (2.19). The last term vanishes for the same reason. So collecting the V 0 -channel contributions to i,j M ij M ji conn , we have

Calculation of i,j M ii M ij
An important observation here is that the sum over j projects the representation V D to the trivial irrep V 0 , which follows from the formula for C 0i in (2.17). This means that when we expand M ii and M ij into S variables as in the first line of (2.54), we only need to keep the term S H0 or S 00 from the expansion of M ij .

3.3.1
Contribution from V 2 , V 3 channels From the above observation, and since V 2 , V 3 appear only in S HH S HH , we immediately see that

Contribution from V H channel
From the above observation, the only non-zero contributions in the V H channel come from S HH→H S H0 , S H0 S H0 and S 0H S H0 . These are Note that Λ H is a symmetric 3 × 3 matrix and Λ −1 H ) 23 = (Λ −1 H ) 32 . In the penultimate step, we have used the normalization equation (C.9). These add up to (3.37)

3.3.3
Contribution from V 0 channel The non-zero contributions come from S 00 S 00 and S HH→0 S 00 . They are

Summing all channels
The disconnected piece is with the first term given by (3.39)

Calculation of i,j M ii M ji
We can write down the answer from inspection of (3.39) The reasoning is as follows. The sum over j projects to V 0 . This means that the only non-zero contributions are, from the V 0 channel, S 00 S 00 conn S HH→0 S 00 conn (3.43) and from the V H channel This identifies the contributing entries of Λ −1 V 0 , Λ −1 H using the indexing in (2.51) and (2.52). given the similarity between the expectation value in section 3.3, we have contributions of the same form, up to taking care of the right indices on Λ −1 V 0 , Λ −1 H . Given the symmetry of F (i, j) under exchange of i, j, the disconnected piece is the same as above The sums over j, k project to V 0 . The non-vanishing contributions are S 00 S 00 conn and S H0 S H0 conn . They add up to The disconnected part is Now we are projecting to V 0 on the first index of both M's. This means that the contributing terms are S 00 S 00 and S 0H S 0H . Repeat the same steps as above in (3.46) to get The only difference is that we are picking up the (1, 1) matrix element of (Λ −1 H ) instead of the (2, 2) element, since we defined S 0H = S {V 0 ;1} and S H0 = S {V 0 ;2} .
Adding the disconnected piece, which is the same as (3.47), we have We are now projecting to V 0 on first index of one of the matrices and second index of the other. Hence the contributing terms are S 00 S 00 and S 0H S H0 . The result is Here we project to V 0 on all four indices, so Adding the disconnected piece we have The contribution from the V 0 channel is given by (3.55)

3.9.2
The V H channel It is convenient to use (3.6) to arrive at (3.56) Useful equations in arriving at the above are the sums which can be obtained by hand or with the help of Mathematica. In the latter case, it is occasionally easier to evaluate for a range of integer D and fit to a form P olynomial(D)/D some power .

The V 2 , V 3 channels
Now calculate the HH → V 2 and HH → V 3 channel.
We used the fact that the projectors for H, V 2 are orthogonal.
Similarly, the contribution from V 3 is zero. Another to arrive at the same answer is to recognise that V 3 is the antisymmetric part, so P H,H→V 3 ab;cd = 1 2 (δ ac δ bd − δ ad δ bc ) . (3.59)

Calculation of i,j M ii M jj
Since i M ii and j M jj are S D invariant, we only have contributions from the V 0 channel. Use the first four terms of (3.5) to get (3.62)

3.11
Calculation of i,j,k M ii M jk Here we get contributions from S HH→0 S 00 conn and S 00 S 00 conn . Adding these up from (3.5) The disconnected part is

Summary of results for quadratic expectation values in a large D limit
It is interesting to collect the results for the connected quadratic expectation values and consider the large D limit. Let us assume that all the Λ V 0 , Λ H , Λ V 2 , Λ V 3 scale in the same way as D → ∞ and consider the sums normalized by the appropriate factor of D The dominant expectation values in this limit are the first, second and ninth. These are the quadratic expressions which enter the simplified 5-parameter model considered in [8] (see Equation (1.1)). It will be interesting to systematically explore the different large D scalings of the parameters in real world data, e.g. the computational linguistics setting of [8] or in any other situation where permutation invariant matrix Gaussian matrix distributions can be argued to be appropriate.

A selection of cubic expectation values
In this section we use Wick's theorem from Appendix A to express expectation values of cubic functions of matrix variables in terms of linear and quadratic expectation values. The permutation invariance condition requires sums of indices over the range {1, · · · , D}. This leads to non-trivial sums of products of the natural-to-hook projector F (i, j). The invariants at cubic order are 52 in number (Appendix B of [8]) and correspond to graphs with up to 6 nodes.

1-node case
Since this is independent of i, we can use (3.60) to get Using A.6 we have Calculating this requires doing a few sums, which can be done by hand or with Mathematica (the function KoneckerDelta is handy).
Using (3.3), we find for the second term in (4.4) For the first term on the RHS of (4.4) The first term in (4.7) can be expressed as a function of the parameters of the Gaussian model using (3.18). The second term is calculated by specialising the fundamental quadratic moments (3.6) and doing the resulting sums over the F -factors. Consider the V 0 contributions to the second term above. The term proportional to (Λ −1 V 0 ) 11 vanishes due to the first of (4.5). The (22) contribution, using the third of (4.5) is The (12) contributions, using the second of (4.5) is Now consider the V H contribution to the second term in (4.7). The ( The sum of products of six F 's is readily done with Mathematica to give Contributions from the (1, 2) matrix element of symmetric matrix Λ H give This uses the second of (4.5). From (1, 3) and (3, 1) we have From (2, 3) and (3, 2), we have (4.14) Now consider the contribution from V 2 . It is convenient to use (3.5) This means that When the expression (4.17) is substituted in (4.15) the first term on the RHS of (4.17) does not contribute because of (4.18). The second term gives The third term gives (4.20) The fourth term gives Collecting terms from (4.19), (4.20) and (4.21) we get Multiplying the factor 3 µ 2 √ D−1 from (4.7) to get a contribution to i,j M 3 ij , we get The contribution from V 3 is Now use the fact that (P H,H→H ) aa;bb = 0 a 1 ,a 2 (P H,H→H ) a 1 a 2 ;a 2 a 1 = a 1 ,a 2 (P H,H→H ) a 1 a 2 ;a 1 a 2 = (D − 1) (4.26) to find Collecting all the contributions we have . (4.29) The first term is along with the first and second of (4.5) we can show that Consider the remaining three terms. Focus on the first of these : We already know the first term from (3.51). So let us consider the second. An easy calculation using (3.5) (or equivalently using (3.6)) shows that the contribution from the V 0 channel is From the V H channel, the contributions are From the V 2 channel, we get By relabelling indices, it is easy to see that The sums over i 1 , · · · , i 6 project to the V 0 representations. As a result, using (A.6), along with (2.59), we have (4.40)

A selection of quartic expectation values
The methods we have used to calculate the cubic expectation values, which were explained in detail above, extend straightforwardly to quartic expectation values. The first step is to use Wick's theorem A.7. Then we use the formulae for quadratic and linear expectation values from Sections 2 and 3. In order to arrive at the final result as a function of D, µ 1 , µ 2 , Λ V 0 , Λ H , Λ V 2 , Λ V 3 we have to do certain sums over products of the natural-to-Hook projector F (i, j). We will give some formulae below to illustrate these steps for the quartic case, without producing detailed formulae as in previous sections.

A 2-node quartic expectation value
The quadratic average is Using this and we can work out the formula for i,j M 4 ij as a function of the 13 Gaussian model parameters. Mathematica would be handy in doing the sums over products of F (i, j) which arise.

5.2
A 5-node quartic expectation value i,j,k,p,q M ij M jk M kp M pq From A.7, we have All the summands on the RHS can be evaluated using 3.5 or 3.6 in terms of F (i, j). The sums can be done with the help of Mathematica to obtain expressions in terms of D, µ 1 , µ 2 , Λ V .

Summary and Outlook
We have used the representation theory of symmetric groups S D in order to define a 13-parameter permutation invariant Gaussian matrix model, to compute the expectation values of all the graph-basis permutation invariant quadratic functions of the random matrix, and a selection of cubic and quartic invariants. In [8] analogous computations with a 5-parameter model were compared with matrix data constructed from a corpus of the English language. A natural direction is to extend that discussion of the English language, or indeed other languages, to the present 13-parameter model. Combining the experimental methods employed in [8] with machine learning methods such those used in [19], in the investigation of the 13-parameter model, would also be interesting to explore. As a theoretical extension of the present work, it will be useful to generalise the representation theoretic parametrisation of the Gaussian models to perturbations of the Gaussian model, where we add cubic and quartic terms to the Gaussian action. Identifying parameter spaces of these deformations which allow well-defined convergent partition functions and expectation values will be useful for eventual comparison to data. If we ignore the convergence constraints, the general perturbed model at cubic and quartic order has 348 parameters, since there are 52 cubic invariants and 296 quartic invariants (Appendix A of [8]). As in the Gaussian case, we can expect that representation theory methods will be useful in handling this more general problem. Further techniques involving partition algebras underlying the representation theory of tensor products of the natural representation will likely play a role (see e.g. [20] for recent work in these directions).
It is worth noting that permutation invariant random matrix distributions have been approached from a different perspective, based on non-commutative probability theory [21,22,23]. The approach of the present paper and [8] is based on the connection between statistical physics and zero dimensional quantum field theory (QFT). It would seem that the approach of the present paper can complement the theory developed in these papers [21,22,23] by producing integral representations (Gaussians or perturbed Gaussians) of random matrix distributions having finite expectation values for permutation invariant polynomial functions of matrices. The results on the central limit theorem from the above references would be very interesting to interpret from the present QFT perspective.
The computation of expectation values in Gaussian matrix models admits generalization to higher tensors. Indeed the motivating framework in computational linguistics discussed in [8] involves matrices as well as higher tensors in a natural way. Generalizations of the present work on representation theoretic parametrisation of Gaussian models and computation of graph-basis observables to the tensor case is an interesting avenue for future research.
In this paper, we have focused on the explicit computation of permutation invariant correlators for general D. Some simplifications at large D were discussed in section 3.12. For traditional matrix models having U(D) ( or SO(D)/Sp(D) symmetries), there is a rich geometry of two dimensional surfaces and maps in the large D expansions which allows these expansions of matrix quantum field theories to have deep connections to string theory [28,29]. It will be interesting to explore the possibility of two dimensional geometrical interpretations of the large D expansion in permutation invariant matrix models.

A Multi-dimensional Gaussian Integrals and Wick's theorem
Consider the multi-variable integral with a Gaussian integrand x ∈ R N . A ∈ C N ×N is a real symmetric positive definite matrix. s ∈ R N is an arbitrary complex vector (see for example [17], [18], Appendix A, Equations (8) and (9) of [26]). One can also consider A more generally to be complex with positive definite real part, but to keep a probabilistic interpretation we keep A real symmetric. Expectation values of functions f (x) are defined by These expectation values can be calculated by taking derivatives with respect to s i on both sides of (A.1). For the x variables Application of this equation, along with the formula for dM in terms of the representation theoretic S-variables (2.67) leads to (2.74),(2.76). For expectation values of quadratic monomials we have We define the connected part as The expressions (2.77) and (2.4) follow from these. For cubic expressions For quartic expressions x i x j x k x l = x i x j conn x k x l conn + x i x k conn x j x l conn + x i x l conn x j x k conn + x i x j conn x k x l + x i x k conn x j x l + x i x l conn x j x k + x j x k conn x i x l + x j x l conn x i x k + x k x l conn x i x j + x i x j x k x l . (A.7) These illustrate a general fact (known as Wick's theorem in the quantum field theory context and Isserlis' theorem in probablity theory [27]) about Gaussian expectation values. Higher order expectation values can be expressed in terms of linear and quadratic expectation values. When applied to permutation invariant matrix models, we still have non-trivial sums left to do, after Wick's theorem has been applied. This is illustrated in the calculations of section 4 and section 5.

B
Rep theory of V H and its tensor products Some basics of rep theory of V H can be presented in a self-contained way, assuming only knowledge of linear algebra and index notation. Alternatively, we can observe that the matrices in V H are the same as Young's orthogonal basis. If we just follow the self-contained route, we define We have the orthogonality property : The homomorphism property It is useful to define C σ(a),i = b D H ba (σ)C bi so the above can be expressed as an equivariance property C σ(a),i = C a,σ(i) . (B.7) which is an equivariance condition for the map V H → V nat given by the coefficients C a,i . This map intertwines the S n action on the V H and V nat . Now define C a,b,c We show that this is an invariant tensor.
showing that the transformation is indeed by the matrix D H . These vectors E a are orthogonal. It is useful to calculate the inner product (B.14) which will be useful in the next section. We also know that (C.9) We can therefore identify