Thermodynamics of the BMN matrix model at strong coupling

We construct the black hole geometry dual to the deconfined phase of the BMN matrix model at strong 't Hooft coupling. We approach this solution from the limit of large temperature where it is approximately that of the non-extremal D0-brane geometry with a spherical $S^8$ horizon. This geometry preserves the $SO(9)$ symmetry of the matrix model trivial vacuum. As the temperature decreases the horizon becomes deformed and breaks the $SO(9)$ to the $SO(6)\times SO(3)$ symmetry of the matrix model. When the black hole free energy crosses zero the system undergoes a phase transition to the confined phase described by a Lin-Maldacena geometry. We determine this critical temperature, whose computation is also within reach of Monte Carlo simulations of the matrix model.


Introduction
Some quantum mechanical systems admit a parametric limit in which they are well described by a classical gravitational theory. Such systems are examples of quantum theories of gravity.
It is not easy to find a system with this property but the gauge/gravity duality offers several cases [1,2].
Perhaps, the most striking example is (0 + 1)-dimensional SU (N ) Super Yang-Mills (SYM) theory. This theory contains a finite number of bosonic and fermionic degrees of freedom which are naturally organized in N by N traceless hermitian matrices X i and Ψ α , respectively. This model is often termed BFSS [3], with action given by where D t = ∂ t − i[A, ] is the covariant derivative and summation over spatial indices i, j = 1, . . . , 9 and spinor indices α, β = 1, . . . , 16 is implicit. By dimensional analysis, one concludes that the 't Hooft coupling λ has units of energy cubed. Therefore, the thermodynamics of this system is controlled by two dimensionless parameters: N and τ = T /λ 1 3 , where T is the temperature. According to the gauge/gravity duality, at large N and small dimensionless temperature τ this theory is dual to 11-dimensional supergravity in the following black hole geometry 1 where the c i are numerical coefficients. The leading term follows from the classical black hole thermodynamics of (2), which gives c 1 = − 1 21 120π 2 49 7/5 . The coefficient c 2 is not known analytically because it follows from unknown α 3 corrections to type IIA effective action [4]. The 1/N 2 terms correspond to quantum corrections associated to string loops in the 10 dimensional picture. Notably, the coefficient c 3 has been computed recently using the quartic curvature corrections to 11D supergravity [5]. It is an outstanding challenge to reproduce analitically these predictions directly from the matrix quantum mechanics (1). In fact, the state of the art is a scaling hypothesis for the several terms of the moduli effective action that correctly predicts the leading low-temperature dependence τ 14 5 but it is unable to fix the coefficient c 1 and any of the subleading terms [6,7]. The mean field approximation of [8,9,10] claimed partial success in reproducing the gravity prediction but their numerical method breaks down for sufficiently low temperature [11]. 1 This follows from the decoupling limit of N coincident D0-branes in type IIA supergravity [2], which leads to a charged spherically symmetric black hole in ten dimensions. From our point of view, it is more convenient to work with the solution uplifted to eleven dimensions, where it is purely geometric and describes a black string with horizon topology S 1 × S 8 . In our conventions, the 11-dimensional Newton constant is given by 16πG N = (2π) 8 g 3 s 9 s and the periodic coordinate z obeys z ∼ z + 2πg s s . 2 We are assuming τ N − 5 9 . For lower temperatures, the black string suffers from the Gregory-Laflamme instability and the stable black hole should have S 9 horizon topology. For even lower temperatures τ ∼ N − 5 6 , the curvature at the horizon of (2) reaches the Planck scale.
In a remarkable series of papers [12,13,14,4,15,16,17], the authors performed Monte-Carlo simulations of the matrix quantum mechanics (1) at finite temperature. In [4], they studied the planar limit (N → ∞) at low temperature and obtained the first two terms in equation (4). Their results agree with all available analytical results from the gravity dual and provide a prediction for c 2 . More recently [17], a study of 1/N 2 effects confirmed the gravitational prediction for the coefficient c 3 . This is among the most impressive tests of the gauge/gravity duality we are aware of. Notice that this includes quantum gravity loop effects and probes the regime of chaotic dynamics where supersymmetry is completely broken and integrability is absent.
In this paper, we study the thermodynamics of a massive deformation of the matrix quantum mechanics (1). This model goes by the name of PWMM (plane wave matrix model) or BMN after the authors of [18]. Its action reads where the indices i, j, k run over 1, 2, 3 and the index a runs over 4, . . . , 9. The matrix γ 123 is equal to 1 6 ijk γ i γ j γ k , with ijk the standard 3 dimensional -tensor. This means that the mass parameter µ breaks the SO(9) global symmetry of (1) down to SO(6) × SO (3). The deformation also retains maximal supersymmetry [18].
The BMN model has three significant advantages over the BFSS model. The first is that it has a discrete energy spectrum and a well defined canonical ensemble. Notice that, strictly speaking, the canonical ensemble of the matrix quantum mechanics (1) does not exist. 3 The reason for this is that the eigenvalues of commuting matrices X i can be made arbitrarily large without energy cost. In fact, the Monte-Carlo simulations work because there is a meta-stable thermal equilibrium with a decay rate that is very small at large N . The second advantage is that the BMN model has a dimensionless coupling constant g ≡ λ/µ 3 that, together with the dimensionless temperature T /µ, parametrize a two-dimensional phase diagram. This means that we can use the dual gravitational description at large N and strong coupling g 1, to predict many observables as functions of the dimensionless temperature T /µ. Finally, the third advantage is that the BMN model is expected to have a phase transition whose critical temperature should be easy to measure in Monte-Carlo simulations. 4 In figure 1 we depict the phase diagram of the theory in the planar limit N → ∞. In the weak coupling regime g 1, the dynamics of the system can be studied using perturbation theory. One starts by expanding the fields around one of the minima of the potential Since this is a sum of squares, the minima are given by X a = 0 and . This means that the minima are in one-to-one correspondence with integer partitions of N because we can form an N × N block diagonal matrix by adding many blocks with SU (2) irreducible representations. In the large N limit, tunnelling between different vacua is suppressed and it is possible to study the thermodynamics associated to each minimum [20]. 5 In this paper, we will focus on the trivial vacuum X a = X i = 0. The excitations above this vacuum are gapped and weakly coupled if g 1. For energies much greater than µ and much smaller than µN 2 the density of states grows exponentially with energy. This leads to a Hagedorn phase transition at T = µ 12 log 3 for g = 0. At weak coupling g, this becomes a first order phase transition whose critical temperature can be computed in perturbation theory [21,22] T c (g) = µ 12 log 3 1 + We call the high temperature phase the deconfined phase because the free energy scales as N 2 . For T < T c the system is in the confined phase where the free energy scales as N 0 .
The main goal of the present work is to determine the fate of this phase transition at strong coupling g 1. It is instructive to compare figure 1 with the phase diagram of N = 4 SYM on S 3 . In this comparison 1/µ plays the role of the radius of S 3 and g plays the role of the 4-dimensional 't Hooft coupling. 6 The 4-dimensional theory also has a first order phase transition that starts with a Hagedorn transition of the free theory [24]. At strong coupling, this transition corresponds to the Hawking-Page transition in the dual AdS 5 gravitational description [25,26]. We will argue that the PWMM has a very similar phase diagram.
In particular, we will find a Hawking-Page like phase transition in the dual gravitational description of the PWMM and predict the strong coupling limit of the critical temperature, 5 More precisely, this works for vacua that are associated with a reducible representation of SU (2) that contains many copies (of order N ) of a few irreducible representations of fixed dimension. If the dimensions of the irreducible representations scale with N and there are a fixed number of them (membrane states) then the free energy is of order 1 and there is no phase transition. In fact, the fluctuations around these vacua become free in the 't Hooft limit [21]. 6 In fact, the action (5) can be obtained from the action of SYM on S 3 by truncating the 4-dimensional fields to their zero modes (more precisely, projecting to SU (2) L invariant modes) [23].

T µ
Deconfined phase 0.076 0 g 0.106 Figure 1: The phase diagram of the PWMM. At high temperature, the system is in a deconfined phase where the free energy scales like N 2 . As we lower the temperature, the system undergoes a first order phase transition to a confined phase where the free energy scales as N 0 . The critical temperature T c depends on the dimensionless coupling g and it can be computed in perturbation theory for g 1. In this paper we determine T c at strong coupling from the study of the black hole dual to the deconfined phase of the PWMM.
It would be remarkable to confirm this prediction with Monte-Carlo simulations of the PWMM at strong coupling. We believe this to be accessible with the methods of [12,13,14,4,17].
The dual geometries to each vacuum of the PWMM were constructed in [27,28]. These SUSY vacuum geometries, including the one dual to the trivial vacuum, are surprisingly complicated (see appendix A) [29]. Nevertheless, they share an important feature in that they asymptote to the plane wave solution of M-theory Fortunately, we will not need the detailed form of these vacuum geometries. Our strategy will be to start from very high temperature (T µ) and gradually decrease it. This means that our starting solution is the uplifted 11D SUGRA solution (2) for which the 4form field strength vanishes. This geometry has the same SO(9) symmetry of the trivial vacuum X a = X i = 0. We will then continuously deform this solution by turning on a nonnormalizable mode of dC that corresponds to the relevant deformation that takes the BFSS to the BMN model. This deformation breaks the SO(9) symmetry of (2) to SO(6) × SO(3), making the field equations analytically intractable. In the next section, we explain how this is done in detail, including the numerical methods to solve the relevant Einstein equations.
In section 3, we determine the free energy of the black hole constructed in section 2 and the strong coupling limit of the critical temperature T c (g). We also calculate thermal expectation values of several operators in the high temperature deconfined phase. We conclude in section 4 with a discussion and comments about open questions.

Deformed Black Hole
Let us start by fixing our conventions for the bosonic piece of the 11-dimensional SUGRA action where η is the space-time volume form, R is the Ricci scalar and C is a 3-form gauge potential. Any stationary solution compatible with the SO(6) × SO(3) global symmetry and invariant under translations along the eleventh direction can be written as where ζ ∼ ζ +2π is the periodicity of the 11-dimensional circle and the functions A, B, F , Ω, T 1 , T 2 , T 3 , T 4 , M, L depend on the radial coordinate 0 ≤ y < 1 and on the angular coordinate 0 ≤ x ≤ 1. We shall see that y = 1 corresponds to the black hole horizon and y = 0 to the asymptotic region, which matches the plane wave geometry (9). The angular coordinate x was introduced to break the SO(9) symmetry of an eight sphere to SO(6) × SO(3). We can think of x = 0 as the S 5 equator and x = 1 as the S 2 pole, of the 8-dimensional surface dη = dζ = dy = 0. This form of the solution is tailored to the numerical methods we will use. In particular, all quantities are dimensionless and the domain of the unknown functions is the unit square. The physical solution can then be obtained by using the scalings of the 11D SUGRA action under the following transformations More concretely, the physical solution will be obtained from (11) by multiplying the metric by r 2 0 and the 3-form C by r 3 0 , and by changing the period of the non-contractible M-theory circle according to Both operations are symmetries of the equations of motion, but change the value of the on-shell action to where we defined the dimensionless action I to be the 11D SUGRA action (10) evaluated on the Ansatz (11) and stripped of the overall factor of 1/(16πG N ).
In the last equality of (14), we used the relations (3) between the gravitational parameters and the variables of the dual matrix quantum mechanics. When computing the action of a solution, care must be taken by adding boundary terms that renormalise the on-shell action. In what follows we shall assume that such counter-terms preserve both of the scaling operations described above. Similarly, the physical Bekenstein-Hawking entropy becomes where S is the dimensionless horizon area computed with the metric (11), explicitly given by where h is the determinant of the induced metric on the horizon, which has S 8 ×S 1 topology.
To see how this works in practice for a simple case, consider the exact solution given by A = B = Ω = T 1 = T 2 = T 3 = T 4 = 1 and F = M = L = 0. Changing coordinates, and multiplying the metric (11) by r 2 0 one recovers the 11-dimensional uplift of the nonextremal D0-brane solution (2). Notice that after the Wick rotation η → iη, the Euclidean time circle of (11) must have period 4π/7 in order to avoid a conical singularity. Using (17) this fixes the periodicity of the dimensionfull Euclidean time, which is consistent with the relations (3) between the temperature and the parameter r 0 . Moreover, using the dimensionless area of the horizon S = 2π Vol(S 8 ), we obtain 7 This exact solution describes the high temperature limit T /µ → ∞ of the PWMM. To lower the temperature, we need to appropriately turn on the 3-form potential C. This is implemented in the Ansatz (11) by requiring the function M = M (x, y) to have the following asymptotic behaviour To find out the physical meaning of the parameter µ, we compute the asymptotic behaviour of the physical field strength dC, determined after multiplying by r 3 0 and changing coordinates as in (17), Identifying rx √ 2 − x 2 as the radial coordinate on the 3-plane that contains the 2-sphere, and comparing with the M-theory plane wave solution (9), we conclude that In section 2.2 below, we will explain the precise boundary conditions that uniquely fix the solution. However, the intuition is clear: we require regularity at the axes of symmetry x = 0, x = 1 and y = 1. In particular, the Euclidean period of the η coordinate is always 4π/7 because we impose A = B at the horizon y = 1. 8 At infinity (y → 0), we impose that A, B, Ω, T 1 , T 2 , T 3 , T 4 → 1, that F, L → 0, and (19). In this way, we obtain a one parameter family of (dimensionless) solutions parametrized by µ. The physical entropy of the system, for example, is then computed using (15). Notice that this agrees precisely with the free energy scaling predicted in [6] from the assumption that the tree level and 1-loop contributions for the moduli effective action are of the same order in the strongly coupled regime. It is also clear that thermal expectation values that are non-zero at µ = 0 (i.e. in the non-extremal D0-brane) get multiplied by a function of µ, again in agreement with [6]. 7 Notice that this is compatible with the first term of (4) and the first law of thermodynamics ∂F ∂T = −S. 8 Recall that in Euclidean signature the horizon is the fixed point of time translation symmetry.

Harmonic Einstein equations -DeTurck method
Our approach to solving Einstein's equations, the so called DeTurck method, was first introduced in [30] and studied in great detail in [31]. Its generalization to finding stationary solutions of the form discussed in this manuscript was first detailed in [32].
We first note that in both the line element and gauge field ansätze (11), we have partially gauge fixed coordinate invariance and completely gauge fixed the gauge redundancy associated with C → C + dΛ, where Λ is a two-form 9 . However, the line element (11) still exhibits full diffeomorphism invariance for arbitrary reparametrizations of x and y, which we will fix using the DeTurck method.
The method can be best understood if we write the 11-dimensional Einstein's equations in the trace reversed form where is the Levi-Civita connection associated with a metric g and g is a reference metric. The reference metric g is chosen to have the same asymptotic and conformal structure as the metric we want to determine, i.e. the metric g. Here, our reference metric is just given by the line element (11) with A = B = T 1 = T 2 = T 3 = T 4 = Ω = 1 and F = 0. One can show that the resulting system of equations obtained via E ab = 0 and , Ω, M, L} and their first derivatives along x and y, and γ ab is a two-dimensional positive symmetric matrix (γ ab is the inverse of the metric tensor (11) restricted to a x, y plane with all the other coordinates fixed). This means that, under the appropriate boundary conditions, which we shall discuss below, E ab = 0 forms a system of Elliptic partial differential equations.
It is clear that any solution to E ab = 0 with ξ = 0 is a solution to E ab = 0, however, the converse is not necessarily true. Under some special circumstances, and for certain types of matter fields, one can show that solutions with ξ = 0, coined Ricci solitons, cannot exist [31]. However, the case under consideration is not under this special class. Fortunately, if the system of partial differential equations E ab = 0 is Elliptic, it can be solved as a boundary value problem for well-posed boundary conditions and the solutions should be locally unique.
This means that an Einstein solution cannot be arbitrarily close to a soliton solution and one should easily be able to distinguish the Einstein solutions of interest from solitons by monitoring ξ. This is the approach we are going to undertake here. However, we first need to show that our boundary conditions give rise to an Elliptic problem, and that they are consistent with ξ = 0.

Boundary conditions
Our solution naturally lives on a square grid, with x = 0 denoting the fixed points of the SO(3) symmetry, x = 1 the fixed points of the SO(6), y = 1 denoting the horizon location and y = 0 the conformal boundary. In an abuse of language, we shall refer to x = 0 as the S 5 equator, and x = 1 as the S 2 pole, of the 8-dimensional surface dη = dζ = dy = 0.
We will be interested in measuring certain quantities near the conformal boundary located at y = 0. The leading order term in the metric near the boundary is just that of the D0brane, since the matrix model massive deformation is irrelevant in the UV. Thus we will This guarantees that, asymptotically, the 8-dimensional surface dη = dζ = dy = 0 becomes a round S 8 and that the total D0-brane charge is fixed. The leading term for the functions M and L is fixed by requiring that the non-normalizable mode dual to the massive deformation of the matrix model is turned on. At leading order in y it suffices to consider linear perturbations of the 3-form potential C, which then fix the leading behaviour in y of the functions M and L. These perturbations are naturally expanded in a basis of harmonics of the asymptotic S 8 . For the particular Ansatz (11) for the 3-form potential, the decomposition in S 8 harmonics only includes the harmonic 2-forms with l odd and where 2 F 1 is a hypergeometric function. These harmonic 2-forms ω l satisfy where 8 is the Hodge dual on S 8 . Note that for odd l, ω l is invariant under the action of the SO(6) × SO(3) subgroup of the SO(9) isometry of the S 8 . Finally, analysing the perturbations of the 3-form potential C, one concludes that the required non-normalizable mode has non-vanishing functions M and L, and can be written in terms of the harmonic form ω l with l = 1. Asymptotically, this fixes their leading behaviour to be given by To better understand the boundary conditions at the conformal boundary we actually need to consider in more detail the asymptotic expansion of the fields. We now turn to this problem.

Asymptotic expansion at conformal boundary
In general, each function will have an asymptotic expansion in powers of y. For example, The equations of motion yield second order coupled differential equations in the variable x for all the coefficient functions like A n (x). These can be easily solved assuming that only smooth solutions on the S 8 of the boundary are allowed, that is to say, all the coefficient functions admit an expansion in harmonics on the S 8 of the boundary. These harmonics can be of scalar, vector or tensor type and must be invariant under the unbroken SO(3) × SO (6) symmetry.
Let us first consider the functions M and L that are associated to the 2-forms on S 8 introduced above. The expansion in powers of y can be seen to arise from the normalizable and non-normalizable modes that are excited, plus their back-reaction. At the linear level there are two independent field perturbations associated to M and L, which are called v 1 and v 2 in the perturbation analysis of [33]. We can drop the perturbation v 2 because we impose that its non-normalizable modes vanish, and its normalizable modes start at a power of y beyond what we consider in this paper. Thus, we have where we denote non-normalizable modes with a tilde and normalizable without. These non-normalizable modes behave near the boundary as f We set them all to zero but the mode l = 1. This is the content of the boundary condition (28), which sets α 1 = µ, and defines the type of relevant deformation we decided to turn on.
Of course we are not free to set the normalizable modes to zero. Their form can only be obtained once the solution is known everywhere, i.e. once regularity deep in the bulk and at the axis is imposed. These modes behave as f l (y) ∼ y 1+l , near the boundary. Notice that the normalizable modes of the perturbations v 2 , which we dropped in (30), have f l (y) ∼ y 8+l . Figure 2 summarizes these facts. In (30), we called back-reaction to all terms that are non-linear in the modes. At each order in the expansion at y = 0, these can also be expanded in harmonic 2-forms on S 8 . In this paper, we consider the first 8 terms in the expansion, where the coefficients β i , δ and γ appear in the expansion of scalar perturbations that we discuss in a moment. Note that the first line in the expansions (31) and (32) contains the terms linear in the modes, while the remaining terms arise from the back reaction of the fields.
The modes of the 2-form perturbation v 1 have a leading behaviour near the boundary of the form y 1−l for non-normalizable modes and y 1+l for normalizable modes. Figure 2 shows the power of y as a function of the spin l. It also includes the powers for the other perturbation v 2 . The dashed horizontal lines cover the region considered in this paper, up to order y 7 .
There are five S 8 scalars in our Ansatz. They are the functions A, B, T 4 , Ω and the trace of the S 8 metric fluctuations that we define by As in the previous case these fields can be expressed in terms of non-normalizable and normalizable modes, as well as the back reaction of all modes. In general, at the linear level there are three independent scalar perturbations, called s 1 , s 2 and s 3 in [33]. The perturbation s 2 vanishes for our Ansatz (this follows because our geometry is static after reducing to type IIA supergravity). We also drop the perturbation s 1 because our boundary condition impose the vanishing of its non-normalizable modes and its normalizable modes only start at order y 14 in the expansion near the boundary. Thus, for the purposes of this paper, we have and similarly for B, T 4 , Ω and Q. The functions S l are S 8 scalar harmonics, that are SO(3) × SO(6) invariant, and have the form with l even and P m n the associated Legendre polynomial of degree n and azimuthal number m. These functions are the usual eigenfunctions of the Laplace operator on S 8 with We also impose that the non-normalizable modes of the perturbation s 3 vanish. As shown in figure 3 this includes the three modes β 2 , β 4 , and β 6 that appear in the expansion near the boundary at order y 5 , y 3 , and y, respectively. The leading behaviour of the normalisable modes is also shown in the figure. The scalar perturbations also contain two zero-modes that we denote by γ and δ. These modes appear first at order y 7 and are constant on the S 8 . We call them zero-modes because they are not the zero frequency limit of any time dependent perturbation. Up to order y 7 the scalar perturbations have the form B(x, y) = 1 + δy 7 S 0 (x) , Next let us consider tensor perturbations. These arise from the fields T 1 , T 2 and T 3 , which we can write as where θ i denotes coordinates on S 8 and the symmetric tensor T ij is traceless with respect to the S 8 metric h ij , i.e.
The trace part of these tensor perturbations is given by the function Q already considered in the scalar perturbations above. The modes that appear in the symmetric traceless tensor T ij can be divided in their divergence and divergence-less parts. The divergence part is obtained by acting on the scalar harmonics with the differential operator where ∆ = ∇ i ∇ i is the S 8 Laplacian. These are, however, the same modes described above for scalar perturbations. In fact, their appearance in the tensor perturbations can be gauge away by imposing the gauge condition ∇ i T ij = 0. However, here they will be present in the tensor perturbations, since we do not have such freedom, because in the DeTurck method a given gauge choice is imposed on us. 10 The divergence-less part of these tensor perturbations comes again with non-normalizable and normalizable modes. We drop the non-normalizable modes and, for present purposes, we can neglect the normalizable modes because at linear level they appear first at order y 9 . Working up to order y 7 , the expansion of tensor perturbations reads where the (T l ) ij are S 8 harmonic tensors that satisfy with l ≥ 2 and with l even to guarantee invariance under the SO(3) × SO(6) subgroup.
Explicitly these harmonics are given by where Notice that, although an independent tensorial normalizable mode of spin l appears first at order y 7+l , these tensor perturbations already make their appearance at lower orders through the back-reaction.
Finally let us consider vector perturbations. For our Ansatz, F is the single SO(9) vector.
It turns out that there are no divergence-less vectors on S 8 that are SO(3)×SO(6) invariant.
Thus the expansion of this field will only contain derivatives of the scalar perturbations, which is indeed confirmed by the expansion 10 Thus, a gauge transformation is necessary to make the precise map between our expansion and that of [33].
All constants in the above expansions that remain to be determined correspond to expectation values of dual operators in the matrix model. Up to order y 7 in the above expansions, these are the constants α 1 , α 3 , α 5 and β 2 , β 4 , β 6 and γ, δ. More normalizable modes show up at higher order, but we decided to only present results for these.
For a more accurate numerical extraction of the remaining normalizable modes, we do a final change of variables that will ease the numerical procedure, namely we define Our numerical procedure aims to solve for all ten Q i (x, y). We impose the following Neumann and Dirichlet boundary conditions at y = 0 These boundary conditions guarantee that all non-normalizable modes (except µ) are set to zero. In particular, the mode β 2 is the hardest to exclude because it only appears at order y 5 in the asymptotic expansion. For example, a non-zero β 2 would give rise to Therefore, the boundary conditions (53) force β 2 = 0.

Symmetry axes
The boundary conditions at the equator x = 0 are just those obtained via smoothness of the solutions. This implies that all Q i should be even functions of x, except Q 3 , which should be odd under x → −x. Moreover, we must have Q 4 = Q 5 at x = 0 to avoid a conical deficit.
In practice, we just impose Similarly, at the x = 1 pole, we require that F is odd and A, B, T 4 , Ω, M, L, are even under reflection around x = 1. Moreover, we avoid conical deficits by imposing These conditions imply ξ x = 0 at x = 0 and x = 1.
At the horizon, which in the Euclidean setting is also a symmetry axis, regularity is easier to impose after changing to a new radial coordinate via 1 − y = (1 − y) 2 . In the y coordinate, the conditions for regularity are that F is odd and all other functions are even under reflection around y = 1. Moreover, we impose A = B at y = 1 to avoid conical deficits with the periodicity ∆η = 4π/7. In practice, we use the boundary conditions These boundary conditions imply ξ y = 0.
It is a relatively easy exercise to show that the boundary conditions detailed above, together with the Einstein-DeTurck equations, form a well posed Elliptic problem [34,35]. Furthermore, at the fictitious boundaries x = 0, x = 1, and y = 1, the boundary conditions induced on ξ are the relevant ones to admit ξ = 0, i.e. Einstein solutions, everywhere in the bulk [31]. We are thus ready to present our results and to detail the numerical method we used to solve the Einstein-DeTurck equations.

Smarr formulae
Having Smarr formulae is particularly important in situations in which the solution is presented numerically. Let us then recall how we can construct such formulae in the present case. For every Killing vector v of the solution (11) we can define an antisymmetric conserved Conservation of this tensor follows from the equations of motion (22), ∇ a G abcd = 0, and from the identities In the language of differential forms, this means that we have a closed 9-form where Integrating d( K v ) over a 10-dimensional surface Σ 12 of constant time with y 1 < y < y 2 , we conclude that where we used the fact that the boundary of Σ 12 has two disjoint components Γ(y 2 ) and Γ(y 1 ) (with opposite orientations). This shows that the integral Choosing v = ∂ ∂η to be the generator of time translations, we obtain To compute I v (1) we used the fact that the horizon is a Killing horizon of ∂ ∂η with surface gravity equal to 7 2 . To compute I v (0) we used the asymptotic expansion of the fields.
Choosing v = ∂ ∂ζ to be the generator of translations along the M-theory circle, we obtain where the integration measure, given in (16), is defined by the horizon metric. This integral measures the momentum along the M-theory circle (or D0-brane charge in the type IIA picture) which is constant as we vary µ. We can think of the first term in (66) as the momentum carried by the black string, and the second term, which is also positive, is the momentum carried by the matter fields outside the horizon. In Fig. 4 we plot the momentum carried by the black string. As µ increases the momentum carried by the fields outside the horizon increases.
It is also useful to integrate the d( K v ) over the 10-dimensional surface of constant ζ. By a similar argument as the one above, we conclude that the following integral is independent of y where Γ(y) is the 9-dimensional surface of constant y and ζ. Choosing v = ∂ ∂ζ we obtain I v (1) = 0 from the behaviour of the solution at the horizon. Thus, using the behaviour as  (66)). As µ increases, the black string and the fields outside the horizon carry less and more momentum, respectively, keeping the total momentum of the geometry fixed.
y → 0, we deduce the following identity relating the parameters of the asymptotic expansion of the fields.
For v = ∂ ∂η we also obtain I v (1) = 0. However, I v (0) depends on higher orders of y in the asymptotic expansion of the fields than those considered above.

Numerical solution
We used a standard pseudospectral collocation in x and y, and solved the resulting system of non-linear algebraic equations with a damped Newton-Raphson method. The dependence in x and y of all the functions was represented using tensor products of two Chebyshev collocation grids, each of which living on the unit interval (0, 1). Our integration domain is thus a square (x, y) ∈ (0, 1) × (0, 1).
In expanding the functions Q i around the relevant boundaries, we have found no sign of non-smoothness. This means that a priori we expect the convergence of our method to be exponential in the number of grid points N and that no patching procedure is required.
The only delicate numerical problem associated with these equations is that we need to accurately extract third and fourth derivatives off of the conformal boundary, in order to read the several constants corresponding to normalizable modes. For this reason, we decided to work with octuple precision and no less than 51 grid points on each integration domain. In addition, due to the very bad condition numbers of the matrices we have to invert, we found useful to use up to twelve patches close to the boundary (depending on the values of µ and how steep our functions behave). These are conforming patches, which are patches that only coincide along a line, and have no overlapping regions. Since we are interested in accurately extracting asymptotic quantities, our patches coincide with lines of constant y and cluster close to y = 0.
In order to monitor the convergence of our numerical method, we monitored χ = ξ a ∞ as a function of the number of grid points N , as well as where S N denotes the entropy computed with N grid points in both directions. Both plots are displayed in Fig. 5, where a linear-logarithmic scale is used and we set µ = 1. The results are consistent with exponential convergence, as dictated by pseudospectral collocation methods.  A perhaps more striking test of our numerics comes from the identity (68). We have checked that this relation is obeyed by our numerical data, never exhibiting a violation above 10 −6 %. Similarly, we checked that the Smarr formulas (65) and (66) are verified by our numerical solutions with an accuracy of 10 −6 %. The Smarr formulae provide a very non-trivial validation of our numerical results because they relate quantities measured at the horizon (y = 1) to quantities measured at infinity (y = 0). This gives us full confidence that our numerical procedure is accurate enough for the physics we want to extract.
In Fig. 6 we plot a typical run of our numerical method. It shows the behavior of Q 1 , Q 9 and Q 10 as a function of x and y. Note that these are all gauge invariant. From these plots we can easily see why we needed octuple precision, namely there is a large hierarchy between the functions. For instance, Q 9 evaluated on the horizon appears to be larger than all the remaining functions. This problem becomes worse as we increase µ. Figure 6: From left to right: three-dimensional plots of Q 1 , Q 9 and Q 10 as a function of x and y for fixed µ = 1.
We now turn to more physical quantities. In particular, we would like to see how the horizon shape is changing as we change µ. It is clear that the geometry will slowly move from having a round S 8 with SO(9) symmetry to a deformed S 8 with a manifest SO(3) × SO(6).
To explicitly quantify how deformed the horizon is from full spherical symmetry, we measure the radius of the S 2 at the pole and the radius of the S 5 at the equator. If the ratio between these quantities is very small, the horizon is highly distorted from spherical symmetry. We plot this quantity in Fig. 7a.
The fact that this ratio reaches such small values might be worrying and suggestive of a Gregory-Laflamme type instability along the S 5 directions. In order to settle this, one would need to perturb this solution, and check its dynamical stability. We are currently undertaking this study, but have no results to report. Finally, we can also plot the normalized area of the horizon as a function of µ, which we will need to reconstruct the free energy. This is done in Fig. 7b, where we see the horizon area decreasing with increasing µ.
We finalize this section by presenting, in Fig. 8, the several extracted expectation values as a function of µ. We obtain these expectation values by computing the first few y-derivatives at y = 0 of the functions Q i (x, y) and fitting them to the asymptotic expansions discussed in section 2.2.1. A detailing of the predictions determined by perturbations around the µ = 0 background (dashed red lines) is given in Appendix B. Any Monte Carlo simulation of the PWMM (5) should hope to reproduce these results.

Thermodynamics
Our numerical solution, expressed in terms of the functions in the Ansatz (11), depends on a single dimensionless parameter µ that determines the asymptotic behaviour of the 3form potential C through the boundary condition (19). Thus, the corresponding on-shell dimensionless action I and entropy S, respectively defined in (14) and (15), are functions of this single parameter. The boundary conditions imposed at the horizon fixed the periodicity of the Euclidean time circle to 4π/7, independently of µ.
Next, to obtain physical solutions from the above single-parameter family of solutions, we scaled the metric by r 2 0 and the 3-form C by r 3 0 , and changed the period of the Mtheory circle according to (13). The new family of solutions, parametrized by µ and r 0 , has the same leading asymptotics of the non-extremal D0-brane solution (2) with an additional 3-form potential C with asymptotic behaviour (20). It is then convenient to parametrize this new family of solutions by the temperature T and the mass deformation µ, which are related to the original single parameter by µ = 7 12π µ T as derived in (21). Moreover, the  on-shell action and entropy of the two-parameter and single-parameter families of solutions are simply related by (14) and (15), which we can rewrite in the form for a known (dimensionfull) constant c 0 . In the particular case of zero mass deformation µ = 0 we recover the scaling with temperature as predicted directly from the matrix quantum mechanics in [6,7]. It is then clear that both the free energy and entropy are restricted to satisfy F (T, µ) F (T, 0) = I( µ) where, by definition, f (0) = s(0) = 1.
The behaviour of the free energy and entropy (71), together with the scaling of the free energy as T 14 5 at zero mass deformation µ, can be used in the first law to relate the functions f ( µ) and s( µ). This leads to the following equation which can easily be integrated where we wrote explicitly the integration constant C. Notice that the boundary condition f (0) = 1 does not determine the constant C. However, assuming that both s( µ) and f ( µ) are analytic around µ = 0, and therefore have a regular Taylor series expansion, removes all ambiguity, Since from computing the horizon area we know the function s( µ) numerically, we can do a polynomial fit to determine the first coefficients s n , and then use it to plot f ( µ) in Fig.  9. The most important feature of this plot is that f vanishes for µ = µ c ≈ 1.7532672. This means that, for µ > µ c , the free energy of the deconfined phase of the PWMM is positive and of order N 2 . Therefore, the confined phase that has a free energy of order N 0 will be  Figure 9: The free energy ratio f ( µ) obtained numerically using (75).
Let us now consider thermodynamical stability. The specific heat of the system is given by From (70) and (71) we may also express the specific heat in terms of the function s( µ) as Since in the range the black hole geometry is thermodynamically favoured, s( µ) is a decreasing function, as shown in Fig. 7b, we conclude that the specific heat is always positive and therefore our solution is thermodynamically stable in this range.

Discussion
Our main result is the construction of the black hole geometry dual to the deconfined phase of the PWMM. This allowed us to determine the value of the critical temperature at strong coupling as depicted in the phase diagram 1. In addition, we determine the thermal expectation values of several observables in the deconfined phase (see Fig. 8).
At this point we would like to discuss an important caveat that we disregarded in the main text. There should be many black hole geometries with different horizon topologies and the same asymptotics as the solution we constructed. One can think of these as the finite temperature and backreacted versions of the many ways to distribute spherical probe M5 and M2 branes in equilibrium in the M-theory plane wave [18,36]. These solutions are in one-to-one correspondence with the many vacua of the PWMM [37]. Our expectation is that the solution with lower free energy in the high temperature limit is the one we found because it has the simplest horizon topology. However, as we decrease the temperature it is possible that other black hole solutions start to dominate the thermal ensemble. 12 Therefore, what we really determined was an upper bound for the critical temperature for the deconfinement transition. Notice that it is sufficient to find one black hole solution with negative free energy at a given temperature to conclude that the system must be in the deconfined phase at that temperature. Even if this solution is dynamically unstable it must decay to another solution with lower free energy, thus the system remains in the deconfined phase.
We hope our results motivate others to start a systematic exploration of the phase diagram of the PWMM by direct simulation of the matrix quantum mechanics, e.g. using the Monte-Carlo methods of [12,13,14,4,17]. Our work provides concrete predictions for the behaviour of several thermodynamic quantities at strong coupling in the deconfined phase.
We also provide predictions for thermal expectation values of several operators. However, the precise map between the gravitational parameters shown Fig. 8 and operators of the PWMM is still missing. This map is known [33] in the limit µ → 0 but its extension to finite µ requires the development of holographic renormalization with plane-wave asymptotics.
In fact, there has been a preliminary Monte-Carlo simulation of the PWMM [38]. In this work, the authors simulate the PWMM at fixed temperature (T /µ = 1/3 in our conventions) and as they vary the coupling, they observe a first order phase transition for 0.03 g 0.045.
This result is not in direct contradiction with our results but it implies a non-monotonic behaviour of the critical temperature as a function of the coupling g, complicating the phase diagram 1. It would be great if this result could be confirmed by a more systematic Monte-Carlo simulation of the PWMM.
It would also be very interesting to perform a Multicanonical Monte-Carlo simulation [39,40,41] of the PWMM that could measure the density of states of the system. This would provide a window into the thermodynamics of the system in the microcanonical ensemble, which is expected to have a richer structure, including a Hagedorn phase. 13 Our black hole solution was constructed starting from the limit µ/T = 0. It would be interesting to understand our solution in the opposite limit µ/T → ∞. It is hard to address this question using our numerical methods because the black hole becomes very deformed and requires a much finer discretization grid. In any case, Fig. 7a suggests that when µ T the black hole looks like a pancake (more precisely, a large 6D ball with a small thickness in the transverse 3 directions, times the M-theory circle). It should be possible to study this limit analytically using the blackfold approach of [42]. The large deformation of the horizon also suggests that the system might be unstable to a topology change to a ring-like horizon with S 5 × S 3 × S 1 topology. It should also be possible to study the low temperature regime of such black holes using the blackfold approach. We leave these ideas for the future.

A Vacuum geometries
The supergravity solutions dual to the vacua of the BMN model were constructed in [27,28].
They are given by with asymptotic behaviour V ≈ ρ 2 z − 2 3 z 3 and boundary condition V (z = 0) = 0. A given vacuum of the BMN matrix quantum mechanics corresponds to a certain distribution of charged disks sourcing the potential V as explained in [28].
For our purposes, it is sufficient to consider the asymptotic behaviour of the potential at large z ∼ ρ, where the P l is the Legendre polynomial and a l characterize each specific vacuum because they are multipoles of the charge distribution that sources the potential V . Inserting this expansion of the potential in the solution (79) we obtain an asymptotic expansion that can be compared with the asymptotic expansion of our Ansatz (11) discussed in section 2.2.1.
More precisely, we perform the following change of coordinates in our Ansatz (11) where the dots denote terms suppressed by powers of ρ 2 + z 2 that are determined so that our Ansatz (11) has the same type of asymptotic expansion as the vacuum solutions (79).
This comparison leads to the following relations µ = 2 √ a 1 , β 2 = a 3 a 1 , β 4 = a 5 a 1 , β 6 = a 7 a 1 , In other words, the vacuum geometries of [27,28] have an asymptotic expansion of the form of section 2.2.1 with all the parameters given in terms the multipoles a l and an arbitrary constant φ 0 that represents the freedom to shift the potential associated to the D0-brane charge (from the 10-dimensional point of view). In addition, the parameter β 2 , that appears for example in (54), vanishes in all vacuum solutions. This confirms our interpretation of β 2 as a state dependent response and β 2 as a source that deforms the theory.

B Perturbations Around the Background Solution
We begin by expanding the functions using the spherical harmonics outlined in section 2, where Q and T ij were introduced in (42). We plug this Ansatz into the equations of motion and expand them up to O( µ 2 ). In particular, the equation of motion d dC = 0 gives rise to linear ODEs for the functions q 9,l (y) and q 10,l (y). Moreover, the boundary conditions discussed in section 2.2 imply that the only non-zero modes are q 9,1 (y) and q 10,1 (y). We find these functions using a single variable version of the Chebyshev method described in the main text to reduce the two linear ODEs to a set of linear algebraic equations, which can easily be solved using Newton's method (we solve the equations in the y coordinates because this simplifies the boundary conditions on the horizon).
The harmonic Einstein equations E ab = 0, expanded up to O( µ 2 ), can also be decomposed into spherical harmonics. The equations of scalar type are E τ τ , E yy , E zz , and E τ z . For example E τ τ (x, y) = l S l (x)f 1,l (y).
The vector equation E yi (where i runs over the S 8 coordinates) can be decomposed as follows Finally, we can use the tracelessness of the tensor harmonics as well as the divergence-less nature of T ij to write the components of the Harmonic Einstein equations corresponding to the S 8 Written in this way, the equations E ab = 0 can easily be projected onto a basis of the scalar harmonics, resulting in a system of seven ODEs, These ODEs are linear in the functions q 1,l (y), q 2,l (y), q 7,l (y), q 8,l (y), q 3,l (y), q l (y), q l (y) and quadratic in the functions q 9,l (y) and q 10,l (y) that can be previously determined from the gauge field equation of motion d dC = 0. 14 The terms quadratic in the known functions q 9,l (y) and q 10,l (y) can be thought of as sources in the linear equations for the other 7 functions. This gives rise to 7 linear non-homogeneous ODEs which can easily be solved using spectral methods. To this order in µ, the sources are only non-zero for l = 0 and l = 2.
The normalizable modes can be extracted from the solutions by comparing their behavior with the asymptotic expansions of the fields. In particular,  14 Note that the equation f l (y) = 0 is automatically satisfied setting q l (y) = 0.
These values satisfy the Komar identities (65, 68) and provide a non-trivial check of the numerics as shown in Figure 8.