Nonperturbative test of the Maldacena-Milekhin conjecture for the BMN matrix model

We test a conjecture by Maldacena and Milekhin for the ungauged version of the Berenstein-Maldacena-Nastase (BMN) matrix model by lattice Monte Carlo simulation. The numerical results reproduce the perturbative and gravity results in the limit of large and small flux parameter, respectively, and are consistent with the conjecture.


Introduction
The holographic principle claims that quantum gravity can be described by a dual nongravitational theory. AdS/CFT duality [1] provides us with concrete realizations of the holographic principle. But in fact, neither AdS nor CFT is crucial. AdS/CFT duality is a special case of gauge/gravity duality that admits non-AdS/non-CFT duality. As yet other examples of the duality emerge, it became natural to ask if the gauge-singlet constraint is crucial on the QFT side. Indeed, this is something that gained more popularity since the appearance of the Sachdev-Ye-Kitaev (SYK) model [2,3]. The SYK model is described by quantum mechanics of N fermions with random coupling between them and it does not have a gauge symmetry. Despite the lack of gauge symmetry, fermion bilinears that look like "singlets" appear to have natural dual descriptions in an emergent two-dimensional spacetime [4]. The same story holds for the Gurau-Witten tensor models [4,5] in which the gauge-singlet constraint can be either imposed or not.
Maldacena and Milekhin [6] discussed whether the gauge-singlet constraint is necessary for more traditional types of gauge/gravity duality. They considered the D0-matrix model [7,8] and its ungauged version, and conjectured that the difference between the gauged model (with singlet constraint) and the ungauged model (without singlet constraint) is exponentially small at low temperature. In particular, they conjectured that the same gravity dual describes the low-energy dynamics of the gauged and ungauged models. Numerical simulation of the D0-matrix model [9] provided results consistent with the conjecture.
There is a natural generalization of the D0-brane matrix model keeping maximal supersymmetry, called the BMN matrix model [10]. Maldacena and Milekhin considered the ungauged version of the BMN matrix model as well. In this paper, we test their conjecture for the BMN matrix model by employing numerical methods. In addition to the original motivation coming from the Gurau-Witten tensor model, yet another motivation, in this case, comes from quantum simulations [11]. The gauge-invariant Hilbert space consisting only of singlet states is rather complicated, and merely writing down the orthonormal basis is already a difficult task. This problem can be avoided by introducing the extended Hilbert space that contains non-singlet states. A potential worry arises if the additional non-singlet degrees of freedom introduce other unexpected technical issues. Specifically, if many light modes emerge, they can easily be excited and lead to a large error. If the conjecture by Maldacena and Milekhin is correct, such light modes do not exist and hence the use of the extended Hilbert space can be rather straightforward.
The BMN matrix model contains the 't Hooft coupling constant λ, the deformation parameter µ, and the effective dimensionless coupling constant is g eff = λ/µ 3 . This, together with temperature, spans the phase diagram of the model, describing different regimes. For works regarding the phase diagram see e.g [12][13][14]. In particular, we will be interested in the strong coupling regime where a dual description by weakly-coupled gravity exists.
The paper is organized as follows: in Sec. 2.1 we are defining the models, while in Sec. 2.2 we describe the gravity duals. In Sec. 2.3 we give more details for the gauged and ungauged versions of the models and we discuss the conjecture. The numerical analysis and the main results for the BMN model [10] follow in Sec. 3. We discuss observables obtained from the partition function and their difference between the gauged and ungauged versions. In Sec. 4 we conclude, while supplementary results are contained in the appendices.

Maldacena-Milekhin conjecture
In this section, we review the conjecture made by Maldacena and Milekhin [6]. Firstly, we specify the theories under consideration in Sec. 2.1. Then, in Sec. 2.2, we introduce the dual gravity description. The details of the conjecture are described in Sec. 2.3.
Both X M and ψ α are in the adjoint representation of the U (N ) gauge group, and the covariant derivative D t acts on them as D t X M = ∂ t X M − i[A t , X M ] and D t ψ α = ∂ t ψ α − i[A t , ψ α ]. The action is given by (2.1) The equations of motion for the gauge field A t give rise to the Gauss constraint In the operator formalism, the gauge-singlet constraint on physical states emerges due to the integration over A t . See Sec. 2.1.4.

BMN Matrix Model
The plane-wave deformed theory, which is called the BMN matrix model [10], is given by 1 The extra terms appearing in the action are mass terms for bosons, fermions and interaction terms. We can take γ 123 as 2 which simplifies the fermionic mass term to (see also Appendix D.1) In addition, these new terms result in a new class of vacua labelled by representations of the SU (2) group. In other words, matrices that minimize the potential in addition to the trivial ones (i.e, X i = 0 = X a = ψ α ) can be written in the form ψ α = 0, X a = 0 for a = 4, · · · , 9, X i = µJ i for i = 1, 2, 3, (2.7) where J i are the generators of SU(2). In the limit µ → 0, the deformation terms vanish and one expects the above model to converge to the BFSS model. This however assumes that there is no phase transition between the models and indeed evidence until now supports this idea. Note also that the singlet constraint (2.2) is not affected by the deformation. One can construct an effective, dimensionless coupling constant via we can also define another dimensionless quantity These g eff control the regimes of validity of the two descriptions which can be put into the following language • Supergravity region: g Note that this r appearing in (2.9) controls the regimes of the duality in a similar spirit with the holographic spatial dimension appearing in the usual AdS/CF T (see Fig. 1). In particular, as we will see in the gravity analysis this r corresponds to energy scales in supergravity but let us briefly clarify this point. It will turn out (see e.g Eq. (2.16)) that r controls the size of the eight-sphere S 8 in the 't Hooft large N limit. Small r means small energies deep in the bulk (even inside the black hole) and very large r means near the "boundary" where the matrix model perturbation is performed (see again Fig. 1). The trustable region for supergravity is discussed below Eq. (2.21).
Being in the canonical ensemble, we study (thermodynamic) features of the models at a specific temperature, which results in specific energy of the system via E = E(T ). In this paper, we are mainly interested in the strong coupling limit of the model, which according to (2.8), (2.9), and (2.20) results in low energies and via the canonical ensemble at low temperatures. The temperature in the matrix models sets the circumference (β = 1/T ) of the Euclidean circle on which we put our models and is also connected with the Hawking temperature on the gravity side. We will comment more on this when we discuss the gravity dual picture in (2.18). In both theories, since conformal symmetry is absent, the 't Hooft coupling λ = g 2 Y M N can be set to 1 by proper rescaling of time t and matrices. In other words, all dimensionful quantities can be made dimensionless by multiplying with appropriate powers of λ (e.g, τ = T /λ 1/3 ). In simulations we set λ = 1. Figure 1: A pictorial representation indicating the gravity and matrix model side. The parameter that controls the different regimes is the effective, dimensionless couplings g (µ) eff = λ µ 3 and g (r) eff = λ r 3 . Near the boundary the supergravity description is not valid while one can apply matrix perturbation theory with expansion parameters g The supergravity solution can be trusted up to r λ 1/3 . This figure is based on [6].

Ungauged matrix model
If we turn-off the gauge field A t in the BFSS or BMN matrix model, we obtain the ungauged matrix model. Specifically, we define the ungauged action S ungauged as (2.10) The thermal partition functions of gauged and ungauged models are defined by Here, the time direction is Wick-rotated to the Euclidean signature and compactified to circumference β = T −1 . We impose the periodic boundary condition for bosonic fields X and A t , and antiperiodic boundary condition for fermionic fields ψ.

Gauge-singlet constraint
Let us discuss the meaning of the ungauging in terms of quantum states in the Hilbert space. By using the gauge-invariant Hilbert space H inv and the extended Hilbert space H ext , thermal partition functions are written as (see e.g., Appendix A.2 in Ref. [16]) dgTr Hext ĝe −βĤ (2.13) = Tr H inv e −βĤ (2.14) and Z ungauged (β) = Tr Hext e −βĤ .
(2. 15) In (2.13), V SU(N ) is the volume of the SU(N ) gauge group, SU(N ) dg is the integral with the Haar measure, andĝ is the operator acting on the Hilbert space as the SU(N ) transformation corresponding to the group element g ∈ SU(N ). By construction, 1 V SU(N ) SU(N ) dgĝ acts as the projector from H ext to H inv . In fact, (2.13) is directly related to the path-integral formulation with gauge field A t described by the canonical partition function (2.11). The projection operator tells us that we should count gaugeequivalent states only once, and hence, we obtain (2.14) from (2.13). The operatorĝ is the counterpart of the Polyakov loop in the path-integral formalism. Integration over SU(N ) is the remnant of the path integral with respect to A t .

The gravity duals
Gravity dual of gauged matrix models Let us start with the 't Hooft large-N limit of the gauged D0-brane matrix model without deformation. We take the effective dimensionless coupling constant λT −3 to be large (which is equivalent to low temperature). The BFSS gravity dual at strong coupling is believed [17] to be the charged black zero brane in type IIA supergravity formed by N coincident D0-branes. The metric is given as There is a killing horizon at g tt = 0 that sets the temperature of the black IIA brane to Knowing the temperature one can pursue a thermodynamic analysis and compare with the relevant quantities of the matrix model [17][18][19][20]. In particular the energy of the system is E 7.41N 2 λ −3/5 T For the BFSS model the dimensionless effective coupling is constructed as 20) in addition to g (r) eff , such that the supergravity region is given at g (E) eff 1, e.g at low energies and via (2.18) at low temperatures. One could use g (T ) eff = λ T 3 as well. Comparison between the BFSS matrix model and the black zero-brane was explored numerically using Monte-Carlo simulations for the internal energy of the theory accessing in this way in a non-perturbative fashion gauge/gravity duality and introducing α corrections to eq. (2.18) (see Refs. [21,22] for the first simulations and Ref. [23] for large-N and continuum limit). The trustable supergravity region demands the effective curvature of the eight-sphere to be small, which is given by the inverse sphere radius (R) as As long as the effective curvature κ is small, we can trust the supergravity description (2.16) which is the case for r λ 1/3 (see Fig. 1)). In addition, we have to ensure that the dilaton on the horizon is small N 10/7 . Considering always the large N limit the latter condition is always satisfied and combining all the conditions we have (2.23)

Gravity dual of ungauged matrix models
The idea behind gauging a symmetry or not in this particular example might become more intuitive if we present it in the language of string theory. The gauge invariant (physical) states are singlets. These states correspond to closed strings and are constructed by acting with combinations of matrix operators on the vacuum state upon taking the trace over the gauge group indices 3 , i.e |physical = Tr X I 1X I 2 · · ·X I L |vacuum .

(2.24)
On the other hand, the ungauged model allows some room for non-singlets without a trace |physical =X I 1X I 2 · · ·X I L |vacuum . (2.25) Even though the singlet sectors of both theories are the same, in the ungauged model there is a new sector that hosts non-singlets. The latter can intuitively be realized as an arbitrary long open string made out of L bits whose endpoints can reach the boundary. Let us now discuss the continuum picture, e.g L → ∞. In the bulk, this configuration of long strings is described by the gravity dual of non-supersymmetric Wilson loops and the difference with the usual supersymmetric Wilson loop is the fact that in the latter case the string obeys Neumann boundary conditions while in the former Dirichlet boundary conditions [6]. This means that the tip of the non-singlet string can freely 3 Here the gauge group is SU (N ) and by taking the trace we mean summing over the SU (N ) indices This is different from the trace over the Hilbert space, Tr Hinv or Tr Hext . move in the bulk and reach even to the boundary. To be able to compare with gravity, we should ask what would be a natural cut-off such that we could approximate the energy of this massive string with supergravity. This is dictated by the validity of (2.16). Recalling that [r] = (energy) 1 we may demand that E min ∼ λ 1/3 . This indeed would be a natural cut-off because when we have a pair of a string and an anti-string, we can arbitrarily lower their energy (length) by placing both its endpoints on the boundary. On the other hand, we cannot approximate its energy from the supergravity side in this case because the latter is not valid at the boundary. Hence, a natural non-zero value of minimal energy calculable on the gravity side is with C adj being a number of order one. The idea then is that the fate of non-singlet adjoint strings is to end on the boundary minimizing in this way its energy. Therefore, contributions from the energies of non-singlet states should be negligible in the large N and low-temperature limit and the gravitational dual of the non-singlet strings in this regime is the same as that of the singlets in (2.16). This is the conjecture made in [6] which we focus on in what follows.

Gauged vs ungauged: the conjecture
The temporal component of gauge field A t is not dynamical. Its role is to impose the gauge-singlet constraint. Matrix models do not have spatial gauge field components A x , A y , · · · , and hence, the only effect of a gauge symmetry is to impose the gauge-singlet constraint. The ungauged version of the model does not obey the singlet constraint. The SU (N ) symmetry is treated as a global symmetry of the system. In [6], Maldacena and Milekhin considered the BFSS matrix model (µ = 0) and BMN matrix model (µ > 0), and claimed that gauged and ungauged versions are essentially the same at large N and strong coupling, in the sense that the contribution of the non-singlet sector in the partition function is negligible in the large N limit. The partition functions of the two models are given as (2.11), (2.13), (2.14) for the gauged model, and (2.12), (2.15) for the ungauged model.
Let us denote the difference of the gauged and ungauged free energies in the large N limit as is a function that depends on the parameter regimes of the model 4 . The factor N 2 comes from the 't Hooft counting.
Ref. [6] discussed the free energy difference in two different regimes and focused on the strong coupling regime. The BMN model has two parameters that control the regimes of the system, that is µ and T .
The high temperature and weak coupling regime In this limit, the difference of the free energies is This is the weak coupling regime at high temperatures and large µ. At these temperatures, we do not have a bulk dual so we will not be interested in this regime.
The low-temperature and strong coupling regime Under the presence of a gravitational dual the free energy difference is This limit is a bit subtle. Note that for the BMN model it is not enough to consider just small temperatures, but one also has to consider the µ → 0 limit to compare with gravity. The reason is that even though for finite µ there is a gravitational dual description it is given by a deformed geometry [15]. From (2.29) we get the difference of energies E = ∂(βF ) to lowest order in temperature and µ = 0. The consensus built in [6] and [9] is that n adj is the degeneracy of the lightest mode and C adj its energy. This result seems now to be understood for the µ = 0 case. Indeed at µ = 0, the factor n adj is an O(1) integer, and C adj is an O(1) positive number. Results of numerical simulation at µ = 0 [9] are consistent with n adj = 2 and C adj 1.

The low-temperature and weak coupling regime
The BMN model can be also weakly coupled at low temperatures, in contrast to the BFSS model. This is due to the effective coupling (2.8). In the large µ region the BMN model admits a perturbative analysis [24][25][26] with perturbation parameter (2.8). The spectrum is discussed in Appendix D.1 and the lightest mode is the one created by B † a acting on |0 . Whether we take the trace or not corresponds to the decision between the gauged or ungauged theory. In the latter case, there is no trace and the lightest mode is created by B † a |0 with energy (D. 13). This is what we called C adj above and we have six of them because we have six harmonic oscillators (one for each direction in the SO(6) part) which gives the degeneracy. Also, in the low-temperature region, we can still use an exponential ansatz and we expect the perturbative result of the lightest mode to be given by This is what we expect in the perturbative, low-temperature regime for the lightest mode of the theory. In fact, if we wish to study heavier modes or if we wish to be more precise we should consider the full U (1) sector given by Table 4) but as we discuss later on and we show in Fig. 6 at very low temperatures, the difference between the contribution of the lightest mode and that of the full U (1) sector is negligible (at finite µ).

Numerical analysis
In this section, we summarize the numerical analysis. We extrapolate the results to large N and continuum limit to eliminate lattice artifacts and finite N corrections.

Lattice regularization
Below, we explain the details of the lattice regularization. The action is the same as the one used in Ref. [23], except that also the deformation terms are added (see also [12]).

Gauge fixing
The action of the gauged BMN matrix model given in (2.4) is invariant under the SU(N ) gauge transformation. For numerical efficiency, we took the static diagonal gauge, Associated with this gauge fixing, we added the Faddeev-Popov term to the action.

Lattice action
We regularized the gauge-fixed continuum theory by introducing a lattice with L sites and spacing a. The time parameter takes the values t = a, 2a, · · · , La. Breaking the action into the bosonic part (S b ), the fermionic part (S f ), the Faddeev-Popov term S F.P. and the mass deformation parts (∆S b and ∆S f ), the respective lattice action is (3.7) Here, U = diag(e iα 1 /L , e iα 2 /L · · · , e iα N /L ), −π ≤ α i < π. The Faddeev-Popov term S F.P. is given in (3.2). For the ungauged theory, we set U to identity and remove the Faddeev-Popov term.

Simulation strategy
The BMN matrix model has a few nice features that make the numerical simulation easier than the BFSS matrix model. A technical challenge for the latter is the existence of flat directions. To tame them, we have to take N very large. In the BMN matrix model, the flat directions are lifted due to the mass term in the flux deformation. Therefore, we can study the BMN matrix model at relatively small values of N . Furthermore, the condition number of the Dirac operator decreases as µ becomes large. This makes simulations more tractable in this regime.
Although the flat directions are absent in the BMN matrix model, there is a somewhat related issue associated with the existence of multiple vacua. In this paper, we study physics around the trivial background (X i = 0 = X a , ψ α = 0). For small µ and T , tunnelings between the trivial background and fuzzy spheres (X a = 0 = ψ α , X i = µJ i ) can take place frequently. The potential barrier between them depends on µ and N . More precisely, by investigating the stability of a minimum of the potential ∂V ∂X = 0 and relating X ∼ r one can find that the barrier between the two backgrounds scales as ∼ µ 4 N 4 (see Appendix A). Therefore, for very small µ and fixed N , the trivial background configuration can tunnel to a fuzzy sphere background and a distinction between the two is not possible. On the other hand, at relatively large values of N , the simulation remains in the trivial background.
The fluctuation of each matrix entry is roughly given by 5 0.6N −1/2 . Demanding that the radius of the maximum fuzzy sphere given by r max = µ s(s + 1) with spin s = (N − 1)/2 is less than the matrix fluctuations gives a constraint on values of µ for which we cannot distinguish a fuzzy sphere configuration and matrix entry fluctuations at finite N . For the values of µ we have used this constraint is always satisfied so we do not have to worry about it. It is furthermore essential in our simulations to consider a large enough number of lattice points L to avoid lattice artefacts. Hence, the strategy we used is the following: • For fixed µ we perform a series of simulations at different temperatures as well as varying N and L checking that we always stay in the trivial background of the confined phase.
• We extrapolate to the large N and continuum limit of the theory for our observables.
• We calculate the difference between gauged and ungauged observables such as the energy and those defined in (3.12) and (3.13), for different temperatures at fixed µ and then fit exponentials of the form (2.29), (3.14), (3.15) respectively.
• Finally, we take the limit µ → 0 to cross-check the results with Ref. [9].
In the past, there were two possible kinds of conjectured phase diagrams of the gauged BMN matrix model at N = ∞ (Fig. 2). The large-µ region admits perturbative 5 The value 0.6N −1/2 is the standard deviation. The order-one factor can be determined numerically, from the expectation value of 9 I=1 TrX 2 I , which is approximately 3.5N [23]. There are 9N 2 matrix entries and hence 3.5N calculations [27,28] and the transition is found to be of first order. The small-µ region has been studied by using the dual gravity description [15] and recently via Monte-Carlo simulations [12]. A first-order transition has been established both analytically and numerically, while in addition a surprising possibility to study aspects of M-theory, like the Schwarzschild black hole has been established. Numerical simulation applies to the intermediate-µ region as well, and the results are consistent with a first-order phase transition [12]. Therefore, the left diagram of Fig. 2 is most likely the correct one. In this paper, we do not need to know the order of the phase transition as it does not affect our argument.
According to the Maldacena-Milekhin conjecture, the gauged and ungauged theories are exponentially close at low temperatures. Therefore, we will focus on the confined phase. In it, observables like e.g the energy, are independent of temperature up to 1/N corrections. We can determine these temperature-independent values by taking the large-N and continuum limit at some fixed value of T . To determine these values reliably, we will study different values of T at each µ. To extrapolate to large N and continuum regime, we use the fit ansatz This ansatz contains 1/L corrections to the continuum limit and 1/N 2 corrections to the large N limit. In practice we have considered only terms up to ε 0,1 to avoid overfitting and considering the large L values of our simulations. The extrapolations of other observables are of the same form.
The values of the flux that have been used are in the range µ = 2, 3, 4, 5 while the range of temperatures is adjusted for each µ with T ∈ [0.2, 1.1]. In addition, the size of the matrices and the lattice spacing runs over N = 8, 12, 16 and L = 12, 24, 48, 96, respectively. We extrapolate to the large N continuum limit of the theory as explained above. This procedure is repeated both for the gauged and ungauged data.

Energy of the system
In the BFSS limit µ = 0, Maldacena and Milekhin conjectured (2.30), which we repeat here: Numerical simulation in Ref. [9] suggests n adj = 2 and C adj 1. We studied the behavior of ∆E at finite µ and low temperature. The results are given in Fig. 3 with perturbative results valid at large µ and two kinds of numerical fits. The punchline at this stage is that there is an exponential decay with 1/T , which verifies the predicted functional form (2.29). We have done fits with a different number of free parameters to check the reliability of our fit ansatz. The red curves in Fig. 3 are fits of the two free parameters D E,2 (µ) and C E,2 (µ) according to The fit results are shown in Table 1. We have extrapolated C E,2 (µ) to µ = 0 assuming the quadratic dependence C E (µ) = C E (0) + Aµ 2 . As long as µ is small enough we can assume that an expansion of this form is valid. 6 The result we get is C E (0) = 0.860(68), which is consistent with the value 0.83(21) obtained in Ref. [9] and the other parameter is estimated as A = 0.087(04). On the other hand, it is difficult to extrapolate D E (µ) to µ = 0; the quadratic ansatz gives 2.755 ± 2.153, which has a very large error within numerical accuracy (see Fig. 4). In addition, we can consider the ratio between C E and D E , and as we can see in Fig. 5 it shows an increasing trend with µ. The ratio n E ≡ D E /C E changes as we vary µ. At µ = 0, simulations in Ref. [9] estimated that the numerical value is n E | µ=0 = 1.91(78) and our findings for the BMN matrix model is consistent with this.  Table 1: We fit ∆E by using the two-parameter ansatz (3.10). The µ = 0 value is the extrapolation of the fits to the µ → 0 limit.

Validity of perturbative results
The perturbative prediction is expected to hold in the limit µ → ∞. Our nonperturbative methods allow us to check how much their range of validity extends towards lower µ. The perturbative estimates are shown in Fig. 3 in comparison to the numerical data. As expected, one observes a convergence of perturbative and nonperturbative data in the low temperature limit. For large µ this range of convergence should extend up to higher temperatures.  It is instructive to discuss in more detail the perturbative predictions. We will focus in this discussion on the largest µ (µ = 5), which should provide the best convergence to the numerical data. In the large µ limit, the Hamiltonian of BMN decouples in two parts (see Appendix D.1). In the energy plots in Fig. 3 the full U (1) sector (2.32) is taken into account, but we can also investigate the relevance of the different contributions and compare the full U (1) sector with the lightest mode (2.31).
At very low temperatures and finite µ the contribution of the full U (1) sector almost coincides with the lightest mode (2.31) and we cannot distinguish them practically.
We observe that at finite µ the perturbative result of the lightest mode and the exponential fits from the numerical data cross at finite values of temperature as shown in Fig. 3. There is, consequently, no asymptotic convergence in the small temperature limit. When we include the one-loop correction given by [6,29] we do not expect things to change by much because the one-loop corrections in the gauge and ungauged models cancel each other [6]. Indeed, we checked that this is the case and at small temperatures, the difference between the two is suppressed. Fig. 6 is an additional illustration of these findings.  There is a crossing of the two-parameter data with the lightest mode (2.31) and its one-loop correction (3.11) but no crossing for the full U (1) perturbative sector (2.32).
[Right]: The difference of the perturbative result (2.32) and the two-parameter fit with respect to temperature (as the black curve) including the error bars of the data points for µ = 5. We observe an indistinguishability between the actual data and the perturbative result using the full U (1) sector within the error bars of the data. At smaller temperatures the difference approaches zero.
However, in view of the overall behavior we still conclude that the full U (1) sector (2.32) provides a cleaner description of the data since there is an asymptotic convergence instead of a crossing with the fitted non-perturbative result. Only a small shift smaller than our statistical uncertainty remains in this case. In the parameter region we studied, considering only the lightest mode shows a large deviation from the numerical data at lager temperatures. Except for the region where the fit and the lightest mode contribution cross, the full U (1) contribution is closer to the data and it seems to capture the asymptotic behaviour better. Therefore we decided to show the contribution of the full U (1) sector in Fig. 3. If we insist that the U (1) sector is protected from finite µ corrections [24][25][26]30] we could carry these results to finite µ and the slight shift from the actual data could be related with the SU (N ) sector where we do not precisely know what to expect at finite µ. This would require, however, a further detailed analysis. Such investigation can reveal more about the fine print of the perturbative calculations. Overall we can, however, conclude that our results reproduce correctly the perturbative predictions.

Other observables
The same analysis can be done for other gauge-invariant quantities such as the sum of traces of the squared matrices defined via Tr (X I ) 2 , (3.12) and the commutator Tr [X I , X J ] 2 . (3.13) Similar to the energy, we may take the difference between the gauged and ungauged results which we expect to be exponentially close to each other in the limit of (2.29). The justification for this comes from the partition functions, which are close to each other. Therefore we expect the relevant observables to scale as 14) with a priori unknown parameters D R , D F , C R and C F . We report the results in Table 2 while the respective diagrammatic fits are shown in Figs. 13 and 8. We have not shown the exponential fit for the observable F 2 because it is similar as can be seen from the logarithmic plots in appendix B. It is important to note that there is no prediction whatsoever for the parameters C R,F and D R,F . In addition, there are no constraints between C R , D R and C F , D F so we always use a two-parameter fit. Due to finite µ deformations, we expect the situation to be different than for the BFSS model. However, in the µ → 0 limit one should recover the BFSS values. Since there are only data available for the F 2 term from [9], we can compare the commutator, and indeed for C F it converges to its BFSS value as µ → 0 (right panel of Fig. 9). On the other hand, for D F we observe a discrepancy with the BFSS value. It is likely that we need to study smaller values of µ to observe convergence with the BFSS results regarding D F . In addition, the small discrepancy at intermediate µ for C E,R,F and D E,R,F could probably be an effect caused by finite µ contributions. The important result is that C E,R,F agree nicely when µ = 0 as we can see from the left panel of Fig. 10. Since, no correlation is assumed between D R,F and C R,F we can not claim the same for D R,F . We can expect that this should be the case when we are studying particularly small µ values, but on the other hand, it is quite challenging to do a large N , continuum, gauged and ungauged study and remain in the confined phase, given the current resources.  action with respect to µ in the positive-r axis. This potential is symmetric with respect to the y-axis. As we increase µ we see that when we are in the trivial background r = 0 the simulation becomes more and more stable, such that we observe a decrease of F 2 as we increase µ. The same scaling behaviour holds also for R 2 .
In Fig. 7 we show that our results are consistent with theoretical expectations. Finite µ effects could be observed even though we are putting the system in the trivial background X I = 0, I = 1, · · · , 9. This can be already observed in Monte Carlo histories (see for example Fig. 14) and by studying the SO(3) potential of the action. We remind that on top of the classical trivial background there are quantum fluctuations that give actually some non-zero expectation values for the matrices. The classical potential of the bosonic SO(3) part of the action is found to be (see Appendix A and specifically eq. On top of that, we have quantum fluctuations on the r value resulting in fluctuations on V (r). For the case of µ = 0, we immediately get the BFSS result V (r) ∼ r 4 which, in the matrix model language is given by the observable F 2 since considering quantum fluctuations the matrices do not strictly commute. We may now ask how this observable behaves for µ = 0 and the answer is clear from (3.16). Bigger values of µ result in smaller values of V (r) because the effect of finite µ is to confine the simulation more into the classical vacuum r = 0. For finite µ, the matrices (to be more precise their eigenvalues) are grouped around r = 0 while for µ = 0 they can spread in the whole range of the potential barrier and this makes (some of) the eigenvalues escape to infinity resulting in the flat direction problem. Indeed, this is the behaviour we observe in the left panel of Fig. 7 both for gauged and ungauged models, and the theoretical explanation is because for bigger values of µ the fluctuations of the potential become smaller due to a bigger potential wall created by finite µ. This is shown in the right panel of Fig. 7 and it affects both the observables {R 2 , F 2 } such that they decrease as we increase µ. At the same time, the observable F 2 for the gauged and ungauged models studied at µ = 0 in [9] is of order O(15) always while here we see clearly a decreasing trend as we increase µ.
Another point of view that arrives at the same conclusion is to consider the large µ limit. In this case, all the bosonic matrices can be written in terms of harmonic oscillators, and they all scale as 1/µ (for the SO(3) part) and 2/µ (for the SO(6) part) as it can be seen from eq.'s (D.4) and (D.5).
An interesting puzzle to understand is the role of the unstable solution of the SO(3) potential given by r = µ/2 in our conventions (see again Appendix A). It is not known in the literature how this term appears in the simulations and how it affects them. Thus, we do not exclude the possibility that the simulation reaches frequently this solution altering the results non-trivially at finite µ.

BFSS values for two-parameter fits
The dependence of the exponential parameters (C R , C F ) on µ is given in Fig. 9 and for C F it converges to the BFSS result [9]. This provides additional evidence that the limit µ → 0 is smooth and consistent but nonetheless considering smaller values of µ will be of much importance such that also D F will approach its known BFSS value. [Right]: the fit for C F is given by the equation C F = 0.825(59) + 0.074(4)µ 2 . The left most point is the BFSS point taken from [9] (see also Table: 3).

Conclusions and Discussions
We have studied the gauged and ungauged BMN model at finite flux µ values and relatively small temperatures with high statistics. We confirmed that the difference between the gauged and ungauged partition function is exponentially small, being proportional to e −C adj /T also for finite µ. This is our main result. This happens also in the case that the limit µ/T 1 is not satisfied, namely at intermediate temperatures and finite µ. Of particular importance is the exponential decay of energies in Fig. 3 supporting the two different limits, namely the gravity limit and the perturbative regime as we vary µ. In particular, at higher µ values and relatively small temperatures, the system behaves in such a way that it converges to the perturbative result (2.32). On the other hand, when we gradually decrease µ towards µ → 0 the system seems to approach the gravitational results C adj = 1 and n E → n adj = 2 obtained in [9]. To verify this numerically for n E turned out to be difficult at this level due to fitting problems, namely if we insist on an even power expansion for C E and D E then the ratio is a rational function that does not have a rapidly converging expansion in (even powers of) µ for a large range of µ. : Values of C E , C F , C R labelled as E, F, R respectively. We expect these values to be approximately the same when the partition functions of both models are exponentially close to each other. We see that this happens as µ → 0 and indeed we expect this to be the case in this limit as we can recall from eq. (2.29).
[Right]: Values of D E , D F , D R . We observe that D for Tr X 2 I and Tr [X I , X J ] 2 observables do not change with respect to µ but for energy does. This is probably due to different degeneracies of the energy eigenstate of the system.
On the other hand, we observed that the degeneracy of the energy states can change with µ as it is shown in Fig. 5. It seems that by turning on µ there could be some massive modes that affect the simulation in a non-trivial way. In particular, in Fig. 5 we observe a consistent decrease of the ratio n E = D E /C E as we decrease µ. This ratio appears to be the degeneracy of the eigenstate of the system [6] and it seems that it changes with µ. This is an interesting and puzzling issue since the respective change is not known theoretically. One possibility is that, at finite µ, the contribution from second-lightest or even higher modes is not yet negligible in this parameter region. We tried to fit ∆E by using two different excitations as ∆E = 2N 2 C E (µ)e −C E (µ)/T + n N 2 C E (µ)e −C E (µ)/T with a few different values of n (including n = 6, which is the number of lightest adjoint modes in the limit of µ → ∞), but we were not able to obtain reasonable fits.
Indeed, it seems that we have a factor behaving like an effective degeneracy. The cause of this behaviour is most likely finite µ effects whose form in the interacting gauged sector has not been studied extensively. It may be the case that several low-energy excited modes with slightly different energies are contributing. The precise details of this remain unknown to us but we hope that this will initialise a more systematic study of the precise finite µ contributions to the partition function of the model, perturbatively and (if possible) also non-perturbatively.
That non-singlet modes are exponentially suppressed can make quantum simulation based on the extended Hilbert space simpler. Suppose that Hamiltonian time evolution is performed on a quantum computer, as |Φ → e −iĤt |Φ . If |Φ is gauge-invariant, e −iĤt |Φ is also gauge-invariant, as long as the time evolution is exactly realized. If |Φ is a specific gauge fixed state, e.g., a fuzzy sphere with a certain representation, then the gauge fixing will not be spoiled via the Hamiltonian time evolution. However, if there were light non-singlet modes, small simulation errors could easily excite nonsinglet modes and lead to a large deviation from the exact result. Our findings in this paper suggests that we do not have to worry about such a possibility.

Data management
No additional research data beyond the data presented and cited in this work are needed to validate the research findings in this work. Simulation data will be publicly available after publication.

A SO(3) potential profile
We are interested in finding the profile of the SO(3) classical potential in the BMN model. This feature also allows finding the probability for the trivial vacuum configuration to transition to a fuzzy sphere configuration.
Let us then concentrate on the bosonic action of the BMN model which has the following form To find extrema of this potential let us consider an ansatz of the form X i (t) = ρ(t)J i and substitute it into the Lagrangian. Using identities of ijk and that for irreducible (considering actually the maximal) SU (2) representations we have Tr (J i ) 2 = N 3 (N 2 − 1), we find with potential We are considering the case where we have just one big fuzzy sphere, which is a configuration that can be distinguished by the simulation and specifically by the Myers observable (the cubic term in (2.4)).
Asking stability of this potential, we differentiate with respect to ρ and we obtain three solutions • ρ = 0 which is stable, • ρ = µ which is stable, • ρ = µ 2 which is unstable. For ρ = 0 we have the trivial background solution, while for ρ = µ we have a fuzzy sphere of radius ρ. This solution is known as the fuzzy sphere background and in fact, ρ gives us the size of the SO(3) subpart of the spacetime. On the other hand, for ρ = µ 2 we have an unstable solution whose interpretation we are not aware of. The maximum of the potential is given by this solution In the large N limit, this potential barrier behaves like ∼ µ 4 N 4 , which means that the combination µ, N is what matters such that we keep under control the system. Specifically, in simulations, it is of much importance to keep this combination large to avoid undesirable tunnellings from the trivial backgrounds to fuzzy sphere backgrounds.

B Details for log plots
In this Appendix we also report the logarithmic plots for all observables E, R 2 , F 2 in Figs. 11, 12, 13 respectively. where

C Monte Carlo histories
with Π M = δL δẊ M being the conjugate momenta for bosonic matrices. The H 0 terms construct the free U (1) sector of the model while H int denotes the interactive SU (N ) part. In the large µ limit, the interactive sector can be treated perturbatively while the free sector is claimed to be protected from contributions to all orders in µ [24][25][26]30]. In addition we can introduce harmonic oscillators defined by the operators which obey canonical commutation relations Then, the bosonic part of the free Hamiltonian results in Complexifying the real spinor matrices ψ α using yields the anticommutation relations We may now use the chirality property of the complexified fermions (iγ 123 )ψ ± = ±ψ ± accompanied with (C ± ) 2 = C ± and C + C − = 0. Recalling the splitting ψ = ψ + + ψ − the fermionic part results in Tr iψ α γ 123 ψ α = 3µ 2 Tr ψ +α ψ −α . (D.10) Summing both we get the U (1) free part of the Hamiltonian written in terms of bosonic and fermionic harmonic oscillators In our conventions, the SO(3) and SO(6) sectors have mass µ and µ/2 respectively while the fermions have mass 3µ/4 7 . The zero-energy ground state of this free Hamiltonian is denoted as |0 and is annihilated by The U (1) free sector of the above Hamiltonian is spanned by excitations of operators of word-length one, e.g Tr [A † i ], Tr [B † a ], Tr [ψ α ] while the SU (N ) free sector is spanned by operators of word-length two and larger. The lightest mode is given by the SO(6) part of the free Hamiltonian with lowest energy For the gauged model this is the first excited state created by Tr B † a |0 , while the remaining excitations of the free sector are shown in Table 4. For the ungauged model, one can simply act with B † a on |0 which results in the same energy. Six oscillators are yielding a six-fold degeneracy of this sector. In addition there are n = 3 oscillators for the SO(3) part with energy E SO(3) = µ and n = 8 fermionic oscillators with energy E fermions = 3µ/4. The ground state energy vanishes since due to supersymmetry we have 14) The lowest adjoint mode of the gauged theory in the perturbative limit µ → ∞ is (D.13) and perturbation theory shows that it is protected at least to first order in µ [6,24,26]. The spectrum of the free Hamiltonian has been studied in perturbation theory (with perturbative parameter ∼ 1/µ, e.g (2.8)) in [24,26,29] and their energy, representations and degeneracy are given in Table 4. It was conjectured that the free spectrum does not receive any perturbative corrections to all orders in µ −1 . On the other hand, short representations of the SU (N ) sector can, in principle, combine and form multiplets, and indeed they may receive perturbative corrections. On top of that we may also note that there could be non-perturbative corrections [26]. A precise analysis of the form of non-perturbative and perturbative corrections for the construction of multiplets has not been done and, therefore, the energy correction can not be estimated precisely. Whether or not this is something that can be justified analytically we do not know, because also the non-perturbative corrections can not be estimated at all. state SO(6) × SO(3) reps. energy degeneracy |0 (1,1) 0 1 Tr B † a |0 (6,1) (1,1)+(20,1) µ 1+20 Table 4: Lowest energy states for the trivial background X = 0, their representations and degeneracy. The first five lines correspond to the U (1) part of the model which is free.

D.2 Representation algebra of the BMN model
Finite µ corrections to n E are not known. We observed that they change with µ and here we can ask whether or not the works of [30] and [26] could provide some insight. On the other hand, finite µ corrections to C adj depend on the mixing of states, presumably coming from the interacting sector of the model whose precise form is also unknown. Let us highlight some of these results in the literature. The classification of the superalgebra in the plane wave limit of M-theory has been considered in [30]. Indeed, in this work and building on the results from Kac [31], it was shown that the complexification of the special unitary Lie superalgebras su(2|4; 2, 0) (for µ > 0) or su(2|4; 2, 4) (for µ < 0) results in A(1, 3) ∼ = sl (2,5). However, here we are always in a scenario of positive µ, and hence the bosonic part of the algebra is given as su(2, 0)⊕su(0, 4) ∼ = so(3)⊕so (6).
The peculiar supersymmetry of the BMN model manifests itself in the time dependence of the supersymmetry transformations. In the large µ limit the Hamiltonian splits into a free and interacting part (see Appendix D.1). The commutator between the interacting Hamiltonian and the supersymmetric charge is proportional to the supersymmetric charge. This implies that the boson and fermion masses differ. In particular, the energy level difference given by the commutator between the supercharges Q α and the Hamiltonian in the interacting SU (N ) sector is µ 4 in our conventions whenever the Gauss constraint is satisfied. The application of the supersymmetry charges Q ± α , where one performs a chirality split as Q ± = C ± Q , C ± = 1 2 (1 ± iγ 123 ) , (D. 16) can be applied at most eight times for Q + or eight times for Q − changing the energy level by Q ± |state = ∓ µ 4 |state . This results in a multiplet consisting of 2 8 = 256 states with irreducible representations at each level. Thus, one can change irreducible representations by acting with Q ± resulting in a shift of the energy at each level. This energy shift coincides with the difference between bosonic and fermionic masses µ − 3µ 4 = µ 4 = 3µ 4 − µ 2 in the SO(3) and SO(6) sectors respectively. The successful application of a supercharge Q α on a state leads to changing one u(1) ⊕ su(2) ⊕ su(4) irreducible representation to another [30]. The multiplets (energy levels) split into a specific positive integer given by m. The energy of each multiplet differs from the vacuum (E 0 ) as [30] E = E 0 + µm 4 , with m = 0, 1, · · · , 8. accompanied by the condition (γ 1 + 1) = 0. This is a supersymmetric deformation of the ungauged theory and in [6] it was shown that this shift compensates the gauge condition in the superalgebra resulting again in [H new , Q α ] = µ 4 Q β γ 123 βα .
(D. 19) We expect that in the ungauged model, extra degrees of freedom living on the boundary of the theory and interpreted as open strings [6] (which in principle can reach deep into the bulk) may contribute in the partition function. These degrees of freedom transform in the adjoint representation. Therefore, one is led to ask what is the adjoint representation of the superalgebra of the BMN model. This question has been answered in [30] and the result is given in Table 5 Energy Representations +µ/4 (2,4) 0 (1,1) ⊕ (3,1) ⊕ (1,15) −µ/4 (2,4)