The BFSS model on the lattice

We study the maximally supersymmetric BFSS model at finite temperature and its bosonic relative. For the bosonic model in p+1 dimensions, we find that it effectively reduces to a system of gauged Gaussian matrix models. The effective model captures the low temperature regime of the model including one of its two phase transitions. The mass becomes p1/3λ1/3 for large p, with λ the ’t Hooft coupling. Simulations of the bosonic-BFSS model with p = 9 give m = (1.965±.007)λ1/3, which is also the mass gap of the Hamiltonian. We argue that there is no ‘sign’ problem in the maximally supersymmetric BFSS model and perform detailed simulations of several observables finding excellent agreement with AdS/CFT predictions when 1/α′ corrections are included.


Introduction
The leading candidate for a non-perturbative formulation of M-theory is expected to be the infinite matrix size limit of a matrix model of some kind. One such proposal is the BFSS model [1][2][3] which was conjectured to capture the entire dynamics of M-theory. 1 Relatives of this model such as the BMN model [4] or models derived from the ABJM model [5] are also considered possible viable candidates for such a non-perturbative formulation. All of these conjectured formulations of M-theory are regularised versions of the supermembrane and are matrix quantum mechanical systems. They are based on the matrix regularisation of membranes introduced by Hoppe [6] and extended to the supermembrane in [1] and [7].
In the second half of the paper we focus on the BFSS model. However, we dedicate the earlier sections of the paper to a careful study of the Hoppe regulated bosonic membrane which is the bosonic part of the BFSS model. This model is also of separate interest since it is the high temperature limit of maximally supersymmetric 1+1 dimensional SU(N ) Yang-Mills (or matrix string theory) compactified on a circle when the fermions decouple due to their anti-periodic boundary conditions at finite temperature. The model has been studied already both theoretically [8] (see also ref. [9]) and numerically [10,11] and our results are in accord with these studies. The zero temperature supersymmetric gauge theory has been studied in a hamiltonian light cone context in ref. [12].

JHEP05(2016)167
We then continue our study with the BFSS model. This model first emerged as N = 16 supersymmetric Yang-Mills quantum mechanics [13][14][15], later it arose as the 11-dimensional supermembrane in Hoppe's regularization and most recently as the BFSS model [2], a candidate for a non-perturbative formulation of M-theory, and it also describes the low energy effective theory of a system of D0-branes [16]. A lattice version of the model which preserves eight of the sixteen supersymmetries was presented in [17].
Our lattice regularisation of the model follows Catterall and Wiseman [18]. In this formulation it is clear that the Pfaffian obtained when the Fermions are integrated out can in general be complex. However, we find that the phase of the Pfaffian is restricted to a narrow band so that cos(Θ pf ) ∼ 1. There is therefore no real phase or 'sign problem' as far as simulations of the model are concerned. We simulate the model using a rational hybrid monte carlo algorithm (RHMC) and find excellent agreement with results reported in [18][19][20][21] though our values for the energy are slightly higher than those of [19,21] but in excellent agreement with predictions of AdS/CFT when leading 1/α corrections are included. Those interested primarily in the supersymmetric model can skip to section 4 for our discussion of the model and results, consulting section 2 for the continuum model and our notation.
The principal results of this paper are: • We perform monte carlo simulations of the bosonic BFSS model and measure the two point correlation function, the mass gap and the eigenvalue distribution of each matrix. All fit with Gaussian matrix quantum mechanics with the same mass as that found from the correlation function.
• We derive an effective description of the bosonic model using a 1/p expansion where p is the number of matrices. The description is in terms of p massive scalar fields that are gauged under the adjoint action of SU(N ) but are otherwise free scalar fields 2 • The effective model has a phase transition in the same region as the full model (see [22]). 3 • We study the maximally supersymmetric BFSS model and present arguments showing that though the lattice model in general has a complex phase, it is only the cosine of this phase that enters in simulations and the model is such that the angle is restricted to regions where the cosine is positive hence there is no sign problem in the full model.
• We simulate the model using a Fourier accelerated rational hybrid monte carlo algorithm confirming results found earlier by other groups and find excellent agreement with predictions of AdS/CFT when subleading 1/α corrections are included. 2 The relevant results here were previously obtained in [8] which we became aware of these after writing this article. 3 The model can also be interpreted as the high temperature description of the maximally supersymmetric 1 + 1 dimensional Yang-Mills compactified on a circle. In this setting its two transitions are the high temperature limit of the black hole black string transition in the dual gravity model [24].

JHEP05(2016)167
The simulations we report provide a useful test of our code as we proceed to examine further observables and the inclusion of longitudinal M5-branes or equivalently D4-branes. Our goal being the Berkooz-Douglas matrix model holographically dual to the backreacted D0/D4 brane system. Such studies can provide a solid test of the AdS/CFT correspondence with flavour degrees of freedom, which is a widely used tool for non-perturbative analysis in flavour dynamics.
The layout of the paper is as follows: in section 2 we introduce the BFSS model and continue in section 3 with a study of the bosonic part of the model describing our lattice discretisation and hybrid monte carlo algorithm. We then present our numerical results for this bosonic model and continue in section 3.5 to develop an expansion for the model in terms of the inverse of the number of matrices. This leads us to introduce our effective Gaussian model and compare it with the full model. In section 4 we return to the supersymmetric model and present our lattice discretisation of this model. We give a discussion of the phase of the Pfaffian in the model and present our results. The paper ends with a discussion of our results in section 5.

The BFSS model
The BFSS matrix model is a one dimensional supersymmetric Yang-Mills theory naturally arising in type IIA superstring theory as a low energy effective description of D0-branes. It is conjectured that in the limit of a large number of D0-branes, N , the model is equivalent to uncompactified eleven dimensional M -theory [2] while for finite N it corresponds to Mtheory compactified on a light-like circle [3]. The easiest way to obtain the BFSS matrix model is via dimensional reduction of ten dimensional supersymmetric Yang-Mills theory down to one dimension. The resulting reduced ten dimensional action is given by [23]: where Ψ is a thirty two component Majorana-Weyl spinor, Γ µ are ten dimensional gamma matrices and C 10 is the charge conjugation matrix satisfying C 10 Γ µ C −1 10 = −Γ µT . We take a representation for Γ µ in terms of nine dimensional (euclidian) gamma matrices γ i : where C 9 is the charge conjugation matrix in nine dimensions satisfying C 9 γ i C −1 9 = γ i T and σ i are the Pauli matrices. With the following choice for the Majorana-Weyl spinor: where ψ is a sixteen component Spin(9) Majorana fermion and Wick rotating to Euclidean time, we obtain the Euclidean action:

JHEP05(2016)167
which, as is manifest, is Spin(9) invariant. Note: we have not imposed any restriction on the nine dimensional spinor basis. For example if we choose γ i to be in the Majorana representation (where the nine γ i are taken to be real and symmetric) then C 9 = 1 16 and we arrive at the more standard form for the action (2.4). However, we are interested in a basis in which the discrete theory is free of fermion doubling, which can be achieved by taking a basis [18] in which C 9 = 1 8 ⊗ σ 1 . We will return to this in section 4, where we consider the discretisation of the full BFSS matrix model.

Bosonic BFSS on the lattice
In this section we focus on the bosonic part of the action (2.4) given by: The zero temperature model was introduced by Hoppe [6] as a gauge fixed and regulated description of membranes. The model has also been extensively studied in the literature. 4 It has been simulated for a first time in refs. [24,25], where certain aspects of the full model were described in terms of simple Gaussian model. Its phase structure at finite temperature has been explored in [8,10,22], where the authors found that as the temperature is decreased the model first undergoes what is probably a 3rd order deconfining-confining phase transition into a phase with non-uniform but gapless distribution for the holonomy. As the temperature is further decreased there is a 2nd order transition to a gapped holonomy with a quadratic decrease in the internal energy to a constant value for lower temperatures. The high temperature expansion of the model was developed in [26]. In what follows we will reproduce the main results of [8,22] and elaborate on the properties of the theory at zero temperature. In particular we will provide evidence that the low temperature phase of the model has an effective description in terms of a free massive scalar which captures many of the finite temperature features of the model including one of its two phase transitions.

Discretisation
We discretise time to Λ sites t n = a n, (n = 0, . . . , Λ − 1), where the lattice spacing is a = β/Λ, and impose periodic boundary conditions identifying the point t Λ = Λa = β with the point 0. To discretise the covariant derivative D t we define the transporter fields: where P denotes a path ordered product. Let us consider for a moment the pure derivative part of D t . On the lattice we have: The model is also the high temperature limit of 1 + 1 dimensional N = 8 supersymmetric Yang-Mills on R × S 1 where β is now the period of the S 1 and the fermions drop out due to their anti-periodic boundary conditions at finite temperature.

JHEP05(2016)167
To make the above expression gauge covariant we have to transport back the field at t n+1 to t n . For the discrete version of the covariant derivative, we obtain: where U n+1,n = U † n,n+1 . Using equation (3.4) for the discrete bosonic action we obtain: where without loss of generality we have taken g = 1 √ N . 5 The action S b can be written in a much simpler form by using the U(n) Λ gauge symmetry of the model. Indeed, at each lattice site we have a local U(N ) symmetry. Using that symmetry we can perform the transformation: for the bosonic action (3.5) we obtain: (3.7) The unitary matrix W has the decomposition W = V DV † , where D = diag{e iθ 1 , . . . , e iθ N } is a diagonal unitary matrix and V is a unitary. But the action (3.7) has the residual symmetry X i n → V X i n V † , which we can use to diagonalise W. Furthermore, it has an additional symmetry X i n → h n X i n h † n , where h n is a diagonal unitary matrix, which we can use to "distribute" the diagonal matrix D among all of the hop terms. Indeed, defining the matrix D Λ = diag{e iθ 1 /Λ , . . . , e iθ N /Λ }, which satisfies (D Λ ) Λ = D, one can verify that under the transformation: the action (3.7) transforms into:

JHEP05(2016)167
Now let us discuss the measure of the transporter fields U n,n+1 . The measure can be written as: where we have used that U 0,1 = W (U 1,2 . . . U Λ−2,Λ−1 ) † and the invariance of the measure. But the action (3.9) depends only on the matrix W (in fact only on the eigenvalues of W). Therefore the integration over the measure of the transporter fields reduces to: (3.12)

Hybrid Monte Carlo
To implement the Hybrid Monte Carlo algorithm it is convenient to work in a gauge in which the holonomy matrix is non-trivial at only one link (we choose the one connecting sites zero and Λ − 1). To this end we omit the diagonal matrices h n in the transformation (3.8). The action (3.9) is then given by: (3.13) The corresponding Hamiltonian for the molecular dynamics part of the HMC algorithm is then: where P i n are the canonical momenta corresponding to the hermitian matricesX i n , and p l d are the canonical momenta associated to the angles θ l . Hamilton's equations read:

JHEP05(2016)167
where dots denote derivatives are with respect to the Monte Carlo time. Using equation (3.13), for the corresponding forces we obtain: which we implement in the Hybrid Monte Carlo.

Phase structure
In this section we reproduce the studies of the phase structure of the bosonic model originally obtained in [22]. The main observables that we focus on are the internal energy of the system E, the "extent of space" R 2 and the Polyakov loop P defined via: where the holonomy matrix, U , is the continuum limit of the link variable U 0,Λ defined in equation (3.2). The expectation value of the Polyakov loop |P | plays the rôle of an order parameter for the confining-deconfining phase transition discussed in [22]. In figure 1 we presented a plot of this order parameter as a function of the temperature. The plot is for N = 16 and lattice spacing a ≈ 0.05. One can see that near temperature T ≈ 0.90 there is a phase transition. The change of the slope of the curves and the fluctuations in the simulations near T ≈ 0.9 is consistent with the existence of a second order phase transition.
In figure 2 we present plots of the energy and "extent of space" as functions of the temperature, for N = 16 and lattice spacing a ≈ 0.05. The dashed curves represent the analytical high temperature behaviour obtained in [26]. Our results agree very well with the corresponding studies in [8,10,22].
A more detailed analysis of the temperature range close to the transition revealed that there are in fact two transitions. To uncover more detail on the nature of the phase transition the authors of [22] analysed the distribution of the holonomy matrix near the  Figure 2. Plots of the scaled energy E/N 2 and the "extent of space" R 2 as functions of the temperature. The dashed curves correspond to the high temperature behaviour obtained in [26].
One can see that near T ≈ 0.9 the plots suggest the existence of a second order phase transition. The energy and temperature in the plots are in units of λ 1/3 . phase transition and uncovered behaviour consistent with the Gross-Witten model [27] and concluded that the holonomy undergoes a transition from a uniform distribution at T c1 0.8758(9) to a gapped distribution at T c2 = 0.905(2). In figure 3 we present our plots of the distribution of the holonomy around the phase transition. The dashed curves in the plot represent fits with the gapped and ungapped forms of the Gross-Witten distribution which are in excellent agreement with those of [8,10,22]. We have not attempted to refine their results, rather our purpose is to uncover a hidden Gaussian structure in the model.

Gap and eigenvalue distribution
In this section we investigate the eigenvalue distribution of the scalar fields. We also perform studies of the mass of the theory at zero temperature. Our results suggest that at JHEP05(2016)167 all temperatures the eigenvalue distribution of any one of the X i is given by a Wigner semicircle, with a radius R λ following the temperature behaviour of the observable R 2 . 6 Therefore, we conclude that while the radius of the distribution experiences a phase transition the shape of the distribution remains unchanged. We believe that the main reason for this behaviour is that for nine scalar fields the commutator squared term can be replaced and approximated by a quadratic mass term in the spirit of [28], where an expansion at large number of scalar fields has been developed. Generalising these techniques, we are able to obtain an estimate of the mass, which agrees very well with both the gap measured from correlation functions and the radius of the distribution which are obtained from Monte Carlo simulations.
In the limit of high temperature the model reduces to the 10-dimensional Yang-Mills matrix model considered in [28]. The obtained behaviour of the radius is in agreement with the large temperature expansion performed in [26] and provides an analytic approximation to the leading coefficients in this expansion.
We begin by considering the model at zero temperature. In this case the holonomy can be completely gauged away and the model simplifies. Furthermore, at zero temperature the correlator: captures the gap m = E 1 − E 0 of the theory. To calculate the gap in the discrete theory, we periodically identify the time direction with period β: Note that although formally β is the same parameter that we have at finite temperature, since we set the holonomy to zero here its meaning is just a periodic coordinate as opposed to inverse temperature. Our result for the correlator for N = 30, β = 10 and lattice spacing a = 0.25 is presented on the left in figure 4. The fitting curve is given by equation (3.22) and when we perform a two parameter fit we obtain A ≈ 7.50 ± 0.2 and m ≈ (1.90 ± .01) λ 1/3 . However, for Gaussian scalar fields of mass m we have A = N 2m(1−e −βm ) and performing a one parameter fit for m yields m = 1.965 ± 0.007 and A = 7.63 ± 0.03. On the right we have presented a plot of the eigenvalue distribution of one of the matrices for the same parameters. The fitting curve represents a Wigner semicircle of radius R λ ≈ 1.01. The fact that the theory is gapped and that the eigenvalue distribution is a semicircle suggests that that at low temperate the model has an effective action: for each of the matrices X i . It is well known [29] that for the action (3.23) the eigenvalue distribution of X is given by a Wigner semicircle of radius:

1/D expansion
Adapting the techniques developed in [28] to the time dependent case that we are considering we can obtain a theoretical estimate of the radius at low temperature. 7 Let us consider again the action of the bosonic model (3.1): where we have rescaled so that β = λ 1/3 T (effectively we set g = N −1/2 ). The commutator square term can be written as: where t a are SU(N ) generators normalised to Tr t a t b = δ ab and the tensor λ abcd is given by: It is convenient also to define the inverse kernel of λ abcd satisfying: We will also use the identities: The action (3.26) can then be written as: Our next step is to add to the action the term: 8 the action S b = S b + ∆S then becomes: Next we define: k ij,lm ≡ t a ij k ab t b lm , (3.34) 7 The next order corrections for the current model were computed in [8] but we were unaware of this work at the time of writing. 8 Note that we can always add this term since Dk e −∆S = const.

JHEP05(2016)167
From the definition of k ij,lm it follows that it is traceless with respect to the first and the second pair of indices and we can invert: k ab = t a ij k ji,ml t b lm . Substituting in the action (3.33), Fourier transforming (via X = n e i 2πn β tX n ) and assuming k ab is time independent we obtain: where we have defined the projector on traceless matrices: and assumed that k is constant. Defining also the double index matrix: we can integrate out the X's: The effective action for the field k then becomes: We now notice that the first term in the expression for the matrix W (3.37) commutes with the projector P . It is natural then to consider an ansatz for k which also has this property. Thus we consider: k ij,lm ≡ k ij P ij,lm = P ij,lm k lm (3.40) The last equality is possible only if all diagonal components of k ij are the same (namely k ii = k jj for all i, j) we also choose k ij to be symmetric. Then: and we have for all powers r that (P.W (n).P ) r ij,lm = P ij,lm (∆ ij (n)) r . Therefore, The corresponding saddle point equation S eff [k] = 0 becomes:

JHEP05(2016)167
We can now sum the series to obtain: In principle we can solve for k ij in terms of θ i − θ j , but we will restrict ourselves to extracting the low temperature dependence. The first term in equation (3.44) has the following expansion: One can see that the effect of the holonomy is exponentially suppressed at low temperature.
To leading order we can then consider a symmetric ansatz k ij = k 0 , which also implies k ab = k 0 δ ab . We obtain: Substituting into the action (3.33) to leading order we obtain: and the corresponding radius of the eigenvalue distribution is: This result agrees within a few percent with R λ ≈ 1.01 obtained by fitting the distribution in figure 4. The exponential suppression of the holonomy corrections to the low temperature saddle point for k ab suggest that a lot of the physics of the model (at least at low temperature) can be captured by the action: where we restored the gauge field in the covariant derivatives. Surprisingly this model knows lots about the phase transition of the full model. An analysis similar to that above shows that the critical temperature for p adjoint gauged Gaussian matrices as in [9] with mass m occurs at T Gaussian c = m ln p which for m = 1.965 yields T c = 0.8943 and for m = p 1/3 yields T c = 0.9467. In figure 5 we have presented our results for the energy and the Polyakov loop. Note that in this approximation the scaled energy E/N 2 is equal to the extent of space R 2 . The plots are for m = 9 1/3 , N = 32 and lattice spacing a ≈ 0.05. One can see that both the energy and the Polyakov loop exhibit the same behaviour as for the full model. There seems to be a phase transition taking place at T ≈ .95 appearing to be of a second order. It is slightly shifted towards high temperatures relative to the critical temperature for the full model. If instead of mass m = p 1/3 λ 1/3 we use the value m (1.90 ± .01)λ 1/3 the phase transition is shifted in the opposite direction (just bellow T = 0.9). This indicates that if one fits the mass parameter one can improve even further the agreement of the gauged gaussian model and the full one. The dashed curve in the second plot is the theoretical prediction for the high-temperature behaviour of the Polyakov loop, again one can observe a very good agreement. The high temperature behaviour of the energy on the other side disagrees with the corresponding behaviour of the full model. This is not surprising since we derived the effective action at low temperature and the behaviour at high temperature is dominated by the highest power of the potential which has been changed from quartic to quadratic. One can also see that at low temperature the energy remains constant. In figure 6 we have presented a plot of the energy versus the temperature zoomed in near the phase transition. The dashed curve represents a fit with: with parameter T 0 ≈ 0.896. Note that although this seems to indicate the presence of a separate third order phase transition, there is not enough evidence to support that. It is also possible that two apparent phase transitions signal the presence of a single first order phase transition, numerically it is difficult to exclude such a possibility. One can perform the large p analysis (see [28]) in the high temperature limit where the model now has p + 1 matrices (the holonomy becomes the additional matrix) and predict that the model in this limit again becomes Gaussian but now the field k ab includes the holonomy and the saddle becomes k ab = √ p + 1δ ab which again predicts a Gaussian matrix model with a high temperature effective action S eff = √ p+1 2 T r((X i ) 2 ) and consequently a Wigner semicircle distribution for the eigenvalues of X i . We conclude this section by presenting results for the eigenvalue distribution of X i for the gauged gaussian model at finite temperature. In figure 7   show a Wigner semicircle. One can see that the shape of the eigenvalue distribution does not change with temperature.

The supersymmetric model on the lattice
In this section we consider the full supersymmetric BFSS model on the lattice. The model has been simulated using a non-lattice approach in [19] and using lattice discretisation in [18] and [21]. The main goal of these studies has been to compare the low temperature regime of the model to the holographically dual black hole geometry. The authors of refs. [19] and [21] also compared the high temperature regime of the model with the high temperature expansion performed in [26]. Our goal is to reproduce some of these studies and calibrate our code.

JHEP05(2016)167
A naive discretisation of the action (2.4) would result in a fermion doubling. This can be avoided [18] if the charge conjugation matrix is taken to be C 9 = 1 8 ⊗ σ 1 . 9 Constructing a basis for which C 9 is of this form is relatively straightforward. For example one can tensor up the Majorana basis in seven euclidean dimensionsγ a E : γ a = −γ a E ⊗ σ 3 , for a = 1, . . . , 7 , γ 8 = 1 8 ⊗ σ 2 , and verify that indeed C 9 is of the desired form (it also satisfies C 9 = γ 9 ). We next proceed in discretising the action (2.4). Since the bosonic part of the action is identical to the one considered in section 3, we will consider only the fermionic part of the action: We begin by splitting the fermions into two eight component fermions: ψ = (ψ 1 , ψ 2 ) and defining the forward and backward derivatives D ± : One can show that the discretised kinetic term then becomes: where the plus/minus sign in the last term corresponds to periodic/anti-periodic boundary conditions for the fermions. 10 Using the gauge from the previous subsection when the holonomy is concentrated on a single link we can write S kin f as: Since all fields transform in the adjoint of SU(N ) instead of dealing with matrices we can use the corresponding real components: X a = tr(t a X) and ψ a = tr(t a ψ), where t a are the standard basis of SU(N ) (introduced earlier) normalised as tr t a t b = δ ab . S kin f can then be written as: K ab mn = (δ m+1,n − δ m,n )δ ab ± δ m,Λ−1 δ n,0 d ab (4.7)

JHEP05(2016)167
where the plus/minus sign corresponds to periodic/ant-periodic boundary conditions. The kinetic term can also be written as: Finally, defining: M mn,αβ,ab = 1 2g 2 We can write: Observe that M is the sum of an anti-hermitian kinetic term and a hermitian potential and M † (X) = −M(−X). Also because the Bosonic action is symmetric in X the Pfaffian in the partition function can be replaced by |Pf(M)| cos(Θ P f ). Now as long as − π 2 < Θ P f < π 2 the cosine is positive and the effective action defines a true probability distribution given by In figure 8 we present a plot of the phase of the pfaffian of the fermionic matrix for N = 3 and four lattice points. 11 One can see that the cosine remains positive. We believe that the drop in the cos θ curve at very low temperatures is a lattice effect and is not present in the continuum limit. Our results show a very good agreement with the earlier studies of ref. [30]. 11 Note that to control the flat directions at low temperature we have added a small mass for the bosonic field. One can see that the phase remains small for all temperature, but drops at very low temperatures. We believe that this is a lattice artifact and is not present in the continuum limit. The flat directions are controlled with a small mass regulator at low temperature.

RHMC and fermionic forces
The next step is to apply the RHMC method [31] to the model. To this end we need the so called fermionic forces. Let us summarise briefly the philosophy.
As we have shown above the model does not suffer from a severe sign problem and so we ignore the phase of the Pfaffian and use that: Here ξ is a 16(N 2 c − 1)Λ dimensional vector consisting of the pseudo-fermionic fields. The idea of the RHMC is to approximate the rational exponent of the matrix M † M with a partial sum: where the parameters α 0 , α i , β i and # depend on the rational exponent δ, the spectral range of the matrix M † M and the desired accuracy. We will need two rational exponents.
To update the pseudo fermions we use that the field η ≡ (M † M) −1/8 ξ has a gaussian distribution and solve for ξ = (M † M) 1/8 η using a multi-shift solver. Therefore, δ = 1/8 is one of the rational exponents that we need. To calculate the fermionic forces and the contribution to the hamiltonian we need to invert (M † M) −1/4 and the second exponent is δ = −1/4.

JHEP05(2016)167
Let us elaborate on the computation of the fermionic forces. We have two type of forces: gradients with respect to X n ij and gradients with respect to the phases of the links θ i . Using the partial expansion (4.20), one can easily derive an expression for the derivatives of the fermionic action: where h i satisfy (M † M + β i ) h i = ξ i and are obtained from the multi-solver.

Simulation results
In this subsection we provide our results from the Monte Carlo simulation of the model. We focus on the same observables that we analysed for the bosonic model in section 3.3. The definitions of the extent of space R 2 and the expectation value of the Polyakov loop P remain the same. The expression for the internal energy is [18]: 12 We have simulated the following configurations. For temperatures T > 2 we have used N = 8 and Λ = 8. For the region 1 ≤ T ≤ 2 we have used N = 8 and two or three different sizes of the lattice (for each point) in the range 8 ≤ Λ ≤ 16. (For T = 1 we also went to Λ = 32). For temperatures lower than one we have used N = 10, 12, 14 and two lattice sizes per point Λ = 8, 16. For all temperatures the Polyakov loop is largely independent on the lattice spacing. The extent of space also experiences very weak lattice effects. However, this is not the case for the internal energy and even for temperatures as high as T = 2 lattice effects can be a factor. In figure 9 we present our results for the energy at T = .9, 1.0 for different lattice spacing. One can see that the lattice effects die out linearly, which allows us to extrapolate the energy to zero lattice spacing. Our results for the internal energy are summarised in figure 10. The dashed curve at high temperatures is the theoretical curve obtained in the high temperature expansion. The dashed curve at low temperature represents the prediction of gauge/gravity. 13 The model becomes unstable for small matrix sizes N , an effect which has been related to Hawking radiation in the dual gravitational theory [30,32]. To compare with the gauge/gravity predictions one needs to consider large matrices. Simulations with large matrix sizes is computationally expensive and as a result our low temperature data is still preliminary. However, even at this point we have excellent overall agreement with the studies of [19] and [21]. 12 Note that this expression is also valid for the bosonic model. The result can be obtained by rescaling the fields such that the kinetic term is temperature independent and removes any temperature dependence from the measure (the Van Vleck Morette determinant generically depends on temperature). Then differentiating with respect to temperature and using the Ward identities associated with the total number of degrees of freedom yields this result. 13 We have used the α corrected expression 1  Figure 10. Results for the internal energy for 8 ≤ N ≤ 14 and 8 ≤ Λ ≤ 16. The dashed curve at high temperature corresponds to the theoretical results of [26], while the low temperature curve represents the prediction for the internal energy from the AdS/CFT correspondence.
Finally, in figure 11 we present our results for the Polyakov loop and the extent of space. Again the dashed curve represents the high temperature theoretical result obtained in [26]. One can see the excellent agreement at high temperatures. Our result for these observables agree with the previous studies preformed in [18][19][20] and [21].  Figure 11. Plots of the expectation value of the Polyakov loop |P | and the extent of space R 2 as functions of temperature. The dashed curves represent the predictions of the high temperature expansion.

Discussion
In this paper we have analysed both the purely Bosonic and the supersymmetric BFSS models. These are "Hoppe" regulated membranes and one would expect that in the large N limit when the regulator is removed that they describe the full quantum dynamics of these membranes. Surprisingly we found that the bosonic model, for sufficiently large embedding dimension reduces to a system of p-massive free bosons with the mass given by m ∼ p 1/3 . For p = 9 we performed detailed simulations of the model evaluating both the fall off of the correlation function and the eigenvalue distribution of the X i fit with a Wigner semi-circle both of which give a consistent mass m (1.965 ± .007)λ 1/3 . This is a fundamental non-perturbative result and gives the mass gap in the full Hamiltonian of the model.
The correspondence of the full model and our gauged Gaussian model is excellent for a wide range of temperatures. Somewhat surprisingly the phase transition region of the full model is faithfully reproduced by the effective model with the two transitions of the full model merged into one. Since the finite temperature version of the model is also the high temperature limit of 1 + 1 dimensional maximally supersymmetric Yang-Mills [24] compactified on a circle, we have established that this latter model should also reduce to a system of free massive scalars in its large radius high temperature phase.
We then study the full supersymmetric BFSS model using a rational hybrid monte carlo simulation with Fourier acceleration to evaluate observables of the model. After describing our lattice discretisation of the model we investigated the phase of the Pfaffian obtained on integrating out the Fermions. The Pfaffian is generically complex, however, its phase is in fact not a problem for simulations. What enters simulations is the cosine of this phase and in the regularisation used in our work this phase is in fact restricted to a region where the cosine is positive once the lattice spacing is sufficiently small. Direct measurements confirm that the phase is indeed small.
Though our results for this part of the paper do not yet go beyond those of [19] or cover as low a temperature as those of [21] they are more precise than those of Catterall

JHEP05(2016)167
and Wiseman [18] who use a similar lattice simulation. We have taken several lattice spacings and then performed an extrapolation to the limit of zero lattice spacing. We find good agreement with earlier results and excellent agreement with the predictions of gauge/gravity once 1/α corrections are included. Our results appear to approach the predictions of gauge/gravity a little more closely than those of [19] but the difference is broadly within the errors. The principal difficulty in simulating the model at lower temperatures is due to critical slowing down and though Fourier acceleration helped for T 0.5, simulations become increasingly more difficult at lower temperatures. An additional difficulty is the instability due to flat directions, which require increasingly larger matrix sizes as the temperature is reduced.
One of the principal aims of this work was to check the claims of previous work and in particular those on the absence of a complex phase problem. We were also interested in calibrating our code as we extend it to include systems with D4-branes. The extension to such systems will allow us to perform more extensive tests of the gauge/gravity correspondence.