Thermal phase transition in Yang-Mills matrix model

We study the bosonic matrix model obtained as the high-temperature limit of two-dimensional maximally supersymmetric SU($N$) Yang-Mills theory. So far, no consensus about the order of the deconfinement transition in this theory has been reached and this hinders progress in understanding the nature of the black hole/black string topology change from the gauge/gravity duality perspective. On the one hand, previous works considered the deconfinement transition consistent with two transitions which are of second and third order. On the other hand, evidence for a first order transition was put forward more recently. We perform high-statistics lattice Monte Carlo simulations at large $N$ and small lattice spacing to establish that the transition is really of first order. Our findings flag a warning that the required large-$N$ and continuum limit might not have been reached in earlier publications, and that was the source of the discrepancy. Moreover, our detailed results confirm the existence of a new partially deconfined phase which describes non-uniform black strings via the gauge/gravity duality. This phase exhibits universal features already predicted in quantum field theory.


Introduction
Bosonic matrix models in one dimension have been studied in various contexts. Despite their simple structure, they display a rich non-trivial phase diagram that can be accessed by analytical methods only in certain limiting cases. An important motivation to study these theories arises from their connections to supersymmetric Yang-Mills theories (SYM) in one and two dimensions, which have various applications to quantum gravity via the gauge/gravity duality [1].
The Euclidean action of the gauged bosonic U(N ) matrix model is 1 where λ = g 2 YM N is the 't Hooft coupling, β is the inverse temperature, I, J = 1, 2, · · · , d with d = D − 1, and X I are N × N hermitian matrices. The covariant derivative D t is defined by D t X I = ∂ t X I − i[A t , X I ], where A t is the gauge field. We will mainly consider the case D = 10, since it is the bosonic version of the BFSS model [2]. In addition to gauge symmetry, this theory has the U(1) center symmetry. An order parameter associated with the center symmetry is the Polyakov loop, where P denotes the path ordering. The Polyakov loop transforms as P → e iθ P under the U(1) transformation. Another important symmetry is SO(d) which rotates the X I scalars as d-dimensional vectors. In this paper, we will focus on the breaking of the center symmetry.
While the full supersymmetric BFSS model in the large-N limit has a dual gravity description [3], no weakly-curved gravity dual is known for its bosonic part. Still, this model can provide insights into gravity, as we will see shortly. For comparison and to study the large-D limit, we will also investigate the case of D = 26.
This theory is deconfined at high temperature, and confined at low temperature for any D ≥ 3. Based on large-D and large-N analytical techniques [4] and, partially, on numerical Monte Carlo results at finite N [5], it was initially believed that the deconfinement transition is not a first order one but rather there are two phase transitions, one of second order and one of third order, in close proximity. More recently however, evidence was presented that there may be only one transition of first order [6]. We study numerically the phase diagram of this model in the large-N limit in order to confront the numerical results with analytical predictions, in particular taking into account the continuum limit which was not previously investigated. In this paper, we will provide robust numerical evidence that for D = 10 there is only one deconfinement transition and it is of first order in the large-N limit.
There are several reasons to be interested in the order of the phase transitions for this model. Firstly, let us point out the connection to the topology change between a black hole and a black string [7,8]. In order to understand it, note that the bosonic matrix model in Eq. (1.1) is the high-temperature limit of two-dimensional maximally supersymmetric Yang-Mills (SYM) theory compactified on the spatial circle S 1 . The bosonic matrix model is obtained by shrinking the temporal circle of the two-dimensional SYM theory to a point. Consequently, what is called the "temporal circle" in the bosonic matrix model actually corresponds to the "spatial circle" in the compactified two-dimensional SYM theory. The phase transition we study in this paper is the remnant of center symmetry breaking/restoration along the spatial S 1 in this higher-dimensional theory, which can be regarded as the black hole/black string topology change in the R 1,8 × S 1 spacetime of the dual gravitational description [9].
In the studies of the black hole/black string transition based on general relativity, a reliable analysis of the topology change is out of reach at the moment due to the unavoidable curvature singularity. This problem can be avoided by using the dual gauge theory, namely the 2d maximal SYM theory mentioned above. The 2d maximal SYM theory contains information about the stringy corrections and can teach us how the singularity can be resolved. At low temperatures and large N , where the stringy corrections are small, the dual gravity description predicts a first order transition when the radius of S 1 is small [9]. However, the details of how the topology change takes place is out of reach in the gravitational approach. At intermediate and high temperatures, the only practical tool to gain insights into the dynamics of the theory is a numerical approach.
Lattice approaches (for example see e.g. Refs. [10][11][12][13] for recent reviews) are applicable to 2d maximal SYM [14][15][16][17]. Numerical approaches are computationally expensive and it is hard to take N sufficiently large at this stage. A theory that is more computationally tractable is the one-dimensional bosonic matrix model Eq. (1.1). Depending on the nature of the transition in the bosonic matrix model, the dual gravity prediction of a first order transition may survive at high temperature (left of figure 1), or it may fail and the transition will split into two, a third order and a second order one as depicted in the right panel of figure 1. In the former case, we might be able to obtain valuable intuition into the lowtemperature region, which is not easy to access with currently available computational resources, from the high-temperature region, which is relatively easier to study. In the latter case, the "non-uniform string" becomes stable due to stringy effects.  [9,18]. The vertical axis is the inverse temperature, β = T −1 and the horizontal axis is the radius R of the spatial circle. The phase transition is related to the breaking of center symmetry along the spatial circle, analogous to confinement and deconfinement. Black hole, uniform black string and non-uniform black string phases on the gravity side correspond to the "completely" deconfined, confined and "partially" deconfined phases, respectively. (See Sec. 2.1 for the meaning of "complete" and "partial" deconfinement.) Our numerical simulations tell us that the left figure is more likely to be true.
Similar problems are of interest in the context of the application of holography to QCD. In fact, this strategy was also used by Witten [19] to study 4d Yang-Mills (YM) as the high-temperature limit of 5d maximal SYM, for which the dual gravity analysis is tractable. Depending on whether the qualitative features of the phase transition change or not, the scope of the holographic approximation can change. (See e.g. Ref. [20] for detailed discussions regarding this point.) Another example is SYM on R 3 ×S 1 deformed by the gaugino mass [21]. The deconfinement transition on a small three-sphere was studied in the same manner [22,23].
The nature of the deconfinement transitions in various gauge theories fits to the framework of partial deconfinement [24] that we will explain in Sec. 2. The matrix model provides us with the simplest setup to study it in a non-perturbative way. A good understanding of these transitions might shed light on the deconfinement transition in gauge theories or the microscopic nature of the QCD crossover from the hadronic phase to the quark-gluonplasma phase.
A related motivation is provided by the black hole information problem. When the phase transition is of first order, there is an unstable phase [24] which is analogous to the Schwarzschild black hole with negative specific heat [25] in the standard setup of the AdS/CFT duality [19]. Therefore, the identification of simple models exhibiting such behavior is the first step towards investigating this issue with a detailed numerical study.
The theory we consider in this paper (Eq. (1.1)) admits an analytic treatment in terms of a large-D expansion [4], which predicts two transitions, one of second and one of third order. Numerical Monte Carlo studies at D = 10 [5], up to N = 32, looked consistent with this large-D analysis, while more recent studies [6] for the same value of N seemed to contradict with that analysis.
In the following we investigate some reasons that motivate our investigation. Firstly, finite-N corrections become more relevant near the critical point, in analogy to finitevolume effects of a statistical system near criticality. The numerical data of Ref. [5] did not include a dedicated study of finite-N corrections, showing only two values of N . In particular, the transition was studied at N = 32 and we will show that this value is not sufficiently large to reveal the order of the transition by looking at |P | (T ). Other papers [26,27] which numerically backed up the observation in Ref. [5] were not able to study the large-N limit either. Evidence for a first order transition was presented in Ref. [6], along with a discussion of why the large-D analysis may fail to correctly predict the order of the transition for small D. The first order signal shown in that study was however not completely clear. Moreover, the study did not consider discretisation effects which may affect the nature of the transition by shuffling the order of temperatures if there are multiple transitions as expected from the large-D analysis (see figure 3). This was due to the limited computational resources available to deal with the considerable increase in numerical cost for larger N and smaller lattice spacing.
Secondly, D = 10 may not be sufficiently large to make contact to the large-D expansion and therefore to trust the analytical expectation. In order to get a rough intuition, let us consider an analogous gravity problem, the black hole/black string transition in general relativity on R D−1 × S 1 . In this case, the transition is first order at D ≤ 13 [28] and large enough D means D > 13. This suggests that it might be dangerous to trust the large-D analysis at D = 10. Note that the order of the transition is particularly sensitive to the value of D. The large-D analysis seems to be more reliable all the way down to small D for other properties of the theory like the approximate location of the transition [29], see also the discussion in [6]. It is therefore important to study in detail the large-N limit at fixed D = 10.
In this paper, we present numerical evidence that the transition in the D = 10 model is of first order. Therefore, in figure 1, the left diagrams are more likely to be true. (Although, strictly speaking, our findings also allow for the possibility that the transition is of higher order at intermediate values of R.) We study D = 26 as well. Somewhat surprisingly, we have observed signals consistent with a first order transition, which suggests the large-D approximation is not precise even at D = 26. The organization of this paper is as follows. We start with explaining theoretical expectations in Sec. 2, in order to define the strategy of our numerical simulations. Then, we will show the results of the simulations and their implications in Sec. 3. Sec. 5 is devoted to discussions and conclusions.

Theoretical expectations
In this Section, we summarize the theoretically expected features of the theory, which will be compared to numerical data and used to determine the order of the phase transition. We introduce three different patterns characterizing the deconfinement transition [24] and explain how they can be distinguished using numerical simulations. These patterns can be explained by introducing the concept of partial deconfinement [24,30,31]. Since our main goal is to identify the order of the transition, we only include relevant details about partial deconfinement and we refer the interested reader to the papers cited above for a comprehensive review. The word 'partial deconfinement' [30] is used to characterize a phase where a subset of the U(N ) group, which we denote as U(M ) in the following, is deconfined. A matrix representation of this concept is illustrated in figure 2. The 'completely deconfined' and 'confined' phases correspond to, as expected, M = N and M = 0, respectively. Intuitively, the reason why partial deconfinement takes place is simpler to understand when working in the large-N limit: complete deconfinement requires an energy of order N 2 , and hence, if the energy is much smaller (say E ∼ N 2 ), only a part of the color degrees of freedom can be excited, U(M ) with M ∼ √ N . The initial motivation of partial deconfinement was to understand the gauge theory description of the Schwarzschild black hole with negative specific heat via the gauge/gravity duality, closely following a very similar mechanism based on partial Higgs-ing [32,33] in gauge theories. There have been various consistency checks for several theories, both at weak coupling and strong coupling [24,30], and explicit demonstrations based on state counting are available for several theories [34].

Partial deconfinement and possible phase structures
In the large-N limit, quite generally, the deconfinement transition in gauge theories was classified in three types [24] according to the picture of canonical ensemble and the concept of partial deconfinement introduced above. Given that there are three phases, we will introduce two temperatures T 1 and T 2 , the first separating the 'completely confined' and the 'partially deconfined' regions, and the second separating the partially deconfined and the completely deconfined regions. The phase diagram is represented by its absolute value |P | as a function of the temperature. In figure 3, the three rows represent the following three distinct possibilities for the order of the phase transition: • First order with hysteresis. There is a local maximum of the free energy corresponding to the partially deconfined phase separating two minima, the completely deconfined phase and the confined phase. A hysteresis sets in at T 2 ≤ T ≤ T 1 . In the microcanonical ensemble, when the volume is sufficiently large, the confined and completely deconfined phases can occupy most of the space, and the partially deconfined phase appears at the interface. In a matrix model, the partially deconfined phase is stable in the microcanonical ensemble because there are no spatial dimensions and a separation in volume cannot take place.
• First order without hysteresis. At the transition temperature T = T 1 = T 2 , there is a Hagedorn string with degenerate free energy [22,23]. It corresponds to the partially deconfined phase denoted by the orange line. When the volume is sufficiently large, different vacua can appear at different locations, but this is only realized in systems with spatial dimensions.
• Two transitions of second and third orders. There are three stable phases with T 1 ≤ T 2 : even the partially deconfined phase is stable, both in canonical and microcanonical ensemble.
A convenient way to distinguish the three different phases (confined, partially deconfined and completely deconfined) is to look at the distribution of the phases of the Polyakov loop ρ(θ) [24]. Note that this is just one of the characterizations of the phases. In fact, partial deconfinement does not necessarily require center symmetry, and the Polyakov loop is not necessarily an order parameter for some theories with partial deconfinement [34]. By definition, θ is distributed between +π and −π. In the confining phase, the distribution is uniform at large N : In order, the three rows correspond to three deconfinement transitions: first order with hysteresis (T 1 > T 2 ), first order without hysteresis (T 1 = T 2 ), and two transitions of second and third orders (T 1 < T 2 ). Dashed lines represent unstable phases, while solid lines represent stable or metastable phases. The blue, orange and red lines correspond to the confined, partially deconfined and completely deconfined phases. Right: Corresponding free energies for the three types of scenarios. See text for details. Here we use |P | = P .
Here, the subscript c stands for 'confined'. We will use p and d for 'partially deconfined' and 'deconfined', respectively. The transition between the partially and completely deconfined phases is the Gross-Witten-Wadia (GWW) transition [35,36], as found in [24]. Namely, in the partially deconfined phase, the distribution is not uniform, but also not gapped, i.e. ρ p (θ) > 0 everywhere in −π ≤ θ ≤ π, while in the completely deconfined phase ρ d is gapped. It is natural to expect [24] where ρ GWW (θ) is the distribution at the GWW transition point. In many cases, the distribution takes a simple form: and Here, we have fixed the U(1) phase factor using P = |P |. In the partially deconfined phase, A = 0 at T = T 1 and A = 1 at T = T 2 (the GWW transition point), and 0 < A < 1 otherwise. In the transition region, we can interpret A = M N . In the completely deconfined phase, A ≥ 1, and M = N regardless of the value of A.
For the bosonic matrix model, the form of ρ d (θ) in Eq. (2.4) has been confirmed by previous studies in a wide region of parameter space [5,24]. In Ref. [5], the form of ρ p (θ) in Eq. (2.3) has also been observed and used as the evidence for the absence of a first order transition. However, as we will see, this is an artifact of the finite-N correction described below in Sec. 2.3.
For our purposes it is important to note, for example from figure 3, that in the case of the first order transition with hysteresis the deconfined phase becomes unstable below P = 1 2 . In other words, if the value of P is stable at 0 < P < 1 2 , the transition cannot be of first order.

Large-D, large-N analysis
An analytic approach to understanding the thermal phase transition in Yang-Mills matrix models with action Eq. (1.1) has been developed in Ref. [4]. The key technical tool employed was an expansion of the functional integral around a non-trivial saddle point in the limit of a large number of matrices d (where d = D − 1). Around the saddle point and in this limit, fluctuations are suppressed by powers of 1/d. A priori, it is not obvious what value of d is large enough to justify this expansion, although hints may be obtained from gravity computations, as noted in Sec. 1. Our numerical results suggest that d = 9 is definitely too small, and even at d = 25 our simulations show a qualitatively different behavior. The analytic agreement with numerical studies [5] mentioned in Ref. [4] might be attributed to fact that N ≤ 32 is too small to reveal the nature of the large-N transition. On the other hand, simulations at N = 32 up to d = 15 were interpreted as consistent with a first order transition [6], in disagreement with the analytical approach.
The results of Ref. [4] can be summarized as follows. Throughout the analysis, a gauge is adopted where A t is time-independent and diagonal. At low temperatures, the eigenvalue distribution ρ(θ) of the Polyakov loop P becomes constant as N → ∞. As a consequence, the Polyakov loop vanishes. This behavior persists up to a temperature T 1 , where a second order phase transition happens. The large N eigenvalue distribution is now given by 2 and |P | continuously increases form 0 to 1/2 as T is increased to a second critical temper- , this is the same distribution as ρ p given by Eq. (2.3). At T = T 2 , the third order Gross-Witten-Wadia type [35,36] phase transition occurs after which |P | increases further but with a smaller slope. The eigenvalue distribution becomes gapped at T = T 2 and eventually approaches a single delta function at very high temperatures.
The predictions for the critical temperatures including the first 1/d corrections at large N are given by [4] We note that for d → ∞, ∆T c := T 2 − T 1 → 0, i.e. the two transitions occur in a very narrow temperature regime, making quantitative numerical checks difficult. For the cases considered in this paper, one obtains the values in Tab. 1.

Finite-N effects
By definition, there is no spatial extent in a matrix model. Hence, the thermodynamic limit has to be realized as the large-N limit. The finite-N corrections can obscure the phase transitions, and it is important to know what kind of corrections are expected, in order to determine the nature of the phase transition numerically. The situation is easier to understand when the large-N transition is not of first order. Because the large-N limit is the thermodynamic limit, the transition becomes sharper gradually as N increases.
Some caution is required when the large-N transition is of first order. The free energy is of order N 2 also in this case, and the fluctuations about the minima are 1/N -suppressed. Therefore, at sufficiently large N , the distributions of the observables such as |P | and E/N 2 should have a two-peak structure since tunneling between them is suppressed as e −const.×N 2 . However, when N is not sufficiently large, the tunneling probability might be so large that the two peaks merge and become one single wide peak. In this way, the first order nature of the transition gets completely hidden. Even worse, because the confining and deconfining phases give ρ c (θ) = 1 2π and ρ d (θ) = 1 2π (1 + cos θ), the mixture of two phases -say of the confining phase with probability 1 − A and of the deconfining phase with probability Agives ρ(θ) = 1 2π (1 + A cos θ), which is exactly the same as ρ p of the partially deconfined phase. Therefore, observing a distribution compatible with ρ(θ) = 1 2π (1 + A cos θ) is not enough to claim that the transition is not of first order.

The order of the phase transition and the large-D limit
In this Section we investigate numerically the phase transition of the bosonic matrix model. We are in particular interested in the D = 10 case, where the large-D approximation might no longer be applicable. The smooth behavior of the order parameter observed at small N turns into a signal for a transition only in the large-N limit. In this limit, two possible scenarios can be discriminated: 1. Two distinct transitions become visible: The lower one will be of second order and the higher one will be indicated by an expectation value of the Polyakov loop |P | = 1 2 . Due to the continuous behavior of the Polyakov loop and the small difference of the transition temperatures, the two transitions can only be distinguished at large enough N .
2. One first order transition appears: In case of the first order transition, the signal will be quite similar to the one of a second order transition. Starting from the low temperature confined phase, there will be a broadening of the minimum of the constraint effective potential and an increase of the susceptibility. However, before the actual second order transition (Hagedorn transition) occurs, a second minimum of the effective potential induces a first order transition.
Because of the nature of this phenomenon, it is difficult to decide about the order of the transition based on the susceptibility, but a two peak signal in the histogram of the order parameter or a hysteresis effect allows a clear distinction of the two scenarios. Therefore, we concentrate on the appearance of these signals at large N in the following investigations.
Our lattice action includes the bosonic part and the gauge fixing part of the BFSS lattice action defined in Ref. [37], which was also used in the numerical study of Ref. [38]. We consider scans in T at a fixed number of lattice points L and matrix size N . The action is parameterized in such a way that the lattice spacing scales as 1/L and quantities like the temperature are all provided in units of the 't Hooft coupling λ, which is set to unity in the numerical simulations. If we assume the scenario of two separate transitions, the increase of the Polyakov loop would have a finite width ∆T c , with a rise starting at T 1 and stopping at T 2 . The slope of |P | at T 1 would increase with N , but the point with |P | = 1 2 at T 2 would remain at a finite distance ∆T c from T 1 . Based on this assumptions, we can deduce a rough estimate of ∆T c and its large-N extrapolation from the width of the transition. At N = 32, which has been the maximal N in previous numerical studies, the obtained width is still compatible with the large-D prediction ∆T c = 0.022. However, the extrapolation towards the large-N limit does not support a finite ∆T c required by the scenario of two separate transitions. In order to substantiate these findings, a more detailed analysis is necessary.
As pointed out above, the best way to discriminate the two scenarios is the two-state signal or a hysteresis of the order parameter. The two-state signal can be deduced from a two-peak structure in the histograms of the order parameter that persists in the large-N limit. In addition, we consider possible effects of the finite lattice spacing by comparing histograms from simulations with a different number of lattice points L. In case of a first order transition, the separation of the two peaks becomes more pronounced at large N since the tunneling rate between the two states is exponentially suppressed with N 2 in this limit. The tunneling is also visible in Monte-Carlo time, but this effect is not unambiguous since it has an algorithm dependence.
For N ≤ 32 we can not observe a clear two-peak signal in the histogram (although hints of a two-peak structure are visible, see the appendix), but at N = 48 a two-state signal can be observed that becomes more pronounced at N = 64, see figure 5. There is evidence for the existence of two phases, one with small |P | and the other with |P | ≈ 0.5. Consequently, a hysteresis of the order parameter is found, see figure 6. We have investigated three different lattice spacings at N = 64 in order to show that the effect persists in the continuum limit. The complete Monte Carlo history for N = 64, L = 24 is shown in figure 22 in the appendix for the transition temperature T = 0.885. It displays repeated tunneling events between the two phases.
Overall, we conclude from these data that there is a first order transition at D = 10 and thus a considerable qualitative deviation from the large-D expansion. However, as figure 7 shows, a continuum extrapolation of T c (which is largely insensitive to N ≥ 48) would lie within T 1 and T 2 obtained from the next to leading order expansion in large D. This indicates that at least the location of the transition is correctly estimated by the analytic formulae in Sec. 2.2.
Consistent with earlier numerical studies [5], we have seen that the signal at smaller N might indicate a different scenario. We have also found that the peak of the low temperature phase at small |P | gets significantly broader around the transition temperature. This is consistent with the expected Hagedorn instability at T = T 1 . Due to this phenomenon, we observe an increase of the susceptibility towards a peak when approaching the critical temperature from below. Consequently, the first order transition occurs just before a second order transition manifests itself, as anticipated in Sec. 2.1. One would expect that the transition changes from first to second order at larger D and the picture becomes more consistent with the large-D analysis.

D = 26 model
The analysis of the previous Section revealed that the large-D expansion fails to describe the order of the phase transition for D = 10. In this Section, we repeat the above investigation for D = 26 which might be large enough for the expansion to be qualitatively applicable, but also small enough to sufficiently limit the computational costs. It turns out that the simulation results are consistent with a first order transition. Figure 8 shows the dependence of the order parameter |P |(T ) near the transition. A naive large-N extrapolation with fixed lattice size L = 24 locates a possible transition window to be between T = 0.873 and T = 0.874. This separation is much smaller than

First order signals from other observables
It is useful to study other quantities in order to provide further evidence for the existence or absence of a first order transition. Let us consider E If the transition is of first order, the two-peak signal should be visible for these observables as well.
In figure 11 and 12 we plot the distribution of E N 2 and R 2 for D = 10. Like for the order parameter, we observe a first order signal from a clear two-peak structure. Note that the locations of the two peaks do not move significantly as a function of N . This is consistent with the theoretical expectation that the maximum of the free energy (partially deconfined phase), rather than the minima (completely deconfined or confined phases), moves.
The first order signal from E/N 2 and R 2 is less pronounced at D = 26, as shown in figure 13 for N = 64. We have only presented R 2 since the histograms of E/N 2 are similar. For N = 48, no clear signal can be observed, but for N = 64, a two-peak structure is visible. It is the clearest for L = 24, where we collected higher statistics. The results  for L = 48 suggest that the two-peak signal persists in the continuum limit. This is also supported by the Monte Carlo history shown in figure 23 of the Appendix. The separation of the phases at D = 26 is less pronounced compared to D = 10, which might be consistent with the transition developing towards a combination of higher order transitions at some critical D > 26.

Partial deconfinement
In the previous Section, we have determined the order of the deconfinement phase transition. Let us go one step further and study details of the phase transition. In particular, we use our numerical data to test the partial deconfinement [24,30,31], reviewed in Sec. 2.
Partial deconfinement is the proposal that the deconfinement transition happens gradually so that at an intermediate stage only M < N color degrees of freedom are deconfined (figure 2). We will now derive some consequences of this assumption for correlations between different observables and the eigenvalue distribution of the Polyakov loop, and test them against our numerical data.
As we have seen, the phase transition is of first order. Then the size of the deconfined sector M should change with temperature as visualized in figure 14. At finite N , there are non-negligible fluctuations around the saddles, and hence, the size of the deconfined sector M fluctuates configuration-by-configuration during the Monte Carlo simulation. Below, we will relate the value of M to observables (Polyakov loop P , energy E and the extent of space R 2 ). This leads to nontrivial relations between observables, which can be used for the consistency check of the partial deconfinement proposal. Then we will confirm those relations numerically.
Let us assume that Eq.   [22,23] seems to hold (within errors) also in the strongly coupled regime. Since Eq. (4.1) holds, it is reasonable to assume partial deconfinement with M N = 2|P |. Consequently the partial deconfinement proposal provides further nontrivial predictions: 1. At fixed temperature, the energy per excited degree of freedom is fixed, and hence the energy above the ground state is proportional to M 2 . Therefore, as a function of |P | and T , we expect where ε 0 represents the zero-point energy. We expect that this relation is precise at large N , where the fluctuation about the planar limit is suppressed. Our data confirm these relations with good precision. We provide the data for D = 10 since the results for D = 26 are similar.
As a first check, we have plotted the density distribution of the configurations in (E, |P |)-plane in figure 16. We have introduced bins in (E, |P |)-plane, counted the number of configurations in each bin, and repeated the same for R 2 as well. As expected, the distribution becomes sharper as N increases due to suppression of fluctuations in the planar limit. This confirms, at least at a qualitative level, that (4.2) and (4.3) are valid.
In order to confirm (4.2) and (4.3) quantitatively, we separate the configurations into bins according to their |P | value. In that way we obtain expectation values E (|P |, T ) of the energy E for fixed |P | and in the same way R 2 (|P |, T ). Figure 17 shows the plot of In Ref. [37], the correlation between energy and Polyakov loop in the D0-brane matrix model has been studied in a parameter region with |P | 0.7, i.e. in the completely deconfined phase. It was found that the energy and Polyakov loop are not correlated at all. In contrast, we find that even for T = 1.5 where |P | ≈ 0.9, E/N 2 and R 2 are correlated, see figure 18. We also studied the low temperature phase at T = 0.5, where there appears to be no correlation.

Conclusion and future directions
In this paper, we have studied the nature of deconfinement in the bosonic Yang-Mills matrix model (1.1). We have concluded that the transition is of first order for D = 10. By interpreting this model as the high-temperature limit of 2d maximal super Yang-Mills (SYM) it is natural to conclude that the phase diagram of 2d maximal SYM is like the left panel of figure 1. Via the gauge/gravity duality, we can rephrase this finding [9] to say that the α -corrections (which become more important at higher temperatures) do not alter the order of the phase transition between black hole and black string in canonical ensemble. We have also observed a first order transition for the D = 26 case, which provides further insights about the validity of the large-D expansion. In comparison with the results of Ref. [4] and in agreement with [6], we conjecture that there is a critical dimension D c > 26 with first order transition at all D < D c , see figure 21. The large-D picture with two separate transitions at T 1 and T 2 applies only for D ≥ D c . The validity of this conjecture is not easy to confirm since the temperature separation ∆T decreases with increasing D and higher accuracy is needed to determine the existence of separate transitions. There are various directions for future follow-up studies. The Berenstein-Maldacena-Nastase (BMN) matrix model [39], which is a one-parameter deformation of the D0-brane matrix model [2,40,41] with the flux parameter µ, would provide us with a numerically tractable setup with a weakly-coupled gravity dual [42]. There are two possible phase diagrams for the BMN matrix model, as shown in figure 19, and in particular the phase structure at µ = 0 is still unclear. For an unambiguous determination of the phase structure, simulations at large enough N and sufficiently low temperature are essential. Despite the extensive numerical studies performed for the D0-brane matrix model, see Refs. [43][44][45] and related studies, this has not been achieved so far, and the dual gravity interpretation [2,3,40] of the results might be affected. The bosonic version of this model is an instructive exercise to test the simulation methods. Before we had obtained our numerical data, we had two possibilities for its phase diagram in mind, which are illustrated in figure 20. Our numerical simulation indicate that the picture on the left without a stable intermediate phase is realized at µ = 0.
The supersymmetric matrix models [2,3,39,40] are dual to black zero-brane in type IIA string theory. In order to understand how gauge theory degrees of freedom describe quantum gravity, it is important to determine which of the two possibilities shown in figure 19 is realized in this case. This requires an investigation of the very low temperature region, as previously studied in Ref. [50], but also at sufficiently large N . So far, the largest value of N in the literature is N = 32 [37], and given the lesson from this paper, it might still be too small at low temperatures.
In this paper, we have not studied the details of the unstable partially deconfined phase, or equivalently the maximum of the free energy. This phase, which is not important in the importance sampling for the canonical ensemble, actually contains very important T μ deconfine confine μ deconfine confine 'partially' deconfine T Figure 19: Two kinds of conjectured phase diagrams of the BMN matrix model in the canonical ensemble. The vertical axis is the temperature T . The large-µ region admits perturbative calculation [46,47] and the transition is found to be of first order. The smallµ region has been studied by using the dual gravity description [42], but the order of the transition has not been established. See Refs. [48,49] for lattice simulations along this direction. information of the theory, because this is the phase connecting the completely deconfined phase and the confined phase in the microcanonical ensemble. It should be possible to study the property of this phase by determining the location of the 'dip' between the two peaks, and picking up the configurations from there. A detailed numerical study of this phase would be useful for understanding black hole evaporation in the context of holography. When the bosonic matrix model is interpreted as the high-temperature limit of 2d maximal SYM, a comparison with the dual gravity calculation at low temperature would be very interesting, because the low and high temperature regions resemble each other at least at the qualitative level as we have shown in this paper. The P -vs-T plot for various d, from large (left) to small (right). Blue, orange and red lines are confined, partially deconfined and completely deconfined phases, respectively. At d = ∞, T 1 = T 2 , and at large but finite d, T 1 < T 2 [4]. We have observed a first order transition, and hence T 1 > T 2 for D ≤ 26.

B Data for N = 32
To compare with the results of Ref. [6], we provide histograms for |P | and R 2 at N = 32, L = 24 in figure 24. Similarly to Ref. [6], the histogram of |P | features a slight shoulder around |P | = 0.25 − 0.3 for the transition temperature T = 0.885. This may be taken as an indication for a first order transition. However, as emphasized in the main text, a detailed analysis including larger N and a continuum extrapolation is needed for a definite conclusion.