Gauged And Ungauged: A Nonperturbative Test

We study the thermodynamics of the `ungauged' D0-brane matrix model by Monte Carlo simulation. Our results appear to be consistent with the conjecture by Maldacena and Milekhin.


Introduction
The meaning of gauge symmetry in holography is not immediately clear. Although the most well understood version is indeed gauge/gravity duality [1], the dual field theory may not have to be a gauge theory.
Gauge symmetry is a source of various headaches at the technical level as well. For example: • In the Hamiltonian formulation of quantum field theory, the algebra of gauge invariant operators has a complicated structure, and even state counting is difficult. This is one of the reasons that the Hamiltonian formulation of lattice gauge theory [2] is not very practical.
• The temporal component of the gauge field is often related to the sign problem.
• It is not easy to keep both gauge symmetry and supersymmetry on lattice.
For (0 + 1)-dimensional theories, the ungauged counterparts do not seem pathological. In the Hamiltonian language, the gauge singlet condition (Gauss's law) is omitted. In the pathintegral language, the gauge field A t is turned off. It does not lead to possible pathologies in higher dimensions-the breakdown of Lorentz symmetry, for example-either. If such an "ungauged theory" makes sense, it would avoid the problems listed above. A priori, however, it is not clear how much the gauged and ungauged theories differ. It is often said that the gauge symmetry is not important in the deconfining phase. a While it is qualitatively true, it is more subtle at the quantitative level; at least at high temperature, where the coupling constant is small, free energies of gauged and ungauged theories are different. b For example, nothing is known for the strongly coupled region near the transition to the confined phase; the two theories may or may not resemble one another. Below the transition temperature, the two theories are clearly different, in that the ungauged theory cannot be confining by definition. Still, the contribution from the non-singlet sector might be so small that the physics in the singlet sector is practically unaffected.
Recently Maldacena and Milekhin [7] considered this problem by taking the D0-brane matrix model [3,[8][9][10] as a concrete example. In terms of string theory, the gauge singlet sector describes closed strings. Gauge non-singlets naturally appear when open strings are allowed. As a string theory, such a setup seems to be fine. Then it would make sense to consider the duality between the ungauged theory and string theory. They proposed a reasonable dual gravity prescription, and made a striking conjecture: the difference between gauged and ungauged theories is exponentially suppressed at low temperature. For example the difference of the energy should scale as d adj N 2 C adj e −C adj /T , where d adj is a positive integer and C adj is an order one positive number which corresponds to the energy of the adjoint excitation. a The explanation is that "color degrees of freedom can be seen directly". This is somewhat misleading, because physical states are still gauge invariant. See e .g. [6] regarding this point. b In the case of the D0-brane matrix model which we will study in this paper, E ∼ 6N 2 T and E ∼ 6.75N 2 T for the gauged and ungauged theories, respectively. See Appendix B.
By doing so, the circumference of the Euclidean circle β is the inverse temperature: β = 1/T . The gamma matrices γ M αβ (M = 1, 2, · · · , 9) are the 16 × 16 left-handed part of the gamma matrices in (9 + 1)-dimensions. ψ α (α = 1, 2, · · · , 16) are N × N real fermionic matrices. This theory is the dimensional reduction of 4D N = 4 super Yang-Mills theory to (0 + 1)dimensions. We often set the 't Hooft coupling λ = g 2 Y M N to one, without losing generality. Equivalently, all dimensionful quantities are measured in units of the 't Hooft coupling; for example the temperature T actually refers to the dimensionless combinationT ≡ λ −1/3 T . It also means the energy scale is related to the strength of the interaction: low temperature (small T ) and strong coupling (large λ) are equivalent, in the sense thatT is small. In the same manner, long distance is strong coupling.
The action and partition function are given by The 'ungauged' theory is defined simply by dropping the gauge field A t , as In the Hamiltonian language, the ungauging procedure we just described is equivalent to removing the gauge singlet constraint.

Dual gravity descriptions
The dual gravity description of the gauged D0-brane matrix model near the 't Hooft large-N limit was proposed in Ref. [10]. c The dual is a black zero-brane consisting of N D0-branes and open strings connecting them. The internal energy E = ∂ ∂β (βF), where βF = − log Z(β), is identified with the energy of the black hole above extremality. At low temperature (strong coupling), the dual gravity calculation d provides us with where A 7.41 is an analytically calculable number, up to the α and g s corrections (higher order inT and 1/N 2 ). This conjectured duality has been confirmed numerically with good precision, by comparing direct numerical Monte Carlo results of the gauged D0-brane matrix model and the dual gravity calculations, including stringy corrections. e The duality has been tested for other quantities as well: the supersymmetric Polyakov loop [20] agrees with dual gravity calculations based on the minimal surface prescription [21], and correlation functions [22] agree with the calculation based on the generalized conformal symmetry [23] analogous to the GKPW relation [24,25].
In the Hamiltonian formulation, the Hilbert space of the gauged theory consists of gauge singlets, which are obtained by acting traces of products of scalars on the vacuum, such as Tr(X M 1X M 2X M 3X M 4 )|Vac . Such single trace operators are similar to Wilson loops in QCD, which are identified with QCD flux strings, and naturally leads to an intuitive interpretation as closed strings. It motivates us to identify microstates of black zero-brane with closed string states, as was speculated before the duality had been found [26,27].
Gauge non-singlets admit products of scalars without trace, introducing open strings. For example, (X M 1X M 2X M 3 · · ·X M L )|Vac (no trace!) can describe an open string consisting of L bits. In the dual gravity description, because there is no special point in the bulk where open strings can end, it is natural to assume the open strings ends at the boundary of the space, as the left picture in Fig. 1 [7].
Such a long open string can interact with itself at self-intersecting point; on the gauge theory side, the kinetic (electric) term describes such an interaction. (It is essentially the same as the self-interaction of the long closed string explained in Ref. [6]. The crucial point is that, although the interaction at each intersection is 1/N -suppressed, the growth of intersection points overcome such suppression.) Correspondingly, on the gravity side, a long open string can split into a long closed string and a simpler open string, as in the middle picture in Fig. 1. But then the open string can shrink toward the boundary, and the bight of the closed string can also shrink, like in the right picture of Fig. 1, so that the energy (which is roughly proportional to the length of the string) is minimized. In this way, Maldacena and Milekhin [7] conjectured that the dual gravitational system is a black zero-brane plus short open strings localized near the boundary. They have estimated the contribution of such short open strings, and concluded that the free energy F should increase by c It has also been proposed that this model can describe M-theory, in a region with much stronger coupling [8][9][10]. d reviewed, for example, in Ref. [16] e The large-N , continuum result is available in Ref. [11]. Numerical studies by several independent groups [12][13][14][15][16][17][18][19] obtained consistent results. where d adj is a positive integer and C adj is the energy of the adjoint excitation. The dots represent terms negligible at large N and small T such as those from higher representations.
if the conjecture is correct. See Ref. [7] for more precise arguments. Our goal in this paper is to test this conjecture by studying the gauge theory side.

Lattice regularization
We use the simulation code developed for Monte Carlo String/M-theory Collaboration, which is freely available to anybody [28]. The ungauged version is obtained simply by turning off the gauge field (and the associated Faddeev-Popov term) from the code for the gauged theory. f In order to make the lattice regularization simpler, we use a slightly different (though equivalent) form of the fermionic action, Here γ M (M = 1, · · · , 10), which are 16 × 16 upper-right block of the 10d gamma matrices Γ M . For the later convenience, we take g f The ungauged version is also available upon request to M. H. g Others are taken as follows: This model is obtained by dimensionally reducing the ten-dimensional N = 1 super Yang-Mills theory to one dimension. The index α of the fermionic matrices ψ α corresponds to the spinor index in ten dimensions, and ψ α is Majorana-Weyl in ten-dimensional sense.
For numerical efficiency, we take the static diagonal gauge, and add the associated Faddeev-Popov term, to the action. The numerical simulations are performed in the phase-quenched limit. The reader can refer to Ref. [11] for more details on this aspect.

'Naive' Regularization
We regularize the theory by introducing a lattice with N t sites. Our lattice action is where U = diag(e iα 1 /Nt , e iα 2 /Nt · · · , e iα N /Nt ), −π ≤ α i < π, and D ± acts on ψ as The action on the X M is the same: D + X M (t) = U X M (t + a)U † − X M (t). Other than the gauge fixing, this action is the same as the one used in [13].

Tree-level Improved Lattice Regularization
The forward and backward derivatives used in the naive lattice discretization are related to the covariant derivative in the continuum theory by We can reduce the discretization error by using The improved action is obtained by replacing D ± with D (imp.) ± . The Faddeev-Popov term remains unchanged.

How to measure the internal energy
In this study we calculate the internal energy in the gauge theory side for both the ungauged and the original matrix model. This quantity can be calculated by using [29] Note that the Faddeev-Popov term is not included in the right hand side. The derivation of (19) is as follows. By rescaling the fields and time, we can take the 'dimensionless coupling' to be λ = λβ 3 . In the lattice action, with this normalization, the βdependence appears only through λ , as an overall factor β −3 in front of the action. From this we obtain E = − 3 β S b + S f , up to an additive constant. Here S f = Tr 1 2 D · D −1 = const. (D is the Dirac operator defined by S f =ψDψ), regardless of the details of the action, as long as the fermionic part is a fermion bilinear. So E = − 3 β S b up to an additive constant. We take the constant so that E = 0 at T = 0 if supersymmetry is not broken. This constant can be obtained just by counting the fermionic and gauge modes. The "−1" in (19) is the constant mode of U(1), which is removed by hand in our calculation.
We also use (19) for the ungauged theory. Then, if the ground state is in the singlet sector, the energy should vanish at T = 0.

Numerical exercise: bosonic theory
In this section we study the bosonic analogue of the D0-brane matrix model, simply obtained by neglecting the fermionic degrees of freedom and their interactions. The model is described by the usual matrix model once the term S f in Eq. (8) is dropped. The conjecture by Maldacena and Milekhin is developed in the context of gauge theories with a gravity dual. Although the bosonic matrix model in this section has no known dual, we want to exemplify our numerical procedure in this numerically easier case, where complications due to fermionic contributions are absent.
The gauged bosonic theory has been studied as the high-temperature limit of D1-brane theory [30][31][32][33]. Unlike the supersymmetric theory, the gauged bosonic theory is confined at low temperature. For properly normalized quantities, such as E/N 2 or R 2 ≡ 1 N β dtTrX 2 M , the temperature independence, up to 1/N corrections, in the confining phase, expected from the Eguchi-Kawai equivalence [34], has been confirmed numerically [30][31][32]. The internal energy for this bosonic theory is related to 4λ F 2 is the potential energy, and (kinetic energy) = 2 × (potential energy) holds due to the virial theorem.
In the ungauged version of the bosonic theory, we expect an exponentially suppressed contribution from the nonsinglet sector, with Boltzmann weights ∼ e −C adj /T , on top of the constant part coming from the ground state. Properly accounting for these new Boltzmann weights when computing the expectation value of different observables, like the energy or R 2 , results in an exponentially small difference.
In order to verify the presence of the contribution from the nonsinglet sector, and that it is suppressed at small temperature, we numerically calculate observables for the regularized ungauged bosonic model at nine different temperatures 0.2 ≤T ≤ 0.6 with N = 6, 8 and 12 and L = 12, 16, 24, 32, 48 and 64. We take the continuum limit (L → ∞) and the large-N limit (N → ∞). This is the first time such limits have been studied systematically for this bosonic model and numerical details for the analysis can be found in Appendix A, together with summary tables. In the following we will show the rescaled quantities measured on the lattice asẼ/N 2 From Ref. [32] we know that the transition temperature isT ≈ 0.9 and we also know that the constant values forẼ/N 2 andR 2 are 0 = 6.695(5) and r 2 0 = 2.291(1) from N = 32 simulations at finite lattice spacing. In Fig. 2, we plotẼ/N 2 andR 2 as a function of 1/T for the ungauged theory at low temperatureT ≤ 0.6-a regime where the gauged theory is confining. Both observables O(T ) are shown together with a fit of the form where C corresponds toC adj ≡ λ −1/3 C adj , while A, shown as a dashed black line in the plots, represents theT = 0 value for the observable. ForẼ/N 2 we obtain A = 6.7118(29) which represents the internal energy of the ungauged theory at zero temperature. ForR 2 we get A = 2.29244(55). We can compare these results with the those obtained for the gauged theory in Ref. [32], despite the latter not being extrapolated to the continuum large-N limit. The agreement is well within the statistical accuracy of the data and it supports the conjecture by Maldacena and Milekhin that the difference of the two theories vanishes at T = 0. Furthermore, the exponent C in our fit is identified withC adj and hence should be independent of the observable used for the fit. We can confirm this expectation within our statistical accuracy: C = 2.043(76) forẼ/N 2 and C = 1.936(71) forR 2 . The value of B for E/N 2 is 20.0(2.9), whose error bar is large, but it could be an integer times C, as predicted by the conjecture. In Fig. 3 we plot the data forẼ/N 2 andR 2 after we have subtracted the corresponding values of A obtained from the fits. In turn, this amounts to calculate the difference between the ungauged and the gauged bosonic theory for the two aforementioned observables, and one can visually check the exponential falloff in the plot. The plots in Fig. 3 show the same fits as the ones reported in Fig. 2. Again, we remind the reader that there is no known gravity dual for the bosonic model we studied in this section, but the results are encouraging and show that a test of the conjecture with good precision is possible with lattice Monte Carlo methods.
the ungauged bosonic matrix model in the low temperature region, where the gauged counterpart is in the confining phase. Both behave as A+Be −C/T , where C 2.0. In the gauged theory, bothẼ/N 2 andR 2 are constant in the confining phase, and the values obtained in Ref. [32] agree with A from the ungauged theory with good precisions.
the ungauged bosonic matrix model in the low temperature region, after subtracting theT = 0 fitted value A (see text and Fig. 2), which represents the energy of the gauged theory in the confining phase. Both differences behave as a falling exponential Be −C/T , where C 2.0. The values of A extracted from the fit agree with the values obtained in Ref. [32] for the gauged theory.
We now move to test the conjecture in the full D0-brane model with fermions, where it is developed.

Taming the flat direction
The biggest obstacle in the simulation of the D0-brane matrix model is the flat directions of the eigenvalues of scalars X M [35,36] [12]. h Namely, the bound state of eigenvalues, which is dual to a black zero-brane, is merely metastable, and the eigenvalues (D0-branes) can be emitted and fly away. In order to measure the energy of the black zero-brane, we have to tame the flat direction, or in other words, simulate the system in the local minimum of the action corresponding to the metastable state. For that purpose we adopt the simplest procedure, which does not require any modification to the action: take the matrix size N sufficiently large [11,12,17,19,37]. The emission rate of the eigenvalues decreases as N becomes large and, for sufficiently large values of N , we can collect sufficiently many configurations of the metastable state before the onset of any instability.
However, the severeness of the instability is regularization dependent. In fact, bosons and fermions contribute to the attraction and the repulsion, respectively. Therefore, if the bosonic contributions dominate due to a regularization artifact (which would disappear only in the continuum limit), the system is more stable, and the instability sets in only as one gets closer to the continuum limit. This is the case for the regularization known as "momentum cutoff method" [12,14,16,19,20,22]. The opposite is true for our lattice regularization: the instability becomes milder as one approaches the continuum limit. With the regularization method we have used, we need to take the lattice size L to be sufficiently large-the exact value depends on N . We have also noticed, a posteriori, that the instability is more severe compared to the gauged theory, which has been studied in detail before [11].
The choices ofT , N and L for the ungauged BFSS theory are summarized in Appendix A. We have performed a simultaneous large-N and continuum extrapolation of the ungauged theory observables using different fit functions and different datasets to assess possible sources of systematic errors. All details are summarized in Appendix A.

Comparison between gauged and ungauged theories
In this section we show the simulation results for E/N 2 , F 2 and R 2 to test the conjecture by Maldacena and Milekhin. The results for E/N 2 , F 2 are numerically under better control, and appear to be consistent with the conjecture. It is hard to make a sharp statement on R 2 because the N -dependence is large, and this observable is more subject to large fluctuations related to the flat directions.

Energy
The difference of the energy should behave as where d adj is integer, if the conjecture is correct, and C adj is different from the one of the bosonic theory in Sec. 3 Table 1: Data for the internal energy of the ungauged and gauged theories in the large-N continuum limit.
For the ungauged theory we have data atT = 0.45, 0.5, 0.6, 0.7, 0.8, 0.9 and 1.0. At each temperature we have simulated various values of L and N and then performed a simultaneous large-N and continuum fit, including 1/N 2 corrections. For the gauged theory, on the other hand, we know from the results in Ref. [11] that there is a good description of the energy as a function of the temperature based on the gauge/gravity duality:Ẽ(T ) = a 0T 14/5 + a 1T 23/5 + a 2T 29/5 . We fit the parameters a 0 , a 1 and a 2 to the continuum large-N data in Ref. [11] (considering onlyT < 1.0 data) and obtain the values ofẼ gauged at the same temperatures ofẼ ungauged . The uncertainties are propagated from the full covariance matrix of the fitted parameters. Table 1 summarizes the data used in the following analysis. The energyẼ ungauged is obtained from a large-N continuum fit excluding L = 8 points and including 1/N 2 corrections (see Appendix A for more details on the fits).
In Fig. 4, we plot both the fitted functional form forẼ gauged and raw continuum large-N data forẼ ungauged in the regionT ≤ 1.0. We also plot ∆Ẽ ≡Ẽ ungauged −Ẽ gauged in the same Fig. 4. This energy difference is consistent with ∆E = d adj C adj e −C adj /T , where d adj = 2 and C adj = λ −1/3 C adj 1.0. We perform several fits on the energy difference data atT ≤ 0.8. First we fit the functional form expected from the conjecture with d adj = 2 and obtain: where A is compatible withC adj 1.0 and the quality of the fit is very good. We also use the same functional form, but fixC adj = 1.0, leading to: where we now get the amplitude compatible with d adj = 2. These fits look promising, but they are not conclusive. In fact, if we fit the same data with a fitting ansatz where the amplitude and the exponent are free parameters we obtain: where the best fit parameters are compatible with the previous results, but with larger uncertainties. These fits are compared in Fig. 5.  Figure 4: [Upper] Large-N and continuumẼ gauged ≡ λ −1/3 E gauged andẼ ungauged ≡ λ −1/3 E ungauged as functions ofT ≡ λ −1/3 T . The former is a fit to data in Ref. [11], and it is compared with the leading order supergravity solution, shown as a dashed black line.
[Lower] Large-N and continuum ∆Ẽ ≡Ẽ ungauged −Ẽ gauged vsT . Two ungauged energy are shown to asses systematic errors: one which includes L = 8 points in the continuum extrapolation and one that does not. AtT = 0.45, the ungauged simulations have not been performed at L = 8 because of numerical instabilities related to the flat direction. We also show some representative functional forms to guide the eye.
The values of such quantities can be different in the singlet and nonsinglet sectors. i However, according to the Maldacena-Milekhin conjecture, the contribution from the nonsinglet sector is suppressed by a Boltzmann factor e −C adj /T relative to the singlet sector, due to the energy gap C adj . Hence the expectation value should be dominated by the singlet sector; we expect and whereF 2 ≡ λ −4/3 F 2 andR 2 ≡ λ −2/3 R 2 , respectively (In Sec. 3, we have numerically observed the same phenomenon happening in the purely bosonic theory, with a different value for C adj , which depends on the theory). Unlike the case of ∆E, there is no particular prediction for the overall coefficients c F 2 and c R 2 .
The details of the analyses forF 2 andR 2 are given in Appendix A. To summarize, we choose to look at two different sets of data, reported in Tab. 2 and Tab. 3 for the ungauged and the gauged theories. For the gauged theory, we use the raw data results reported in Ref. [11] and we extrapolate them to the continuum and large-N limit. One dataset is extrapolated to the N → ∞ limit, while the other one is not, using instead the largest N at each temperature as a proxy to estimate what systematic errors we would make by neglecting finite-N corrections.
We have extrapolated results in both theories for six temperaturesT ≥ 0.5 j at which we compute the differences that enter Eq. (25) and Eq. (26).  Table 2: Data for the potential energy termF 2 and the extent of spaceR 2 of the ungauged and gauged theories in the large-N continuum limit.
In the upper panel of Fig. 6 and Fig. 7, we have plotted ∆F 2 and ∆R 2 , respectively, as a function of the inverse temperature. We have also plotted the relative difference ∆F 2 /F 2 , whereF 2 is just the average value ofF 2 between the ungauged and the gauged theory (similarly forR 2 ). Despite large error bars at 1/T ≥ 2.0, the trend that ∆R 2 and ∆F 2 go to i The contribution from short open strings hovering at the boundary would be small but not parametrically suppressed.
j Unlike E/N 2 , we do not know reasonable fit ansätze for F 2 and R 2 . Therefore we do not try to obtain the values atT = 0.45 by fitting the temperature dependence from other values ofT . We estimate the difference only atT ≥ 0.5.  Table 3: Data for the potential energy termF 2 and the extent of spaceR 2 of the ungauged and gauged theories at the largest N available for eachT and in the continuum limit. N = 32 for all temperatures in the ungauged theory, and for the lowest temperature in the gauged theory, but N = 24 forT ≥ 0.6. zero asT → 0 can be seen. The relative difference is on the order of a few percent everywhere atT ≤ 1.0, which suggests typical matrix configurations important in the path integral are very close. Note that the relative difference of the energy, ∆Ẽ/Ē, is larger. However, it does not seem to be an immediate problem, becauseĒ, the average internal energy between the ungauged and gauge theory, is zero at zero temperature due to the cancellation between bosonic and fermionic contributions. If we take the ratio of ∆Ẽ and another natural energy scale, e.g. the zero-point energy in the bosonic theory, it can be equally small. Moreover, it is interesting to stress that, at very high temperatures, the ratios above are identical to that of the bosonic theory, which is ∆Ẽ/Ē = ∆F 2 /F 2 ∼ 12% (see Appendix B), and that at our simulated temperatures we are already seeing a dramatic reduction, with ratios down to the 2% level.
For ∆F 2 we can attempt to do a fit to the data, similarly to what we did for the internal energy, while, unfortunately, it is hard to make any quantitative statement for ∆R 2 , where the error bars are larger. The exponential decay of ∆F 2 is consistent with e −C adj /T withC adj 1.0, as we can see in Fig. 8. The fits we performed on ∆F 2 are similar to the ones we did for the internal energy in Eq. (22) to Eq. (24) and they all have reduced χ 2 close to one. Fig. 8 shows all three fits on the aforementioned two different datasets independently. Although we don't find any theoretical reason that the overall coefficient is integer, numerically it appears to be rather close to 3, while it was closer to 2 for the internal energy. The exponential decay of ∆Ẽ has the same scale:C adj 1.0. We have seen this observable-independence in the bosonic case as well, and we remark that this is expected from the conjecture.  [Upper] ∆R 2 as functions of 1/T together with some representative functional forms to guide the eye (as done for ∆F 2 , but with ten times smaller coefficients). The datasets shown are in the large-N continuum limit or at our largest N value.
[Lower] ∆R 2 /R 2 as functions of 1/T , whereR 2 is the average value between the ungauged and the gauged theory. There is a clear decreasing trend as the temperature decreases, but the statistical uncertainties become too large at the smallest temperatures to say anything conclusive about theT = 0 limit. [Lower] The same exponential fits are shown using a different dataset where the largest-N results are used at each temperature (reported in Tab. 3). In both cases, only five rightmost data points are included in the fits, corresponding toT ≤ 0.9. The fits are all compatible with the data and among each other, with the two-parameter fit having the largest uncertainties. Although we don't find any theoretical reason that the overall coefficient is integer, numerically it appears to be rather close to 3.

Conclusion and Discussion
In this paper we have tested a recent conjecture by Maldacena and Milekhin [7] by studying the ungauged D0-brane matrix model with numerical lattice Monte Carlo simulations. Our results shown in Sec. 4 appear to be consistent with the conjecture, given our statistical accuracy. More detailed tests of this conjecture -with higher precision, at lower temperatures, and for more observables -are straightforward in principle and should be performed in the future to increase confidence in the conjecture. While analytic methods, such as supersymmetric localization [38], may work k for certain problems, Monte Carlo simulation can be applied to more generic situations. It is important to pursue various approaches which can lead to complementary results. If the conjectured duality between gravity and the ungauged theory is correct, there are several interesting directions. On the gravity side of the story, in addition to the issues raised already in Ref. [7], whether the ungauged theory fits into the supermembrane interpretation [9], in which the gauge transformation is identified with the area-preserving diffeomorphism, would be an important question to address. It would also be interesting to see how the basic results like D0-brane scattering [39][40][41] may or may not be modified. Another interesting question is whether such ungauging procedure can make sense for other theories which do not have dual gravity interpretations; hopefully the ungauging procedure can solve some technical issues already mentioned in the introduction.

A Monte Carlo data analysis
This appendix contains some technical details about the numerical analysis of our lattice Monte Carlo data. We describe our data analysis workflow and we report various details about the extrapolations to the continuum limit and to the large-N limit. Summary tables for the observables measured in the ungauged bosonic matrix model and the ungauged full matrix model are reported in this section.

A.1 Bosonic theory
We have studied the ungauged version of the bosonic matrix model by using the 'naive' regularization of the continuum action: The simulation parametersT , N and L are summarized in Tab. 7 together with the accumulated Monte Carlo samples and the average value of two observables, the internal energỹ E/N 2 and the extent of spaceR 2 . The Monte Carlo samples are retained after 5000 samples of thermalization, and are binned in independent bins that are five times larger than the autocorrelation time of each observable, also reported in Tab. 7. Different observable can have different autocorrelation functions and this is observed in the different number of independent bins that are available forẼ/N 2 andR 2 : the latter observables tend to have less bins because it has longer autocorrelation times. For the simultaneous continuum and large-N extrapolations of both our observables, we have used the fitting ansatz defined in Ref. [11], where O(T , N, L) is the expectation value of an observable directly measured on the lattice at fixedT , N and L. We performed a 4-parameter fit with o 00 , o 01 , o 02 and o 10 while the other coefficients are set to zero. The fits are done on the observablesF 2 andR 2 and representative results are shown in Fig. 9 and Fig. 10. The plots show both the fitted data points and the fitted curve-the different curves are constant-N slices of the best-fit surface as a function of 1/L, while the black curve is the large-N limit. The final results for the internal energỹ E/N 2 and the extent of spaceR 2 are summarized in Tab. 4 and Tab. 5.

A.2 Full theory with fermions
For the gauged theory, we take the data from Ref. [11]. The interested reader can find the tables in Appendix B of Ref. [11]. For the observableẼ/N 2 , the large-N and continuum values for the gauged theory were presented in Ref. [11] and were fitted tõ E/N 2 (T ) = a 0T 14/5 + a 1T 23/5 + a 2T 29/5 .  Figure 9: Continuum and large-N limit extrapolation ofF 2 = 4 3Ẽ /N 2 for the bosonic ungauged matrix model atT = 0.45. The title of the plot summarizes the large-N continuum limit value ofF 2 , the χ 2 with corresponding degrees of freedom, and the Q value of the fit.  Figure 10: Continuum and large-N limit extrapolation ofR 2 for the bosonic ungauged matrix model atT = 0.2. The title of the plot summarizes the large-N continuum limit value ofR 2 , the χ 2 with corresponding degrees of freedom, and the Q value of the fit.   We perform the same fit on the data atT ≤ 0.9 and reproduce the results of Ref. [11] within statistical accuracy: The full covariance matrix for the fitted parameters is reported in Tab. 6 and it is used to obtain accurate error bars onẼ/N 2 (T ) for the same values ofT used in the ungauged theory.
For the other observables,F 2 andR 2 , we perform simultaneous large-N and continuum fits using the data in Ref. [11]. However, a functional dependence on the temperature can not be fitted for these observables, since there is no corresponding expectation from the gauge/gravity duality.
to the data in Tab. 8.
The analysis of the Monte Carlo data for the observables of the ungauged theory proceeds in a similar manner to the one done in the gauged theory [11]: • for each N , L andT value, we look at the Monte Carlo history of the observables and remove the configurations where we find instabilities due to the flat direction; this rarely happen because we simulate large enough N at eachT • for each N , L andT value, we remove 5000 configurations from the beginning of the simulation to reduce thermalization effects • for each N , L andT value, we compute the integrated autocorrelation time τ using the Madras-Sokal algorithm and we make sure that the sample size is sufficient to include at least 10 times this autocorrelation. If this is not possible, we discard the whole ensemble because of limited statistics.
• for each remaining N , L andT value, we compute the average and the standard deviation coming from binned measurements of size 5 × τ .
• for eachT , we extrapolate the results from the previous step to {N, L} = {∞, ∞} with a simultaneous large-N and continuum limit.
• for eachT , we also extrapolate the results from the previous step to L = ∞ at fixed N , for comparison.
Similar challenges to the gauged theory are present in the ungauged theory, like the treatment of long autocorrelations in the Monte Carlo samples and the sensitivity of the continuum and large-N extrapolations to the fitting ansatz and data set. However, there are also some differences that we highlight in the following.
For several parameters {N, L,T }, we have simulated multiple streams of Markov Chain Monte Carlo. This is a common procedure in order to increase the statistics with multiple independent copies of Monte Carlo samples (sometimes called 'replicas').
The first step in the analysis of the multiple streams amounts to combining them as if they were separate independent 'experiments', using a weighted average (or constant fit). Although this is a straightforward procedure, it can happen that some streams have central values that are not compatible with the weighted average within two standard deviations. This is reflected by a poor 'p-value' (or equivalently a large χ 2 in this case).
An example of this happens at N = 16, L = 64 andT = 0.6, where one of the five independent streams is ∼ 2.5 standard deviations from the average of the streams when   looking at the energy observable. This is shown in the upper panel of Fig. 11 and it is actually the only occurrence we encountered in the analysis. We looked for possible sources of systematic errors in the numerical simulations of these particular streams, like a small acceptance rate, a limited number of samples, unusually large fluctuations, etc... We did not find anything suspicious other than the possibility of long autocorrelations, much longer than the currently accumulated samples, which could result in underestimating the uncertainty. Therefore, we keep this point at N = 16, L = 64 andT = 0.6 in the analysis, and we interpret it as a statistical fluctuation. We have also checked that a fit to the continuum at N = 16 andT = 0.6 with a quadratic function in 1/L is not dramatically influenced by including or removing this L = 64 point, as can be seen in the lower panel of Fig. 11. Another issue is that the N -dependence can be very mild and often the data is not sufficient to resolve the 1/N 2 corrections unambiguously. Let us focus on the energy observable where this is more dramatic. Similarly to the gauged theory case, we use the fitting ansatz where we only include {e 00 , e 01 , e 02 , e 10 } terms, where e 00 is the continuum large-N contribution, e 10 is the leading continuum large-N correction, and the other coefficients characterize discretization artifacts and finite-N corrections. The L −j terms are needed to model the data at different lattice sizes L up to the coarsest lattices (L = 8). The dependence on N is only modeled by the N −2 term, but the corresponding coefficient e 10 is often not constrained by the data at N = 16, 24 and 32. When this is the case, a simple fit without 1/N 2 corrections, or equivalently with an average of the datasets at different N values, can be used to represent the data. In other words, we can use results at fixed N as if they were already in the N = ∞ limit. Fig. 12 shows the results for the coefficient e 00 , the energy in the large-N and continuum limitẼ, as a function of temperature, for four different extrapolations. It compares fits with and without the leading 1/N 2 corrections when all the datapoints are included (all values of L and N ) and fits to data where L = 8 points have been removed. In all cases we note that the quality of the fit, measured by the reduced χ 2 (χ 2 /dof), is always around or below 2, being a little worse for the higher temperatures. At the same time, including or excluding 1/N 2 corrections gives compatible results. Removing points with the coarsest lattice spacing helps estimating additional systematic uncertainties of our extrapolation models. For the fitting models considered above, there is a discrepancy of a bit more than one standard deviation when two different datasets are used at higher temperaturesT ≥ 0.8. Our most conservative approach is to consider the results with the largest statistical errors, the ones obtained by fitting the datasets without L = 8 and with the 1/N 2 corrections. They are compatible within two standard deviations with the results from all the other fits.
We also noticed that observables likeF 2 andR 2 have relative statistical errors that are orders of magnitude smaller than the ones for theẼ/N 2 observable. This turns out to be an issue when doing least-square extrapolations because highly-precise data can artificially drive the reduced χ 2 of the fit to values much larger than one: it does not mean that the functional form used in the fit is necessarily bad in reproducing the data, but only that the final errors on the fit parameters is under-estimated. Of course, this latter scenario has to be tested, for example by manually increasing the uncertainty on the individual data points and tracking how the fit parameters change. In Fig. 13 we show the continuum extrapolation of F 2 data at N = 16 andT = 0.7 in the ungauged model. The three panels show data points with statistical uncertainties increased by a factor of 1, 2 and 4 (clockwise), but they are still too small to be seen on the plots. At the same time, each panel shows the continuum extrapolation with a second order polynomial in 1/L performed on the full data set and on the data set without the L = 8 point. The two fits in each panel are always compatible with each other and the functional form goes through the data points quite nicely, indicating that the fit can be used to reliably extrapolate the points to the continuum. However, the reduced χ 2 is always much larger than one, unless the individual point error bars are increased by a factor of 4. The four-fold error bar increase is reflected by a corresponding four-fold increase in the fit parameters error bars. We adopt this error-inflating procedure forF 2 andR 2 for both the ungauged and the gauged theory, but for the latter we only need a factor of 2 increase. The simultaneous large-N and continuum extrapolated results forF 2 at eachT are shown in Fig. 14 for the ungauged theory and in Fig. 15 for the gauged theory. Fits to the datasets with and without L = 8 are performed with and without 1/N 2 corrections. We notice that for the ungauged theories, the four different fits give compatible results, although fits with L = 8 have a somewhat worse fit quality. On the other hand, for the gauged theory, fits with L = 8 give significantly different results at low temperaturesT ≤ 0.7. Given the worse quality of the fits in those cases, we therefore only consider fits without L = 8 when we construct differences with the ungauged theory at each temperature.
A different approach to the N = ∞ limit that has been used in past studies is to consider continuum results at fixed N , and if N is large enough one can neglect 1/N 2 corrections. We plot the fixed-N continuum results forF 2 in the ungauged and gauged theory at each  Figure 13: Continuum limit at fixed N = 16 andT = 0.7 for the ungaugedF 2 with and without L = 8 points. Each panel, starting from the top, shows fits to points with inflated error bars by a factor of 1, 2 and 4, respectively. The fits in each panel are consistent, showing that a quadratic extrapolation function works well, even for L = 8, but the quality of the fits is bad. This indicates that our errors are underestimated. temperature in Fig. 16. When multiple values of N are present at the same temperature in the continuum, we see negligible finite-N corrections, as all points are statistically compatible. This justifies the approach of considering continuum results for the largest N at each temperature as a proxy for the N = ∞ result.
In the case ofR 2 , however, the N dependence is larger.R 2 is, by definition, very sensitive to the flat direction and, even when simulations do not show samples affected by the instability, the value ofR 2 tends to be larger at smaller N . We can see this trend in Fig. 17, where for some temperatures, the largest value of N can still be different than the value obtained by extrapolating at N = ∞, and has a systematic shift to larger values. In order to assess the error one is making by not considering a full simultaneous large-N and continuum limit, we keep two data sets in the analysis of the conjecture: the largest N dataset and the dataset obtained by including 1/N 2 corrections and without the coarsest lattice spacing L = 8.

B High temperature behavior
At high temperature, the D0-brane matrix model and its bosonic analogue both reduce to the classical matrix model. The kinetic and potential energies are related by the virial theorem as (kinetic energy) = 2 × (potential energy), and hence the total energy E satisfies E = 3 2 × (kinetic energy) = 3 × (potential energy). The kinetic energy in the classical theory   simply counts the number of degrees of freedom: (Strictly speaking, we need to take into account the conservation laws coming from the constant shift of eigenvalues and SO(9) rotations. They affect subleading terms with respect to 1/N 2 .) Hence E/N 2 = 6T and 6.75T for gauged and ungauged theories, respectively. E and F 2 are related by E/N 2 = 3 4λ F 2 , because (potential energy) = N 2 4λ F 2 .