Symmetry, Integrability and Geometry: Methods and Applications The Multitrace Matrix Model of Scalar Field Theory on Fuzzy CP n⋆

We perform a high-temperature expansion of scalar quantum field theory on fuzzy CP^n to third order in the inverse temperature. Using group theoretical methods, we rewrite the result as a multitrace matrix model. The partition function of this matrix model is evaluated via the saddle point method and the phase diagram is analyzed for various n. Our results confirm the findings of a previous numerical study of this phase diagram for CP^1.


Introduction
Fuzzy spaces are noncommutative geometries which arise from quantizing certain compact Kähler manifolds. The most prominent such space is the fuzzy sphere, which was first constructed by Berezin [1]. In the original construction, the aim was the same as that of geometric quantization, i.e. to provide a general quantization prescription for a particle whose phase space is an arbitrary Poisson manifold. Today, fuzzy spaces attract most interest for different reasons: First, fuzzy spaces appear quite naturally in various contexts in string theory where they replace parts of the classical geometry of the target space with an approximate quantum geometry. Closely related is the observation that fuzzy geometries seem to emerge from the dynamics of matrix models and thus they could be crucial in background independent formulations of theories of gravity. And finally one can regulate quantum field theories on Kähler manifolds by putting the theory on the corresponding Berezin-quantized or fuzzy manifold.
The idea of using fuzzy spaces as regulators for quantum field theories goes back to the early 1990's [2,3]. This approach is very appealing, as the definition of scalar quantum field theories on fuzzy spaces is under complete control: All functional integrals are automatically well-defined because the algebra of functions on a fuzzy space is finite dimensional. Taking the large volume limit of the fuzzy space, we can even regulate scalar quantum field theories on flat spaces and thus try to compete with the lattice approach. The main advantage of fuzzy regularization over the latter is that all the isometries of the original Kähler manifold survive the quantization procedure.
Particularly nice spaces to use in a fuzzy regularization are the complex projective spaces, as they are the Berezin-quantizable manifolds with the largest possible symmetry groups. Fur-thermore, their quantization is straightforward and can be done completely in terms of group theory. As usual in a "good" quantization, real functions are mapped to hermitian operators on a Hilbert space, which is finite dimensional in Berezin quantization. Real scalar field theories are therefore simply hermitian matrix models.
The most prominent hermitian matrix models are given by a potential consisting of a trace over a polynomial in the matrix variable. One can therefore switch directly to an eigenvalue formulation. In the case of scalar field theories on fuzzy spaces, this is not possible because the kinetic term yields a coupling to a number of fixed "external" matrices.
A first attempt at gaining an analytical handle on fuzzy scalar field theories was made in [4]. A new method to overcome the problem of external matrices was then proposed in [5]. Here, a high-temperature expansion of the kinetic term in the partition function was performed and the resulting expressions could be evaluated analytically via group theoretic methods. It was shown that the resulting partition function can be rewritten as the partition function of a multitrace matrix model. This partition function can then be computed analytically for both finite and infinite matrix sizes using, e.g., orthogonal polynomials or the saddle point approximation. For P 1 , this computation was performed to second order in the inverse temperature β in [5]. In this paper, we continue this work and generalize the results to third order in β and to arbitrary P n .
One of the motivations for this work is to explain the phase diagram for scalar field theory on fuzzy P 1 which has been obtained via numerical methods in [6], see also [7] for a more detailed study as well as [8] for a review and further numerical results. The numerical results suggest that the phase diagram is invariant under a particular multiscaling. We can therefore restrict ourselves to the limit of infinite matrix size, in which we can use the saddle point approximation to compute the partition function of our model.
Further reasons to compute the multitrace matrix model to third order in β are the possibility to use this result in a similar study of scalar field theory on Ê × P 1 as well as our intent to discuss the link to (deformed) integrable hierarchies in future work.
In the analysis of the phase diagram, we will focus our attention on the three lowestdimensional fuzzy spaces P 1 F , P 2 F and P 3 F . In the first case, the goal will be to compare the resulting phase diagram with the numerically obtained one. The quantum field theory on the second space corresponds in the large volume limit to a scalar quantum field theory of φ 4 -type on Ê 4 . While admittedly it is not clear what the Lagrangian of the field theory on Ê 4 being regularized actually is, this presents an example of both a well-defined and renormalizable fourdimensional noncommutatively deformed φ 4 -theory. The theory on P 3 could be interpreted as a regularization of a non-renormalizable field theory, and one might hope for signs of this in the matrix model.
The paper is structured as follows. In Section 2, we review the construction of fuzzy P n and scalar field theory on this noncommutative space. Section 3 describes the high-temperature expansion in detail and the results are presented to order β 3 . In Section 4, we analyze the thus obtained multitrace matrix model for the three lowest-dimensional fuzzy P n and we conclude in Section 5. Conventions, rather technical details and helpful intermediate results are given in the appendix.
2 Scalar f ield theory on fuzzy P n The general mathematical framework containing the quantization of complex projective space which is referred to as fuzzy P n in the physics literature is known as Berezin-Toeplitz quantization, see e.g. [9] and references therein for a detailed discussion. In the case of P n , there is a shortcut to the general constructions of Berezin-Toeplitz quantization which originates from the fact that P n is the coset space U(n + 1)/U(1) × U(n). We will use this group theoretic approach here, as it has the additional advantage of allowing for simple computations of quantities like spectra of quadratic Casimirs and their eigenspaces, which we will need for our further discussion.

Berezin quantization of P n
The Hilbert space H ℓ which we use in Berezin quantizing P n is the space of global holomorphic sections of the line bundle O(ℓ) over P n with ℓ ≥ 0. As a vector space, H ℓ is spanned by the homogeneous polynomials of degree ℓ in the homogeneous coordinates z 0 , . . . , z n on P n . Recall that P n ∼ = SU(n + 1)/S(U(1) × U(n)), and H ℓ forms a representation of SU(n + 1) which is given by the totally symmetrized tensor product of ℓ fundamental representations. In terms of Dynkin labels, this representation reads as (ℓ, 0, . . . , 0) and has dimension N n,ℓ := dim(H ℓ ) = dim(ℓ, 0, . . . , 0) = (n + ℓ)! n!ℓ! .
The generators of su(n + 1) are represented on End (H ℓ ) by the adjoint action of hermitian matrices L i , and we introduce the quadratic Casimir operator according to The eigenvalues of C 2 are positive and given on the irreducible subspace with Dynkin labels (m, 0, . . . , 0, m) by 1 2m(m + n). The degeneracy of each of these eigenspaces is given by Because of C 2 (σ −1 ℓ (f )) = σ −1 ℓ (∆f ), where f ∈ Σ ℓ and ∆ is the Laplace operator on P n , it is justified to identify C 2 with the Laplace operator on fuzzy P n .
The matrices L i represent the generators of su(n + 1) and thus satisfy the algebra where the f ijk are the structure constants of su(n + 1). We choose the L i such that In the adjoint representation R = (1, 0, . . . , 0, 1), they satisfy the Fierz identity from which we conclude that With the above relation, one readily verifies the following identity for the structure constants: Using the overcompleteness relation for the Perelomov coherent states, dµ |z, ℓ z, ℓ| = vol( P n )½, where dµ = ω n n! is the Liouville measure obtained from the Kähler form ω yielding the Fubini-Study metric, one readily deduces a formula for integration: Given a function f ∈ Σ ℓ , the integral can be written as a trace over the quantized function σ −1 ℓ (f ) ∈ End (H ℓ ): 2.2 Quantum scalar field theory on P n F As we are interested in matrix models, it is convenient to switch from the label ℓ of our representations to the label N n,ℓ and drop the subscript. One should, however, keep in mind that only for P 1 , there is an ℓ for every value of N . In the following, we will represent elements of End (H ℓ ) by hermitian matrices Φ of dimension N × N .
In the previous section, we collected all the necessary results for writing down a scalar field theory on fuzzy P n . Putting everything together, we arrive at the following action functional 2 on End (H ℓ ): As we work with hermitian generators L i , the quadratic Casimir operator C 2 has positive eigenvalues and for r ∈ Ê and g > 0, the action is therefore bounded from below. This, together with the finite dimensionality of End (H ℓ ), enables us to introduce the well-defined functional integral where dµ D (Φ) is the Dyson measure on the set of hermitian matrices of dimension N × N . Recall that we can diagonalize a hermitian matrix Φ according to Φ = ΩΛΩ † , where Ω ∈ U(N ) and Λ = diag(λ 1 , . . . , λ N ) is the diagonal matrix of eigenvalues of Φ. Under this decomposition, the Dyson measure splits into an eigenvalue part and an "angular" integration over U(N ): where dµ H (Ω) is the Haar measure 3 and ∆(Λ) is the Vandermonde determinant In the case of simple hermitian matrix models consisting of traces (and multitraces) over polynomials in Φ, the angular integration is trivial, because tr(Φ n ) = tr(Λ n ), and reduces to a constant volume factor. The remaining integral over the eigenvalues can then be computed by standard methods as e.g. the saddle point approximation or orthogonal polynomials. Here, however, the kinetic term contains the fixed external matrices L i which obstruct a straightforward translation to the eigenvalue picture.

2.3
The toy models N = n + 1 on P n In the case N = n+1, i.e. when ℓ = 1 and End (H 1 ) forms the adjoint representation of su(n+1), the kinetic term of our model (3) can be evaluated explicitly by using the Fierz identity (1). We find here that where tr(K) stands for the sum over the eigenvalues of C 2 on End (H 1 ). Note that, as necessary, the kinetic term vanishes for Φ ∼ ½. We will use this class of toy models for consistency checks of our computations below.

The high-temperature expansion
As it does not seem possible to compute the partition function (4) analytically, we perform a high-temperature expansion as suggested in [5]. That is, we separate out the kinetic term in the functional integral and Taylor-expand its exponential, assuming β to be small. As β is usually inversely proportional to the temperature in statistical mechanics models, this expansion is also known as a a high-temperature expansion in the literature. For each of the terms appearing in this expansion, the integral over the angular part of the Dyson measure can be performed -in principle straightforwardly -using group theoretic methods. The results can be rewritten in terms of multitrace terms, and, after putting them back into the exponential of the functional integral, one ends up with a multitrace matrix model.

Setup of the expansion
Let us consider our model (3) on fuzzy P n F with the dimension of the quantum Hilbert space H ℓ being N . The space of quantized functions End (H ℓ ) is spanned by the generators 4 τ µ , µ = 1, . . . , N 2 of u(N ). We start by rewriting the kinetic term of the action in the following way: The expansion of the kinetic term in the action now reads as and we will restrict our attention in the following to the terms up to order O(β 3 ). We want to perform the integral over the U(N ) part of the Dyson measure, i.e. to integrate out the angular degrees of freedom in Φ. For this, we decompose the hermitian matrix Φ according to Φ = ΩΛΩ † , where Ω ∈ U(N ) and Λ = diag(λ 1 , . . . , λ N ). The integrals we have to evaluate at order O(β k ) are thus of the form where the essential part in index notation is given by Various algorithms have been proposed in the literature to compute integrals of the type (7), cf. e.g. [10,11,12]. The most involved integral of the form (7) which we are interested in is the one for k = 3, which is already very difficult to handle by the suggested methods. Fortunately, the integrals (6) allow for a further simplification [5], which is then accessible via group theoretic methods. Using tr(A) tr(B) = tr(A ⊗ B) and AB ⊗ CD = (A ⊗ C)(B ⊗ D), we rewrite (6) according to The idea presented in [5] is now to use the orthogonality relation of the Haar measure (12) to evaluate these integrals. We thus have where the sum is taken over the irreducible representations contained in the tensor product of 2k fundamental representations of SU(N ). The traces tr ρ are taken in the representation ρ, and we have where χ ρ (Λ) denotes the character of Λ in the representation ρ. Characters of representations of SU(N ) can easily be calculated using e.g. the formulas in [13]. The remaining challenge is therefore to evaluate tr ρ (τ m 1 ⊗ τ n 1 ⊗ · · · ⊗ τ m k ⊗ τ n k ).

The restricted traces tr ρ (·)
Consider again the generators τ m of su(N ), and denote their matrix components by τ αβ m , α, β = 1, . . . , N . The full trace over the tensor products of matrices in index notation is given by To evaluate the restricted traces tr ρ (·), we need to project onto the irreducible representations which we do using projectors P (i,j) 2k constructed from Young symmetrizers. The technical details of the construction of these projectors are given in Appendix B. Explicitly, we let the projector P (i,j) 2k act onto the indices β appearing in the Kronecker deltas to restrict to a representation ρ (i,j) : The completeness relation (14) for the projectors P 2k translates into the following completeness relation for the restricted traces: which can serve as a first consistency check of the correctness of the calculated projectors P A second test is to verify that each individual restricted trace indeed reduces to the character if all the τ m are equal: Let us now compute the combined sums of the restricted traces for each type of Young tableaux and contract the τ m with the K mn to simplify the results. That is, we compute the following expressions: For k = 1 and k = 2, corresponding to the contributions at orders O(β) and O(β 2 ), these sums have already been calculated in [5]. They are : : The lengthy result for k = 3 is given in Appendix C.
In calculating these results, we used many identities which we will briefly comment on now: First of all, using the Fierz identity for the generators of u(N ) as well as the relations for the L i , we compute that for arbitrary A, B ∈ u(N ), Applying these relations to tr(K) := K µν tr(τ µ τ ν ) yields The identities (8) allow us to successively rewrite expressions involving K µν in terms of traces over products of the L i , which in turn can be reduced using L 2 i = c L ½ N and the identity for the structure constants (2). Some useful intermediate results are collected in Appendix C.

The multitrace matrix model
Combining the reduced traces with the characters in the various representations, we arrive at the following expressions for the I k : The result for I 3 is lengthy and because it can be easily calculated from the list of restricted traces given in Appendix C, we refrain from presenting it here. Note that all the above integrals pass the first consistency check: To rephrase the perturbative expansion in terms of an effective action, we re-exponentiate the terms. That is, we write where we demand that the S i are polynomials in the eigenvalues of order 2i. As the same holds by definition for the I i , we can match both sides order by order and arrive at We can now perform the second consistency check of our result and compare the re-exponentiated action with the toy model N = n + 1 from Section 2.3. In the representation N = n + 1, we have Plugging this into the expressions for S 1 , S 2 and S 3 , we find that S 1 indeed reduces to the kinetic term of the toy model (5). Since the terms S 2 and S 3 vanish as required, our results pass this consistency check as well. By re-inserting the integration over the Haar measure, we can return from the eigenvalues to the full hermitian matrices Φ. We thus obtain the multitrace matrix model with action S = S 1 + S 2 + S 3 with etc. The involved expressions for S 2 and S 3 are again lengthy but easily calculated from the results given above. Altogether, we obtained a multitrace matrix model whose partition function approximates the partition function of fuzzy scalar field theory on complex projective space up to order O(β 3 ). This approximation should be valid in particular for large values of the couplings r and g.

Large N solutions of the model
The partition function of the multitrace matrix model we obtained in the previous section can now be evaluated analytically for finite N using the methods of orthogonal polynomials. As we are mainly interested in the phase diagram, we consider instead the large N limit and use the saddle point method to determine the partition function here. We try to be self-contained and present the involved steps in detail.

The large N limit
The phase diagram determined numerically in [6] is invariant under the multiscaling limit where N → ∞ and N 2 βg as well as N 3/2 βr are kept fixed. This justifies to solve our model in the large N limit to compare it with the phase diagram. In this limit, the discrete set of eigenvalues goes over into a continuous function: We rescale λ i → λ(i/N ) =: λ(x), with 0 < x ≤ 1. The traces turn correspondingly into integrals: tr The formulas for general n turn out to be very lengthy and difficult to handle. Therefore we will restrict our attention in the following to the three projective spaces P 1 , P 2 and P 3 . The first case n = 1 is interesting, as we would like to compare the resulting phase diagram to the one numerically obtained in [6]. The second case n = 2 is a well-defined, four-dimensional quantum field theory with quartic potential. The third case n = 3 is interesting as the corresponding scalar field theory on Ê 6 is not renormalizable.
The eigenvalues of the quadratic Casimir, the degeneracy of the corresponding eigenspaces and the dimension of the representation (ℓ, 0, . . . , 0) ⊗ (ℓ, 0, . . . , 0) for the cases n = 1, 2, 3 are listed in the following table: 2ℓ(ℓ + 2) 2ℓ(ℓ + 3) degeneracy of eigenspaces 1 + 2ℓ (1 + ℓ) 3 As the function N n,ℓ is only surjective for n = 1, and thus we cannot find an ℓ for every value of N , we will rewrite the multitrace matrix model in terms of ℓ. From the table above, we easily evaluate the various traces over K appearing in the action of the multitrace matrix model. We have for P 1 : for P 2 : and for P 3 : In the limit ℓ → ∞, the expressions for the various matrix models simplify. Switching to the eigenvalue description and using the moments we can write them down explicitly. On the three fuzzy P n s, we have the models where the repulsive log-term arises as usual from exponentiating the Vandermonde determinant.
The ℓ-dependence of the log-term is due to the factor N 2 n,ℓ , which in turn originated from rewriting the double sum as a double integral. In the above expressions, subleading terms in ℓ have been suppressed in each summand.
As a next step, we have to find the appropriate multi-scaling behavior of the constants β, r, g and the continuous eigenvalues λ(x). The coefficients of the log-terms determines the desired scaling behavior of the total action. We fix the remaining scalings by demanding that the whole action scales homogeneously and that βg scales with N 2 , as for an ordinary hermitian matrix model. We thus find the following rescalings: Note that the scalings for P 1 indeed agree with the ones numerically determined in [6] as well as the ones calculated in [5]. As a final simplification, we note that our theory is invariant under Φ → −Φ, as the potential is even. We expect the eigenvalues to respect this symmetry 5 , and therefore we put all the odd moments c 2n+1 , n ∈ AE to zero. Moreover, we replace the integral over x by an integral over the eigenvalue density ρ(λ) := dx dλ . We thus eventually arrive at the following three models, which we wish to solve:

Solving the models
We will now calculate the partition functions of our models using the saddle point method: In the large ℓ limit, the path integral localizes on classical solutions, or saddle points, of the action 6 . These solutions, which are valid only for a restricted range of the coupling constants, can be easily obtained using standard methods in random matrix theory. We start from the action 7 where I is the union of open intervals on the real line over which ρ(λ) has support. The saddle point equation is obtained by varying the above equation with respect to ρ(λ): Note that our potentials satisfy V (λ) = 0 at λ = 0 and we can therefore determine the Lagrange multiplier ξ by solving the saddle point equation at this special point if 0 ∈ I: otherwise, one has to choose a different value of λ to obtain ξ. We define the free energy F as F := − log(Z), where Z is the partition function of our model. In the saddle point approximation, this reduces to F = βS[ρ(λ)], which we can evaluate using (9):

C. Sämann
To find the eigenvalue density ρ(λ), it is convenient to replace (9) with its derivative 8 with respect to λ: This is a singular integral equation, and its general solution can be found e.g. in [14], see also [15]. First of all, one introduces the resolvent W (λ), which is an analytic function on \I, defined according to Note that for large λ, we have W (λ) ∼ 1 λ . The resolvent is related to the eigenvalue density ρ(λ) and the Cauchy principal value appearing in the equation of motion through the Plemelj formula, and we arrive at The first equation determines ρ(λ) in terms of the resolvent, and the second equation is a much simpler equation than (10), which fixes the resolvent 9 and thus the eigenvalue density. One can show that the resolvent satisfies the Schwinger-Dyson equation The solution to the above equation reads as where ω(λ) describes the part of W (λ) containing the branch cuts. Explicit solutions are now obtained by making assumptions about the support I of the eigenvalue density ρ(λ). The simplest assumption is that I consists of a single interval. This is expected if the potential either consists of one deep well or if the eigenvalue filling is such that all the local minima of the potential are more than filled up. In this case, the resolvent has to have a branch cut over I := (δ 1 , δ 2 ) and the corresponding solution is therefore known as a single-cut solution. The resolvent's singular part has to contain exactly two roots, cf. e.g. [16]: One can now make a general ansatz for the polynomials M (λ) and R(λ). Together with the self-consistency condition that all the moments c n satisfy their defining relation c n := dλ ρ(λ)λ n , 8 When doing this, one obviously has to vary each moment: δc 2 2 = 2c2δc2, etc. 9 Strictly speaking, it fixes the resolvent only up to regular terms, which, however, are absent as can be seen from the large λ behavior W (λ) ∼ 1 λ .
we can solve for all unknowns and determine ρ(λ). Note that the normalization condition on the eigenvalue density c 0 = 1 is equivalent to the less involved condition that the asymptotic behavior of the resolvent is W (λ) = 1 λ + O( 1 λ 2 ). When having a double well potential, we also expect solutions where I is given by the union of two disjoint intervals I = (δ 1 , δ 2 )∪(ε 1 , ε 2 ). Correspondingly, the singular part of the potential contains four roots and we make the ansatz:

This solution is known as a double-cut solution.
It is important to stress that in general, all solutions will be valid only on a subset of the full parameter space of the model under consideration. This subset is characterized by the condition ρ(λ) ≥ 0 (and therefore c 2n ≥ 0) as well as the condition that I is of the assumed form. It is called the existence domain of a solution, and its boundary in parameter space can correspond to a phase transition.
In the following, we will present all the solutions for the various models together with their existence domains. For P 1 , we also give the explicit expressions for the free energies. We will consider three kinds of solutions: the symmetric single-cut solution, the symmetric double-cut solution and the asymmetric single-cut solution. In the latter case, we should strictly speaking include all the odd moments c 2n−1 , n ∈ AE, which we dropped in our actions. However, the full action would be very difficult to handle analytically, and we hope to make at least qualitative statements with our truncation.
The solutions for closely related models, in which r is kept fixed, the coefficient of c 2 2 is a parameter and the coefficient of c 3 2 vanishes, have been computed in [17] for the symmetric single-cut type and [18] for the two other types.

Solutions of the model on P 1
For the symmetric single-cut case, we assume that the eigenvalue density ρ(λ) has support on the interval I = (−d, +d). The solution we obtain from the procedure described above together with the conditions W (z) ∼ 1 z and I dλ λ 2 ρ(λ) = c 2 reads as The free energy of this solution is given by The solution exists, if both ρ(λ) and c 2 are nowhere negative. The condition ρ(λ) ≥ 0 amounts to while c 2 ≥ 0 is always satisfied in the existence domain of the solution (11). Next, we assume a symmetric double-cut support for ρ(λ) on The solution here reads as

C. Sämann
To evaluate the free energy, we compute ξ at λ = √ s and use the relation The result is The symmetric double-cut solution exists, if s = c 2 > 0, i.e. if r < −1 and 3g > β or r > −1 and 3g < β and if s > d. The latter condition yields and we see that the boundary of the existence domain of the symmetric double-cut solution matches that of the symmetric single-cut solution. We therefore expect a phase transition at this boundary. At the point (r 0 , g 0 ) = (−1, β 3 ), something interesting happens: Here, the equation for s becomes trivially satisfied and s is unconstrained. The two cuts can thus be arbitrarily far apart. At this point, the action reduces to and we have a competition of the single-trace and the multitrace potential term. The asymmetric single-cut solution with support on I = (s − d, s + d) with s = 0 is given by To evaluate the free energy, we determine the Lagrange multiplier ξ at λ = s. Our definitions then yield The asymmetric single-cut solution exits if ρ ≥ 0 and if d is real and positive. The first condition implies and c 2 and s are automatically positive. The second condition amounts to g > β 3 .

Phase structure on P 1
If existence domains do not overlap, we expect a phase transition at the boundary. If, however, two solutions exist for the same parameters, then the solution with the lowest free energy will be adopted. The resulting phase diagram is depicted in Fig. 1.
Here, the symmetric single-cut, the double-cut and the asymmetric single-cut solutions are labelled as I, II and III. The boundary of the existence domain between I and II describes the usual second order phase transition of the hermitian matrix model. The existence domain of III is fully contained in the existence domain of II. There is indeed a region of the parameter space, where the asymmetric filling yields a lower value for the free energy than the symmetric filling. This is particularly interesting, as it was very difficult to extract the II-III phase transition from the numerical data in [6]. We can thus confirm the numerical findings. Furthermore, there is the forbidden region in the parameter space for g < β 3 . We expect that higher-order corrections in β would deform this boundary.
Altogether, we have obtained the general features of the phase diagram found in [6]: We have three distinct phases, which come together at the point (r 0 , g 0 ) = (−1, 1/6), which corresponds to (b, c) = (−0.5, 1/12) in the conventions of [6]. This compares to numerically found values of (b, c) = (−0.8 ± 0.08, 0.15 ± 0.05). The discrepancy is due to the fact that the triple point is in a region of the parameter space, where the kinetic term is not small compared to the potential terms.
The effects due to the asymmetric single-cut region have to be considered only qualitatively, because we have dropped all the odd moments from the action using symmetry arguments. This was done to keep the solutions under analytical control. The critical line found in [6] corresponds here to the dashed curve. Including the odd momenta would presumably straightened this curve.
The discrepancies compared to [5] arise from the fact that there, a contribution to the action, labelled K K, was neglected in the large N limit while we included it here. This yielded a different model with the opposite sign of the c 2 2 term.

Solutions of the model on P 2
Let us now be brief in repeating the analysis for P 2 : The single-cut solution with support on I = (−d, d) is given by .
Note that the left and the right bound touch at a single point in the r-g-plane.
The asymmetric single-cut solution has support I = (s − d, s + d) and is given by This solution is valid for The moment c 2 and the center of the cut s are automatically positive. Contrary to the case of P 1 , there is no upper bound on r for this solution.
The expressions for the free energies on P 2 are not presented as they are lengthy but can be calculated quite straightforwardly.  We list now the three solutions for the model on P 3 . This model differs from that on P 2 only in the magnitude of the coefficients in the action. The supports are chosen in the same way as on P 1 and P 2 .
Asymmetric single-cut solution: The boundaries for the existence domains are presented in Fig. 3. The meaning of the lines is the same as for P 2 . Not surprisingly, the phase diagram is essentially identical to that for P 2 . Unfortunately, there is no feature hinting at the non-renormalizability of φ 4 -theory on Ê 6 .

Conclusions
In this paper, we computed the partition function of scalar quantum field theory on fuzzy P n to third order in the inverse temperature β, generalizing the results of [5]. As this theory can be interpreted as a noncommutative deformation of scalar quantum field theory on Ê 2n in the large N limit, we also demonstrated the existence of a nontrivial such theory.
We started by expanding the exponential of the kinetic term in the partition function. We then used group theoretic methods to integrate out the zero modes of the action and obtained a multitrace matrix model. This model was then solved via the saddle point approximation in the large N limit. In principle, however, the partition function of the matrix model could have been computed as well at finite N using orthogonal polynomials.
We presented the explicit classical solutions on which the partition function localizes in the large N limit and discussed the arising phase diagrams for P 1 , P 2 and P 3 . We confirmed the findings of the numerical analysis of [6,7] for P 1 and reproduced qualitatively -and partly quantitatively -the phase diagram found numerically. That is, we confirmed the existence of three distinct phases and confirmed analytically their properties suggested by the numerical studies. We also found a triple point which agrees to an acceptable degree with the one found numerically.
Here, it was particularly interesting that we found a large region of the parameter space in which an asymmetric single-cut solution was energetically favorable to a symmetric double-cut solution, even though our potential was symmetric. Such situations have been studied in the past, see e.g. [19]. Physically, the existence of this spontaneous symmetry breaking in the large N limit can be explained as follows. Consider the matrix model at finite N . Here, we have to introduce explicitly a symmetry breaking term for an asymmetric phase to exist, as tunnelling of the eigenvalues would otherwise restore symmetry in the eigenvalue filling. After taking N to infinity, there is no more tunnelling, and the symmetry breaking term can be safely switched off, preserving the asymmetric configuration.
It would be interesting to push the analysis of the phase diagrams further and, for example, to include all odd moments and examine the possibility of a smooth transition of the filling fractions between phases II and III. We then might be able to reproduce the slope for the linear phase boundary found in [6]. It might also be interesting to compare our results to the findings of [20], where the phase structure of noncommutative field theories on Moyal space was analyzed and [21], where questions arising from [20] were discussed on the fuzzy sphere. Moreover, we intend to use the results found here to study scalar field theory on Ê × S 2 as well as the relation of our multitrace matrix model to (deformed) integrable hierarchies in the future. Finally, recall that multitrace matrix models had been proposed as candidates for conformal field theories with c > 1 coupled to gravity [17]. One might be able to make sense of our models in this context, as well.
then label them by their type and their position in the ordering: The tableau λ (2,3) , for example, is given by the third tableau of the second type: λ (2,3) = 1 2 5 3 4 6 .
These expressions can be easily checked for n = 1 and N = 2, 3, 4 by comparing them to the results of a direct computation.