Finite-temperature phase diagram of the BMN matrix model on the lattice

The idea of gauge/gravity duality allows us to investigate aspects of quantum gravity by studying the dual gauge theory and vice versa. For example, the thermodynamic properties of maximally supersymmetric Yang–Mills (SYM) theories in $(p+1)$ dimensions with large SU($N$) gauge groups and strong ’t Hooft couplings capture the semi-classical thermodynamics of non-extremal black D$p$ branes Itzhaki et al. (1998). The first concrete gauge/gravity duality proposal Maldacena (1998a) considered $p=3$, conjecturing an equivalence between conformal $\mathcal{N}=4$ SYM and Type IIB supergravity on AdS${}_{5}\times S^{5}$. In lower dimensions $p<3$, the situation is more complicated because the supersymmetric gauge theories are not conformal. This motivates the use of non-perturbative lattice field theory calculations to analyze these strongly coupled gauge theories to gain insights into quantum aspects of the dual gravitational systems.

In recent years, there has been progress in extracting $p<2$ black hole physics from numerical Monte Carlo analyses of the dual strongly coupled gauge theories at finite temperatures and large $N$ Catterall and Wiseman (2007); Anagnostopoulos et al. (2008); Catterall and Wiseman (2008); Catterall (2009); Hanada et al. (2009a, b); Catterall and Wiseman (2010); Nishimura (2009); Catterall and van Anders (2010); Catterall et al. (2010, 2012); Kadoh and Kamata (2012); Hanada et al. (2014); Giguère and Kadoh (2015); Kadoh and Kamata (2015); Filev and O’Connor (2016); Hanada et al. (2016); Berkowitz et al. (2016a, b); Kadoh (2017); Rinaldi et al. (2018); Catterall et al. (2018); Jha et al. (2018); Berkowitz et al. (2018); Asano et al. (2018); Schaich et al. (2020); Bergner et al. (2022a); Schaich et al. (2022); Pateloudis et al. (2022, 2023), which is currently being extended to $p=2$ Catterall et al. (2020); Sherletov and Schaich (2022, 2023); Schaich and Sherletov (2025, in preparation). See Ref. Schaich (2023) for a recent review of these developments, including further references. In the context of lattice calculations, the lower dimensionality helps both to reduce the computational costs as well as to simplify the process of taking the continuum limit: As discussed in Refs. Golterman and Petcher (1989); Giedt et al. (2004); Bergner et al. (2008); Catterall and Wiseman (2007); Giedt et al. (2018), little to no fine-tuning is required to recover the supersymmetric target theory in the continuum limit.

When $p=0$, maximal SYM reduces to a quantum-mechanical theory de Wit et al. (1988) best known due to a conjecture by Banks, Fischler, Shenker and Susskind (BFSS) Banks et al. (1997) that the large-$N$ limit of this model describes the strong-coupling (M-theory) limit of Type IIA string theory in the infinite-momentum frame. Several groups have numerically studied this ‘BFSS model’, using either a lattice discretization of (Euclidean) time Catterall and Wiseman (2007, 2008, 2010); Kadoh and Kamata (2012, 2015); Filev and O’Connor (2016); Berkowitz et al. (2016a, b); Rinaldi et al. (2018); Berkowitz et al. (2018); Pateloudis et al. (2023) or a ‘non-lattice’ momentum-space approach Anagnostopoulos et al. (2008); Hanada et al. (2009a, b); Nishimura (2009); Hanada et al. (2014, 2016). These investigations have reached sufficiently low temperatures and large $N$ to provide confidence that the observed thermodynamic behavior approaches the leading-order prediction of the dual Type IIA supergravity. Fits to these numerical Monte Carlo results have also been used to predict $\alpha^{\prime}$ and $g_{s}$ corrections to the classical supergravity solutions from the dual SYM quantum mechanics Hanada et al. (2009b); Kadoh and Kamata (2015); Hanada et al. (2016); Berkowitz et al. (2016a, b); Pateloudis et al. (2023).

In the BFSS model, the thermal partition function is not well defined since it includes integration over a non-compact moduli space. This instability cannot be seen directly in the black D0-brane thermodynamics described by large-$N$ supergravity, as it is a $1/N$ effect. Even so, it means that Monte Carlo calculations of the BFSS model are unstable for any finite $N$ Catterall and Wiseman (2010); as $N$ increases, the system may spend longer fluctuating around a metastable vacuum, but this will eventually decay given sufficient time. The standard way to address this issue and stabilize numerical calculations is to deform the theory by adding a regulator that gives a small mass to the scalars. This explicitly breaks supersymmetry, and the extrapolation necessary to remove this deformation can significantly increase the complexity and costs of the calculations.

This motivates the alternative of studying the Berenstein–Maldacena–Nastase (BMN) deformation of the BFSS model Berenstein et al. (2002), which describes the discrete light-cone quantization compactification of M-theory on the maximally supersymmetric “pp-wave” background of 11d supergravity. The ‘BMN model’ preserves maximal supersymmetry and has a well-defined thermal partition function because of the absence of flat directions, making it a better-behaved system to study using Monte Carlo methods. In the $N\to\infty$ planar limit, it can be thought of as a compactification of four-dimensional $\mathcal{N}=4$ SYM over a spatial sphere with large R-charge Ishiki et al. (2009). In addition, the BMN model has no $\alpha^{\prime}$ corrections in the regime where the dual supergravity description is valid Horowitz and Steif (1990), in contrast to the BFSS model.

Finally, the BMN model possesses a non-trivial finite-temperature phase diagram with a high-temperature deconfined phase and a low-temperature confined phase.¹¹1While the BFSS model has historically been expected to be deconfined for all temperatures, Refs. Bergner et al. (2022a); Pateloudis et al. (2023) recently argued that a confined phase persists in the BFSS limit in which the BMN deformation is removed. In this sense, it resembles higher-dimensional theories Catterall et al. (2018); Jha et al. (2018); Catterall et al. (2020); Sherletov and Schaich (2022, 2023); Schaich and Sherletov (2025, in preparation), while involving more modest computational costs in numerical analyses. However, the gravity dual of the BMN model is more complicated than that of the BFSS model and, therefore, is less explored Lin and Maldacena (2006); Costa et al. (2015).

In this work, we study the phase diagram of the BMN model on the lattice, determining critical deconfinement transition temperatures at a range of couplings spanning three orders of magnitude from the perturbative regime to stronger couplings that approach the supergravity solution. We carry out calculations with multiple lattice sizes and numbers of colors $N$ while also checking discretization artifacts for our lattice action and evaluating the Pfaffian of the fermion operator to ensure we do not encounter a sign problem. We begin in the next section by reviewing the BMN model and previous lattice studies of it. In Sec. III, we discuss the simple lattice formulation that we employ in this work, also addressing discretization artifacts and the Pfaffian. Our numerical results for the phase diagram are presented in Sec. IV, including comparisons to predictions from perturbation theory and the dual supergravity, as well as analyses of the order of the transition. The data leading to these results are available through Ref. Jha et al. (2024). We conclude in Sec. V with a discussion of our plans for future work.

Our starting point is the BFSS model in Euclidean time $\tau$. With an anti-Hermitian basis for the SU($N$) generators, $\text{Tr}\left[T^{a}T^{b}\right]=-\delta^{ab}$, the continuum action takes the form

	$\displaystyle S_{\text{BFSS}}=\frac{N}{4\lambda}\int d\tau\ \mbox{Tr}\Bigg{[}$	$\displaystyle-\left(D_{\tau}X_{i}\right)^{2}-\frac{1}{2}\sum_{i<j}\left[X_{i},X_{j}\right]^{2}$		(1)
		$\displaystyle+\Psi_{\alpha}^{T}\gamma_{\tau}^{\alpha\beta}D_{\tau}\Psi_{\beta}$		(2)
		$\displaystyle+\frac{1}{\sqrt{2}}\Psi_{\alpha}^{T}\gamma_{i}^{\alpha\beta}\left[X_{i},\Psi_{\beta}\right]\Bigg{]},$

where $D_{\tau}=\partial_{\tau}+[A,\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ ]$ is the covariant derivative corresponding to the gauge field $A$, $X_{i}$ are the nine scalars of the theory, and $\Psi_{\alpha}$ is a sixteen-component spinor. The indices $i,j=1,\cdots,9$ while $\alpha,\beta=1,\cdots,16$; the latter will generally be suppressed, and we have also suppressed the $\tau$ dependence of all the fields. $\gamma_{i}$ and $\gamma_{\tau}$ are $16\!\times\!16$ Euclidean gamma matrices; in the next section, we will specify the representation we employ for them. All fields transform in the adjoint representation of the SU($N$) gauge group. In this ($0+1$)-dimensional setting, both the ’t Hooft coupling $\lambda\equiv g_{\text{YM}}^{2}N$ and the Yang–Mills coupling $g_{\text{YM}}^{2}$ are dimensionful.

The theory defined by Eq. (2) has flat directions when the scalar fields $X_{i}$ commute with each other. Supersymmetry preserves some of these flat directions at the quantum level, which results in an ill-defined thermal partition function as discussed in Catterall and Wiseman (2010). The BMN model changes this by adding the following terms to the action:

	$\displaystyle S_{\mu}=-\frac{N}{4\lambda}\int d\tau\ \mbox{Tr}\Bigg{[}$	$\displaystyle\left(\frac{\mu}{3}X_{I}\right)^{2}+\left(\frac{\mu}{6}X_{M}\right)^{2}$		(3)
		$\displaystyle+\frac{\mu}{4}\Psi_{\alpha}^{T}\gamma_{123}^{\alpha\beta}\Psi_{\beta}$		(4)
		$\displaystyle+\frac{\sqrt{2}\mu}{3}\epsilon_{IJK}X_{I}X_{J}X_{K}\Bigg{]},$

with dimensionful deformation parameter $\mu$. The indices $i,j$ divide into the two sets where $I,J,K$ take values in $1,2,3$ while $M=4,\cdots,9$. We also define

$$\gamma_{123}\equiv\frac{1}{3!}\epsilon_{IJK}\gamma_{I}\gamma_{J}\gamma_{K}=\gamma_{1}\gamma_{2}\gamma_{3}.$$

(5)

For non-zero $\mu$, the terms in Eq. (4) explicitly break the SO(9) global symmetry of Eq. (2) down to $\text{SO(}6\text{)}\!\times\!\text{SO(}3\text{)}$ while preserving all 16 supersymmetries of the theory. Integrating over the fermions produces the Pfaffian of the fermion operator, $\text{pf}\,(\mathcal{M})$. The Euclidean path integral in the continuum then takes the form

$$Z=\int[\mathscr{D}A][\mathscr{D}X]\ \text{pf}\,(\mathcal{M})\ e^{-S_{B}},$$

(6)

where the bosonic action is

	$\displaystyle S_{B}=-\frac{N}{4\lambda}\int d\tau\ \mbox{Tr}\Bigg{[}$	$\displaystyle\left(D_{\tau}X_{i}\right)^{2}+\frac{1}{2}\sum_{i<j}\left[X_{i},X_{j}\right]^{2}$		(7)
		$\displaystyle+\left(\frac{\mu}{3}X_{I}\right)^{2}+\left(\frac{\mu}{6}X_{M}\right)^{2}$		(8)
		$\displaystyle+\frac{\sqrt{2}\mu}{3}\epsilon_{IJK}X_{I}X_{J}X_{K}\Bigg{]},$

and the fermion operator is

$$\mathcal{M}(A,X)=\gamma_{\tau}D_{\tau}+\gamma_{i}[X_{i},\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ ]-\frac{\mu}{4}\gamma_{123}.$$

(9)

We can study this theory at finite temperature $T$ by formulating it on a temporal circle of circumference $\beta=1/T$. The non-zero temperature breaks supersymmetry, and according to gauge/gravity duality, the finite-temperature theory corresponds to a black-hole geometry in the dual supergravity. To make sure that the $\alpha^{\prime}$ corrections near the horizon of the dual black hole are small, the system must be in the regime $T/\lambda^{1/3}\ll 1$ and $\mu/\lambda^{1/3}\ll 1$. It is convenient to define the dimensionless coupling $g\equiv\lambda/\mu^{3}$, in terms of which the latter constraint is $g\gg 1$. In our numerical calculations, we will identify the critical temperature of the deconfinement transition by fixing $g$ and scanning in the dimensionless ratio $T/\mu$, which we will call just “the temperature”.

There are predictions for the critical temperature $\left(T/\mu\right)_{\text{crit.}}$ in both the weak- and strong-coupling limits. In weak-coupling limit $g\to 0$, perturbative calculations Furuuchi et al. (2003); Spradlin et al. (2004); Hadizadeh et al. (2005) predict a first-order deconfinement transition with critical temperature

$$\lim_{g\to 0}\left.\frac{T}{\mu}\right|_{\text{crit.}}=\frac{1}{12\ln 3}\approx 0.076.$$

(10)

In the $N\to\infty$ planar limit, Refs. Spradlin et al. (2004); Hadizadeh et al. (2005) also calculate the $\mathcal{O}(g)$ and $\mathcal{O}(g^{2})$ corrections,

$$\left.\frac{T}{\mu}\right|_{\text{crit.}}=\frac{1}{12\ln 3}\left[1+C_{\text{NLO}}\left(27g\right)-C_{\text{NNLO}}\left(27g\right)^{2}+\cdots\right]$$

(11)

where the perturbative expansion parameter is $27g$, while

	$\displaystyle C_{\text{NLO}}$	$\displaystyle=\frac{2^{6}\cdot 5}{3^{4}}\approx 4$
	$\displaystyle\text{and\leavevmode\nobreak\ }C_{\text{NNLO}}$	$\displaystyle=\frac{23\cdot 19\,927}{2^{2}\cdot 3^{7}}+\frac{1\,765\,769\ln 3}{2^{4}\cdot 3^{8}}\approx 71.$

For strong coupling $g\to\infty$ with $T\ll\lambda^{1/3}$, the numerical construction of the dual supergravity solutions at finite temperature by Ref. Costa et al. (2015) predicts a larger

$$\lim_{g\to\infty}\left.\frac{T}{\mu}\right|_{\text{crit.}}\approx 0.106,$$

(12)

with the transition still first order.

Our ambition is to use non-perturbative lattice calculations to track $\left(T/\mu\right)_{\text{crit.}}$ as the system interpolates between the weak-coupling limit of Eq. (10) and the strong-coupling limit of Eq. (12), similar to what we did for the bosonic sector of the BMN model in Refs. Dhindsa et al. (2022, 2024). The work we present here finalizes the investigations reported in preliminary form by our recent conference proceedings Schaich et al. (2020, 2022). Other groups have also studied the phase diagram of the BMN model on the lattice Catterall and van Anders (2010); Asano et al. (2018); Bergner et al. (2022a).²²2More recent lattice studies of the BMN model have focused on ‘ungauging’ Pateloudis et al. (2022) or the energy of the dual black D$0$-branes Pateloudis et al. (2023) as opposed to the phase diagram of interest here. All of these three previous lattice studies take different approaches to explore the phase diagram of the theory, which also differ from our approach presented below.

The earliest study, Ref. Catterall and van Anders (2010), was limited to a small number of lattice sites ($N_{\tau}=5$) and small numbers of colors ($N=3$ and 5). Apart from constant rescalings of fields and parameters, this work uses the same lattice action as we describe in the next section. Converting to our conventions, Ref. Catterall and van Anders (2010) fixes $T/\mu=1/3$ and scans in the coupling $g$ — essentially the opposite of our approach — finding a deconfinement transition with critical coupling $g_{\text{crit.}}\simeq 0.035$. While $T/\mu=1/3$ is larger than both the $g\to 0$ and $g\to\infty$ limits discussed above, the small values of $N_{\tau}$ and $N$ considered by Ref. Catterall and van Anders (2010) may make those limits inapplicable.

In addition to significantly increasing the number of lattice sites (focusing on $N_{\tau}=24$) and colors (focusing on $N=8$), Ref. Asano et al. (2018) also employs a second-order discretization of the covariant derivative, which introduces four new parameters whose values were chosen to minimize lattice artifacts in the inverse Dirac operator. This work fixes $g$ by considering fixed values of $\mu/\lambda^{1/3}=1/g^{1/3}$, focusing on $9\geq\mu/\lambda^{1/3}\geq 1$ that correspond to $0.00137\leq g\leq 1$. For each fixed $\mu/\lambda^{1/3}$, Ref. Asano et al. (2018) scans in a dimensionless temperature $T/\lambda^{1/3}$, which also differs from our approach described below. For larger $\mu/\lambda^{1/3}\gtrsim 3$ (weaker couplings $g\lesssim 0.037$), this work observes distinct deconfinement and SO(9)-symmetry-breaking transitions, all first order, which merge for larger $g$.

Finally, Ref. Bergner et al. (2022a) employs gauge fixing to the static diagonal gauge, $A(\tau)=\mbox{diag}(\alpha_{1},\cdots,\alpha_{N})/\beta$, which reduces the number of degrees of freedom and allows investigations of lattice sizes up to $N_{\tau}=48$ and up to $N=48$ colors. This work again fixes $\mu/\lambda^{1/3}$ and scans in $T/\lambda^{1/3}$, focusing on smaller $5/3\geq\mu/\lambda^{1/3}\leq 0.1$ that imply stronger couplings $0.216\leq g\leq 1000$, with particular interest in the $\mu\to 0$ limit where the BMN model reduces to the BFSS model. Ref. Bergner et al. (2022a) confirms the first-order nature of the transition by presenting two-state signals in the Polyakov loop for $2/3\gtrsim\mu/\lambda^{1/3}\geq 0.17$ corresponding to $3.4\lesssim g\leq 216$. For $\mu/\lambda^{1/3}\approx 1/3\to g\approx 27$, the resulting $\left(T/\mu\right)_{\text{crit.}}$ agrees with the strong-coupling limit in Eq. (12). Larger critical temperatures at stronger couplings, with fixed $N=12$, are attributed to finite-$N$ corrections.

Like Ref. Catterall and van Anders (2010), and unlike Refs. Asano et al. (2018); Bergner et al. (2022a), we use a simple gauge-invariant lattice discretization to formulate the BMN matrix model on a one-dimensional Euclidean lattice with thermal boundary conditions (periodic for the bosons and anti-periodic for the fermions). The temporal extent of the lattice, $\beta=aN_{\tau}$, is divided among $N_{\tau}$ sites separated by lattice spacing ‘$a$’. The corresponding temperature is $T=1/(aN_{\tau})$. Our numerical calculations involve the dimensionless lattice parameters $\mu_{\text{lat}}=a\mu$ and $\lambda_{\text{lat}}=a^{3}\lambda$. Rather than considering the ratio $T/\lambda^{1/3}$ used by Refs. Asano et al. (2018); Bergner et al. (2022a), we focus on the dimensionless parameters introduced in the previous section, which remain consistent in both the lattice and continuum theories:

$\displaystyle\frac{T}{\mu}$

$\displaystyle=\frac{1}{N_{\tau}\mu_{\text{lat}}}$

$\displaystyle g$

$\displaystyle\equiv\frac{\lambda}{\mu^{3}}=\frac{\lambda_{\text{lat}}}{\mu_{\text{lat}}^{3}}.$

(13)

We use the following lattice discretization of the covariant derivative:

$$\mathcal{D}_{\tau}=\begin{pmatrix}0&\mathcal{D}_{\tau}^{+}\\ \mathcal{D}_{\tau}^{-}&0\end{pmatrix},$$

(14)

where $\mathcal{D}_{\tau}^{-}$ is the adjoint of $\mathcal{D}_{\tau}^{+}$. The finite-difference operator $\mathcal{D}_{\tau}^{+}$ acts on the fermions at lattice site $n$ as

$$\mathcal{D}_{\tau}^{+}\Psi_{n}=U(n)\Psi_{n+1}U^{{\dagger}}(n)-\Psi_{n},$$

(15)

where $U(n)$ is the Wilson gauge link connecting site $n$ and site $n+1$. That is, under an SU($N$) lattice gauge transformation, $U(n)\to G(n)U(n)G^{{\dagger}}(n+1)$. $U^{{\dagger}}(n)$ is the adjoint link with the opposite orientation. Our choice of $\mathcal{D}_{\tau}$ provides the correct number of fermions, free from extraneous ‘doublers’, which is readily seen by setting the gauge links to unit matrices and computing $\det\mathcal{D}_{\tau}=\det\left(\Delta^{+}\Delta^{-}\right)$, the determinant of the scalar Laplacian.

The resulting bosonic action of the lattice theory is

	$\displaystyle S_{B}=-\frac{N}{4\lambda_{\text{lat}}}\sum_{n=0}^{N_{\tau}-1}\mbox{Tr}$	$\displaystyle\Bigg{[}\left(\mathcal{D}_{\tau}X_{n}^{i}\right)^{2}+\frac{1}{2}\sum_{i<j}\left[X_{n}^{i},X_{n}^{j}\right]^{2}$
		$\displaystyle+\left(\frac{\mu_{\text{lat}}}{3}X_{n}^{I}\right)^{2}+\left(\frac{\mu_{\text{lat}}}{6}X_{n}^{M}\right)^{2}$		(16)
		$\displaystyle+\frac{\sqrt{2}\mu_{\text{lat}}}{3}\epsilon_{IJK}X_{n}^{I}X_{n}^{J}X_{n}^{K}\Bigg{]},$

while the fermion operator $\mathcal{M}(U,X)$ matches Eq. (9) with the covariant derivative $D_{\tau}$ replaced by the finite-difference operator $\mathcal{D}_{\tau}$. For the gamma matrices we use a representation where $\gamma_{8}$ and $\gamma_{9}=-i\mathbb{I}_{16}$ are diagonal while

$$\gamma_{\tau}=\begin{pmatrix}0&\mathbb{I}_{8}\\ \mathbb{I}_{8}&0\end{pmatrix}=\sigma_{1}\otimes\mathbb{I}_{8}$$

(17)

and the other $\gamma$ matrices are all real, with $i\sigma_{2}$ in the outermost Kronecker product — see App. A for further details or Ref. Bergner et al. (2022b) for the explicit construction.

The lattice action as a whole is finite in lattice perturbation theory. The only divergences that can occur arise from a one-loop fermion tadpole, which vanishes in a non-abelian theory because of gauge invariance. This is shown for the BFSS model in Ref. Catterall and Wiseman (2007), which remains applicable in the BMN case. Hence, at the classical level, the system will flow to the correct continuum supersymmetric target theory without fine-tuning as the lattice spacing is reduced. Thus, the continuum limit is obtained by extrapolating $N_{\tau}\to\infty$ with $g$ and $T/\mu$ fixed, while the thermodynamic limit corresponds to increasing the number of colors, $N\to\infty$.

In numerical calculations, we use the standard rational hybrid Monte Carlo (RHMC) algorithm Clark and Kennedy (2007), which we have implemented in the publicly available parallel software package for lattice supersymmetry Bergner et al. (2022b) presented by Ref. Schaich and DeGrand (2015). In contrast to the higher-dimensional theories that Ref. Schaich and DeGrand (2015) focuses on, for the BMN model, we employ the gauge-invariant Haar measure in the lattice path integral

$$Z=\int[\mathscr{D}U][\mathscr{D}X]\ \text{pf}\,(\mathcal{M})\ e^{-S_{B}},$$

(18)

and do not preserve any exact supersymmetries at non-zero lattice spacing $a>0$. The RHMC algorithm can encounter instabilities when the coupling $g$ is too strong, in a way that depends on both the lattice action as well as the lattice spacing set by $N_{\tau}$. For the simple lattice action above, and considering our coarsest lattice spacing with $N_{\tau}=8$, we encounter instabilities for $g\gtrsim 0.025$, significantly smaller than the values reached by Refs. Asano et al. (2018); Bergner et al. (2022a). This may motivate switching to a more complicated improved action in future work.

Lattice discretization introduces $a$-dependent artifacts in numerical calculations, which also depend on the lattice action and typically increase with $g$. In particular, we have observed that for our simple lattice action, discretization artifacts break the expected SO(6) symmetry Schaich et al. (2020, 2022). Specifically, the six gauge-invariant $\text{Tr}\left[X_{M}^{2}\right]$ split into a set of two with larger values and a set of four with smaller values — still significantly larger than the three $\text{Tr}\left[X_{I}^{2}\right]$. To quantify this effect, we consider the ratio

$$R_{\text{SO(6)}}\equiv\frac{\left\langle\text{Tr}\left[X_{(2)}^{2}\right]\right\rangle-\left\langle\text{Tr}\left[X_{(4)}^{2}\right]\right\rangle}{\left\langle\text{Tr}\left[X_{(6)}^{2}\right]\right\rangle},$$

(19)

where $\left\langle\text{Tr}\left[X_{(2)}^{2}\right]\right\rangle$, $\left\langle\text{Tr}\left[X_{(4)}^{2}\right]\right\rangle$, and $\left\langle\text{Tr}\left[X_{(6)}^{2}\right]\right\rangle$ average over the two larger traces, the four smaller traces and all six of them, respectively. In Fig. 1, we plot this ratio for three different couplings $g=0.001$, $0.002$, and $0.01$, considering a small SU(4) gauge group to access large lattice sizes up to $N_{\tau}=128$ close to the continuum limit. We find that the SO(6) breaking vanishes linearly in the $N_{\tau}\to\infty$ continuum limit, confirming that it is merely a discretization artifact, which may be reduced or removed by employing an improved lattice action in future work.

One other complication is that the RHMC algorithm treats the factor of $\text{pf}\,(\mathcal{M})e^{-S_{B}}$ in the path integral Eq. (18) as a real, positive Boltzmann weight. However, our lattice discretization of the BMN model allows the Pfaffian to be complex, $\text{pf}\,(\mathcal{M})=|\text{pf}\,(\mathcal{M})|e^{i\phi}$. We proceed by ‘quenching’ the phase $e^{i\phi}\to 1$ Schaich and DeGrand (2015), which in principle requires reweighting in order to recover the true expectation values $\left\langle\mathcal{O}\right\rangle$ from phase-quenched (‘${}_{\text{pq}}$’) calculations. That is,

	$\displaystyle\left\langle\mathcal{O}\right\rangle=\frac{\int[\mathscr{D}U][\mathscr{D}X]\ \mathcal{O}e^{-S_{B}}\ \text{pf}\,(\mathcal{M})}{\int[\mathscr{D}U][\mathscr{D}X]\ e^{-S_{B}}\ \text{pf}\,(\mathcal{M})}=\frac{\left\langle\mathcal{O}e^{i\phi}\right\rangle_{\text{pq}}}{\left\langle e^{i\phi}\right\rangle_{\text{pq}}},$		(20)
	$\displaystyle\mbox{where}\leavevmode\nobreak\ \left\langle\mathcal{O}\right\rangle_{\text{pq}}=\frac{\int[\mathscr{D}U][\mathscr{D}X]\ \mathcal{O}e^{-S_{B}}\ |\text{pf}\,(\mathcal{M})|}{\int[\mathscr{D}U][\mathscr{D}X]\ e^{-S_{B}}\ |\text{pf}\,(\mathcal{M})|}.$		(21)

Reweighting requires computationally expensive measurements of the Pfaffian phase $\left\langle e^{i\phi}\right\rangle_{\text{pq}}$ and fails if this expectation value is consistent with zero (a ‘sign problem’).

In Figs. 2 and 3, we present some checks of the phase of the Pfaffian $\left\langle e^{i\phi}\right\rangle_{\text{pq}}$, which show that our phase-quenched numerical results are not significantly affected by phase reweighting. In Fig. 2 we again consider gauge group SU(4), plotting the real part of $\left\langle e^{i\phi}\right\rangle_{\text{pq}}$ against $\lambda_{\text{lat}}=g/[(T/\mu)N_{\tau}]^{3}$. The right-most point for each data set corresponds to $N_{\tau}=8$, while the largest systems we consider here have $N_{\tau}=24$. For $N_{\tau}=24$, the Pfaffian phase measurement for a single field configuration is roughly 300 times more computationally expensive than generating a single molecular dynamics time unit (MDTU) with the RHMC algorithm. The plot confirms that the Pfaffian becomes real and positive in the $N_{\tau}\to\infty$ continuum limit, with no sign problem for sufficiently large lattices.

Finally, Fig. 3 shows the same quantity for some SU(8) systems similar to (but smaller than) those we analyze to determine the phase diagram in the next section. Here $N_{\tau}=8$ and the computational cost for each Pfaffian phase measurement increases to roughly 3000 times the cost of generating an MDTU with the RHMC algorithm. We consider three values of $g=10^{-5}$, $10^{-4}$, and $0.001$ (from top to bottom) and plot the results against values of $T/\mu$ ranging from the confined phase to the deconfined phase. We can clearly see that the Pfaffian phase fluctuations increase as the coupling increases and as the temperature decreases. Because the right-most green point for $g=0.001$ matches the corresponding point in Fig. 2 apart from the different gauge group, we can also check that the fluctuations increase from $1-\left\langle e^{i\phi}\right\rangle_{\text{pq}}=0.0113(14)$ for $N=4$ to $0.0501(84)$ for $N=8$. These Pfaffian phase fluctuations are likely related to the instabilities we observe for $N_{\tau}=8$ and $g\gtrsim 0.025$, motivating us to impose a minimum $N_{\tau}\geq 16$ for strong couplings $g\geq 0.001$ in our main calculations, to which we now turn.

We organize our lattice investigations of the deconfinement transition in terms of the dimensionless parameters $g$ and $T/\mu$ discussed above, Eq. (13). We consider four fixed values of the coupling that span three orders of magnitude: $g=10^{-5}$, $10^{-4}$, 0.001, and 0.01. In each case, we determine the transition temperature by scanning in the ratio $T/\mu$, generating $\mathcal{O}(10)$ ensembles that range from the confined phase to the deconfined phase. As discussed at the end of Sec. II, this differs from the approaches taken by prior studies Catterall and van Anders (2010); Asano et al. (2018); Bergner et al. (2022a), most of which also employ different lattice actions. To explore the thermodynamic limit where the number of colors $N\to\infty$, for each $g$ we repeat this procedure for up to three values of $N=8$, 12, and 16. To explore the continuum limit where the number of lattice sites $N_{\tau}\to\infty$, for each $(g,N)$ we consider two lattice sizes, either $N_{\tau}=8$ and $16$ for the two weaker couplings $g\leq 10^{-4}$ or $N_{\tau}=16$ and $24$ for the two stronger couplings $g\geq 0.001$. Table 1 in App. B summarizes these details. In total, these calculations involve 261 ensembles. Full information about them and the 35 additional SU(4) and SU(8) ensembles considered in the previous section is available in our open data release Jha et al. (2024).

To study the deconfinement transition of the BMN model, we focus on the Polyakov loop

$$\left\langle|PL|\right\rangle=\frac{1}{N}\left\langle\left|\text{Tr}\left[\prod_{n=0}^{N_{\tau}-1}U_{\tau}(n)\right]\right|\right\rangle\equiv\frac{1}{N}\left\langle|\text{Tr}\left[\mathbb{P}\right]|\right\rangle.$$

(22)

That is, by “Polyakov loop” we refer specifically to the magnitude of the trace of the holonomy along the temporal circle, with that holonomy itself corresponding to the $N\!\times\!N$ matrix $\mathbb{P}$. In the large-$N$ limit, $\left\langle|PL|\right\rangle$ is an order parameter for the spontaneous breaking of the $Z_{N}$ center symmetry, which vanishes in the confined phase while remaining non-zero in the deconfined phase.

This behavior is illustrated by the left-hand plots in Figs. 4 and 5, which show the Polyakov loop increasing as $T/\mu$ increases from the confined to the deconfined phase. Figure 4 compares all four values of the coupling $g$ for fixed $N_{\tau}=16$ and gauge group SU(12), while Fig. 5 compares the three gauge groups for fixed $N_{\tau}=24$ and $g=0.01$. (There is very little dependence on $N_{\tau}$ with $g$ and $N$ fixed, which can be seen from Table 1 in App. B.) The right-hand plots in these figures present the corresponding Polyakov loop susceptibility

$$\chi_{{}_{|PL|}}\equiv\left\langle|PL|^{2}\right\rangle-\left\langle|PL|\right\rangle^{2},$$

(23)

which exhibits a peak at the critical temperature $\left(T/\mu\right)_{\text{crit.}}$ of the deconfinement transition. By eye, Fig. 4 clearly shows the expected behavior, with $\left(T/\mu\right)_{\text{crit.}}$ moving from around the weak-coupling limit Eq. (10) towards the strong-coupling supergravity prediction Eq. (12) (marked by the vertical red line) as the coupling increases. Recall that the latter value corresponds to the large-$N$, $g\to\infty$ limit of the continuum theory and is not a strict bound on our finite-$N$ calculations.

We can obtain simple estimates for $\left(T/\mu\right)_{\text{crit.}}$ from the locations of the susceptibility peaks. These estimates are collected in Table 1, but are limited by the discrete values of $T/\mu$ we have analyzed. In principle, this can be improved by using multi-ensemble reweighting to interpolate between these values Kuramashi et al. (2020); Ferrenberg and Swendsen (1988). Here, we instead take a simpler approach of interpolating our results for the Polyakov loop itself. Specifically, we fit $\left\langle|PL|\right\rangle$ to the four-parameter sigmoid ansatz Schaich et al. (2020, 2022)

$$\Sigma=A-\frac{B}{1+\exp\left[C(T/\mu-D)\right]}.$$

(24)

All four parameters are necessary since our finite-$N$ Polyakov loop results are non-zero for all temperatures. We also need to omit the highest- and lowest-temperature points from the fits in order to keep the $\chi^{2}/\text{d.o.f.}$ under control. These fits correspond to the solid lines in the left-hand plots of Figs. 4 and 5, which extend only as far as the subset of points that are included.

The transition temperature $\left(T/\mu\right)_{\text{crit.}}$ is given by the fit parameter $D$ that corresponds to the inflection point where the sigmoid’s slope is maximized:

	$\displaystyle\frac{d^{2}\Sigma}{d(T/\mu)^{2}}\propto 1-\frac{2\exp\left[C(T/\mu-D)\right]}{1+\exp\left[C(T/\mu-D)\right]}$	$\displaystyle=0$
	$\displaystyle\implies C(T/\mu-D)$	$\displaystyle=0.$

As shown by Table 1 in App. B, our fits produce values for $\left(T/\mu\right)_{\text{crit.}}$ with small uncertainties, which are fully consistent with the susceptibility peaks. However, Table 1 also shows that the fits typically involve uncomfortably large $\chi^{2}/\text{d.o.f.}$, suggesting that these small uncertainties are likely underestimated. We therefore proceed by using the central values from the sigmoid fits but the more conservative uncertainties from the susceptibility peaks. This approach should account for systematic effects (e.g., from the choices of ansatz and fit range) in addition to purely statistical uncertainties.

The results we obtain in this way are used to produce the $\left(T/\mu\right)_{\text{crit.}}$ vs. $g$ phase diagram shown in Fig. 6. While we have only two lattice sizes $N_{\tau}$ for each $\{N,g\}$, from Table 1 in App. B we can see that all these pairs of $\left(T/\mu\right)_{\text{crit.}}$ results agree within uncertainties. We therefore average each pair by fitting to a constant, which leaves a single degree of freedom and produces very small $0.005\lesssim\mbox{$\chi^{2}/\text{d.o.f.}$}\lesssim 0.3$. Figure 6 shows these averages for each $\{N,g\}$, making it clear that there is also no visible $N$-dependence in our results.

At weak coupling $g\lesssim 10^{-4}$ we find critical temperatures in excellent agreement with NNLO perturbation theory in the $N\to\infty$ planar limit, Eq. (11). Because we use the perturbative expansion parameter $27g$ on the horizontal axis, these points appear at $27g\lesssim 0.003$. Perturbation theory breaks down around $27g\approx 0.05$ due to the relatively large coefficient $C_{\text{NNLO}}\approx 71$ in Eq. (11). Our non-perturbative numerical results monotonically increase as $g$ increases, approaching the large-$N$, strong-coupling supergravity prediction from Eq. (12), which is marked by a horizontal red line in Fig. 6. In fact, it is remarkable how close we get to this $g\to\infty$ limit given our modest $g\leq 0.01$. If we were able to reach stronger couplings, we might observe our finite-$N$ lattice calculations exceeding the large-$N$ supergravity prediction, as in Ref. Bergner et al. (2022a), but this will require future work to implement an improved lattice action. In parallel, our results also provide motivation for dual-supergravity calculations to attempt to determine the functional dependence of $\left(T/\mu\right)_{\text{crit.}}$ on $g$, to compare with the numerical results we have obtained.

Although we have now presented our main results for the BMN model phase diagram in Fig. 6, there are two more analyses we can carry out based on the Polyakov loop. First, we confirm that we are, in fact, dealing with the expected transition between confined and deconfined phases by considering the eigenvalues of the $N\!\times\!N$ matrix $\mathbb{P}$, Eq. (22). In the high-temperature deconfined phase, these $N$ eigenvalues for a given field configuration should all be aligned around one of the $Z_{N}$ vacua, with a localized distribution of complex phases. In the low-temperature confined phase, we should have instead a uniform distribution of eigenvalue phases around the unit circle Maldacena (1998b); Witten (1998); Aharony et al. (2004a). While it is possible to estimate $\left(T/\mu\right)_{\text{crit.}}$ (in the large-$N$ limit) by determining the temperature at which the eigenvalue phase distribution first spreads out to cover $[-\pi,\pi)$ with no gap remaining Aharony et al. (2004b, a), here we simply check systems on either side of the transition. Figure 7 shows a representative check, for $N_{\tau}=16$ and $g=0.01$, comparing SU($N$) gauge groups with $N=8$, $12$ and $16$. For $T/\mu=0.094$ below the transition, we see that the distributions become broader as $N$ increases, consistent with the uniform distribution expected for the confined phase in the large-$N$ limit. For $T/\mu=0.104$ above the transition, the distributions become more localized as $N$ increases, as expected for the deconfined phase.

Second and finally, we study the scaling of the Polyakov loop susceptibility peak $\chi_{\text{max}}$ with the number of degrees of freedom $\propto N^{2}$, which can distinguish between first-order and continuous transitions. Already in Fig. 5, we can see that the peak becomes higher as $N$ increases towards the $N\to\infty$ thermodynamic limit. In Fig. 8, we plot $\chi_{\text{max}}$ against $N^{2}$ on log–log axes for couplings $g=0.01$ (top) and $g=0.001$ (bottom), and fit the data to a power law,

$$\chi_{\text{max}}=CN^{2b},$$

(25)

where $C$ and $b$ are fit parameters. For the strongest coupling $g=0.01$, the critical exponent is $b=0.866(29)$ for $N_{\tau}=16$, increasing to $0.909(26)$ for $N_{\tau}=24$ (with only statistical uncertainties). These results are consistent with the first-order value $b=1$ Imry (1980); Fisher and Berker (1982); Binder and Landau (1984); Challa et al. (1986); Fukugita et al. (1989) in the continuum limit $N_{\tau}\to\infty$. The $g=0.001$ results are also close to one, although with no dependence on $N_{\tau}$ at all: $b=0.816(32)$ and $0.815(34)$ for $N_{\tau}=16$ and $24$, respectively. These fits leave a single degree of freedom and feature reasonable $0.2\lesssim\mbox{$\chi^{2}/\text{d.o.f.}$}\lesssim 1.7$.

Using perturbative calculations to compute an effective action, Ref. Hadizadeh et al. (2005) suggests that the BMN phase transition remains first order in the weakly coupled regime. For $g\leq 10^{-4}$, we can report only rougher results due to the smaller $N_{\tau}=8$ and $16$ we consider (cf. Table 1), each with only two gauge groups SU(8) and SU(12). Still, these $\left(T/\mu\right)_{\text{crit.}}$ allow us to estimate the critical exponent (neglecting uncertainties), and we observe that it now becomes smaller as $N_{\tau}$ increases towards the continuum limit. Specifically, for $g=10^{-4}$ we find that $b$ decreases from $0.85$ to $0.61$ as $N_{\tau}$ increases from $8$ to $16$. Similarly, for $g=10^{-5}$, $b\simeq 0.71$ for $N_{\tau}=8$ decreases to $0.59$ for $N_{\tau}=16$. These results suggest that although we have a first-order transition for strong couplings $g\gtrsim 0.001$, the transition appears to be continuous for weaker couplings $g\lesssim 10^{-4}$. Future work is needed to confirm this and search for a possible critical endpoint $g_{\star}$ to the line of strong-coupling first-order transitions.

We have presented results from our numerical Monte Carlo investigations of the BMN matrix model on the lattice, featuring our determination of the critical temperature of the model’s deconfinement phase transition for dimensionless couplings $g$ that span three orders of magnitude from the perturbative regime towards the domain of the dual supergravity, Fig. 6. Considering multiple SU($N$) gauge groups and lattice sizes $N_{\tau}$ for each $g$, we reproduce $\left(T/\mu\right)_{\text{crit.}}$ from perturbation theory for $g\lesssim 10^{-4}$ while approaching the large-$N$ holographic strong-coupling limit from Ref. Costa et al. (2015) for $g\gtrsim 0.001$. We have also analyzed the order of the transition, finding that the first-order transition for strong couplings appears to become continuous for weaker $g\lesssim 10^{-4}$.

While our results appear broadly consistent with previous lattice studies of the BMN model’s phase diagram Catterall and van Anders (2010); Asano et al. (2018); Bergner et al. (2022a), our approach differs by both our choice of lattice discretization as well as our strategy for analyzing the system. In particular, discretization-dependent instabilities in RHMC calculations restrict us to couplings $g\lesssim 0.025$ significantly smaller than those reached by Refs. Asano et al. (2018); Bergner et al. (2022a). We have discussed related discretization artifacts, including the phase of the fermion operator Pfaffian, which we observe to fluctuate significantly at strong coupling and large lattice spacing, Fig. 3. Despite this, a sign problem can be avoided by working with sufficiently small lattice spacings, which make the Pfaffian real and positive, Fig. 2.

Given the relatively small $g\leq 0.01$ we consider, it is remarkable how close our results in Fig. 6 get to the $g\to\infty$ limit from Ref. Costa et al. (2015). As discussed at the end of Sec. II, it is possible that finite-$N$ corrections may allow $\left(T/\mu\right)_{\text{crit.}}$ to exceed this strong-coupling limit Bergner et al. (2022a). This motivates further calculations with stronger couplings as well as larger values of $N$ in order to carry out better-controlled extrapolations to the $N^{2}\to\infty$ thermodynamic limit. From our investigations of discretization artifacts and the Pfaffian in Sec. III, we can appreciate that this will require either (or both) analyzing larger lattice sizes $N_{\tau}$ or switching to an improved lattice action. Larger $N_{\tau}$ will also enable better-controlled continuum extrapolations than were possible in this work. In addition, considering larger SU($N$) gauge groups is also important at weaker couplings to improve the scaling analyses that suggest the deconfinement transition becomes continuous rather than first order. If this observation is confirmed, it will be interesting to search for the critical endpoint $g_{\star}$ to the line of first-order transitions at stronger couplings.

In general, the BMN model remains less explored than the $\mu\to 0$ BFSS limit, and there are further investigations we may pursue using numerical lattice Monte Carlo calculations. These include computations of the internal energy of the system as a function of temperature Pateloudis et al. (2023), as well as analyses of additional ‘fuzzy sphere’ vacua Asano et al. (2018) or the effects of ‘ungauging’ the system Pateloudis et al. (2022).

We thank Simon Catterall and Toby Wiseman for helpful conversations, Denjoe O’Connor for useful discussions at the ECT* workshop “Quantum Gravity meets Lattice QFT”, and Yuhma Asano for sharing numerical data from Ref. Asano et al. (2018). DS was supported by UK Research and Innovation Future Leader Fellowship MR/S015418/1 & MR/X015157/1 and STFC grants ST/T000988/1 & ST/X000699/1. Additional support came from the International Centre for Theoretical Sciences Bangalore (ICTS) during visits to participate in the two editions of the program “Nonperturbative and Numerical Approaches to Quantum Gravity, String Theory and Holography” (ICTS/Prog-NUMSTRINGS2018/01, ICTS/numstrings-2022/8); we thank ICTS for its hospitality. Numerical calculations were carried out at the University of Liverpool, on USQCD facilities at Fermilab funded by the U.S. Department of Energy, and at the San Diego Computing Center and Pittsburgh Supercomputing Center through XSEDE supported by National Science Foundation grant number ACI-1548562.