New Algorithm and Phase Diagram of Noncommutative Phi**4 on the Fuzzy Sphere

We propose a new algorithm for simulating noncommutative phi-four theory on the fuzzy sphere based on, i) coupling the scalar field to a U(1) gauge field, in such a way that in the commutative limit N\longrightarrow \infty, the two modes decouple and we are left with pure scalar phi-four on the sphere, and ii) diagonalizing the scalar field by means of a U(N) unitary matrix, and then integrating out the unitary group from the partition function. The number of degrees of freedom in the scalar sector reduces, therefore, from N^2 to the N eigenvalues of the scalar field, whereas the dynamics of the U(1) gauge field, is given by D=3 Yang-Mills matrix model with a Myers term. As an application, the phase diagram, including the triple point, of noncommutative phi-four theory on the fuzzy sphere, is reconstructed with small values of N up to N=10, and large numbers of statistics.


Introduction
The goal of this article is to reconstruct by means of a (hopefully) novel, and efficient Monte Carlo method the phase diagram of noncommutative phi-four on the fuzzy sphere. This was originally done in [1]. The basic theory is given by the following two-parameter matrix model In this equation L a are the SU (2) generators in the irreducible representation with spin s = (N − 1)/2, Tr H 1 = N , b is the mass parameter, and c is the coupling constant. The parameter a can always be chosen to be equal to 1. There are three known phases in this model. The usual Ising transition between disorder and uniform order. A matrix transition between disorder and a non-uniform ordered phase, and a (very hard to observe) transition between uniform order and non-uniform order. The three phases meet at a triple point [1,2]. The non-uniform phase, in which rotational invariance is spontaneously broken, is simply absent in the commutative theory. The non-uniform phase is the analogue of the stripe phase observed on the Moyal-Weyl spaces [14], whereas the disorder-to-non-uniformorder transition is the generalization of the one-cut-to-two-cut transition, observed in the Hermitian quartic matrix model [15,16], to the fuzzy sphere. This is a highly non-trivial problem, which is due mainly, to the more complicated phase structure of matrix scalar phi-four. It involves transitions between vacuum states, with very low probability distributions, and as a consequence, they are extremely difficult to sample correctly with the Metropolis algorithm. In particular the non-uniform-touniform transition is virtually unobservable in ordinary Metropolis, due to the absence of tunneling between the identity matrix, corresponding to the uniform phase, and the other idempotent matrices, corresponding to the non-uniform phase. This means simply that the Metropolis updating procedure does not sample correctly, and equally, i.e. according to the Boltzmann weight, the entire phase space which includes an infinite number of vacuum states. This was circumvented, in [1,2], by a complicated variant of the Metropolis algorithm, in which detailed balance is broken. This problem was also studied in [3,4,7,8]. The analytic derivation of the phase diagram of noncommutative phi-four on the fuzzy sphere was attempted in [17][18][19].
The related problem of Monte Carlo simulation of noncommutative phi-four on the fuzzy torus, and the fuzzy disc was considered in [5,6], and [9] respectively.
The main strategy employed, in this article, towards a better resolution of this problem, is to reduce the model down to its eigenvalues, without actually altering it. This is achieved by: 1) coupling the scalar field to a U (1) gauge field, in such a way, that in the commutative limit N −→ ∞, the two modes decouple completely, and thus we return to an ordinary phi-four theory, and 2) diagonalizing the scalar field by means of a U (N ) gauge transformation, viz Φ = U ΛU + , and then integrating out the unitary matrices U and U + from the path integral.
In this algorithm, we thus trade off the Monte Carlo simulation of the unitary matrices U and U + , in the original model (1), with the Monte Carlo simulation of a U (1) gauge field on the fuzzy sphere, which we know is very efficient using ordinary Metropolis [10].
The primary interest, of this article, is therefore Monte Carlo simulation of a noncommutative phi-four theory, coupled to a U (1) gauge field on the fuzzy sphere, using the Metropolis algorithm with exact detailed balance. The scalar field transforms in the adjoint representation of the U (1) gauge group, and as a consequence, the scalar and gauge degrees of freedom decouple in the commutative limit N −→ ∞. In other words, this theory becomes an ordinary phi-four theory in the commutative limit. In this theory, the usual scalar kinetic action ∼ −T r[L a , Φ] 2 is replaced with ∼ −T r[X a , Φ] 2 , where X a is itself obtained by Monte Carlo simulation of an appropriate gauge action, which will be centered around ∼ L a , in the so-called fuzzy sphere phase 1 . The pure gauge action is given by D = 3 Yang-Mills action, with a Chern-Simons (Myers) term. For b = c = 0 the full action is in fact D = 4 Yang-Mills action, with a Chern-Simons (Myers) term.
This article is organized as follows. In section 2, we present the detail of the U (1) gauge covariant noncommutative phi-four theory on the fuzzy sphere, and also explain the Metropolis algorithm employed in our Monte Carlo simulations. In section 3, we report our first numerical results, on the phase diagram of noncommutative phi-four on the fuzzy sphere, using our new algorithm. We give independent measurements, of the three transition lines, discussed above, and then derive our estimation of the triple point. These results are obtained with small values of N up to N = 10, and large numbers of statistics. In section 4, we give a construction of a one-parameter family of noncommutative phi-four models on the fuzzy sphere, which define, a regularization of duality covariant noncommutative phi-four on the Moyal-Weyl plane. We conclude in section 5, with a brief summary, and outlook.

The Action
Instead of the basic model (1), which is the primary interest in this article, we consider a four matrix model given by the action The fuzzy sphere phase is given by the background The values m 2 = 2c 2 M and µ = −9β/α 2 , of interest, are (with c 2 = (N 2 − 1)/4 being the Casimir operator) In the remainder we will be interested in the first case.
The first scaled parameter is [10] In the notation of [1] 2 , after replacing with X a = αϕL a in S m , we have a =α 2 ϕ 2 a 0 /2, b = r and c = u. The other scaled parameters are therefore given bỹ The dependence of the model on the coupling constant a 0 is fully taken into account by consideringb andc instead of b and c. The situation with the coupling constantα is more subtle. We expect that for large values ofα the gauge sector S g describes a U (1) gauge field on the fuzzy sphere, and as a consequence, the matter sector S m describes a (real) scalar field in the adjoint representation of the gauge group on the fuzzy sphere. More precisely we have in general X a = αϕ(L a + A a ), where A a is the U (1) gauge field which depends generically onα. For large values ofα, the gauge field is weakly coupled to the scalar field, and in the commutative limit N −→ ∞, the two fields become fully decoupled due to the commutator structure of the interaction. This is one of the main principles underlying our algorithm. Hence, the dependence of the model on the coupling constant α is also fully taken into account, in the limit N −→ ∞, by consideringb andc instead of b and c. The theory S g + S m describes therefore, for large values ofα and large values of N , a scalar phi-four on the fuzzy sphere. In all of the simulations reported in this article, we take a 0 = 1 andα = 10 for concreteness. The choice forα is dictated by the fact that a fuzzy sphere phase, in the model with r = u = 0 (the four dimensional Yang-Mills action), is known to persist only for values ofα given by [10] 3α ≥α * = 2.55 ± 0.1.
Let us discuss the phase structure of the pure potential model V (Φ). The ground state configurations are given by the matrices We compute V [Φ 0 ] = 0 and V [Φ γ ] = −r 2 /4u. The first configuration corresponds to the disordered phase characterized by < Φ >= 0. The second solution makes sense only for r < 0, and it corresponds to the ordered phase characterized by < Φ >= − r 2u U γU + . There is a nonperturbative transition between the two phases which occurs, not at r = 0, but at r = r * = −2 √ N u, which is known as the one-cut-to-two-cut transition 4 . The idempotent γ can always be chosen such that It is not difficult to show that this dimension is maximum at k = N/2 (assuming that N is even), and hence from entropy argument, the most important two-cut solution is the so-called stripe configuration given by γ = diag(1 N/2 , −1 N/2 ).
There are therefore three possible phase transitions, and as a consequence, there exists a triple point. The famous 2nd order Ising phase transition 0 −→ ± −r/2u 1 N . The famous 3rd order matrix phase transition 0 −→ ± −r/2u(1 N/2 , −1 N/2 ). Clearly then, there must exist also a transition between the Ising and matrix configurations, viz 1 −→ γ, which is expected to be a continuation of the Ising line to large values of the coupling constant u, and thus it is expected to be 2nd order.
In the numerical simulations, we will be interested in the values m 2 = µ = 0. As a test of our simulations, we will use the following exact Schwinger-Dyson identity 5 4 In terms ofb andc the critical value occurs atb = −2 √c . If the relation between r * and u were on the other hand linear, viz r * ∼ u, then we would have insteadb/ √ N ∼c. 5 By changing X a to X ′ a = (1 + ǫ)X a and Φ to Φ ′ = (1 + ǫ)Φ, in the partition function, we can derive from the invariance of the path integral this identity.
The operator IDE is given by

Algorithm and Simulation
The path integral we want to simulate is Let us now diagonalize the hermitian N × N matrix Φ by writing the polar decomposition Φ = U + ΛU , Λ = diag(λ 1 , ...., λ N ) for unitary N × N matrices U . The measure becomes In above [dU ] is the Haar measure on the group U (N ), whereas ∆ N (x) is the Vandermonde determinant. By using now gauge invariance of the above path integral, we can reabsorb the unitary matrix U , by changing X a as X a −→ U X a U + , and as a consequence, the integral over U decouples. The path integral becomes then The scalar action is, then, given by We will apply the Metropolis algorithm in which we change the eigenvalues λ i one at a time. Under the change of the eigenvalue λ i (fixed i), i.e. under λ n −→ λ ′ n = λ n + δ ni ǫ, the action S[Λ] changes as The first line is the variation of the kinetic term, the two first terms of the second line provide the variation of the potential, whereas the last term is the variation of the Vandermonde determinant. The variation of the action S[Λ], under the change of the entry (i, j) of one of the matrices X a , say X a −→ X a + ∆X a , is given by We choose (∆X a ) mn = δ ni δ mj ǫ * + δ nj δ mi ǫ.
The variation becomes We remark that for diagonal elements, i.e. i = j, this variation vanishes identically. This is simply due to the fact that the scalar kinetic action does not depend on diagonal elements of the matrices X a . The full variation under the change of the entry (i, j) of one of the matrices X a , which will enter the Metropolis algorithm, will naturally contain contributions coming from the pure gauge action. This part has been used elsewhere with great success [10]. The identity in this case still reads as in (16), with the operator IDE given by The Vandermonde action contributes to the integer 4N 2 , and as a consequence, it does not appear in IDE. It is very hard to generate, in the simulation, a sample of gauge and scalar configurations which satisfy this exact identity, due to the large degree of auto-correlation observed in the fuzzy sphere phase, i.e. for large values ofα. To reduce this undesirable effect, we separate any two successive configurations used in our measurements, by a large number of unused Monte Carlo configurations.
We measure the expectation value of the action < S m >, the total power P T , the power in the zero mode P 0 , the kinetic term < K >, the specific heat C v 6 , the magnetization M and the susceptibility χ. The action has already been defined . The other observables are defined by We use the Metropolis algorithm to update configurations, and we use the jackknife method to estimate error bars. The choice of the initial state is irrelevant. The Metropolis algorithm and the initial state used are discussed below in more detail. Typically, starting from a given/prepared initial state we run the Metropolis algorithm for T T thermalization steps to achieve thermalization, and T MC Monte Carlo steps for the actual Monte Carlo evolution. We record all of the T MC configurations and compute averages over them. Each two successive Monte Carlo steps are separated by T C auto-correlation steps. The value of T C can be chosen to be at least equal to the auto-correlation time, for a given set of parameters, which can be computed using the usual formula.

The Ising Phase Transition
In this case, the Metropolis updating procedure consists in going through the entries of each matrix X a , and through each of the eigenvalues of Λ, sequentially, and then attempting to change them in the usual way.
The initial state is prepared as follows. First, we start from Λ = 0 and X a = αL a , atb = 0, which we know is the true minimum at this point, and then run a Metropolis updating procedure, on this initial state keeping X a fixed, without taking into account the effect of the Vandermonde determinant, which is obviously the hardest part to thermalize, to obtain the actual initial state forb = 0. Using this initial state, we launch the full Metropolis updating procedure.
Next, we start changingb adiabatically (slowly), in such a way that the initial configuration for each new value ofb is the last configuration obtained for the previous value ofb. Each time, we run starting from this initial state, a Metropolis updating procedure, keeping X a fixed, and without the effect of the Vandermonde determinant, to obtain the actual initial state for that particular value ofb, before we launch the full Metropolis updating procedure.
We have checked that the location of the disordered-to-uniform-ordered transition does not depend on the above procedure, and thus it is fully independent of the initial conditions utilized.
A simulation consists typically of 2T T + T C × T MC steps where T MC = T T = 2 13 (N = 4, 6) or T MC = T T = 2 14 (N = 10), and T C = 2 5 . The first T T steps is done at fixed X a = αL a , and without the Vandermonde determinant.
The disordered-to-uniform-ordered transition is shown on figure (1). This transition can appear only for small values ofc. We take for examplec = 0.1. The 2nd order Ising transition (location of the peaks in the specific heat and the susceptibility) occurs  Table 1: The Ising transition points.
Using just these two points, we can determined the boundary between the disordered and the uniform-ordered phases, as a straight line, with slope given by The fit to the uniform-ordered-to-disordered transition line is given by (suppressing error bars because they are quite insignificant in this case) This agrees with [1]. We note that we have dropped out the intercept in the fit equation because it is, within statistical errors, completely negligible. This confirms the general expectation that the Ising line must go through the origin (c,b) = (0, 0). We note finally that this transition can also be obtained using the usual Metropolis algorithm with the ordinary pure scalar action, i.e. with the action (4), with X a fixed given by X a = αL a .

The Uniform-to-Non-Uniform Phase Transition
Thermalization and Tunneling: The non-uniform-ordered-to-uniform-ordered transition can appear only for medium and large values ofc. It is a second order phase transition, which is the continuation of the Ising transition, to larger values ofc.
The non-uniform-ordered phase is the phase associated with spontaneous breaking of rotational/translational symmetry on the fuzzy sphere 7 . This fact lies at the heart of its fundamental importance.
Firstly, we note that this transition is virtually impossible to be observed using the Metropolis algorithm, with the action (4), where X a = αL a .
We can probe the uniform-to-non-uniform phase transition (although still very difficult), using the Metropolis algorithm with the action (21), where X a is obtained itself via the Metropolis algorithm, with the action (3), where M = β = 0.
The initial state, for a fixedb andc, is prepared as follows. We start from a random configuration Λ, and from X a = αL a , and then run a Metropolis updating procedure for T T steps on this initial state, without taking into account the effect of the Vandermonde determinant at fixed X a , to obtain the actual initial state. Starting from this resulting state, we run a full Metropolis updating procedure for T T + T C × T MC steps. This whole process consists a single simulation.
A simulation consists typically of 2T T + T C × T MC steps, where T MC = T T = 2 13 , and T C = 2 5 for N = 6, and T MC = T T = 2 14 and T C = 2 6 for N = 8.
We have studied thermalization in great detail. Typically, we tend to repeat the same simulation T S = 2 7 +1 times, where each simulation is started from the final state obtained in the previous simulation. The goal is to assess tunneling transitions between the different vacua < Φ >∼ 1, γ and γ k .
As pointed out earlier the vacuum state < Φ >∼ 1 has always the smallest energy, whilst the vacuum state < Φ >∼ γ has always the largest energy. The other states are naturally somewhere in between. However, from entropy considerations, it is the state < Φ >∼ γ, which has the largest phase space volume, which can be seen from the size of the Grassmannian manifold U (N )/(U (k) × U (N − k)), given by the dimension At infinite N , we therefore expect that only < Φ >∼ 1 and < Φ >∼ γ are stable vacua and thus must be observed, while for finite N , tunneling transitions to other states are expected and will in fact also be observed.
Some of our results are: • We present, in figure (2) and (3), scatter plots for the kinetic action K and the magnetization M respectively, forc = 2.5, and various values ofb, for N = 6. Each point is a single simulation consisting of 2T T + T C × T MC steps. There are at most T S points. The first simulation has been started off from a random Λ and X a = αL a , whereas each successive simulation is started off from the final state obtained in the previous simulation.
Conversely, the magnetization in the vacuum state < Φ >∼ 1 corresponds to the largest plateau, while the magnetization in the vacuum state < Φ >∼ γ corresponds to the smallest (almost vanishing) plateau. In this case there are clearly four distinct plateaus.
• We observe, in figures (2) As |b| decreases, transitions away from < Φ >∼ 1 become more frequent, and scatter plots start showing various plateaus corresponding to the other vacuum states.
We conjecture that if we repeat the simulation a sufficient number of times T S , then the system will settle into its true minimum. This may take a long time only in the transition region between large and small |b|. It is immediately obvious, from the above discussion, that for large |b| the minimum is < Φ >∼ 1, while for small |b| the minimum is < Φ >∼ γ.
Eigenvalues Distributions: It is quite obvious, that the most revealing order parameter, is the eigenvalue distribution of the scalar field Φ. In our approach, the eigenvalues are precisely the degrees of freedom which we are sampling. We can then use immediately the T MC sets of eigenvalues λ i obtained in the Monte Carlo evolution, for a fixedc and b, to construct appropriate histograms. These are precisely the eigenvalue distributions ρ(λ) of the scalar field Φ.
In figure (5), we plot the eigenvalue distributions for various values ofb, across the uniform-to-non-uniform transition point, for N = 6 andc = 2.5. We observe that we go from the one-cut solution, centered about + −r/2u, to the two-cut solution, centered about ± −r/2u, aroundb = −10.5 ± 0.5, which agrees with our other measurement (see below).
Although in the two-cut solution we know that the eigenvalues are ± −r/2u, we can not tell how many of them are pluses, and how many of them are minuses 9 . In order to determine the distribution of the plus and minus signs, we may then plot, the probability distribution of the values of the magnetization T rΦ. Alternatively, we can directly look at the eigenvalues themselves, to see which matrices are involved. As it turns out, in the transition region between large and small |b|, the vacuum states are not given simply by the pure states 1, γ, γ 1 and γ 2 ,but they are, typically, given by admixture of these pure states.
Critical Values: According to [1], the non-uniform-ordered-to-uniform-ordered transition, should occur at the value ofb, where the susceptibility and the specific heat are peaked, which is something we were not able to reproduce in our scheme in any consistent way.
The determination of the location of the non-uniform-ordered-to-uniform-ordered transition, can also be based, on the location of the "discontinuity/jump" in the expectation value of the kinetic term. This discontinuity is also associated with a discontinuity in the total power, power in the zero mode and magnetization.
As opposed to all other simulations reported in this article, we will attempt in the current case to cross the critical line by holding −b fixed, while varyingc. In this way, we are guaranteed to cross, first, the non-uniform-ordered-to-uniform-ordered transition, as we increasec, at some fixed value of −b. If we fixc instead, and start increasing −b, we will hit the matrix phase transition first (see next subsection), then the non-uniformordered-to-uniform-ordered transition.
The detail of this simulation goes as follows. The initial state, for a fixedb andc, is prepared by starting from a random configuration Λ, and from X a = αL a , and then run a Metropolis updating procedure for T T steps on this initial state, at fixed X a without taking into account the effect of the Vandermonde determinant. We repeat this process for T S = 2 4 steps to get the actual initial state. Starting from this resulting state, we run a full Metropolis updating procedure for T T + T C × T MC .
We work always with T T = T MC = 2 13 , and T C = 2 4 , for N = 6, 8, 10. The constraint on the identity is < IDE > /N 2 = 4.00 ± 0.30 (N = 6), and < IDE > /N 2 = 4.00 ± 0.25 (N = 8, 10). The results are shown on figure (4). In the graphs of the total power, and the power in the zero mode, we can find from the scaling (9), that in the Ising phase P 0 = P T ∼ −b/( √ Nc), which is why the graphs for the powers for different N do not collapse.
We will take, as our measurement of the non-uniform-ordered-to-uniform-ordered transition points, the arithmetic average of the critical points, obtained from the discontinuity/jump in the expectation value of the kinetic term for different N 10 . We drop here the calculation of the error bars which requires much more efforts. Some results are given in table (2).
The fit to the non-uniform-ordered-to-uniform-ordered transition line, as computed from table (2), is given byc = −0.22b + 0.38. The slope is very close to the slope of the Ising transition line given by equation (34). This confirms the general conjecture of [1], that the non-uniform-ordered-to-uniform-ordered transition line, is the continuation, of the Ising transition line, to general values ofc andb. However, the intercept of the fitc = −0.22b + 0.38 seems to be quite large. We claim that this is, only, due to our limited number of data points, and lack of error bars. Clearly, forc = 0, there is no Ising transition, nor a non-uniform-ordered-to-uniform-ordered transition. In other words, the non-uniform-ordered-to-uniform-ordered transition line must go through the origin (c,b) = (0, 0). The fit to the non-uniform-ordered-to-uniform-ordered transition line, as computed from table (2) plus the point (c,b) = (0, 0), is now given bỹ The error in the intercept is found to be 0.1, while the error bar in the slope is negligible. The measured slope, as well as the measured small intercept, are reasonably close to the values measured in [1].

The Matrix Phase Transition
The non-uniform-ordered-to-disordered transition, also called matrix transition, appears for medium and large values ofc. We perform simulations in a similar fashion to the Ising case, with the exception that we start from a random configuration for each value ofb. We takeb in the range [−15, 0], with step equal 0.25, and values ofc in the range [2,25].

The Matrix Transition in the Limit of Large Couplings: It is expected that
for large values of the coupling constantc, the matrix transition in the full model, will be given approximately, by the matrix transition in the pure potential model, i.e. the model without kinetic term. This approach becomes exact in the limitc −→ ∞.
We include in figure (6), the behavior of the magnetization M = |TrΦ|, the zero power (power in the zero modes) N 2 P 0 , the sepcific heat C v /N 2 , and the average action < S m > forc = 16. We plot the pure potential model for comparison.
It is well known that the matrix transition occurs, in the pure potential model, at the point where the specific heat divided by the number of degrees of freedom becomes equal to 1/4, after passing through its minimum as we increase |b|. This corresponds, for any fixed valuec, to the transition pointb * = −2 √c . This transition is anticipated by the intersection point, which is N −independent, seen on the graph of the action < S m >, and by the location of the wide maximum, seen on the graphs of the magnetization M and the zero power P 0 . However, all these estimates, provide only an under estimation of the actual transition point in the pure potential model.
If we take, as our measurement of the matrix transition in the full model, the point where the specific heat becomes equal to 1/4 after passing through its minimum, then we find, as opposed to the pure potential model, an under estimation of the transition point. The intersection point of the action < S m > provides, as before, also an under estimation of the transition point.
In the full model, we have observed that, for sufficiently large values ofc, a reasonable estimation of the matrix transition point, which compares favorably to the theoretical prediction coming from the pure potential model, can be given by the location of the broad maximum, seen on the graphs of the magnetization and the zero power.
We search for this maximum for values ofb much smaller than the discontinuity point relevant for the non-uniform-to-uniform transition.

Eigenvalues Distributions and The Behavior Near the Triple Point:
We have also investigated the matrix transition at the level of the eigenvalues distributions. In principle, the matrix transition occurs where the eigenvalues distributions split into two disjoint supports (cuts). In other words, it occurs at the point, where the distribution goes from a symmetric one centered around 0 (as opposed to being centered around either + −r/2u or − −r/2u in the case of the non-uniform phase), to a distribution with two symmetric cuts centered respectively around −r/2u and − −r/2u. A sample of the eigenvalues distributions, in the full model and in the pure potential model, are shown on figure (7) for N = 6 andc = 6. We have used the eigenvalues distributions of Φ, as the primary set of order parameters, employed in the determination of the matrix transition point, for smaller values of the coupling constantc. Following [1], we have considered the regime [2,3]. This is the regime of interest to the calculation of the triple point (more on this below). We note that, the method employed above (maximum of magnetization and zero power), becomes unpractical in this regime. The results obtained for N = 4, 6, 10 are included in table (4), and compared to the estimation of [1].
We have determined the matrix transition point according to the following (somewhat arbitrary) criterion. We have looked for the value ofb, for which the eigenvalues distribution ρ(λ) at λ = 0, drops below 1. The transition point is taken as the arithmetic average of this value ofb, and the next one, for which, typically, the eigenvalues distribution at λ = 0 becomes distinctly below 1. A similar technique, to determine the matrix transition point, is employed in the recent thesis [11].
A sample of the eigenvalues distributions of Φ is shown on figure (8) forc = 2.5. We also include, a sample of the probability distribution of TrΦ, which may be used to determine the actual content of a given configuration Φ. The number of pluses and minuses, can only be inferred, from the plot of the probability distribution of TrΦ. If Φ is a fluctuation about 0 or γ, then the probability distribution of TrΦ will contain a single symmetric peak around 0. There is also the possibility that Φ is a fluctuation about γ k , then the probability distribution of TrΦ will contain a single symmetric peak around −r/2u(2k − N ). Typically, Φ will fluctuate about a mixed state, and as a consequence, several peaks will be present in the probability distribution of TrΦ . For example, if Φ is a mixture of γ k and γ, then, two peaks centered around −r/2u(2k − N ) and 0 will be present. Some examples are shown on figure (8).
Using the results shown in table (4), we can determine the non-uniform-ordered-todisordered boundary. The fit to the matrix (non-uniform-ordered-to-disordered) transition line is given byc This line is slightly different from the one measured in [1], which may be due to our criterion for determining the matrix transition point. However, we should also recall that their result was obtained using a modification of the Metropolis algorithm which breaks detailed balance.

Triple Point and Phase Diagram
The most reliable estimation of the triple point can be obtained from the intersection point of (34) and (36), because these two lines are the easiest, and the most accurate, to obtain with our gauge fixed Metropolis algorithm, and also with the algorithm of [1].
In fact, they can even be accessed using the plain Metropolis algorithm. We deduce immediately that the triple point is located at The phase diagram is shown on figure (9). See also figures (10) and (11), where a close-up look at the matrix, the Ising, and the non-uniform transition lines is shown.

Comparison of Various Algorithms
The algorithm used in [1] to compute the phase diagram is based on, a very complex variation, of the Metropolis algorithm, which does not preserve detailed balance. In the region of the disordered phase, their algorithm behaves essentially as the usual Metropolis algorithm, with a processing time per configuration, with respect to the matrix size, proportional to N 4 . The new Metropolis algorithm, described in [1], behaves better and better, as we go farther and farther, from the origin, i.e. towards the regions of the uniform and non-uniform phases. The processing time per configuration, with respect to the matrix size, is claimed to be proportional to N 3 , for the values of N between 4 and 64. See graph 9.12 of F.G Flores' doctoral thesis 11 , where we can fit this region of N with a straight line. Also, it is worth noting, that this new algorithm involves, besides the usual optimizable parameters found in the Metropolis algorithm, such as the acceptance rate, a new optimizable parameter p, which controls the compromise between the speed and the accuracy of the algorithm. For p = 0 we have a fast process with considerable relative systematic error, while for p = 1 we have a slow process but a very small relative error. This error is, precisely, due to the lack of detailed balance. Typically we fix this parameter around p = 0.55 − 0.7.
The algorithm of [1] is the only known method, until now, which is successful in mapping the complete phase diagram of noncommutative phi-four on the fuzzy sphere. However we had found it, from our experience, very hard to reproduce this work. 11 Not available on the ArXiv.
Our first original goal was to find an alternative method which is, i) conceptually as simple as the usual Metropolis method, and ii) without systematic errors, and iii) can map the whole phase diagram. This goal was achieved by the algorithm described and used in this article. The processing time per configuration, with respect to the matrix size, in our algorithm, is proportional to N 4 , which is comparable to the usual Metropolis algorithm, but with the virtue that we can access the non-uniform phase. There is no systematic errors in this algorithm, and hence no analogue of the parameter p mention above.
How does our algorithm compares with the algorithm of [1], is a much harder question, since we have no complete understanding of the detail of their algorithm. Their algorithm is faster, but this can not be the only concern. Accuracy of the method, and conceptual simplicity, are also very important virtues, especially, for difficult problems, such as this one, where the physics is extremely interesting, but very hard to attain. Our algorithm satisfies both these two requirements.
Our other goal, in this article, was to compare the results obtained by the two methods for the non-uniform phase. There are still discrepancies between the two methods which is very puzzling. The non-uniform phase is characterized, in this article, by the "discontinuity/jump" in the expectation value of the kinetic term, the total power, power in the zero mode and magnetization. According to [1], this jump is also associated with a peak in the susceptibility and specific heat indicative of a second-order behavior, which is something we were not able to reproduce in our scheme, in any consistent way. This is very troubling, to say the least, because we could not, from what we have and know at this point, ascertain whether this is due to a technical problem, or if it is a genuine discrepancy. 4 The Self-Dual Noncommutative Φ 4 on the Fuzzy Sphere 4.1 Self-Dual Noncommutative Φ 4 We consider, for simplicity, a real scalar field on the noncommutative (Moyal-Weyl) plane [x µ ,x ν ] = iθ µν . The phi-four theory on the noncommutative plane is, a particular limit, of a one-parameter family of phi-four models on the noncommutative plane, obtained by the addition of an extra operator, the harmonic oscillator potential , to the kinetic part of the action. The action reads explicitly We know that derivations on R 2 θ are inner, given by the adjoint action, viẑ Alternatively, the action can be rewritten as This is the Grosse-Wulkenhaar model. The addition of the harmonic oscillator potential to the kinetic action modifies, and thus allows us to control, the IR behavior of the theory. A particular version of this theory was shown to be renormalizable by Grosse and Wulkenhaar in [12]. It was shown in [13], that this action is covariant under a duality transformation which exchanges, among other things, positions and momenta. The value Ω 2 = 1, in particular, gives an action which is invariant under this duality transformation. The theory at Ω 2 = 1 is called the Langmann-Szabo model or the self-dual Grosse-Wulkenhaar model. The usual phi-four theory on the noncommutative plane corresponds to the limit Ω −→ 0. The other interesting limit is Ω −→ 1, which corresponds to the self-dual Grosse-Wulkenhaar model. The main technical simplification, occurring in the limit Ω −→ 1, is the observation that the off-diagonal term in the action drops, and we end up with the action 12 Let us now introduce creation and annihilation operators a + and a satisfying [a, a + ] = θ byx The number operatorN is defined byN = a + a/θ. We can verify, for example, that x 2 µ = 2θN + θ. We will work in the number basis defined bŷ N |n >= n|n > , a + |n >= θ(n + 1)|n + 1 > , a|n >= The components ofφ, in the number basis, are given byφ nm =< n − 1|φ|m − 1 >. In the number basis {|n >} the action S Ω reads explicitly This is a three-parameter model, where the mass parameter r and the quartic coupling u, are given by The other coupling is the harmonic oscillator coupling Ω. 12 After regularization this action becomes the Penner matrix model.

Fuzzy Sphere as a Regulator
In the remainder of this section, we will write down a non-perturbative regularization of this theory on the fuzzy sphere. We only need to consider the kinetic term. Let L a be the generators of SU (2) We will work in the basis {|m >} defined by the usual relations L 3 |m >= m|m >, L ± |m >= l(l + 1) − m(m ± 1)|m ± 1 >, where l = (N − 1)/2 and L ± = L 1 ± iL 2 . We relabel the basis as |m >= |i >, where m = i−l−1. We compute ( Rotating around the x-axis, with an angle π, we have L 1 −→ L A real scalar fieldφ is a hermitian N ×N matrix which will be expanded in the obvious wayφ We start by considering a more general Laplacian, obtained by adding a harmonic oscillator potential to ∆ 0 , in the most obvious way. First, we introduce the coordinates operatorsx a , on the fuzzy sphere, byx a = RL a / √ c 2 , which satisfy [x a ,x b ] = iRǫ abcxc / √ c 2 andx 2 a = R 2 . We define the right-acting coordinate operatorsx R a byx R aφ =φx a , and then introduce the coordinates operators X a by We define the Laplacian In other words, we consider the kinetic term 13 By sitting on the north pole, i.e.x 3 = R1 N , and taking the limit N −→ ∞, and R −→ ∞, keeping R 2 / √ c 2 = θ fixed, the fuzzy sphere reduces to the noncommutative plane.
The normalization 4πR 2 is chosen such that in the commutative limit N −→ ∞ we have ( The last term is not present on the noncommutative plane which is, clearly, an unwanted effect. After some trial and error, we have discovered, that the correct Laplacian on the fuzzy sphere, which reproduces precisely the effect of the harmonic oscillator potential, is given by This will describe a squashed fuzzy sphere, which is more appropriate, for the nonperturbative description of the noncommutative plane. Indeed, we compute We can now include a mass term and a phi-four coupling in a trivial way. The full action, on the fuzzy sphere, will read We scale the field asφ =φ/ √ 2π, and also introduce the parameters The full action can then be rewritten as This is a one-parameter family of phi-four models on the fuzzy sphere which generalizes (1). Coupling to a U (1) gauge field is straightforward, i.e. we make the replacement L a −→ N/2X a . The analogue of (4) is obviously given by

Conclusion and Outlook
In this article, we have proposed a new algorithm for the Monte Carlo simulation of noncommutative phi-four on the fuzzy sphere, and also reported our first numerical results on the corresponding phase diagram, obtained with small values of N up to N = 10, and large numbers of statistics. Basically, the new algorithm employs gauge invariance in order to reduce the scalar sector to the core eigenvalues problem. The phase diagram is complex consisting of three transition lines: the Ising or uniform-to-disorder, the matrix or non-uniform-to-disorder, and the uniform-to-non-uniform transition lines. These lines intersect at a triple point. The measurement of the uniform-to-non-uniform transition line, using our algorithm, remains very demanding but tractable. The measurements, included in this article, are largely consistent with those reported originally in [1].
The first immediate extension of this work is to optimize the algorithm further, and push the calculation of the phase diagram to higher values of N , with reasonably large numbers of statistics, especially in the case of the uniform-to-non-uniform transition line. We note that a major improvement of our algorithm, may be achievable, by replacing the Metropolis updating procedure, for the scalar eigenvalues problem, by the Hybrid Monte Carlo algorithm, whereas we may keep using the very efficient Metropolis for the gauge sector.
Another immediate line of investigation is the calculation of the phase diagram of the self-dual noncommutative phi-four on the fuzzy sphere, constructed in the last section. The main question, here, is what happens to the Ising transition line, as Ω goes from Ω = 0 to Ω = 1, and as a consequence, what is the fate of the triple point.           (36), which must be extrapolated to even smaller values ofc, in order to deduce an estimation of the triple point. We also compare, for small values of c, with the measurement of [1]. The discrepancies between the two measurements, for small c, is stemming from our criterion, based on the eigenvalues distributions, for determining the location of the matrix transition, which is different from the one used in [1].  Figure 11: The Ising (disorder-to-uniform) and the non-uniform-to-uniform transition lines.
The Ising transition appears for small values ofc, while the non-uniform-to-uniform transition appears for large values ofc. The Monte Carlo measurements of the Ising transition is fully consistent: more data points can be included quite easily, error bars are under control, large N extrapolation is straightforward, and result obtained by our algorithm coincides with the measurement of [1]. On the other hand, the two Monte Carlo measurements of the non-uniformto-uniform transition, included in this graph, required much more calculation than their Ising and matrix counterparts put together. We did not attempt, here, to determine their error bars. The measured slope and small intercept, of the resulting non-uniform-to-uniform fit, are reasonably close to the measurements of [1]. Work on this major problem, i.e. a fully consistent determination of the non-uniform-to-uniform transition line, is still in progress. We also plot the matrix line where the intersection points, with the Ising and the non-uniform-to-uniform lines, provide our two estimations of the triple point.