On the Phase Structure of Commuting Matrix Models

We perform a systematic study of commutative $SO(p)$ invariant matrix models with quadratic and quartic potentials in the large $N$ limit. We find that the physics of these systems depends crucially on the number of matrices with a critical r\^ole played by $p=4$. For $p\leq4$ the system undergoes a phase transition accompanied by a topology change transition. For $p>4$ the system is always in the topologically non-trivial phase and the eigenvalue distribution is a Dirac delta function spherical shell. We verify our analytic work with Monte Carlo simulations.


Introduction
Multi-matrix models arise in a wide variety of settings from matrix string theory [1] to M-theory. Although a non-perturbative formulation of M-theory in terms of its fundamental degrees of freedom is still lacking, the best candidate for such a formulation appears to be the infinite matrix size limit of a matrix model of some kind. The leading candidate for such a formulation is the BFFS model [2,3] which was conjectured to capture the entire dynamics of M-theory. Relatives of this model such as the BMN model [4] or models derived from the ABJM model 1 [5,6] are also considered possible viable candidates for such a non-perturbative formulation. An alternative matrix model formulation is based on the IKKT model [7] which is a pure matrix model where even time has disappeared.
All of these conjectured formulations of M-theory are regularised versions of the supermembrane. They are based on the matrix regularisation of membranes introduced by Hoppe [8] and extended to the supermembrane in [3] and [9]. They also arise as dimensionally reduced 4-dimensional or 3-dimensional supersymmetric field theories.
Multi-matrix models further arise in lower dimensional variants of the IKKT model [10], in the low energy dynamics of D-branes [11] and simple models of emergent geometry [12,13] and emergent gravity [14,15]. Many of these models will have regimes where commuting matrices play a rôle.
In [16,17] it was established that the unique rotationally invariant three dimensional joint eigenvalue distribution that corresponds to a parabolic one dimensional distribution is the uniform distribution within a ball of radius R. It was also established that the strong coupling limit of Hoppe's two matrix model [8] which describes the low energy dynamics of D0-branes [11] in N = 1 supersymmetric Yang-Mills in four dimensions was captured by commuting matrices. In part the motivation for the current paper arose from this earlier work coupled with a desire to understand commuting matrices in and of themselves.
To our knowledge no systematic study of commutative matrix models with general potential, has been undertaken prior to the current work. An understanding of commutative matrix models fills a gap in the literature and because of the simplicity of these systems the results may prove useful in a wider context. Such models, of course, also have an intrinsic interest in their own right.
In this paper we show that due to rotational invariance we can recover the full joint eigenvalue distribution from that of the one matrix distribution, but only when the eigenvalue distribution of the full system is topologically trivial. We begin by studying Gaussian distributions (considered previously in refs. [18][19][20][21]) and find that the generalisation of the Wigner distribution for p = 1 becomes the uniform distribution within a disk for p = 2, but for p = 3 the distribution is ρ (3) 3 − x 2 , which is divergent at the boundary but still integrable. We find a special rôle is played by p = 4 as it is the critical dimension where the distribution is a Dirac delta function on the unit sphere: ρ (4) . For all p > 4 only spherical shells occur and the Gaussian distributions of commuting matrices have eigenvalue distributions ρ p ( x) = 2 Ω p−1 δ(1− x 2 ) where Ω p−1 is the volume of the unit p−1-sphere. When considering models with the quartic potential, a| x| 2 + b| x| 4 , we find that for p = 1, 2 and 3 the system has a phase transition at the critical values a c = −2 √ b for p = 1; a c = 0 for p = 2 and the surprising positive value a c = √ 20b 3 for p = 3. There is no transition for p > 4, rather the distribution is concentrated on the sphere irrespective of the potential.
The principal results of this paper are: • We find that there is a special rôle played by p = 4, it is the critical dimension where shell solutions become the energetically preferred eigenvalue configurations.
• The eigenvalue distributions for Gaussian ensembles of p rotationally invariant commuting matrices with p = 2, 3 and 4 can be obtained by lifting the Wigner semi-circle distribution. The distribution for p = 4 is a spherical δ-function shell. The distributions for p > 4 are δ-function shells but cannot be obtained by lifting the Wigner distribution. We derive an analytic technique for the reduction (or lifting) of commuting models with arbitrary rotationally invariant potentials.
• Commuting matrices with quartic potential V ( x) = a| x| 2 + b| x| 4 have phase transitions of 3rd order for p = 1, 6th order for p = 2 and 4th order for p = 3. In these transitions the eigenvalue density undergoes a one-cut to two-cut transition for p = 1, a disk to annulus transition for p = 2 and a ball to shell transition for p = 3. For p ≥ 4 there is a phase transition from a spherical shell to a metastable phase comprising of a mixture of shell and uniform distributions. The metastable phase exists only for negative b and sufficiently large a.
• The critical transitions occur at a c = −2 √ b for p = 1, a c = 0 for p = 2 and a c = √ 20b/3 for p = 3. For p = 4 the metastable shell-mixture transition occurs at b c = 0 with a 2 > |6b|. There is also an instability transition at a 2 = −6b. For all p > 4 and b > 0 the strong eigenvalue repulsion forces all of the eigenvalues onto a shell and there is no transition.
The structure of the paper is as follows: In Section 2 we describe the family of commuting matrix models we consider and obtain the integral equation satisfied by the joint eigenvalue distribution for these systems in the large matrix size limit. We further show how the eigenvalue density, integral kernel and effective action can be reduced to a lower dimensional system and lifted back to the original dimension due to rotational invariance.
In section 3 we study Gaussian systems in different dimensions. We show that for p = 2, 3 and 4 the eigenvalue distribution is simply the rotationally invariant lift of the Wigner semicircle. We further show that a further lift to p = 5 does not yield a normalisable positive distribution, however we establish by studying the effective action that the least action is given by spherical shells. Spherical shells are the preferred distributions for all p > 4. We finish the section by confirming this conclusion with Monte Carlo simulations.
In Section 4 we study the quartic potential V ( x) = a| x| 2 + b| x| 4 and study the phase structure of these systems. We find that the well known 3rd order transition at a c = −2 √ b of the p = 1 model becomes a 6-th order transition for p = 2 and occurs at a c = 0 while for p = 3 the transition occurs at the positive value a c = √ 20b/3 and is fourth order. We conclude the section by showing that for p = 4 there is a phase transition from a spherical shell eigenvalue distribution to a metastable phase comprising of a mixture of shell and uniform distributions. The metastable phase exists only for negative b and sufficiently large a.
The paper finishes with our conclusions and discussion in Section 5.
The results of this paper should have applications wherever an ensemble of commuting matrices form a good approximation.

Commuting matrix model 2.1 The model
We consider a commuting SO(p) invariant p-matrix model with partition function: where X is an array of p, N × N commuting hermitian matrices,D X is the corresponding invariant measure and V ( X) is an SO(p) invariant potential. The set of commuting hermitian matrices X, can be parameterised by a set of real diagonal matrices Λ and an unitary matrix U : The corresponding Jacobian is given by: where θ = U † dU and u lm are coordinates on SU (N ). The partition function (2.1) can be written as: The resulting effective action (we divide by N 2 ) for the eigenvalues λ is: (2.5) At large N the dynamics is dominated by the saddle point. Varying with respect to λ i we obtain: Equation (2.6) determines the eigenvalue distribution in the large N limit and admits rotationally invariant shell solutions. The only shell solution consistent with SO(p) invariance is a p−1 dimensional spherical shell. These solutions have been considered in refs. [18][19][20][21] for gaussian potential, where the authors argued that the radius of the spherical shell is independent on the number of the commuting matrices. One can show that the same holds for any potential. Indeed, it is straightforward to verify that the vector equation (2.6) is satisfied by a homogeneous spherical eigenvalue distribution of radius R, provided the radius satisfies: Equation (2.6) admits also p-dimensional ("fat") rotationally invariant solutions, which may or may not be energetically favoured relative to the shell solution. To explore these solutions we consider a course grained approximation: and extremize the following functional: Upon variation with respect to ρ we obtain the integral equation: differentiating equation (2.10) with resect to x we obtain: which we recognise as the continuous limit of equation (2.6). This of course is not surprising, since as long as we are dealing with p-dimensional distributions it shouldn't matter when we take the continuous limit. It turns out that instead of directly solving equation (2.11) in p dimensions one can use the properties of the logarithmic kernel in equation (2.9) to reduce the problem to a lower dimensional one. Let us study this in more details.

Reducing the effective action
The rotational invariance of the potential V p ( x) suggests that the system settles in a rotationally invariant eigenvalue distribution, which is fully characterised by its radial distribution. In such cases the distribution can be reduced to lower dimensions without any loss of information. Furthermore, the reduced distribution can also be lifted back to higher dimensions. This opens up the possibility to reduce a higher dimensional problem down to one or two dimensions where it can be analysed more easily, the obtained one-dimensional distribution can then be lifted back to higher dimensions. What makes this approach valuable is that the logarithmic kernel in the effective action (2.9) is preserved under such dimensional reduction. Furthermore, for polynomial potential the reduction just alters the coefficients of the polynomial. This suggests (naively) that the saddle point equation of a given problem can be analysed only in one dimension and the solution to the analogous problem in higher dimensions can be obtained by simply lifting the one dimensional distribution. It turns out that this is the case only for distributions with simple topology and the description breaks down if the distribution undergoes a topology change transition (look at section 4). This still leaves a large class of problems for which reducing the distribution can be useful. To describe how this works let us first focus on the reduction from p to p − 1 and p − 2 dimensions, we have: These relations can be inverted by solving the integral equations (2.12) and (2.13). The result is [16]: Our strategy is to describe how the p dimensional action (2.9) reduces to p − 2 dimensions and then to show how a two dimensional action reduces to one dimension. In this way we can reduce both odd and even dimensional actions down to one dimension. Let us begin by reducing the potential term in (2.9). Using equation (2.15) we obtain: where: and Ω p−1 is the volume of the p−1 dimensional sphere. Note that if V p is a polynomial of x of a certain degree, V p−2 is also a polynomial of the same degree, just the coefficients change according to (2.17). Next we focus on the reduction of the logarithmic kernel in (2.9). Using the rotational invariance of the distribution ρ p (x), we can write: where the kernel K p (x, x ) is given by: Substituting equation (2.15) for ρ p (x) into equation (2.18) and integrating by parts for p > 2 we obtain: (2.20) One can show that: where Ω p−3 is the volume of the unite p − 3 sphere. For the reduced logarithmic term we obtain: (2.22) Using (2.15) one can also show that: Equation (2.23) implies that if ρ p is normalised to one so is ρ p−2 , which suggests that the last term on the right-hand site of equation (2.22) is just the constant 2/(p − 2). Finally defining µ p−2 = µ p we can write: and because the equations of motion for µ p and µ p−2 imply that both ρ p and ρ p−2 are normalised to one, the effective actions S p and S p−2 differ only by a constant and describe equivalent physics. If p is odd and p > 2 one can repeat this procedure until one reduces the problem down to one dimension. For the relation between the effective actions one obtains: where H n is the harmonic number. For the saddle point equation for ρ 1 we obtain: Finally, one can use equation (2.14) (see Appendix A) to show that:

Gaussian model
In this section we focus on the properties of commuting matrix models with a quadratic potential: We begin by studying the joint eigenvalue distributions for various number of commuting matrices.

Gaussian model in various dimensions
Using equation (2.17) one can reduce the potential (3.1) to two or one dimensions depending on whether p is even or odd. In even dimensions one can use equation (A.3) to reduce the potential further to one dimension. It is easy to verify that the reduced potential is: Substituting V 1 into equation (2.26) and differentiating with respect to x we obtain the integral equation: whose solution is a Wigner semi-circle, which if normalised to one has a radius R 2 p = 4/p: Therefore we conclude that for gaussian potential the p-dimensional joint eigenvalue distribution is obtained by lifting a Wigner semi-circle distribution using equations (2.14) and (2.15). Let us see how this works in different dimensions.
For p = 1 we trivially obtain a Wigner semi-circle of radius R 1 = 2.
For p = 2 using equation (2.14) we obtain that the joint eigenvalue distribution is a uniform disk of radius R 2 = √ 2: The distribution (3.5) can easily be obtained directly in two dimensions by using the fact that log( x − x ) 2 is proportional to the Green's function of the laplacian in the two dimensions (see for example ref. [18]).
For p = 3 we use equation (2.15) to lift the Wigner semi-circle (3.4). We obtain: The distribution in equation (3.6) diverges at the boundary, however it is still integrable. In the next subsection we will compare this distribution to Monte Carlo simulations at large (but finite) N and we will confirm that it is indeed approached by the physical distribution in the large N limit.
For p = 4 it is convenient to first lift the Wigner semi-circle (3.4) to two dimensions using equation (2.14) and then lift from two to four dimensions using equation (2.15). One easily obtains: Note that the distribution ρ (4) is a shell and is thus three (rather than four) dimensional. In fact this is the spherical shell saddle point that we analysed in the previous section. Indeed, if we substitute the potential (3.1) into equation (2.7) we arrive at unit radius R 4 = 1. It is intriguing that the equation (3.7) which we derived under the assumption of a four dimensional ("fat") distribution agrees with the derivation of the shell saddle point above equation (2.7). In fact as we are going to see this is no longer the case for dimensions higher than four.
For p > 4 we run into troubles. In even dimensions p = 2n (n > 2) using first equation (2.14) and then equation (2.15) one can show that the distribution is a shell proportional to derivatives of a delta function: , which is not a positive function and cannot represent joint eigenvalue distribution. In odd dimensions p = 2n + 1 (n > 1). The distribution is not integrable. Indeed, using equation (2.15) one can show that ρ (2n+1) ( x) ∝ 1/(4/(2n + 1) − x 2 ) (n−1/2) , which is not integrable near the boundary for n > 1. Therefore we conclude that although the saddle point extremising the effective action (2.9) can be constructed mathematically by lifting the Wigner semi-circle distribution (3.4), for dimensions higher than four (p > 4) the mathematical solutions are not physical and cannot be realised as eigenvalue distributions. However the spherical shell saddles derived in the previous section still exist. It is then natural to conclude that for p > 4 the joint eigenvalue distribution is given by [18]: where Ω p−1 is the volume the unit p − 1-sphere.
Overall, we see that the eigenvalue distribution depends crucially on the number of commuting matrices. The different eigenvalue distributions can be split into two classes: The first class is for p ≤ 4, when the joint eigenvalue distributions are obtained by lifting the Wigner semi-circle distribution (3.4). The second class is for p ≥ 4, when the spherical shell saddles are realised and the radius depends only on the shape of the potential but not on the dimension. Interestingly these two classes overlap at p = 4 since the three-sphere shell can be obtained in both approaches. In the next subsection we analyse this behaviour and argue that it follows from the principle of least action, which should be valid in the large N limit.

Least action analysis
As we observed above for the gaussian potential (3.1) the possible eigenvalue distributions split into two classes. In particular, we showed that for p > 4 the joint eigenvalue distribution does not extremise the effective action (2.9) and is given instead by a spherical shell of unit radius. However we could still reduce the spherical shell to one dimension. One can easily show that the spherical distribution (3.8) reduces to: Inspired by equation (3.9) we will assume that in general the reduced distribution is composed of terms of the form (R 2 − x 2 ) α . If we define η = R/x, then for a very broad class of distributions the reduced distribution can be written as: The normalisation condition dηρ(η) = 1 imposes the following constraint on the coefficients c n and R: It turns out that we can impose one more constraint on c n and R without referring to the saddle point equation (see appendix B for derivation for general potential). For the gaussian potential (3.1) it reads: Applying this to the distribution in (3.10) we obtain: Clearly in general the two constraints in equations (3.11) and (3.13) are not sufficient to determine the coefficients c n and the radius R in equation (3.10). However, they can determine these parameters for pure states, that is when only one of the coefficients c n is non-vanishing. If the non vanishing coefficient is c α one easily obtains: (3.14) Let us consider such a pure state: We will show that for a given p the pure state with the lowest α ≥ 1 has the lowest energy (Note also that in general we could take α to be continuos). To compare the energies of the different pure states we have to evaluate the reduced effective action S 1 , however since the pure states are normalised to one and the potential term is fixed by the constraint (3.13) we need just to evaluate the term with the logarithmic kernel, thus we define: The easiest way to evaluate E α for integer α is to uplift the pure state (3.15) to α + 3 dimensions, where it is a spherical shell and use equation (2.25) to evaluate E α (look at appendix C for a derivation). The result is: Let us calculate the derivative of E with respect to α, we obtain: . Note that ∂E ∂α is independent of p. One can also verify that: for α > 1 one has ∂E ∂α > 0, for −2 < α < 1 one has E ∂α < 0 and finally for α = 1 one has E ∂α = 0. This clearly indicates that for α = 1, E has its minimum (look at figure 1), as it should since for α = 1 the pure state is the Wigner semi-circle (3.4), which extremises the effective action S 1 and after uplift S p . However, as we observed in the previous section for p > 4 the joint eigenvalue is a shell of unit radius, which reduces to a pure state with α > 1. The reason is that the uplift ofρ α from equation (3.15) is physical only up to p = α + 3 dimensions (when it is a shell). A further lift would produce either negative shell (derivative of a delta function) or a non-integrable distribution. Therefore, for p > 4 we cannot lift the Wigner semicircle and a pure state with α > 1 should be realised. Furthermore, since E is monotonically increasing function of α, for α > 1 we should always pick the lowest possible value of α. This suggests that for p = 5 we should pick α = 2, but this pure state can be lifted at most to p = α + 3 = 5 and therefore for p = 6 we should pick the next one: α = 3. Following the same argument again, one concludes that in general for p > 4 the pure state with α = p−3 is realised, which is always a shell as equation (3.9) suggests. Furthermore, using equation (3.14) for the radius of the distribution we have that R p−3 = (p − 3 + 3)/p = 1. We arrive at the result that for p ≥ 4 the radius is independent of the dimension and is equal to one, which is the same result that we obtained above using saddle point arguments. Now we have a better understanding why the spherical shell saddles considered above equation (2.7) are not realised for p < 4. It is because the uplifts of the Wigner semi-circle (3.4) are energetically preferred and whenever they are physical (correspond to positive and integrable distribution) they are realised.
So far our analysis involved only pure states. In general we can have a distribution which is a "mixture" of pure states (see equation (3.10)). However, the pure states have different energies and it is plausible to assume that the pure state with the lowest possible energy will have lower energy than any mixed state since this will involve mixing with pure states of higher energy. Generally this is not true for arbitrary potential. However, for a gaussian potential the above considerations suggest that this is the case. We also explicitly verified that pure states are energetically more favoured than mixed states of two and three pure states and believe that it is true for any mixed state.

Monte Carlo simulation of the gaussian model
In this subsection we perform Monte Carlo simulations of the gaussian model with potential given in equation (3.1). To this end we implemented the algorithm of Metropolis into a C++ commuter program. Over all we find excellent agreement with the distributions derived in section 3.1.
For p = 1 the model is just an ordinary one-matrix model with a Wigner semicircle distribution, therefore we will begin with the p = 2 case. In this case the distribution is a uniform disk of radius √ 2. Numerically it is easier to analyse the radial distribution. Using equation (3.5) and that we are in two dimensions we obtain: In the left panel of figure 2 we have presented our numerical results for the radial distribution (3.19). The red dashed curve represents the N → ∞ result (3.19). One can see that the agreement with the Monte Carlo simulations improves as the size of the matrices increases and at N = 8000 it is already excellent. In the right panel of figure 2 we have present a plot of the reduced distribution (the distribution of one component of the eigenvalue). One can observe the excellent agreement of the numerical result for N = 8000 with the Wigner semi-circle distribution from equation (3.4) for p = 2.
Next we consider the p = 3 case. Using equation (3.6) and that we are in three dimensions, for the radial distribution we obtain: Our next focus is the case p ≥ 4. In section 3.1 we showed that for p ≥ 4 the joint eigenvalue distribution is a spherical shell of unit radius. We also learned that the reduced one-dimensional distribution is given by equation (3.9), which for p = 4  agrees with a Wigner semi-circle, but for p > 4 differs significantly. In figure 4 and figure 5 we have presented our numerical results for p = 4, 5 and p = 6, 7, 8. The left panels represent the radial distributions. One can see that as the size of the matrices is increased the radial distributions approach spherical shells of unit radii. In the right panels we have presented the reduced distributions. One can see the excellent agreement with equation (3.9) for p = 4, 5, 6, 7, 8. These results support the analysis of the previous chapters and that of ref. [18]. We experimented with higher values of p > 8 and found the same behaviour confirming that there are only two classes of solutions the Wigner semi-circle family for p ≤ 4 and the spherical shell distributions for p ≥ 4.

Non-Gaussian potentials
In this section we consider non-gaussian potentials. We will focus on potentials of the form: V ( x) = a| x| 2 + b| x| 4 , containing a quartic term. Note that in order for the model to be stable we have to impose the restriction b ≥ 0, where the value b = 0 is allowed only if a is positive 2 .

Quartic potential in one dimension.
It is instructive to review the properties of a one dimensional matrix model with potential of the form (4.1). The one dimensional random matrix versions this model has been extensively studied in the literature [23,24] and it has been shown that as the parameters of the potential are varied, the model undergoes a phase transition. This phase transition is reflected in a change of the topology of the eigenvalue distribution. Let us describe the solution to the one matrix model in some details. We will then discuss the generalisation to our setting of p-commuting matrices. The integral equation determining the eigenvalue distribution is given by: The potential (4.1) is even, which implies that the eigenvalue distribution should also be even. This allows us to rewrite the integral equation (4.3) as: which can be brought to a Cauchy form by the reparametrisation: z = a + 2bx 2 , y(z) = ρ 1 (x(z))/x(z); . We obtain: where c 1 and c 2 are given by: here r = 0 for connected distribution and R > r > 0 for disconnected distributions. Let us first consider the case of connected distribution, in this case the boundary of the eigenvalue distribution is at x = ±R and we seek solution to equation (4.5), which is bounded at c 2 and unbounded at c 1 . The unique such solution is given by: and for the eigenvalue distribution we obtain: The radius can be determined by normalising the distribution to one, we obtain: Note that the distribution (4.8) is well defined only for a certain range of the parameter a. Indeed, it is easy to show that the minimum of the distribution is achieved at x = 0 and then requiring that the distribution is positive at its minimum results in the restriction: At a = −2 √ b we have a "critical" distribution which vanishes at x = 0: a further reduction of a results in a phase transition to a disconnected distribution.
To find the form of the distribution we have to search for solutions of equation (4.5) that are bounded at both ends. In fact we can look for solution symmetric with respect to z = 0. Substituting c 1 = −z 0 and c 2 = z 0 , which implies z 0 = b(R 2 − r 2 ) and a = −b(R 2 + r 2 ) , for the unique such solution we obtain: Going back to variables x and ρ 1 for the eigenvalue distribution we obtain: requiring that ρ 1 is normalised to one and using the relation a = −b(R 2 + r 2 ) for the outer and inner radii we obtain: ; . (4.14) One can see that at the critical distribution, when a = −2 √ b, one has r = 0. This justifies the name "critical" since it belongs to both classes: the connected and the disconnected distributions which are more commonly referred to as the "one-cut" and "two-cut" solutions respectively. One can also see that for b < −2 √ b, which is the regime when the "one-cut" solution (4.8) is inconsistent, both radii of the "twocut" solution are well defined and the system is described by the "two-cut" solution (4.13). The system in fact goes through a 3-rd order phase transition at a = −2 √ b. To show this we have to analyse the behaviour of the specific heat of the model across the phase transition. The easiest way to calculate the heat capacity is to calculate the derivative of the internal energy with respect to the "temperature". To this end we calculate the expectation value of the potential (4.1) for the eigenvalue distributions (4.8) and (4.13) with rescaled couplings a → a/T, b → b/T . The next step is to calculate the derivative with respect to T and then take T = 1. We obtain: where C 1 v and C 2 v are the specific heats of the 'one-cut' and 'two-cut' solutions, respectively. One can easily see that at a = −2 v at this point. This confirms that the phase transition is of a third order.
In figure 6 we have compared the large N analytic expressions (4.8) and (4.13) to Monte Carlo simulations of the model for N = 800 and for definiteness we have set b = 1/2. The figure shows the excellent agreement with the theoretical large N results. Furthermore, one can see that at a = −2 √ b the critical embeddings is realised, which confirms the phase transition is continuous. In figure 7 we have compared the analytic expressions for the specific heat (4.15) to Monte Carlo simulations for N = 100 and N = 400.

Commuting matrix model with quartic potential in two dimensions
There are many possible extensions of the one matrix model to rotationally invariant two matrix models. The most obvious extension would be to consider (4.1) where x are two random matrices, which do not commute. To our knowledge this model has not been solved. An alternative approach is to build a non-hermitian matrix from Φ = X + iY and consider a non-Hermitian model with Hermitian Hamiltonian and quartic potential built from Φ † Φ. Such a system was solved by [24] (see also ref. [25]) using their method of Hermitization. Because the matrices don't commute and Φ † Φ = X 2 + Y 2 + i[X, Y ] this model is significantly different from those we consider but should reduce to our model if the contribution from the commutator is forced to zero. To our knowledge the commuting matrix models described below are new.
Here we perform an analogous investigation with emphasis on the relation between the two-matrix model and the reduced one matrix model. Our starting point is the integral equation: (4.16) Applying the Laplacian on both sides of the equation and using the two dimensional identity , one arrives at: where D is the domain of the distribution. Rotational invariance requires that the domain is either a disk or an annulus or a more exotic configuration of numerous concentric disks. The intuition that we gained from the one dimensional model suggest that for a quartic potential only the disk and the annulus are realised. Indeed, stability of the model requires that b ≥ 0, where b = 0 is allowed only for positive a 3 . For a > 0 the distribution (4.17) is positive and well defined for all x inside a disk of radius R. Normalising the distribution to one, we obtain: The eigenvalue distribution of the disk phase is then: If a = 0 we have critical distribution, which goes to zero at x = 0. For a < 0 the expression in equation (4.17) is negative at x = 0 and vanishes for some x > 0. It is therefore unphysical in this region and we need to modify our expression for the distribution. However, the functional form of the distribution (4.17) is independent on the shape of the domain D, this is a special property of the logarithmic kernel in two dimensions. Because of this property we are free to modify only the range of the distribution. A natural choice is to keep the same outer radius R and choose the inner radius r in such a way that the integral | x|<r d 2 x ρ 2 (x) vanishes. This results in: and we can write the eigenvalue distribution of the annulus phase as: Let us now calculate the reduced distribution: For a > 0 we reduce the disk distribution (4.19) to obtain: which as expected looks like equation (4.8) for the connected distribution in one dimension. In fact, if we reduce the potential according equation (A.3) we obtain: It is easy to convince oneself that to derive the connected one dimensional distribution for the reduced potential (4.24), one has to take a → 2a and b → 8b/3 in equation (4.8). In doing so one arrives at equation (4.23), confirming that indeed the disk phase of the commuting two-matrix model maps to the connected phase of the onematrix model, which is what we expect.
Let us now reduce the annulus phase. Naively we might expect this phase to map to the disconnected phase of the one-matrix model. However, this is not the case. For the reduced distribution of (4.21) we obtain: (4.25) which is profoundly different from the disconnected distribution (4.13). This is an important observation. In all previous examples the "shadow" of the higher dimensional model (namely the reduced distribution) corresponded to the physical distribution of the lower dimensional problem (with the reduced potential). Now we see that this does not hold uniformly. In particular for phases with non-trivial topology, the shadow of the higher dimensional problem does not reduce to the physics of the lower dimensional one. One should not be surprised by this result. Indeed, although the annulus phase has a non-trivial topology, it still corresponds to a connected distribution, this is clearly not the case for the disconnected phase of the one-dimensional model which has two disconnected components and is thus quite different.
Physically, this can be understood, because the quartic potential in one dimension is a double well and thus drives the theory into two disconnected phases associated to the different vacua, i.e. the moduli space of vacua is two points. The rotationally invariant quartic potential in two dimensions corresponds to a Mexican hat and the associated moduli space of vacua is the circle. So all vacua are connected and hence one expects the distribution to remain connected. With this revised intuition we can correct our naïve expectation to anticipate that the topology of the space of eigenvalues undergoes a transition from a disc to an annulus, in accord with the observation above.
These differences between the one-dimensional and two-dimensional models are not manifest when the theories are in the trivial vacuum (at the origin) and both distributions are topologically an interval and a disk, respectively. This is the reason we can map the dynamics of the disk phase of the two-matrix model to the dynamics of the connected phase of the one-matrix model.
One may wonder what happens in the interval a ∈ [− 8b 3 , 0], which still corresponds the one-cut distribution of the one-dimensional model, and whether there is anything special happening at a = − 8b 3 i.e. to the parameter value of the onedimensional transition. It turns out that in the two dimensional model there is no further non-analyticity at this value. The two dimensional transition is shifted to a = 0 and it is at this value that a hole appears in the eigenvalue distribution. What is special about a = − 8b 3 is that the inner radius occurs at the maximum of the reduced distribution, but we find no further non-analyticity at this parameter value.
To emphasise the different physics described by the one-and two-matrix models let us calculate the specific heat and explore its behaviour across the disk-annulus phase transition. Following the same path as in the analysis of the one matrix model we arrive at the following result for the heat capacity: ; for a ≥ 0 ; for a ≤ −0 , (4.26) One can see that . This shows that at a = 0 the heat capacity and its first three derivatives are discontinuous at a = 0 and it is the fourth derivative of the heat capacity which has a finite jump. Since the heat capacity is already a second derivative of the free energy this suggests that the phase transition is of a sixth order making it extremely smooth crossover. The heat capacity has another intriguing property, it is exactly 1/4 at a = 0, just like in the one-matrix model case. This is due to the constraint (B.6). Furthermore, it is odd with respect to the point (0, 1/4) (see figure 10).  eigenvalue distributions for the disk and annulus phases as well as for the critical distribution. One can observe an excellent agreement between the large N theoretical predictions and the numerical results. Figure 9 represents the spread of the eigenvalues for the disk and annulus phases. The middle plots represents a critical disk for a = 0. Finally, in figure 10 we have compared the plot of the large N result for the heat capacity (4.26) with the results for the heat capacity from numerical simulations for N = 100 (blue diamonds) and N = 400 (red diamonds) one can see the excellent agreement with the theoretical large N results.

Quartic potential in three dimensions
In this subsection we investigate the properties of a three-matrix commuting model with quartic potential. Our starting point is the integral equation: (4.27) Applying the operator | x| ∇ 2 x on both sides of the equation and using that ρ 3 is spherically symmetric to perform the angular integrals we obtain: where r = 0 for a "one-cut" solution with the topology of a ball and r > 0 for a "twocut" solution with the topology of an annulus. The easiest way to solve equation (4.28) is to reduce it to an integral equation with a Cauchy kernel. To this end we differentiate with respect x and change variables to z = 3a/π + 30b x 2 /π and y(z) = 2x(z) ρ 3 (x(z)) we obtain: where c 1 = 3a/π + 30b r 2 /π and c 2 = 3a/π + 30b R 2 /π. Our intuition from the gaussian case suggest that we look for a solution to (4.29) which is unbounded at z = c 2 (corresponding to x = R) and bounded at z = c 1 (x = r). The unique such solution is given by: Going back to x and ρ 3 we obtain: The corresponding reduced one dimensional distribution is given by: where Θ(x) is the step function. The explicit form of the distribution for r = 0 can be obtained in terms of elliptic integrals, we will use this solution to compare to numerical simulations. Equation (4.31) is our candidate for the "two-cut" solution. To get the "one-cut" solution we simply take the limit r → 0 in (4.31) obtaining: (4.33) Using equation (4.32) with r = 0 for the corresponding reduced one dimensional distribution we obtain: To obtain the radius of the one-cut solution we normalise it one. For the radius we find: One can check that with this radius the one-cut distribution also satisfies the constraint (B.6).
Obtaining the inner and outer radii of the two-cut solution is more subtle. The normalisation condition for ρ 3 can be used to find the outer radius as a function of a, b and the inter radius r. We obtain: To specify completely R and r we need to use the constraint (B.6), for the two-cut distributions it is given by: (R 2 − r 2 ) 4a 2 (r 2 + 3R 2 ) + 4ab(11(r 2 + R 2 ) 2 + 4R 4 )+ +5b 2 (9(R 2 + r 2 ) 3 − 14r 4 R 2 + 6R 6 = 1 (4.37) Solving equations (4.36) and (4.37) for R, r results in complex algebraic expressions, which we do not write explicitly, but we will keep in mind that in principle R and r are known as functions of a and b.
Note that the one-cut distribution (4.33) achieves its minimum at x = 0 and hence is well defined when ρ 3 (0) ≥ 0, which implies: At a = √ 20b/3 we have a critical solution and for a < √ 20b/3 we expect a phase transition from a ball phase (the one-cut solution) to an annulus phase (the two-cut solution). Let us analyse the heat capacity of the model. To calculate the heat capacity of the model we need the internal energy of the system as a function of both a and b. The internal energy is given by the expectation value of the potential with respect to the eigenvalue distribution ρ 3 . Next we multiply the internal energy by T rescale a → a/T, b → b/T and find the derivative with respect to T setting T → 1 afterwords. Note that this procedure requires knowing the derivatives of R and r with respect to a and b. While we didn't provide an explicit solution for the radii, the derivatives can be easily obtained indirectly by differentiating equations (4.36) and (4.37). Our final expression for the heat capacity is: Next using equations (4.36) and (4.37) we obtain the following expansion for C 1 v and C 2 v near a = √ 20b / 3 : Therefore we conclude that the heat capacity and its first derivative are continuous at a = √ 20b/3, while the second derivative has a finite jump. Therefore, the phase transition is of a fourth order.
In figure 11 we present a plot of the heat capacity as a function of a for b = 1. An interesting property that stands out is that C v = 1/4 at a = 0 (just like in the one-and two-matrix models), which is a consequence of the constraint (B.6). Furthermore, the heat capacity appears odd with respect to the point (0, 1/4). In fact by expanding C 2 v near a = 0 one can show that it is indeed odd with respect to the point (0, 1/4). Remarkably this symmetry persist as an approximate symmetry even across the phase transition for a > √ 20b/3. There is also a striking similarity with the heat capacity of the two-matrix model (look at figure 10). The diamonds in the figure 11 represent the results of Monte Carlo numerical simulations. One can see the good agreement between numerical results and the large N predictions (4.39).
Let us also compare our results for the eigenvalue distribution with numerical simulations. In figure 12 we present plots of the one-eigenvalue distributions for the ball phase, the annulus phase and for the critical distribution (at a = √ 20b/3). While one can see very good agreement in the ball phase (for the one-cut solution), one can see that for the annulus phase (the two-cut solution) the agreement is good only away from the inner radius of the distribution. Near the inner radius numerical simulations imply a sharp fall off, and a probable jump, in the distribution (similar to the one in the two-matrix model), while the analytic expression (4.31) falls gradually. This discrepancy is enhanced as N is increased. At present we don't have a theoretical way of describing such a sharp fall, since the bounded solutions of the Cauchy kernel integral equation (4.29) necessarily vanish at the boundary. A possible way would be to attack numerically the integral equation (4.28), but such studies are beyond the scope of this paper. Furthermore, the very good agreement of the heat capacity of the annulus phase obtained using the two-cut solution implies that it is very close to the real saddle point.

Quartic potential in four dimensions
In four dimensions out starting point is the integral equation: at the following result for the solution for general a and b: where R 2 = a − √ a 2 + 6b −6b .
As one can see the eigenvalue distribution is a mixture of an uniform distribution with density proportional to −b and a spherical shell distribution. One can also see that the distribution is physical only for b < 0 and a 2 > |6b|. However, for b < 0 the potential is unstable. Therefore this solution can be realised, for large N , only as a metastable phase trapped near the local minimum of the potential at x = 0. The absence of tunnelling stabilises this phase in the large N limit. Since increasing a broadens the well of the potential, while lowering the radius of the distribution, for sufficiently large a the eigenvalues spill out of this local well. This transition occurs at the upper bound at a 2 = −6b which represents the critical value and corresponds to the quantum gravity transition of the one dimensional model [22]. We will not investigate this transition further in this paper. For b > 0 the model is stable for any value of a, but the solution (4.44) is unphysical, therefore we expect that the shell saddle (2.7) is realised.
We conclude that if a is sufficiently large one should encounter a phase transition at b = 0 from the spherical shell phase to a mixed phase comprising of a spherical shell distribution and an uniform distribution inside the shell.
In figure 13 we have presented our results of Monte Carlo simulations. The first plot from left to right represents the spherical shell phase for b = 1, a = 1 and N = 400. The vertical dashed line represents the radius of the shell determined by equation (2.7). The second plot represents the mixed phase for b = −1, a = 3/2. The vertical dashed line represents the radius of the shell, while the horizontal dashed line represents the density of the uniform distribution both determined by equation (4.44). One can observe a very good agreement of the numerical results with the large N predictions.
It would be interesting to explore deeper the onset of instability in the mixed phase as a approaches |6b|. It would be also interesting to study the heat capacity of the system and determine the order of the phase transition. We leave these interesting studies for future work.

Conclusions
We have performed a systematic study of commutative SO(p) invariant matrix models with quadratic and quartic potentials. We found that the physics of these systems depends crucially on the number of matrices with a critical rôle played by p = 4. For p ≤ 4 and a quartic potential the system undergoes a phase transition, while for p > 4 the system is always in the low temperature phase. In terms of the joint eigenvalue distribution of the matrices, for p = 2 the transition is from a disc distribution to an annular one at the critical value a c = 0. This is precisely where one would expect the transition in the absence of fluctuations.
The physics here is straightforward: for a > 0 the potential V ( x) has a unique ground state and the resulting eigenvalue distribution is a disc. The precise distribution is given by (4.19), i.e. ρ 2 (x) = 4b and becomes the uniform distribution for b = 0 with R 2 = 1 a . For a < 0 the moduli space of ground states of V ( x) is the circle of radius r = −a 2b . The eigenvalue distribution then spreads into an annulus around this circle. The surprise is that eigenvalue repulsion is sufficiently strong that the annular phase emerges even at a = 0 corresponding to the pure quartic potential. Furthermore in contrast to the one-dimensional model the transition in which the eigenvalue distribution changes from a disk to a shell is in fact sixth order.
For p = 3 the physics is very similar to that for p = 2: One has an eigenvalue ball for large positive a. For negative a the moduli space of vacua of the potential is now a sphere and the eigenvalues spread about this sphere to give a spherical shell distribution. The transition between the two occurs at the surprisingly positive critical value a c = √ 20b/3, so that even a small quartic potential is not sufficient to guarantee some eigenvalues near the origin. Also the transition in this case turns out to be fourth order.
Surprisingly for p = 4 there is no longer a standard ball to annulus phase transition. For positive b, when the quartic potential is stable, the spherical shell phase (given by (3.8)) is the only possible phase, since the effective action's saddle is unphysical. However, for negative 4 b and sufficiently large a (for a 2 > |6b|) there is a mixed metastable phase comprising of a spherical shell with an uniform distribution inside the shell. As a result for a 2 < |6b|, there is a phase transition at b = 0 from the spherical shell phase to the mixed phase. Since the mixed phase contains an uniform ball, this transition can also be viewed as a topology changing phase transition for the eigenvalue distribution.
For p > 4 and quartic potential there are no transitions and one is always in the "broken"-symmetry phase. In fact for p > 4 the joint eigenvalue distribution is the infinitely thin spherical shell given by (3.8) where Ω p−1 is the volume of the unit p − 1-sphere.
There are several generalisations of this work that can be undertaken. One is to consider supersymmetric systems, this should be quite straightforward. A second is to consider the matrix quantum mechanics of commuting matrix models. Further generalisations are to consider non-rotationally invariant systems and more general potentials. We hope to return to these topics in the near future.

A Reducing from two dimensions
In this section we derive equation (2.28) relating the effective action in two and one dimensions. Let us write the effective action in two dimensions: We start with the first term in equation (A.1). Using equation (2.14) and integrating by parts we obtain: where we defined the reduced potential V 1 : Using again equation (2.14) it is easy to show that: which takes care for the last term in equation (A.1). Finally, we focus on the term containing the logarithmic kernel. Defining: and using equation (2.14) and the kernel (2.19) for p = 2 we can write: were S 1 is given by: with µ 1 = µ 2

B General constraint
In this section we derive a general constraint for the model (2.4), which in the special case of gaussian potential and in the N → ∞ limit reduces to equation (3.12). Our starting point is the mathematical identity: where S eff is the action (2.5). The identity can by proven by integrating by parts and using that the integrant vanishes at | λ i | → ∞. Performing the differentiation in (B.1) we obtain: which after dividing by the partition function (2.4) can be written as: Using equation (2.5) it is easy to check that: C Analytic expression for the free energy In this section we obtain an analytic expression for the non-constant part of the free energy E defined in equation (3.16). The idea is to uplift the calculation to α + 3 dimensions, where the pure stateρ α lifts to a spherical shell distribution. It is also convenient to rescale the distributionρ α to the range (−1, 1), To this end we define: ρ α (η) = R αρα (ηR α ) = Γ( α+3 2 ) π 1/2 Γ( α+2 2 ) (1 − η 2 ) α/2 , (C.1) where R α = (α + 2)/p as given in equation (3.14). Next we write E in terms ofρ α : and uplift the calculation of the second term in (C.2) to α + 3 dimensions. The distributionρ α lifts to: where Ω α+2 is the volume of the unit α + 2 dimensional sphere. The crucial step is to use equation (2.22) and the same considerations that lead to equation (2.25) to write: