The non-perturbative phase diagram of the BMN matrix model

We study the maximally supersymmetric plane wave matrix model (the BMN model) at finite temperature, $T$, and locate the high temperature phase boundary in the $(\mu,T)$ plane, where $\mu$ is the mass parameter. We find the first transition, as the system is cooled from high temperatures, is from an approximately $SO(9)$ symmetric phase to one where three matrices expand to form fuzzy spheres. For $\mu>3.0$ there is a second distinct transition at a lower temperature. The two transitions approach one another at smaller $\mu$ and merge in the vicinity of $\mu=3.0$. The resulting single transition curve then approaches the gauge/gravity prediction as $\mu$ is further decreased. We find a rough estimate of the transition, for all $\mu$, is given by a Pad\'e resummation of the large-$\mu$, 3-loop perturbative, predictions. We find evidence that the transition at small $\mu$ is to an M5-brane phase of the theory.


Introduction
The BMN or plane wave matrix model model [1] is the matrix model description of Mtheory in discrete lightcone quantisation (DLCQ) when the M2-brane is propagating on an eleven dimensional plane wave background. This background preserves the maximal, sixteen, supersymmetries and is a massive deformation of the BFSS model [2][3][4].
In contrast to the BFSS model, the BMN matrix model has a discrete energy spectrum. It is conjectured to capture the entire dynamics of M-theory on the plane wave background and provide a non-perturbative definition of M-theory itself. Alternatively, it describes a system of D0-branes of IIA supergravity. It has BPS ground states 1 , where the D0-branes expand to fuzzy spheres [5,6]. In [7] it is argued that the model also provides a regularisation of NS5/M5-branes.
For large µ the model becomes a supersymmetric gauged Gaussian model. Perturbation theory around this limiting model was performed up to three loop order [8,9] and a large µ series for the transition in the Polyakov loop was developed. The dual gravitational theory was studied by perturbing around the BFSS dual geometry to linear order in small µ in [10] and the corresponding dual gravity prediction for the transition was obtained to this order.
The gravity dual has a remarkably rich structure. At zero temperature it is described by bubbling geometries. The dual bubbling geometries are in eleven-dimensional supergravity compactified in one translationally invariant direction on the two-dimensional subspace [11], or equivalently, in type IIA supergravity [12]. The supergravity solutions are labeled by configuration of D2/M2 and NS5/M5 charges, discretely located on a one-dimensional subspace 1 The vacua of the BMN model correspond to 1 2 -BPS states of the eleven dimensional or IIA supergravity. They are invariant under infinitesimal transformations of the SU (2|4) supergroup, which has sixteen supercharges.
of the geometries, and each of them corresponds to a vacuum of the BMN model. In the large N limit there are infinitely many solutions, and interesting special solutions, such as those corresponding to a single stack of D2/M2-branes or NS5/M5-branes, exist. At finite temperature, the system has a Hawking-Page transition-a transition between thermal and black hole spacetimes. The solution examined in [10] has S 8 × S 1 horizon topology, which corresponds to the trivial vacuum in the BMN model, where the S 1 is the M-theory circle. It is surprisingly difficult to obtain general gravity solutions dual to the BMN model at finite temperature; even the solution in [10] required numerical computations to obtain the free energy.
There is some understanding of the dual geometry from the gauge-theory side. Since the corresponding supergravity solutions are supposed to be obtained at low energy, the geometry should be constructed by a low-energy moduli operator in the gauge theory. In [13,14], the BPS operator that is expected to pick up the low-energy moduli was computed by the supersymmetric localisation method, and in the appropriate large-N and strong coupling limit, it was found that its eigenvalue distribution satisfies the same integral equations as those determining the supergravity solution. These equations govern a non-trivial part of the dual metric, which is not determined by the isometry 2 . Although these results successfully reproduced part of the supergravity solution, the dynamics of the emergence of those geometries, such as how the geometries are favoured or superposed as the temperature changes, is beyond the scope of such methods; therefore numerical simulations for the thermal theory provide important information on which geometry emerges naturally. This paper is dedicated to obtaining a non-perturbative estimate of the high temperature phase boundary as a function of the mass parameter, µ.
We find that a Padé resummation of the large µ perturbative series [9], with the assumption that the transition temperature goes to zero linearly with µ, gives a result very close to that of the gravity calculations [10] when expanded to linear order in small µ. We use this interpolating Padé resummed curve as a guide to the location of the transition. In practice, we find it is an excellent guide, giving a reasonable approximate location of the transition for all values of µ.
We study the BMN model non-perturbatively using Rational Hybrid Monte Carlo techniques 3 . We map the phase diagram in the (µ, T ) plane and determine the phase boundary as the system cools from the high temperature phase. We find that for large µ, a sharp transition in the Polyakov loop locates the transition. However, for µ ∼ 7.0 we find that there are two phase boundaries and a new phase appears characterised by the Myers cubic term 2 Related with the isometry part, the spherical shell distributions of S 2 and S 5 can be obtained in the BMN model in the limits where M2-and M5-branes are realised, respectively. The radii of these spheres completely agree with the M2 and M5 radii in the brane picture [15,16]. 3 A preliminary study of the BMN model concerning a phase transition was carried out in [17]. The results differ significantly from ours at lower temperatures due to, we believe, the larger lattice effects inherent in their first order code. In [18] and [19,20], the transition temperature was estimated by using an effective theory. Another numerical simulation of the full BMN model was done around a special vacuum in [21].
in the action acquiring a significant non-zero vacuum expectation value and the approximate SO(9) symmetry no longer holds. We observe that this transition, which we call a Myers transition, is from a phase where the system fluctuates around the trivial configuration to one with fluctuations around expanded fuzzy spheres. In the transition region we observe that the system also fluctuates between different intermediate fuzzy sphere configurations. The transition should be in the universality class of the emergent geometry transitions discussed in [22][23][24]. As the system is further cooled, with fixed µ (e.g. see µ = 6.0 in Figure 2), the Polyakov loop undergoes a separate transition but at a significantly lower temperature. For µ < 3.0 we find the two transitions merge, and the combined transition is visible in both the Myers observable and the Polyakov loop and has a significant jump in the energy, E = H /N 2 , where H is the Hamiltonian (see µ = 2.0 in Figure 3). For µ < 3.0, the transition in the Polyakov loop is rather difficult to observe but one can check that the eigenvalue density of the holonomy undergoes a transition from a gapped to a gapless spectrum as the system is cooled through the Myers transition temperature.
For small µ, in the transition to a stable fuzzy sphere phase the Polyakov loop tends to increase sharply (see Figure 4) to a higher value before decreasing slowly. In turn, for fuzzy sphere backgrounds, the transition in the holonomy to a gapless phase occurs at a considerably lower temperature than the Myers transition; the latter transition coincides with the region where the holonomy associated with the trivial background becomes gapless. These different temperatures would correspond to different Hawking-Page transitions for different solutions in the gravity dual. This result is consistent with the upper bound for the critical temperature discussed in [10].
At large N the system can be trapped in a given vacuum phase and the transition between such vacua, being a quantum effect, is suppressed. This stability was used in [21] to study N = 4 super Yang-Mills from the BMN model. The parameter region they used is located at lower temperatures than the Myers transition we observe, and thus consistent with our findings.
A key ingredient in ensuring that the system was in the preferred thermodynamic state, and not trapped in a false vacuum, was to cool the system adiabatically so that E(T ) remained a monotonic function of temperature.
In the intermediate µ regime the transition in the Myers term occurs at a significantly higher temperature than that predicted by our Padé curve while the transition in the Polyakov loop tracks it more closely. For smaller µ < 2 the Myers transition again approaches the Padé curve and begins to qualitatively agree with the dual gravity prediction [10]. The observed transition also is from small to large Myers term. We interpret this as meaning that on the gravity side a spherical black hole becomes unstable to asymmetric deformations. To our knowledge, the relevant gravity solution describing the black hole corresponding to a fuzzy sphere vacuum has not yet been constructed.
The principal results of this paper are: • A determination of the high temperature phase boundary of the BMN matrix model (see Figure 5).
• A Padé approximant estimate for the phase transition in the trivial vacuum.
• A detailed non-perturbative study of the model for N = 8.
• An extrapolation to large N of the Myers transition for µ = 2.0.
• For µ = 2.0 we find the Myers observable takes a finite, N independent, non-zero value in the range studied. This suggests the transition can be interpreted as one to an M5-brane phase of the theory.
The paper is organised as follows. Section 2 describes the model and the Padé resummation of the large µ transition curve. Section 3 summarises our lattice formulation of the model. Section 4 discusses the observables we measure. Section 5 gives our results and describes the phase diagram we find. We end with some concluding remarks and discussion.
In the large µ limit the model (2.1) reduces to a supersymmetric Gaussian model. This simple Gaussian model undergoes a confinement-deconfinement phase transition as the temperature is lowered and a straightforward calculation gives the critical temperature in this limit as T c = µ 12 ln 3 . This transition has been studied perturbatively in 1/µ [8,9], where in a three loop calculation it was found that This result, while reliable for large µ becomes untrustworthy as µ decreases and passes through zero for µ 13.4. However, if we perform a Padé resummation of (2.2), assuming that T c → 0 linearly with µ as µ → 0, then we can rewrite (2.2) as which to linear order in small µ leads to the prediction T c = µ 12 ln 3 1 + 1638400 5499852 + 1765769 ln 3 0.0925579µ.
This value is surprisingly close to which was obtained from a rather involved dual gravity computation [10]. We use the Padé resummed result (2.3) as a guide to where one might expect the transition in the full model as µ is decreased and study the system using the rational hybrid Monte Carlo algorithm and a novel lattice discretisation described in [25].
The action (2.1) is minimised by the fuzzy sphere configurations, where J r are generators of SU (2) in an arbitrary representation of total dimension N . These are BPS states and are protected ground states of the quantum Hamiltonian.

Lattice Formulation
We use a second order lattice discretisation of the model. In this formulation the time interval is replaced by a periodic lattice with Λ sites. The matrices are located on the lattice site, and the lattice spacing is a = β Λ . The bosonic Laplacian is discretised using the lattice version In the fermionic action the Dirac operator is discretised as is a Wilson term that suppresses fermionic doublers and is a slightly more general discretisation of the derivative (the standard first order lattice derivative is recovered setting r = 0). Here, K a is an anti-symmetric lattice operator while K w is symmetric. The overall Dirac operator CD Lat is an anti-symmetry matrix. The charge conjugation matrix C is symmetric as are Cγ i , while Cγ ij and Cγ ijk are the only anti-symmetric elements of the Clifford algebra. An alternative to using Σ 123 would have been to use Σ 89 = iγ 89 which was the choice used in [17,[26][27][28][29] and our earlier code [30]. However, Σ 89 , in contrast to Σ 123 , does not anti-commute with the mass term. Different components of the spinor then have different dispersion relations and µ-dependent lattice effects which grow with µ. We also found that using the Σ 89 prescription and a first order discretisation (i.e. where r, r b and r 2f are zero) had large lattice effects for moderate values of µ. We discuss details of our code and other technical issues in [25]. The dispersion relations for the two options are illustrated in Figure 1, where one can see two curves for the Σ 89 prescription. As the temperature is lowered, with fixed Λ, the lattice spacing becomes larger and the splitting becomes more extreme.
In practice, the parameter values for our simulations were r = − 1 3 , r 1f = 0, r b = 1 12 , r 2f ≈ 0.148. We chose r = − 1 3 as it eliminates the cubic term in the expansion of K a and we found that choosing r 1f = 0 minimised a lattice artifact in the inverse to D Lat . For comparison the parameters used in [28] were r = −1, r 1f = 0, r 2f = 1 4 , r b = 3 4 and Σ 123 should be replaced with Σ 89 .
For µ = 6.0 and different parameter choices, we performed tests of our Σ 123 formulation against the Σ 12 option. Fixing parameters to those of our current study we found results of the Σ 123 choice with Λ = 24 comparable to the Σ 12 option with Λ = 48.
In this paper we concentrate on Λ = 24 and N = 8; however, preliminary simulations comparing Λ = 24 and Λ = 48 show that the lattice effects are small especially for µ > 4.0 in the range of temperatures we studied. More generally, we found no observable difference in the location of the transitions between the two lattice sizes.

Observables
The BMN model is expected to have many phases and a rather complicated phase structure would be natural at low temperatures due to the multitude of zero energy BPS states. Our goal in this paper is to map the high temperature phase boundary to this region. We therefore study the system as the temperature is lowered, typically beginning our study at T = 2.0 for a fixed µ.
In the path integral formulation used in our study one can show, following the arguments of [31], that the energy is given by For the BFSS model this reduces to and (4.2) together with the Ward identity, which reads was typically used to express the energy without fermionic terms as While this expression provides coding convenience it is not necessary as one can measure these fermionic observables using pseudo fermions. Also, note that (4.4) involves the difference of large numbers and is limited by the precision of the rational approximation to the fermionic contributions. Using this expression to determine the energy in a precision calculation can    become a problem especially when approaching the continuum limit, a limit involving sending Λ to infinity. In practice we find that the direct observable (4.2) behaves better numerically. For the BMN model we implement (4.1) directly using pseudo fermions following the strategy discussed in section 4.2 of [30] and used in [35]. We also, for completeness, observe the BMN Ward identity analogue of (4.3).
The observables we follow are then: E as defined in (4.1), (4.5) The energy and specific heat of the system, which due to supersymmetry, should both decrease to zero as T approaches zero, are given by We see, furthermore, that E must be a monotonic function of T . This proves especially useful in the simulations as tracking the energy as a function of T was a crucial clue in identifying when the system transitioned to a new level 4 . With this strategy and careful simulation as the transition was approached we identified the phase diagram as shown in Figure 5.

Results and Phase Diagram
For large µ, e.g. µ = 18.0, we find that the transition is well tracked by the 3-loop perturbative result (2.2); however, deviations arise as µ is reduced. We concentrate our efforts on µ ≤ 9.0. We find that for µ ≤ 6.5 the system undergoes a phase transition from a small Myers observable to a large one and there is no longer an approximate SO(9) symmetry. This is clear when observing R 2 ii , which at temperatures close to the transition has large fluctuations in the SO(3) components (see Figure 2). A reasonable order parameter for this transition is the Myers observable (4.5), also shown in Figure 2 for µ = 6.0 as a function of temperature. Figure 2 shows there are two separate transitions with the Myers transition occurring at the higher temperature. The figure also shows R 2 ii for each of the nine matrices during a run corresponding to the transition value T c1 = 0.690 marked as the red vertical line in the Myers plot and the higher temperature line in the energy plot. One can see that the SO(3) matrices fluctuate between different fuzzy sphere vacua, while there are only small variations in the SO(6) matrices, which are concentrated at smaller values.
The deconfinement transition is clear in the plot of the Polyakov loop and is measured to occur at T c2 = 0.465. It is the second red line marked in the energy plot. The corresponding eigenvalue distributions of the holonomy for the two transition temperatures are plotted below and one can see that the confinement-deconfinement transition involves the eigenvalue distribution changing from gapped to ungapped. One can also see evidence for both transitions in the plot of E, though the Polyakov transition is less clear due to the relatively small jump in energy across this transition. However, as the temperature is lowered the transition in the energy becomes more pronounced and the two transitions approach, eventually effectively merging in the vicinity of µ = 3.0. Figure 3 shows our observables for µ = 2.0 with N = 8 and Λ = 24. Only one transition was observed and this transition is clear in the energy plot.
Once the system is cooled to a temperature below T c1 the SO(3) matrices settle into the thermodynamically favoured fuzzy sphere configuration, and the SO(6) ones fluctuate around smaller values. Though larger fuzzy sphere vacua exist, further transitions between fuzzy sphere configurations were not observed as the system was cooled further.
We present our overall measured phase diagram in Figure 5. The figure is restricted to N = 8 with Λ = 24 and was determined from the Myers observable, Polyakov loop and energy observables. The simulations were performed beginning at T = 2.0; then T was slowly decreased for fixed µ using the thermalised input obtained at the higher temperature as an initial configuration for the lower T . The strategy of slowly decreasing the temperature for fixed µ proved crucial in determining the transition, as in many cases hysteresis made it difficult for the simulation to find the correct phase at the lower temperature. Which phase was correct was determined by carefully observing that the energy remained a monotonic function of T . The simulations indicate that the observed transitions are all first order, with a significant latent heat for µ = 2.0.
In the second graph in Figure 5 we focus on µ = 2.0 and extrapolate the critical temperature to large N estimating that T c (∞) = 0.35 ± 0.01, where the measured value for N = 8 was T c (8) = 0.345 ± 0.006. From this example and preliminary results for other µ we expect that the large N curve is shifted upwards slightly towards higher temperatures and we infer that the transition curve as measured in Figure 5 is a lower estimate for the large N coexistence curve. Figure 4 shows the energy and Polyakov loop for N = 6, with µ = 2.0. We see that the transition occurs at T c (6) = 0.31 ± 0.02, which is lower than T c (8); both are used in the 1 N extrapolation of Figure 5. Figure 4 also shows the Polyakov loop and we see it makes a sharp transition to a larger value. From this we infer that the fuzzy sphere configurations have typically larger Polyakov loop at a given temperature than the approximately SO(9) symmetric phase. We further find that in this phase the holonomy becomes ungapped at T ∼ 0.15 ± 0.02.

Conclusions
We have mapped an initial phase diagram for the BMN matrix model in the (µ, T ) plane. We have concentrated on relatively small matrix sizes, with most of our data for N = 8. The resulting diagram is presented in Figure 5.
We found that for small enough N and large enough temperature ergodicity is not a problem. However, for large N ergodicity is lost in the simulations, we therefore restricted our study to small N . For N = 6 we tracked the system to relatively low temperatures (down to T ∼ 0.1 with µ = 2.0). A little below the transition temperature it was apparent that maintaining ergodicity was beginning to be a problem as transitions between levels required rather lengthy simulations. For N = 8 and N = 11 there were difficulties with ergodicity, however, cooling the system in sufficiently small temperature intervals allowed us to access the transition. We suspect that going beyond N = 16 would require new techniques to overcome these difficulties.
We have observed two transitions for larger values of µ, which appear to merge for µ ∼ 3.0. The higher temperature transition curve, which shows the Myers transition curve, is the transition from the trivial vacuum to a non-trivial fuzzy sphere vacuum. In this transition the ground state consists of the SO(3) matrices blowing up into fuzzy spheres and the transition has the characteristics of that discussed in [22][23][24]. The Myers transition became very difficult to observe for µ > 6.5 as it occurred in a very narrow temperature interval. We have not successfully tracked it to very large µ and we suspect that, at least in the large µ region of the phase diagram, it is a finite matrix effect.
The second transition is observed in the Polyakov loop and is associated with the eigenvalue density of the gauge field, or holonomy, transitioning from gapped to ungapped (from covering a small interval of the unit circle to covering the entire circle). For µ ∼ 6.0 this occurs in a narrow temperature interval. However, as µ is lowered the Polyakov loop is not a monotonic function of temperature; it jumps upward at the Myers transition, where the system transitions to larger fuzzy sphere configurations. Each BPS ground state has its own typical value for the holonomy. In the Hamiltonian picture the holonomy implements the Gauss law constraint and is therefore sensitive to the density of accessible SU (N ) non-singlet states 5 . It plays an essential rôle when there are many non-singlets around the relevant energy. One would expect that different fuzzy sphere configurations have different densities of non-singlets, which would account for the behaviour of the holonomy, as it is effectively probing these.
We found that as µ was decreased the two transitions merge at about µ ∼ 3.0. This single transition curve then approaches the Padé prediction (2.3) and the gauge/gravity [10] prediction (2.4). We have not tracked the transition to µ less than 1.0. 6 We find strong evidence that the Polyakov loop transition coincides with the Myers transition for µ < 3.0. Also, our results are in qualitative agreement with the gravity predictions of [10]. However, we expect that better agreement over a larger range of µ would be obtained should a black hole solution dual to a general vacuum be available.
Another remarkable observation is that, at strong coupling, µ = 2.0, the Myers term seems to have a finite, non-zero value in the large-N extrapolation. If one considers a typical vacuum as a representation with m copies of the k dimensional representation, it yields at low temperatures Myers ∼ µ 3 3 4 N Tr(J r J r ) ∼ µ 3 (k 2 − 1). Such a configuration provides an example where this observable is non-zero and does not diverge, when k is finite. Thus the easiest way to realise configurations that have finite Myers term at large N is to have many copies of representations that fluctuate around a typical one of relatively small dimension. This should correspond to a state of five-branes [7], and in turn suggests that the Myers tranitions for µ < 3 can be regarded as a transition to a five-brane phase.
Our current study is clearly only the beginning. However, it indicates that an effort to construct such a black hole configuration in the dual gravitational theory, even a numerical one, would be very useful.
In the near future we plan to add D4-brane probes to the BMN model [34] in analogy with our studies [35,36]. This should give an alternative, more detailed probe, of the dual geometry.