The large dimension limit of a small black hole instability in anti-de Sitter space

We study the dynamics of a black hole in an asymptotically AdSd × Sd space-time in the limit of a large number of dimensions, d → ∞. Such a black hole is known to become dynamically unstable below a critical radius. We derive the dispersion relation for the quasinormal mode that governs this instability in an expansion in 1/d. We also provide a full nonlinear analysis of the instability at leading order in 1/d. We find solutions that resemble the lumpy black spots and black belts previously constructed numerically for small d, breaking the SO(d + 1) rotational symmetry of the sphere down to SO(d). We are also able to follow the time evolution of the instability. Due possibly to limitations in our analysis, our time dependent simulations do not settle down to stationary solutions. This work has relevance for strongly interacting gauge theories; through the AdS/CFT correspondence, the special case d = 5 corresponds to maximally supersymmetric Yang-Mills theory on a spatial S3 in the microcanonical ensemble and in a strong coupling and large number of colors limit.


Introduction
The most well-understood example of gauge/gravity duality is the correspondence between type IIB string theory on AdS 5 × S 5 and N = 4 super Yang-Mills with SU(N ) gauge group [1][2][3]. In the large N and large 't Hooft coupling λ limit, the bulk theory reduces to classical gravity. Using general relativity as a tool, gauge/gravity duality has led to better understanding of strongly interacting field theories in general and N = 4 super Yang-Mills in particular. An early success of the correspondence was the understanding of a deconfinement type phase transition in N = 4 Yang-Mills on a sphere as a Hawking-Page phase transition between thermal AdS in global coordinates and an asymptotically AdS geometry containing a black hole [4]. While this thermodynamic phase transition is a well established story in the canonical ensemble, in the micro-canonical ensemble the situation is much less clear. In the dual geometry, there exist small black holes which are both thermodynamically and dynamically unstable. In analogy with the Gregory-Laflamme instability for black strings, refs. [5,6] conjectured the endpoint of the instability to be a black hole with S 8 topology, localized on the S 5 , with a corresponding spontaneous breaking JHEP02(2018)167 of the SO(6) R-symmetry to SO (5). Later numerical analyses [7][8][9] pinned down the onset of a Gregory-Laflamme like instability as a function of horizon radius. The endpoint of the classical instability has as yet been studied only numerically [9,10]. The goal of this work is to better understand both the onset and ultimate endpoint of the dynamical instability by taking a novel limit -the number of dimensions to be large.
While increasing the number of dimensions typically makes Einstein's equations more difficult to solve, simplifications emerge if the number d is taken large enough while at the same time imposing that the solution maintain a high degree of symmetry. In a large dimension limit, the gravitational effects of a black hole horizon die off very quickly with distance. The radial direction can be "integrated out" of Einstein's equations, leaving a simpler set of hydrodynamic like equations to solve that govern the physics of a membranelike horizon. Such an approximation has already proven useful for studying classical black holes [11][12][13][14][15]. The approximation has been successfully applied to black holes with planar topology in the context of applied holography [16][17][18], fluid/gravity duality [19][20][21][22][23], and the membrane paradigm [24][25][26][27]. Relevant for us, the formalism can also be extended to black holes with spherical topology [28][29][30][31][32].
In this paper, we apply the large dimension limit to a study of black holes in an asymptotically AdS d × S d space-time. As we want the d = 5 case to be dual to N = 4 Yang-Mills, we include an anti-symmetric d-form field strength F d in the equations of motion. We start with a linearized analysis of the black hole instability, governed by coupled fluctuations of the metric and d-form. Through a numerical zero mode analysis in d = 5, ref. [7] calculated the onset of the instability as a function of black hole radius. Several years later, ref. [8] provided a more detailed numerical analysis of the associated quasinormal modes, providing plots of the dependence of the mode frequency on horizon radius and revealing that the frequency is purely imaginary. Our large d limit allows us to provide analytic power series expressions for these quantities. We give the threshhold radius to O(d −3 ) and the quasinormal mode dispersion relation at leading order in a large dimension limit. Generalizing the numerical analyses of refs. [7,8] to d > 5, we find our series expressions give good approximations of the numerical results.
Next, we investigate the full non-linear time evolution of these small black holes at leading order in our large d expansion. Assuming that the endpoint of the black hole instability preserves an SO(d) symmetry, we take an ansatz for the metric and d-form that depends only on time t and radial coordinate r of AdS d and a polar angle θ on the S d . In fact, we consider two slightly different large d expansions. In the first, we keep θ as a dynamical variable while in the second, we rescale the polar angle u ∼ √ n(θ − π/2) to zoom in on the dynamics at the equator of the S d . Ultimately, we derive two sets of hydrodynamic-like partial differential equations, one that depends on (t, θ) and another on (t, u).
We find a large class of exact solutions of the (t, θ)-system. At a discrete set of horizon radii, these solutions can be static. The discrete radii correspond to the onset of the initial Gregory-Laflamme like instability and further instabilities at yet smaller radii. In general, however, the solutions we find are strongly time dependent, with a e e t type functional form. The (t, u)-system is more complicated. We can find simple formulae for static solutions JHEP02(2018)167 which exist not just at discrete values of the horizon radius but over a range of values. We can also analyze the time dependence of the system numerically.
The solutions we find share qualitative similarities with numerical studies [9,10] in d = 5, which use the de Turck method to search for stationary solutions of Einstein's equation. These studies have shown evidence for several different types of solutions. There are fully localized spot and belt type solutions on the S 5 which have horizon topology S 8 and S 4 × S 4 respectively. There are also "lumpy" solutions which share the symmetries of the localized solutions but which maintain the original S 3 × S 5 horizon topology. In our large d limit, our solutions are of the "lumpy" variety with S d−2 × S d topology. While strictly speaking, the topology of all of our solutions remains S d−2 ×S d , the solutions can be in some sense "exponentially close" to fully localized solutions with S 2d−2 and S d−1 × S d−1 topology.
In our analysis, energy and entropy are equal at leading order in our large d expansion. Thus, we are not able to distinguish the lumpy and homogeneous solutions based on thermodynamic considerations alone. In the future, we hope to distinguish entropy and energy by extending our analysis to subleading orders. In the meantime, we try to gain some insight into which configurations are preferred by examining the time evolution of perturbations of the homogeneous solutions. We find that at long times, there is run-away behavior, with the horizon radius either growing or shrinking uncontrollably at the poles of the S d . Conceivably subleading terms in the 1/d expansion could stabilize the runaway behavior. Alternately, a logical possibility remains that the evolution does not settle down to an equilibrium configuration [33].
The paper is organized as follows. In section 2, we study linear fluctuations of the small black holes in a large d limit. In section 3, we provide our nonlinear analysis. Section 4 contains a brief discussion and conclusion. Appendices contain various auxiliary technical results.
2 Quasinormal modes of black holes in a large d limit

Linearized equations of motion
We begin with the formulation of the small black hole solutions in an arbitrary number of dimensions. In the special case of ten dimensions, this spacetime can be formed as a solution of the type IIB supergravity equations of motion with the self dual 5-form switched on. In arbitrary dimension, we write the action as The equations of motion take the form

JHEP02(2018)167
In compliance with supergravity, we require F to be self-dual, which means that we are limited to d ≡ p = q, and the norm |F | vanishes. The self-dual field vanishes for d even, so we assume d is odd for the purposes of deriving the equations of motion. (In the final effective equations, d will be a numerical parameter that can take any real value, however.) The canonical stationary solution for this setup is the anti-de Sitter-Schwarzschild metric smeared over a sphere The line element dΩ 2 d is constructed from the usual metric on a d-dimensional sphere of radius one.
We would like to consider linearized fluctuations over this background. We choose a transverse-traceless gauge, in which the fluctuation of the d-form field can be set to zero consistently. For the metric fluctuation g M N = g where ω is the quasinormal mode frequency. This fluctuation is subject to the linearized Einstein equations where ∇ 2 x and ∇ 2 y are the d'Alembertians on AdS d and S d respectively and µ, ν, . . . index the AdS d coordinates. For the zero mode, ω = 0, we can drop χ 3 , and it is straightforward to obtain the equations for the χ i 's. The gauge conditions fix χ 1 and χ 4 , leaving a second order differential equation for χ 2 coming from the rr component of Einstein's equations There are higher order differential equations from the other components of Einstein's equations, but this equation is sufficient to satisfy all the others. We may call it the master equation [36]. We can similarly proceed to derive the master equation for general ω. We

JHEP02(2018)167
use the gauge conditions to fix χ 3 and χ 4 , and then use the rr Einstein equation to fix χ 1 . The rest of Einstein's equations can be manipulated to yield a single second-order differential equation for χ 2 , χ 2 + P χ 2 + Qχ 2 = 0, (2.9) where P and Q are lengthy expressions in terms of d, r, ω 2 , , f , and derivatives of f . We defer them to appendix A. Since (2.8) or more generally (2.9) is the only equation to solve, we now drop the subscript on χ 2 . In the analysis of quasinormal modes, one solves the χ equation (2.9) with physical boundary conditions and obtains the dispersion relation ω = ω(r + , ). This relation contains information about how quickly a perturbation to the black hole (or equivalently the dual nonzero energy state of the field theory) dies off. In particular, the zero mode, ω = 0, will give the threshold radius r * of instability. Since these quantities can be readily computed numerically, we shall obtain some insight into the quality of the large d approximation for our setup. In addition, the analytic expression for ω will serve as a consistency check for the hydrodynamic equations that we derive later.

Threshold radius and dispersion relation
In an asymptotically AdS d spacetime with a S d−2 × R conformal boundary (i.e. global coordinates), it is well-known that there are two branches of black hole solutions, small and large. Using the usual trick of analytic continuation into the Euclidean metric, the Hawking temperature is Fixing T , the resulting quadratic equation can be solved for two possible values of the horizon radius r + (see figure 1). We pass three points of interest on this curve as we shrink the black hole radius down from infinity. The first point we pass, at (r HP , T HP ) = L, d−2 2πL , is the Hawking-Page phase transition. This first order phase transition is to a thermal AdS geometry [39], and is dual, through the AdS/CFT correspondence, to a confinement/deconfinement phase transition in the dual Yang-Mills theory [4]. The second point of interest, is where the small and large black hole branches meet and the heat capacity C v ∼ dr + /dT switches from a positive to a negative quantity. In the limit d → ∞, this critical radius is simply r HC = L, which is also where the Hawking-Page phase transition occurs. What we reconfirm in the analysis is that the dynamical, or Gregory-Laflamme like, instability at r + = r * sets in for slightly smaller black holes. In what follows, we choose L = 1 which sets the distance scale.
To obtain the threshold radius for the dynamical instability, we need to examine the case ω = 0, for which the χ equation takes a considerably simpler form (2.8). Due to the term 1/r d−3 in f , it is appealing to define n ≡ d − 3. The large n limit needs to be taken carefully. The gravitational effects of massive objects fall off very rapidly with distance at Approaching from the large black hole branch, the circle indicates the Hawking-Page phase transition to thermal AdS 5 in the canonical ensemble, the square is the point at which the heat capacity becomes negative, and the triangle is the Gregory-Laflamme type instability [7]. As d gets larger, the points approach each other.
large n. To keep the black holes from completely decoupling, we follow [12] and define a rescaled radial coordinate R = (r/r + ) n . The eigenvalues of the spherical harmonics diverge with n, and need to be rescaled. Thus we defineˆ = ( + n + 2)/n. We impose the usual boundary conditions on χ, i.e. χ is regular near the horizon R → 1, and decays sufficiently fast near the AdS boundary R → ∞. Examining the equation (2.8), there are two possible behaviors near the horizon, χ ∼ 1 and χ ∼ 1/(R − 1). A double power series in (R − 1) and then in 1/n gives for the first behavior At the zeroth order in (R − 1), it is simply constant with no higher order terms in 1/n. This makes it easy to impose the horizon boundary condition, which is to eliminate all the terms of order ∼ O(R − 1) 0 at all higher orders in 1/n. There are also two modes near the AdS boundary, which scale as χ ∼ 1 and χ ∼ 1/R at the leading order in 1/n. The second mode can be expanded as It is not entirely straightforward to match this expression to a series where we have first taken the large n limit and then taken R → ∞. Further expanding each R −2j/n term above in a large n limit will produce log(R) terms that need to be resummed. In what follows, we are able to apply boundary conditions simply by forbidding the χ ∼ 1 behavior and allowing for log(R) j /R terms, where j is a small integer, in the large R limit.

JHEP02(2018)167
Given the boundary conditions, the most straightsforward way to solve the χ equation is to expand χ and r + as (2.14) The differential equations at each order can be solved through DSolve in Mathematica. 1 We can fix r (j) + by imposing the boundary conditions. However, the equations get rapidly complicated, and it is difficult to go beyond first order for r (j) + . Instead, it was found useful to use the Green's method as laid out in [12], which allows us to move up one more order. At each order, we can write the χ equation in the self-adjoint form as and L (j) is the differential operator in the master equation expanded in 1/n so that j n −j L (j) χ = 0. The solutions at each order are given by the Green's method as where u = 1 R and v = 1 1−R . We then have the consistency condition on the sources S (j) , which is derived from the Wronskian of χ (j) [12] lim This condition gives the same result as imposing the boundary condition near R → ∞ on χ (j) , saving the effort to compute χ (j) itself at the corresponding order. We were able to compute S (j) up to the fourth order, from which we find This gives the threshold radius of the Gregory-Laflamme like instability to third order at large n.

JHEP02(2018)167
Dispersion relation for ω and . We can similarly proceed for the ω = 0 case to obtain the dispersion relation, ω = ω( , r + ). The main difference is the appearance of waves near the horizon. Near the horizon at large n, the χ equation (2.9) admits outgoing and infalling modes The infalling mode is conceived as physical, whose horizon expansion takes the following form 2 (2.21) As in the case of ω = 0, at the AdS boundary, we only required that χ is suppressed by 1/R.
Even with the Green's method, however, it is difficult to go beyond the second order for S (j) , because the integral (2.18) becomes hard to solve analytically. The leading order result is This gives the dispersion relation in the large n limit, and yields the same leading order result for the critical radius (2.19). At leading order in n, it makes no difference in the above expression if we writeˆ or . In fact, the higher order corrections appear to be numerically smaller if we use . Furthermore, it is that naturally shows up in the nonlinear analysis we perform later.

Comparison to numerical results
For the purpose of numerical computation, it is useful to transform the original χ equation (2.9) to simplify its asymptotic behavior. We introduce [8] X(y) ≡ (1 − y) 1+iω/4πT y d+ +1 χ(y), y = r + r . (2.23) In this variable and field, the quasinormal boundary conditions amount to The difference between the analytic formula for the threshold radius r * (2.19) and numerical calculation is near the order of expected error ∼ 1/n 3 , shown in figure 2. Further, the two consecutive orders in the formula provide lower and upper bounds, at least for the range of we computed numerically. Both the analytic and numerical results indicate that the threshold of Gregory-Laflamme instability approaches asymptotically to r * ∼ 1/ √ 1 + for large n, reassuring us that the dynamical instability only develops below the thermodynamic threshold r * ∼ 1.
2 There is a subtle relationship between the horizon boundary conditions when ω = 0 and when ω = 0.
Naively, the limit ω → 0 should lead to 1/(R − 1) and log(R − 1)/(R − 1) leading horizon behaviors. In fact, as can be seen from the horizon expansion, the log(R − 1)/(R − 1) term drops out, and the roots of the indicial equation change from −1 ± iω/4πT to 0 and −1. We would like to thank L. Yaffe for discussion on this point.   including the terms up to n −1 and n −2 order, serving as the lower and upper bounds respectively. The dots are the numerical results; the solid line is the n → ∞ limit. The large n approximation of the dispersion relation ω = ω(r + , ) (2.22) also shows an expected deviation ∼ 1/n from the numerical calculation. In the numerical calculation, we assumed that ω is purely imaginary, which is manifest in the large n limit (2.22). As seen in figure 3, although the deviation is not small for the dimension of interest, d = n + 3 = 5, the approximation still captures the qualitative features.
3 Evolution of black holes in the large d limit

Effective equations
In the large d limit, planar black holes are found to behave like a viscous fluid described by effective hydrodynamic-like equations [24]. A similar description is possible for our spherical black hole. To find such a description, we consider a metric and d-form ansatz of JHEP02(2018)167 the following restricted form where g's and f 's are functions of v, θ, and r only. The coordinate v is time-like, and θ is one of the coordinates on S d . A key difference between this nonlinear and the previous linear analysis is that the S d here is allowed to fluctuate only with respect to the polar angle θ. (Our metric is similar to the Kerr black hole, but we use the polar angle instead of the azimuthal angle.) We will restrict F d to be self-dual in accordance with supergravity. Similar to the previous section, we use the scaled radial coordinate R = (r/r c ) n , where r c is some reference scale and n = d − 3. We define the horizon radius R + to satisfy g vv (R + ) = 0. In general, R + will depend on v and θ. However, by the relation r + = r c (R + ) 1/n , the original horizon r + converges to r c in the large n limit. Our large n limit forces us to probe the geometry close to the reference scale r c , which is independent of v and θ.
In solving the Einstein equations (2.3), it is convenient to exploit the fiber bundle structure of our metric ansatz. The Ricci tensor can be expressed in terms of the curvatures of the S d−1 and S d−2 plus additional contributions (see the equations (B.21) in the appendix B). As to the boundary conditions, we require that the spacetime becomes AdS d × S d near the boundary R → ∞ and regular near the horizon R → R + . 3 It is then tedious but straightforward to expand the g's and f 's in 1/n and impose the boundary conditions. Some of the leading terms in the solution (see appendix C) are We chose the field e in such a way that R + = e + O(n −1 ). The fields e and j play the roles of fluid density and fluid velocity. They are subject to the constraints These hydro-like equations describe the evolution of the black hole in the large d limit.
They can be cast into a conservation law for a stress tensor on the hypersurface (with the 3 Note that we also expand the horizon radius as R+ = R

JHEP02(2018)167
normal n = ∂ R ) near the AdS boundary, The covariant derivative on the AdS boundary takes the form .
Then e and j constitute a stress tensor This form is very similar to that in [21]. We are not claiming that this stress tensor is the stress tensor of a dual field theory through the usual AdS/CFT dictionary. In fact, given the internal θ index, such an interpretation is obscure at best.
Zooming in on the equator. One troubling aspect of the conservation equations (3.3) is that ∂ θ T vθ , or equivalently ∂ θ j, does not appear at leading order in the 1/n expansion.
We can restore such a θ derivative by rescaling θ and j: (3.8) In the large n limit, this rescaling zooms in on the equator of the S d , where θ ∼ π/2. In refs. [20,21,24], a similar rescaling of the distance ruler was motivated in part by a desire to keep the sound modes in the fluctuation spectrum. For a conformal fluid, the speed of sound scales as 1/ √ n. By changing the distance ruler, the speed of sound is kept finite rather than infinitesimal, and a spatial gradient of j reappears in the energy conservation equation. In our case, the sound modes are gapped by the curvature of the sphere, but a spatial gradient of the current is restored in a similar way.
Solving the equations of motion order by order with this new scaling, we find Unlike the previous effective equations (3.3), we were not able to recast these equations in the form of a conservation law for a stress tensor.  the virtue of taking a derivative of the momentum-like quantity j + r c ∂ u e. They resemble the equations for the black strings [20], which is perhaps not surprising as stretching the polar angle θ makes the sphere look like a string.

Solving the unscaled system
We were able to find a large class of solutions to the ostensibly nonlinear equations (3.3). In fact, taking a logarithm of e(v, θ) leads to a linearization and an ability to superpose different solutions. The first equation can be solved to give an expression for j(v, θ) in terms of e(v, θ): (3.10) The following expression for e(v, θ) then gives a solution of the coupled system: where a ± and c are integration constants and the frequencies obey the dispersion relations This dispersion relation is precisely what we found earlier (2.22) for the quasinormal modes driving the instability. One choice of sign ω − always leads to a damped response. However, ω + can lead to an exponentially growing response when the radius is below the threshold for the corresponding instability, r 2 c ≤ 1 +1 . In figures 4 and 5, we plot snap-shots of these time evolving solutions for = 1 and 2.
For small a ± , we can expand out the exponential and consider the time dependence perturbatively. The quasinormal modes we considered before were constructed from spherical harmonics Y (Ω d ). The cos θ dependence we find here is in fact a large d limit of one  Figure 5. Snap shots of solutions (3.11) at different times. The almost uniform black holes evolve into the non-uniform black holes, but there is no end-point. We display the density plot on the sphere in insets. First row : We used the = 1 term with a 1+ > 0 (all other a ± = 0), and set the radius r c = 0.7. Second row : we used the = 2 term with a 2+ < 0 (all other a ± = 0), and set the radius r c = 0.5. of these spherical harmonics. As our metric ansatz depends only on the polar angle θ of the S d , we should be careful to select those spherical harmonics which have only such a θ dependence. Such spherical harmonics exist, and in the large d limit, they simplify to cos θ (see appendix D).
Detailed numerical calculations in the d = 5 case [9,10] indicate the existence of static black hole solutions that break the SO(6) symmetry of the S 5 down to SO(5). The black hole horizon radius becomes a function of the polar angle θ on the S 5 . In limiting cases, the black hole radius can even shrink to zero for a range of the θ-parameter, yielding black belt and spot solutions where the topology of the black hole changes to S 4 × S 4 and S 8 respectively. While our solutions (3.11) are in general strongly time dependent, they share some qualitative features with these spot and belt solutions. The odd solutions lead to single spots at the north or south pole, depending on the sign of a ± . The even solutions in contrast lead to either belt solutions or solutions with two spots -one at the north and one at the south pole -depending again on the sign of a ± .
We cannot make e(v, θ) arbitrarily small (or large) and stay within the regime of validity of our large d approximation. We need to keep |log e(v, θ)| n. Thus, we cannot strictly speaking witness a topology change in the shape of the black hole horizon. Nevertheless, the exponential growth and/or suppression in the solution (3.11) is very suggestive that such endpoints may exist, that there exist black spot and belt solutions with topology S 2d−2 and S d−1 × S d−1 respectively. One might conjecture that subleading terms in the JHEP02(2018)167 1/d expansion could halt the exponential growth of the solutions at late times, and also possibly lead to topology changes in the horizon.
At precisely the instability threshold ω + = 0 for integer , the corresponding cos θ mode can lead to a inhomogeneous static solution. Naively, for a general r c there exists a non-integer for which cos θ also leads to a inhomogeneous static solution. The issue with such solutions is that one needs an absolute value sign, |cos θ| , to keep the solution real on the whole sphere, and then there are corresponding non-analyticies in the solution at the equator and distributional terms which do not cancel out of the equations of motion.

Solving the rescaled system
We have not been able to find general time dependent solutions of the rescaled hydrodynamic like equation (3.9). We first investigate steady state solutions of the rescaled system. If we assume no time independence and make the ansatz then static configurations for e(u) must satisfy the following nonlinear equation: where f (u) ≡ e (u)/e(u) and denotes ∂ u . One very simple solution is a Gaussian config- In order for this solution to damp out as |u| gets large, we require r 2 c < 1 3 , which is the stability threshold for the = 2 mode. Because this static solution is symmetric under u → −u, it also exists on a hemisphere where we quotient out by this reflection symmetry about the equator, eliminating all of the odd parity modes, including the most unstable = 1 mode. It is therefore sensible that it is the = 2 mode which determines whether this static solution can exist or not. 5 We can then further look at fluctuations about the Gaussian solution. We look for fluctuations of the form e(v, u) = e s (u) + e −ivω δe(u) , We find solutions 5 We have found a couple of additional static solutions at discrete values of the horizon radius rc. At

JHEP02(2018)167
where and where H q (x) is a Hermite polynomial and c an arbitrary constant. 6 There are zero modes when . (3.22) In between these points, the frequency goes into the upper half plane, and the mode exponentially grows in time. The upper threshold for stability is the same as that for the = 1 mode, r 2 c = 1/2. Interestingly, the lower threshold q 1+2q is always greater than or equal to 1/ √ 3. We need r 2 c < 1/3 to have the Gaussian background solution, and so these fluctuations appear to be stable with respect to the Gaussian background.
We now investigate the full, nonlinear time evolution of these black hole solutions in the large d limit and in the rescaled variables. The global effective equations near the equator (3.9) are simple enough that they can be solved numerically with prepackaged differential equation solvers, such as the NDSolve function in Mathematica. 1 For these equations, provided we begin with a black hole radius r 2 c < 1/3 with a perturbation that preserves the u → −u symmetry, we can evolve to an approximate steady-state that closely matches the Gaussian profile (3.15). See figure 6 for one such numerical evolution. We say approximate because we have observed a different type of instability that kicks in at long times, where the Gaussian profile starts to move toward one pole or the other, spontaneously breaking the parity symmetry. We suspect this instability is a remnant of the = 1 instability we saw in the original θ variable, but we have not been able to find an analytic form for it in the rescaled case. In the regime 1/3 < r 2 c < 1/2, we do not find any approximate steady-state solutions that resemble Gaussians. Indeed, our earlier analysis indicated the Gaussian solution only exists for r 2 c < 1/3. Instead, the numerics suggests evolution toward a spot-like solution at the north or south pole. However, the time evolution crashes before a stationary solution is reached, at least in part because of the exponentially growing nature of the energy density at the poles. For r 2 c > 1/2, perturbations to the black hole solution die out, and the SO(d + 1) symmetric solution appears to be stable.
The rescaled equations (3.9) are second order rather than first order in the spatial derivative, necessitating the introduction of a second pair of boundary conditions. In addition to the Dirichlet condition on the current j = 0 at the poles, we introduce a Neumann condition on the energy density, that ∂ θ e = 0. These boundary conditions imply that the integral of e over the θ direction is a conserved quantity. There are some technical issues with these rescaled coordinates. The points u → ±∞ are singular and we can only solve the differential equations in some range −L < u < L. In practice, the boundary conditions are not applied at the poles but at u = ±L. A further difficulty is that with the Gaussian type solutions we find for e, the ratio j 2 /e in the differential equations can only be evaluated accurately provided L is not too big. A final issue here is that at any fixed and finite n, once u ∼ √ n, we need to include subleading terms in our large n expansion.

Entropy and conservation of energy
We examine energy and entropy constraints. From the first of the evolution equations (3.3), it follows that the following integral over the energy density is time independent, to leading order in n: E tot ≡ π 0 e(v, θ) sin n θ dθ . (3.23) The sin n θ factor is a bit like a Dirac delta function, peaked at θ = π/2. The behavior of e(v, θ) at the poles is thus only very weakly constrained by energy conservation in this large n limit. In this light, it is perhaps less surprising that the time evolution discussed above sometimes leads to diverging values of e(v, θ) near the poles. On the other hand, the value of e(v, π/2) is very strongly constrained. In fact, its time derivative vanishes, as can be seen directly from the equations (3.3).
In the microcanonical ensemble, the preferred black holes are the ones with the most entropy. In the AdS 5 × S 5 case, ref. [10] demonstrated numerically that the spot-like solutions have higher entropy. Unfortunately, we cannot distinguish our black holes using entropy at leading order in a large n expansion. At a fixed energy, all of our solutions have the same entropy because the two quantities are one and the same at leading order in 1/n. The entropy S can be derived from the horizon area of the black hole using the

JHEP02(2018)167
Bekenstein-Hawking formula. In the large n expansion, the leading term is given by S ∝ π 0 R + sin n θ dθ . (3.24) The horizon location is the energy density in our notation, R + = e(v, θ). A similar feature was noticed in a large d limit of planar black holes in ref. [21]. To calculate differences in entropy, we presumably need to calculate higher order terms in our large n expansion. The time evolution itself can indicate which configurations are preferred. In ref. [21], entropy production was encoded in second order derivative terms in the evolution equations. The type of allowed evolution was then constrained by the sign of the coefficients of these second order terms. In our case, if generic initial conditions evolve to spot-like rather than belt-like black hole solutions, we could tentatively conclude that the spot-like solutions are preferred. (In practice, of course, we find no stable endpoint to the time evolution of the equations (3.3) when r 2 c < 1/2.) A possible drawback is that the evolution equations (3.3) are purely first order, and thus may not encode any type of entropy production. The rescaled evolution equations (3.9), on the other hand, have some second derivative terms.
The rescaled equations (3.9) also lead to energy conservation, i.e. that is time independent, applying the boundary conditions that j = 0 and ∂ u e = 0 as |u| → ∞. At leading order in n, the entropy is again proportional to the total energy. Thus all the black holes connected by the evolution equations (3.9) have the same entropy at leading order in n. On the other hand, as these equations have second derivative terms, it may be more reliable to draw conclusions about the preferred endpoint of the instability from the time evolution of generic initial data.

Summary and outlook
We have applied the large d approximation to the study of black hole instability and its evolution in an AdS d × S d background. We derived analytic expressions for the instability threshold and quasinormal spectra order by order in the large d expansion. The results agree with the numerical calculations. We also studied a nonlinear ansatz for the metric, restricted to fluctuate along three of the 2d coordinates: time, the polar angle of the S d , and the radial coordinate of the AdS d . Simple hydrodynamic-like equations were found, one in the polar angle θ and the other zooming in on the equator of the sphere with u = θ − π 2 √ n. We find that the systems admit solutions very like the lumpy black spots and belts numerically constructed in ref. [9]. We provided simple analytic formula for these solutions at leading order in the large d limit. While the solutions in refs. [9,10] were static, our solutions are time dependent. In the stable regime r 2 + > 1/2, fluctuations quickly relax to the homogeneous Schwarzschild solution. In the unstable regime r 2 + < 1/2, we find exponentially growing solutions in the θ variable that do not approach a static endpoint. Despite their strong time dependence, the solutions do resemble at least qualitatively the static belt and spot solutions found in refs. [9,10]. In the rescaled u coordinate, we can JHEP02(2018)167 find a quasi-static solution. At intermediate times, with r 2 + < 1/3 below the threshold for the = 2 instability, the energy density is well fit by the static belt-like solution that we found. However, at long times, numerical evolution suggests this belt-like solution is unstable as well.
We would like to be able to figure out the ultimate fate of these unstable small black holes in AdS d ×S d . The answer to this question has implications both for general relativity and field theory. On the geometric side, it is interesting to understand what types of black hole horizon topologies are possible and how they may or may not arise through time evolution. The behavior at large d may shed light on issues concerning cosmic censorship, in particular on when and whether black holes can break apart. On the field theory side, through the AdS/CFT correspondence in the special case d = 5, the behavior of these small black holes maps out the phase diagram of maximally supersymmetric four dimensional SU(N ) Yang-Mills theory at large N and strong 't Hooft coupling in the micro-canonical ensemble.
Based on an analysis of static solutions in the d = 5 case, refs. [9,10] conclude that the branch of lumpy black holes emanating from the = 1 instability at r + = 0.4402 have lower entropy and higher energy than the homogeneous Schwarzschild solution, and are thus disfavored in the micro-canonical ensemble. However, they find an additional set of black hole solutions where the horizon has S 8 topology. This set has greater entropy than the homogeneous solution, and an energy range that extends both above and below the = 1 instability. Based on an extrapolation of this family, they conjecture that there exists a first order phase transition to these S 8 black holes from the homogeneous solutions at an r + > 0.4402. In fact, one can plausibly connect these S 8 solutions and lumpy black holes at a cusp, forming a classic swallowtail shape for a first order phase transition. See figure 7.
It is important to note that refs. [9,10] construct static solutions with no time dependence. While their conjectured phase transition is economical and reasonable, ideally one should be able to get to these static solutions through time evolution. Without evidence from time evolution, there remain other possibilities for the behavior of these small black holes [33]. There could be other static or stationary solutions that preserve less symmetry on the S 5 or involve other fields from type IIB supergravity. It could be that the instability never reaches any kind of static or even stationary final state, and the solution remains chaotically time dependent, or evolves to the string scale where general relativity ceases to be a useful description.
With regard to the arguments in the previous two paragraphs, our analysis comes up short. We only determined the time evolution at leading order in d. At this leading order, entropy and energy are the same. At fixed energy, we cannot determine which of our homogeneous, spot-like, and belt-like solutions are preferred with respect to entropy, as they all have the same entropy. We anticipate, however, that with some extra work subleading 1/d corrections to the entropy could be worked out and used to distinguish the different solutions.
We did not see dynamical evidence that stable spot-like solutions were preferred. Our failure to see these spot-like solutions, however, should not be construed as evidence that a large d generalization of the conjecture of refs. [9,10] is incorrect. One possibilitythe most likely in our view -is that the run-away behavior we did see could be stabilized by including 1/d corrections to the evolution equations. As several run-away simulations involve the energy density growing uncontrollably at one of the poles of the S d , it is plausible that stabilizing this run-away behavior will lead to a black hole with effectively S 2d−2 horizon topology at the corresponding pole.
If 1/d corrections stabilize these black holes and lead to stationary or static configurations at low energy, one may hope to recover a scenario qualitatively similar to the phase diagram of the black string [37]. In the black string case, the order of the phase transition depends on the number of dimensions. In d < 13 dimensions, the phase transition is first order, while in d > 13 dimensions, the transition becomes second order. Moreover, the critical dimension can be accurately computed in a large d expansion [38]. In our case, figure 7 suggests one may find a first order phase transition in d = 5, and it will be interesting to explore the story in d > 5.
Another possibility is that the 1/d corrections will not stabilize the run-away behavior, and that general relativity will cease to be a useful description once the curvatures get sufficiently large. One thing we can tentatively rule out is fractal behavior at large d, as observed for the black string in [35]. We did not see any evidence for fractal time evolution in our large d limit. The energy density evolved reliably to form one or two lumps on the S d with little or no substructure.
Like most of the literature on this subject of small black holes in anti-de Sitter, our ansatz limits functional dependence to θ, time, and the radius of AdS d . Moreover, we only turn on the metric and a self-dual d-form field strength. We can say very little at this point about solutions which involve more complicated coordinate dependence. We can also say JHEP02(2018)167 very little about solutions which generalize to large d type IIB supergravity solutions with more fields turned on. It should be straightforward and probably worthwhile to include an extra coordinate or two on the S d and AdS d and search for solutions that preserve less symmetry. It could be true, for instance, that a path to a static solution at late times, even if that solution preserves an SO(d) symmetry on the S d , may at intermediate times require the presence of less than SO(d) symmetry. Generalizing type IIB supergravity to large d, on the other hand, seems like a challenging and less well-defined enterprise.
From the preceding discussion, working out 1/d corrections to our results is a key next step. In addition to resolving some of the puzzles about entropy production and runaway behavior in the time evolution, one may hope that with enough subleading terms in the expansion, the results may provide not just a qualitative understanding of the most interesting case d = 5, but a quantitative one as well. A possible caveat is non-perturbative corrections. Ref. [19] argued that the 1/d expansion is not perturbative but asymptotic, that in certain cases corrections of order d −d/2 , or in some cases d −d , may appear.
Another important project for the future is to work out in more detail the entries of the holographic dictionary in this large d limit. For instance, how does the quantity e(v, θ) we called energy density relate to the actual energy density of the field theory through the usual prescription for calculating the stress tensor? In the case d = 5, the expectation values of scalar operators in the Yang-Mills theory can be read from the holographic stresstensor [9,10]. Can something similar be done in this large d limit?
A Equation of motion for χ 2 Given the linearized ansatz (2.6), the fluctuation of the rr component of the metric, χ 2 , follows the differential equation where we have definedḟ ≡ ∂ r f (r) andf ≡ ∂ 2 r f (r).

B Ricci tensor for a metric with two internal spaces
In this appendix, we consider the Ricci tensor for the following type of metric The metricḡ αβ represents the base space, whose coordinates are indexed with α, β, γ, · · · . The metricsg θ 1 θ 2 andg φ 1 φ 2 represent the internal spaces fibered over the base space. Their coordinate indices are θ 1 , θ 2 , · · · and φ 1 , φ 2 , · · · with the numbers of dimensions n A and n B . Capital letters run over the coordinates of the full space. We will derive the expression for the Ricci tensor using the tetrad formalism, in which the computation is simpler than in the coordinate basis. In terms of the vielbeins, we can write the line element as We use the Greek letters for the coordinate indices and the Latin letters for the vielbein indices. The vielbeinsē a ,ẽ k , andẽ l are those of the base and internal spaces, whose elements satisfȳ In terms of these vielbeins for the base and the internal spaces, we may write the full vielbeins as (with the one-form bases explicit)

JHEP02(2018)167
Recall that, in the tetrad formalism, the torsion and the curvature two-form are defined as The torsion free condition T M = 0, along with a bit of guessing, allows us to determine the spin connection ω M N . Because of the symmetry between two internal spaces and the antisymmetry of two indices M and N , we only need to determine These can be inferred from the a and k components of the torsion free condition. Writing out the a component, we have where ω a b is the spin connection on the base space, which satisfies the torsion free condition of its own. We also examine the k component, which is Again, we follow the obvious choice, whereω k 1 k 2 is the spin connection on the one of internal spaces. We are left with Since ω k 1 a is proportional to e k 1 = e k 1 θ 1 dx θ 1 (B.9), we are forced to choose The rest of spin connections can be determined by symmetry. Collecting the non-vanishing spin connections, we have From this result, we can proceed with the second of the defining relations (B.7) to obtain the curvature two-form. There are only four components that need to be determined, which are

JHEP02(2018)167
It is straightforward to show that where R a b and R k 1 k 2 are the curvature two-forms of the base and internal spaces. Writing out the ak 1 component reveals where, in the last step, we used the torsion free condition and the tetrad postulate, As before, the bars indicate objects in the base space. From these results, we can compute the Ricci tensor in a coordinate basis by converting it from the vielbein basis, With a bit of care in reading off the components of differential forms, and noticing that we finally find the Ricci tensor The components R φ 1 φ 2 are the same as R θ 1 θ 2 with simple exchanges, θ ↔ φ and A ↔ B.
Perhaps the most interesting part is the cross term (∂ α A)(∂ α B). This term would have been absent if there had been only one internal space in the fibration, and therefore hard to guess a priori. For the metric ansatz of this paper (3.1), the Ricci tensor and scalar can be written as

C Large n expansion of the metric and field strength
In this appendix, we present all the terms of the metric and d-form field that we have determined for the ansatz (3.1) in the large n expansion. It is convenient to introduce the one form f α and two form f αβ on the base space parametrized by (v, r, θ) to describe the self-dual F d form. The self-duality of F d implies Then the eight equations of motion for the metric can be neatly expressed as where α, β run through (v, r, θ), and i, j and m, n are the indices on the S d−2 and S d−1 .
For the d-form field, we have the four constraints These equations of motion are solved by the following metric and field solutions The e(v, θ) and j(v, θ) satisfy a pair of differential equations (3.3), which appear as the vv and vθ components of the Einstein equations at n 1 and n 0 orders. In addition, the vv component at n 0 order also gives a differential equation for e 1 (v, θ) and j 1 (v, θ). These solutions solve all of the metric and d-form equations of motion up to n 0 order.
Near the equator. We also provide the field and metric solutions in the rescaled polar angle u, which is related to θ as 7 θ = π/2 + u/ √ n. (C.5) We perform the coordinate transformation R = (r/r c ) n explicitly in the ansatz. Our ansatz is then ds 2 = g vv dv 2 +2r c n −1 R 1/n−1 dvdR+2g vu dvdu+g uu du 2 +g A dΩ 2 d−2 +g B dΩ 2 d−1 , (C.6) where α, β run over (v, R, u). The effective equations for e and j are the same whether we express the ansatz in the r or R. However, some of the Einstein equations are pushed to higher orders in 1/n in this setup, leaving fewer equations to solve. The self-duality of F d implies f vR = g −1/2 uu (−f R g vu + n −1 r c R 1/n−1 f u ), (C.7) f vu = g −1/2 uu nr −1 c R 1−1/n f R (−g 2 vu + g vv g uu ) + f u g vu − f v g uu ,  7 We can use an alternative rescaling, e −u 2 /2 = sin n θ, which is at leading order θ = π/2 + u/ √ n + O(n −3/2 ) and leads to the same effective equations for e and j. where The e(v, u) and j(v, u) obey the effective equations of motion (3.9), which are obtained from the vv and vu components of the Einstein equations at the n 1 and n 0 order. These solutions solve most of the equations of motion up to n 0 order. The equations that are not solved in terms of (e, j, e 1 , j 1 ) are the vv, ij, and mn components of the Einstein equations (C.2) at n 0 order. The ij and mn components are the same, which can be solved with standard techniques. We are then left with the vv component of the Einstein equations at n 0 order, which consists of g (2) vv , g (3) uu , g (2) A , g A , g B , and their first and second derivatives with respect to R, along with e, j, e 1 , and j 1 . The expression is too long to solve explicitly, but given the freedom to choose the aforementioned functions, we believe that it can be solved consistently with the boundary conditions. D Spherical harmonics in the large d limit We show that there exist a class of spherical harmonics that depend only on the polar angle θ of S d and that tend to Y ∼ (cos θ) in the large d limit. In general, the spherical harmonics on the sphere S d are given by