Lattice Black Branes: Sphere Packing in General Relativity
LLattice Black Branes: Sphere Packing in General Relativity ´Oscar J. C. Dias, ∗ Jorge E. Santos, † and Benson Way ‡ STAG research centre and Mathematical Sciences, University of Southampton, UK Department of Applied Mathematics and Theoretical Physics,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK Department of Physics and Astronomy, University of British Columbia,6224 Agricultural Road, Vancouver, B.C., V6T 1W9, Canada
We perturbatively construct asymptotically R , × T black branes with multiple inhomogeneousdirections and show that some of them are thermodynamically preferred over uniform branes inboth the microcanonical and canonical ensembles. This demonstrates that, unlike five-dimensionalblack strings, the instability of some unstable black branes has a plausible endpoint that does notrequire a violation of cosmic censorship. Introduction –
Higher-dimensional black holes aremarkedly different from those in four-dimensions. In par-ticular, they can be unstable [1–6], violate black holeuniqueness [2, 3, 7, 8] and also lead to a violation ofweak cosmic censorship [9–11].The prototypical representative of such behaviour isthe uniform five-dimensional black string, which is theproduct space of a Schwarzschild black hole with a cir-cle. The black hole therefore has spatial horizon topology S × S , and is asymptotically R , × S . When the massof the black string is much smaller than the size of thecircle, the S extent of the horizon is much smaller thanthe S length. Such a separation of length scales leads toan instability that was first discovered by Gregory andLaflamme in [1].What is the endpoint of this instability? When firstproposed by [1], there were two candidates: a nonuni-form black string or a localised black hole. Nonuniformblacks strings have the same horizon topology as uniformblack strings, but they no longer respect the symmetriesof the S . The nonuniform solutions branch off fromuniform solutions at a zero mode located at the criti-cal onset of the Gregory-Laflamme (GL) instability. Alocalised black hole has spherical horizon topology S .They were expected to exist since it is possible to placea small spherical black hole within R , × S . The possi-bility of evolution towards localised black holes was espe-cially intriguing because the change in horizon topologyimplies a violation of weak cosmic censorship.Later results [12, 13] ruled out the nonuniform stringsas an endpoint since they have lower entropy (horizonarea) than the uniform strings, while localised black holesremained a possibility since they are entropically pre-ferred [14]. The time evolution of the instability wasfinally performed in [9], providing the first numerical ev-idence of a violation of weak cosmic censorship.The Gregory-Laflamme instability is a generic featureof black objects with extended directions. The instabil-ity occurs in asymptotically flat black holes, as well astheories with compact directions such as those in stringtheory and AdS/CFT [15–18]. However, unlike the blackstring, most of these scenarios contain multiple extended directions. The situation is therefore more complex sincethere is a larger space of unstable perturbations, and con-sequently a larger space of nonuniform solutions. Wewish to shed light on which of these perturbations leadto the most entropically favourable nonuniform solutions,and whether these solutions can be dominant over uni-form solutions.Consider then the six-dimensional black braned s = − (cid:16) − r r (cid:17) d t + d r − r r + r dΩ +d x +d x , (1)where r is the horizon radius. The linear stability cal-culation for this system is identical to the black string[1]. Taking a metric perturbation ansatz g µν = g (0) µν + e − iωt + ik x + ik x h µν , one finds that for fixed r , the fre-quency ω depends only on the norm k + k . In particular,the zero mode ω = 0 occurs at a critical k GL = (cid:112) k + k .Each frequency, including the zero mode, therefore con-tains a circle of degenerate perturbations. Because of thisdegeneracy, it is possible for multiple static solutions toappear from the same zero mode.To restrict some of these possibilities, let us take the x and x directions to be identified on an oblique toruswith a single length scale Lx ∼ x +( n + γn ) L , x ∼ x + n (cid:112) − γ L , (2)for any integers n and n and angle γ ∈ [0 , D -dimensional black strings consisting ofa D − D ≥
14 [19, 20]. Time evolution within a large D ex-pansion indicates that evolution proceeds towards thenonuniform solution, avoiding a violation of cosmic cen-sorship. Black branes on oblique torii were also studied at a r X i v : . [ h e p - t h ] D ec θ FIG. 1: The fundamental torus region is given by the solidblack parallelogram with side lengths L , which is angledat θ = arcsin γ . The arrows show the longest wavelengthwavevectors 2 π(cid:126)k/ | k | that generate parallel lines in the lat-tice. lowest order in large D [21], but unfortunately differentstatic solutions have the same thermodynamic quantitiesat this order.We also mention that related holographic lattice so-lutions in AdS were studied in [22], where a triangularlattice is thermodynamically preferred. Perturbative and numerical methods –
The met-ric perturbations only depend on the norm k + k , butthe torus boundary conditions do not allow arbitrarywavevectors. Each wavevector (cid:126)k = ( k , k ) describes aplane wave which must fit within the torus. The allow-able wavevectors with the longest wavelengths can beviewed in Fig. 1 where the torus has been extended toa periodic lattice spanning R . The parallel lines on thislattice allow us to visualise the permissible plane waveswith the longest wavelengths. We see that the torus gen-erates three sets of parallel lines, corresponding to threedifferent wavevectors. By demanding that the torus hasa single length scale L , two of these wavevectors musthave the same norm. The third wavector also has thisnorm only when γ = 1 / x = L π ( x + γy ) , x = L π (cid:112) − γ y , (3)where the torus identifications (2) imply that x ∼ x + 2 πn x , y ∼ y + 2 πn y . (4)In these coordinates, a static metric perturbation mustbe of the form g µν = g (0) µν + e in x x + in y y h µν , for integers n x and n y . Mapping that back into x and x coordinates,we can identify the wavevector components k and k .The norm of the wavevector is given by | k | = (cid:113) k + k = 1 λ (cid:113) n x + n y − n x n y γ , (5)where we have defined the length scale λ = (cid:112) − γ π L. (6)When L is sufficiently small relative to r , there is noinstability. Now consider increasing L . Since the uniformblack brane is unstable for perturbations with | k | < k GL ,the smallest | k | (longest wavelength) perturbations willbe the first to generate static nonuniform branes. Wetherefore wish to determine which integers n x and n y generate the smallest | k | perturbations. These are n x = ± , n y = 0; and n x = 0 , n y = ± , (7)which have | k | = 1 /λ , and n x = n y = ± , (8)which has | k | = 2(1 − γ ) /λ . Both sets of perturbationshave the same | k | precisely at γ = 1 /
2. Naturally, thereare shorter wavelength perturbations than those in (7)and (8). These perturbations only become unstable forlarger torii, and we do not consider them here.At this point, we would like to distinguish between twotypes of perturbations that generate static nonuniformsolutions: those that have a single wavevector, and thosethat are a linear combination of different wavevectors.Let us call nonuniform branes generated by the former as“semi-nonuniform” black branes, and those by the lateras “fully nonuniform” black branes. Semi-nonuniformblack branes preserve translation invariance in one di-rection and are therefore equivalent to the product of a5-dimensional nonuniform black string and an extra flatdirection. Fully nonuniform black branes have no suchsymmetry.In this work, we focus on the fully nonuniform blackbranes. For γ (cid:54) = 1 /
2, such black branes can be gen-erated by the two perturbations in (7). For γ = 1 / x and y coordinates as well as a newradial coordinate from the transformation r = r − z . (9)Our metric ansatz is given byd s = r (cid:26) − z q d t r + 4 q d z (1 − z ) + q (1 − z ) dΩ S + λ (cid:20) q − γ (cid:18) d x + γq d y + q λ d z (cid:19) + q (cid:18) d y + q λ d z (cid:19) (cid:21)(cid:27) , (10)where q ’s are functions of { x, y, z } , and we have also de-fined λ = λ/r . Note that implicit in this ansatz is ahorizon at z = 0 and asymptotic infinity at z = 1. Theuniform black brane (1) can be obtained in these newcoordinates by setting q i ( z, x, y ) = ¯ q i with¯ q , , , , , = 1 , ¯ q , = 0 . (11)The eight functions q i are obtained by solving the vac-uum Einstein equation R µν = 0, subject to the appro-priate boundary conditions. We will do this perturba-tively about the black brane solution stopping at the or-der where thermodynamic quantities become corrected.We use the Einstein-DeTurck formalism [23–26], whichis valid non-perturbatively as well. This method requiresa choice of reference metric g , which contains the samesymmetries and causal structure as the desired solution.The reference metric we choose is the black brane metricgiven by (11). The DeTurck method modifies the Ein-stein equation R µν = 0 into R µν − ∇ ( µ ξ ν ) = 0 , ξ µ ≡ g αβ [Γ µαβ − Γ µαβ ] , (12)where Γ and Γ define the Levi-Civita connections for g and ¯ g , respectively. Unlike R µν = 0, this equationyields a well-posed elliptic boundary value problem. Itwas proved in [24] that static solutions to (12) necessarilysatisfy ξ µ = 0, and hence are also solutions to R µν = 0.We now discuss the boundary conditions. At z = 1, thesolution must be asymptotically R , × T which requires q i (cid:12)(cid:12) z =1 = ¯ q i . At the horizon z = 0, we require regularityand impose q (cid:12)(cid:12) z =0 = q (cid:12)(cid:12) z =0 , ∂ z q i (cid:12)(cid:12) z =0 = 0 for i = 2 , . . . , q (cid:12)(cid:12) z =0 = 0 = q (cid:12)(cid:12) z =0 . Finally, we impose periodicboundary conditions on x and y according to (4).Since r is just an overall scale, our ansatz (10) isparametrised by λ and γ . However, it’s more convenientin our perturbative calculation to replace λ by | k | . λ , λ or L can be recovered via (5), (6), and λ = λ/r .Now expand the metric functions and | k | in powers of (cid:15) : q i = ¯ q i + ∞ (cid:88) n =1 (cid:15) n q ( n ) i ( z, x, y ) , r | k | = k GL + ∞ (cid:88) n =1 (cid:15) n k ( n ) , (13)where recall that the instability is critical when r | k | = k GL . The periodicity of the torus allows us to express the x and y dependence of q ( n ) i as a sum of Fourier modes. The expansion (13) is placed into the Einstein DeTurckequation (12) and expanded order by order in (cid:15) .At O ( (cid:15) ), we have an eigenvalue problem for { k GL , q (1) i } , where all perturbations with the same r | k | = k GL give the same eigenvalue problem. We there-fore have some freedom to choose the perturbation, whichwill affect results at higher order. We are only interestedin perturbations generated by multiple wavevectors withthe same | k | .For γ (cid:54) = 1 /
2, we take perturbations generated by (7)which takes the form q (1)1 , , = f (1)1 , , ( z ) [cos( x ) − cos( y )] , (14)while for γ = 1 /
2, we take a perturbation generated byboth (7) and (8), which goes as q (1)1 , , = f (1)1 , , ( z ) [cos( x ) − cos( y ) + cos( x + y )] . (15)Up to symmetries, the perturbations above are actuallyquite general. We have used translation invariance to fixthe phases, and the relative amplitudes (including theamplitude of the sin( x + y ) and cos( x + y ) terms in the γ = 1 / ξ (1) µ = 0.As expected, the equations of motion at linear orderfor both of these perturbations are identical, and inde-pendent of γ . They consist of two algebraic relations thatcan be used to determine q (1)1 , ( z ) in terms of q (1)2 ( z ) and q (1) (cid:48) ( z ), as well as the following eigenvalue problem with k appearing as an eigenvalue: q (1) (cid:48)(cid:48) ( z ) + 3 + 4 z − z z − z + 3 z q (1) (cid:48) ( z )+4 8 (cid:0) − z (cid:1) − k (cid:0) − z (cid:1) (1 − z ) (1 − z ) q (1)2 ( z ) = 0 . (16)Together with the boundary conditions that q (1)2 must beregular at the horizon z = 0 and vanish exponentiallyat the asymptotic boundary z = 1, we can solve thisproblem numerically using methods described in detailin [26]. We find that k GL (cid:39) . γ = 1 / γ (cid:54) = 1 / O ( (cid:15) ), theEinstein DeTurck equation is sourced by the linear ordersolution { k GL , q (1) i } . Then the allowed Fourier modes forthe n = 2 functions come from squaring the linear com-bination of modes we had at linear order. For example,in the γ = 1 / q (2) i for i = 1 , . . . , q (2) i ( z, x, y ) = f (2 , i ( z ) + f (2 , i ( z ) cos( x ) (17)+ f (2 , i ( z ) cos( y ) + f (2 , i ( z ) cos( x + y )+ f (2 , i ( z ) cos( x − y ) + f (2 , i ( z ) cos(2 x )+ f (2 , i ( z ) cos(2 y ) + f (2 , i ( z ) cos[2( x + y )]+ f (2 , i ( z ) cos(2 x + y ) + f (2 , i ( z ) cos( x + 2 y ) , and for i = 7 , → sin and f (2 , , ( z ) = 0. There are a total of 78 functions f (2 ,α ) i ( z )which need to determined. Since Fourier coefficient de-couple, each fixed α gives an independent system of dif-ferential equations to be solved numerically using themethods detailed in [26]. Note that three of these sys-tems with α = 1 , , i.e. the three corresponding toFourier coefficients that were also excited at linear or-der) depend upon the correction k (1) . The solution totheir differential equations will determine both f (2 ,α ) i ( z )and k (1) independently. Since k (1) is unique, it is a con-sistency check to verify that all three systems yield thesame k (1) , which they do within machine precision. Nu-merical convergence is shown in the Appendix.The computation continues in a similar manner forhigher orders in (cid:15) . Details can be found in the Appendix.In the end, we stop at O ( (cid:15) ) and find k GL (cid:39) . , k (1) (cid:39) − . ,k (2) (cid:39) . , k (3) (cid:39) . . (18)At O ( (cid:15) ), we can compute the first perturbative correc-tions to various thermodynamic potentials, and we havedecided not continue to higher order.The calculation for γ (cid:54) = 1 / k ( n ) = 0 when n is odd, and that one needs to reach O ( (cid:15) ) in perturba-tion theory to find thermodynamic corrections. To ob-tain our results, the ODEs are solved numerically for anumber of specific values of γ . Thermodynamics –
We now compute the thermo-dynamic quantities of the perturbative solutions. Thetemperature of (10) is T H = 1 / (4 πr ) and the entropy S H is the horizon area divided by 4 G N ( G N is the six-dimensional Newton’s constant). The energy E is com-puted from the asymptotic decay of the gravitational fieldusing the Hamiltonian formalism presented in [27]. TheHelmoltz free energy is then F = E − T H S H . We computethe dimensionless densities in terms of the torus volume V T = r π λ (1 − γ ) − / : τ H = T H V / T , σ H = S H G N /V T , E = EG N /V / T , f = F G N /V / T . (19)As we have mentioned earlier, semi-nonuniform blackbranes in this scenario are never thermodynamically pre-ferred over uniform solutions. Therefore, we only need to FIG. 2: Phase diagram in the microcanonical ensemble for γ = 1 /
2: entropy density difference ∆ σ H between the tri-angular (upper blue curve) or hexagonal (lower black curve)nonuniform black branes and the uniform black brane for agiven energy density E . The green dot locates where the uni-form black brane is critically unstable, with lower energiesbeing fully unstable. compare our fully nonuniform solutions with the uniformones. To compare these solutions in the microcanonicalensemble, we take the entropy difference at the same en-ergy ∆ σ H = σ H ( E ) − σ H ( E ), where σ H and σ H are theentropy densities of the fully nonuniform, and uniformsolutions, respectively. For comparisons in the canonicalensemble, we similarly define the free energy difference∆ f = f ( τ H ) − f ( τ H ).For lattices with γ = 1 /
2, to order O (cid:0) (cid:15) (cid:1) , the thermo-dynamic densities are∆ σ H (cid:39) . (cid:15) ,τ H (cid:39) . . (cid:15) − . (cid:15) − . (cid:15) , E (cid:39) . − . (cid:15) + 0 . (cid:15) +0 . (cid:15) , ∆ f (cid:39) − . (cid:15) . (20)The difference in entropy densities between the γ =1 / E = E GL =3 / k GL / (4 √ π ) (cid:39) . γ = 1 / σ H ∼ (cid:15) in(20). The fact that this power is odd rather than even isequivalent to the fact that δg and − δg are physically FIG. 3: Coefficient c ∆ σ H of the entropy density difference∆ σ H ∼ c ∆ σ H (cid:15) as a function of the lattice angle γ = sin θ for γ (cid:54) = 1 /
2. Perturbation theory breaks down at γ = 1 / distinct perturbations. This distinction can be inter-preted as the difference between triangular and hexago-nal arrangements. We contrast this with the nonuniformstrings, where translation invariance implies that δg and − δg are physically equivalent [12, 13, 19]. Nonuniformstrings therefore only extend towards higher or lower en-ergies (not both), depending on the dimension.Consider now the γ (cid:54) = 1 / σ H = c ∆ σ H (cid:15) + . . . , ∆ f = c ∆ f (cid:15) + . . . , (21)where the coefficients c ∆ σ H and c ∆ f depend on γ . Notethat we conclude that the even power of (cid:15) in (21) impliesthat each γ yields a single branch of lattice solutions (nottwo, like the γ = 1 / c ∆ σ H as a function of γ . Observe that c ∆ σ H < γ . Therefore, fully nonuniform blackbranes with those values of γ do not dominate the micro-canonical ensemble near the zero mode. However, there isa small window 1 / < γ (cid:46) .
538 where c ∆ σ H >
0. In thiswindow, the fully nonuniform branes are preferred overthe uniform phase. Moreover, we also find that withinthis window, the fully nonuniform brane solutions ex-tend from the zero mode towards lower energies E < E GL where the uniform solutions are unstable, but otherwiseextends towards higher energies. Discussion –
To summarise our results thus far, wehave found black branes in R , × T that are fully inho-mogeneous in the torus directions by perturbing aboutthe zero mode. These separate into the case where thetorus contains triangular/hexagonal symmetry γ = 1 / γ (cid:54) = 1 /
2. In the γ = 1 / γ (cid:54) = 1 / γ = 1 / γ > / E < E GL and τ H > τ GL will evolve in time towardsa fully-nonuniform black membrane.However, localised solutions are expected to exist thatmight compete with some of these nonuniform solu-tions. In the space of static solutions, nonuniform blackstrings are known to eventually connect to localised blackholes (which have spatial horizon topology S ) througha topology-changing conical transition [14]. Likewise inour case, we expect spherical S black holes to exist.However, a direct transition from fully nonuniform blackbranes to localised black holes requires a transition wherean entire circle pinches off on the horizon simultaneouslyin moduli space. We find such a scenario unlikely dueto the lack of symmetry. Instead, it is more plausiblethat the horizon pinches off on specific points, leading tomesh-shaped black objects resembling connected blackstrings. Those might later proceed through a secondtopology transition towards localised black holes. Sort-ing out this phase diagram requires constructing thesehypothetical localised solutions, as well as extending ourresults for nonuniform branes to nonperturbative regions.Note that the extra symmetries of the γ = 1 / S radius ofthe black brane is shown in Fig. 4, where it is possibleto extrapolate which arrangement of localised black holesour fully nonuniform black branes are approaching. Wesee that the higher-entropy branch of these black branesapproaches the triangular arrangement (left panel), whilethe lower-entropy branch approaches a hexagonal ar-rangement (right panel). Because a hexagonal arrange-ment of localised black holes would contain two blackholes per torus volume, we expect this situation to beunstable.It is worth considering what happens to these solutionsat higher dimensions, where the Schwarzschild black holeis replaced by Schwarzschild-Tangherlini. In D ≥ γ > /
2, there is an additionalsemi-nonuniform solution coming from the perturbations(8) since they have longer wavelength than those of (7)that generated the fully nonuniform solutions.We also mention that our methods here are not spe-cific to two brane directions, though the situation growsin complexity with increasing brane dimension. It is in-teresting to note that for two brane directions, the mostentropic nonuniform phase moves towards the triangu-lar arrangement of black holes, which is also the densest
FIG. 4: Structure of triangular and hexagonal nonuniform black branes ( γ = 1 / (cid:15) = ± − ). Radius of the S at the horizon(which is proportional to the square root of the entropy density) as a function of the lattice directions x and x for the blackbrane with triangular (left) and hexagonal (right) arrangements. The lighter yellow regions represent the locations of the latticepoints where the entropy density is higher. sphere packing in two dimensions. It is therefore tempt-ing to propose a connection between thermodynamicallypreferred Kaluza-Klein black holes and the mathemati-cal problem of sphere packing. However, we expect thepreferred localised black hole phase to contain only oneblack hole per torus volume. If there is indeed a con-nection to sphere packing, it would be the problem ofpacking spheres on fundamental lattice sites, and not themore sophisticated general sphere packing problem.Our work raises interesting questions for the study ofunstable gravitational systems with multiple extended di-rections. Although our study was restricted to asymp-totically Kaluza-Klein flat spacetimes, we expect similarphysics to be present in the context of asymptoticallylocally anti-de Sitter spacetimes. For instance, the anal- ysis of the so called spinoidal instability performed in[28] and [29] was restricted to co-homogeneity 2 + 1, andour work raises the possibility of finding interesting newphases with fewer symmetries. It would also be inter-esting to pursue the phase structure of the holographicduals of p − branes compactified on T d torus which wasinitiated in [30]. Acknowledgements –
O.J.C.D. is supported bythe STFC Ernest Rutherford grants ST/K005391/1 andST/M004147/1. O.J.C.D. further acknowledges supportfrom the STFC “Particle Physics Grants Panel (PPGP)2016”, grant ref. ST/P000711/1. J.E.S. was supported inpart by STFC grants PHY-1504541 and ST/P000681/1.B.W. is supported by NSERC.
Appendix
Details of the perturbation theory for γ = 1 / Here, we provide details to our perturbative calculation for the γ = 1 / i.e. triangular/hexagonal lattices).Our method is similar to the one used by [12, 13, 19] to explore the existence of asymptotically R , × S non-uniformblack strings.We perturb the metric functions in a power series written as q i ( z, x, y ) = ¯ q i + ∞ (cid:88) n =1 (cid:15) n q ( n ) i ( z, x, y ) , r | k | = k GL + ∞ (cid:88) n =1 (cid:15) n k ( n ) , (22)where ¯ q i describes the uniform black membrane. The torus implies that the perturbation functions q ( n ) i can beexpanded in a Fourier series.But triangular/hexagonal lattice symmetries (as shown in Fig. 4) limit the number of Fourier coefficients that areavailable. As explained in the discussion of of (7)-(8) of the the main text, for γ = 1 / x ) , cos( y ) and cos( x + y ). The equations of motion and boundary conditions at O ( (cid:15) ) imply k GL (cid:39) . q (1)1 ( z, x, y ) = (cid:32) z + 11 − z q (1)2 ( z ) + z (cid:0) − z (cid:1) − z q (1) (cid:48) ( z ) (cid:33) (cid:2) cos( x ) − cos( y ) + cos( x + y ) (cid:3) ,q (1)2 ( z, x, y ) = q (1)2 ( z ) (cid:2) cos( x ) − cos( y ) + cos( x + y ) (cid:3) ,q (1)3 ( z, x, y ) = − (cid:32) z + 11 − z q (1)2 ( z ) + z (cid:0) − z (cid:1) − z ) q (1) (cid:48) ( z ) (cid:33) (cid:2) cos( x ) − cos( y ) + cos( x + y ) (cid:3) ,q (1) i ( z, x, y ) = 0 , for i = 4 , · · · , k GL and q (1)2 ( z ) are the (numerical) solutions of the eigenvalue equation q (1) (cid:48)(cid:48) ( z ) + 3 + 4 z − z z − z + 3 z q (1) (cid:48) ( z ) + 4 8 (cid:0) − z (cid:1) − k (cid:0) − z (cid:1) (1 − z ) (1 − z ) q (1)2 ( z ) = 0 . (24)The linear order solution (23)-(24) now sources the second order solution at O ( (cid:15) ). The Fourier modes that areexcited at this order come from squaring the sum of the three O ( (cid:15) ) modes. In full, they can be written q (2) i ( z, x, y ) = f (2 , i ( z ) + f (2 , i ( z ) cos( x ) + f (2 , i ( z ) cos( y ) + f (2 , i ( z ) cos( x + y )+ f (2 , i ( z ) cos( x − y ) + f (2 , i ( z ) cos(2 x ) + f (2 , i ( z ) cos(2 y ) + f (2 , i ( z ) cos[2( x + y )]+ f (2 , i ( z ) cos(2 x + y ) + f (2 , i ( z ) cos( x + 2 y ) , for i = 1 , · · · , ,q (2) i ( z, x, y ) = 0 + f (2 , i ( z ) sin( x ) + f (2 , i ( z ) sin( y ) + f (2 , i ( z ) sin( x + y )+ f (2 , i ( z ) sin( x − y ) + f (2 , i ( z ) sin(2 x ) + f (2 , i ( z ) sin(2 y ) + f (2 , i ( z ) sin[2( x + y )]+ f (2 , i ( z ) sin(2 x + y ) + f (2 , i ( z ) sin( x + 2 y ) , for i = 7 , , (25)We wish to determine the various f (2 ,α ) i ( z ), with α = 0 , · · · f (2 , , ( z ) = 0). Since eachFourier mode α decouples from the others, we solve an independent coupled ODE system for each α (a second-ordersystem of 6 ODEs for α = 0 and 8 for the others). These ODEs here, and at higher orders in (cid:15) take the form O ij q ( n,α ) j = S ( n,α ) i + k ( n − s ( n,α ) i (26)where O ij is a second order differential operator and S ( n,α ) i , s ( n,α ) i are sourced by the lowest order solution (23)-(24).For α = 1 , , i.e. precisely the modes that were excited at linear order, one has s ( n,α ) i (cid:54) = 0 so we can use each ofthese systems to determine k (1) . We verify for consistency that we obtain the same k (1) − namely the one given in(18) − from each α = 1 , , O ( (cid:15) ), we have a total of 19 Fourier modes, each ofthe form cos( A (3 ,α ) ), where A (3 ,α ) ∈ { , x, y, x + y, x − y, x, y, x + y ) , x + y, x +2 y, x, y, x + y ) , x − y, x − y, x + y, x +3 y, x +2 y, x +3 y } , (27)which are precisely those that follow from taking the cubic power of the sum of the three linear-order modes. Theperturbation functions therefore take the form q (3) i ( z, x, y ) = f (3 , i ( z ) + f (3 , i ( z ) cos( x ) + f (3 , i ( z ) cos( y ) + f (3 , i ( z ) cos( x + y )+ f (3 , i ( z ) cos( x − y ) + f (3 , i ( z ) cos(2 x ) + f (3 , i ( z ) cos(2 y ) + f (3 , i ( z ) cos[2( x + y )]+ f (3 , i ( z ) cos(2 x + y ) + f (3 , i ( z ) cos( x + 2 y ) + f (3 , i ( z ) cos(3 x ) + f (3 , i ( z ) cos(3 y )+ f (3 , i ( z ) cos[3( x + y )] + f (3 , i ( z ) cos(2 x − y ) + f (3 , i ( z ) cos( x − y )+ f (3 , i ( z ) cos(3 x + y ) + f (3 , i ( z ) cos( x + 3 y ) + f (3 , i ( z ) cos(3 x + 2)+ f (3 , i ( z ) cos(2 x + 3 y ) , for i = 1 , · · · , ,q (3) i ( z, x, y ) = 0 + f (3 , i ( z ) sin( x ) + f (3 , i ( z ) sin( y ) + f (3 , i ( z ) sin( x + y )+ f (3 , i ( z ) sin( x − y ) + f (3 , i ( z ) sin(2 x ) + f (3 , i ( z ) sin(2 y ) + f (3 , i ( z ) sin[2( x + y )]+ f (3 , i ( z ) sin(2 x + y ) + f (3 , i ( z ) sin( x + 2 y ) + f (3 , i ( z ) sin(6 πx ) + f (3 , i ( z ) sin(3 y )+ f (3 , i ( z ) sin[3( x + y )] + f (3 , i ( z ) sin(2 x − y ) + f (3 , i ( z ) sin( x − y )+ f (3 , i ( z ) sin(3 x + y ) + f (3 , i ( z ) sin( x + 3 y ) + f (3 , i ( z ) sin(3 x + 2)+ f (3 , i ( z ) cos(2 x + 3 y ) , for i = 7 , , (28)We now have 150 functions f (3 ,α ) i ( z ) to solve for. As before, each Fourier mode α decouples from all others, so wehave 19 independent systems of ODEs to solve (one for each α ). These ODEs again take the form (26). Again, thesystems corresponding to α = 1 , ,
3, and only these, have s (3 ,α ) i (cid:54) = 0, and hence depend on k (2) . Therefore, for eachof these three α ’s we solve the equations of motion to find f (3 ,α ) i ( z ) and k (2) . These three systems of equations areindependent but k (2) is unique, so we must to get the same k (2) in each of them. This k (2) is presented in (18).To find the wavevector correction k (3) , which is required to compute the thermodynamic quantities up to O ( (cid:15) ),we still need to find the metric solutions at order O ( (cid:15) ) since k (3) only appears in the equations of motion at thisorder (in practice it is enough to analyse the cos x, cos y and cos( x + y ) Fourier sectors since k (3) only appears in theequations associated to these modes). At this order, a total of 31 Fourier modes are excited. Schematically, each ofthese is of form cos( A (4 ,α ) ), with A (4 ,α ) ∈ { , x, x, x, x, y, y, y, y, x ± y, x ± y ) , x + y ) , x + y ) , x ± y, x ± y, x + 2 y ) , x + y ) ,x ± y, x ± y, x + y, x + 4 y, x + 2 y, x + 3 y, x + 3 y, x + 4 y } , (29)which are precisely those that follow from taking the fourth power of the sum of the three linear-order modes in (23).Solving the corresponding ODE system we determine the associated Fourier coefficients f (4 ,α ) i ( z ) for α = 0 , · · · , k (3) (cid:54) = 0 given in (18).We stop our analysis at this order O ( (cid:15) ), since this is the first order where the leading contribution to the entropydifference ∆ σ H , and free energy difference ∆ f are obtained. Namely, the entropy, temperature, energy and free energydensities of the γ = 1 / σ H (cid:39) . − . (cid:15) + 0 . (cid:15) + 0 . (cid:15) ,τ H (cid:39) . . (cid:15) − . (cid:15) − . (cid:15) , E (cid:39) . − . (cid:15) + 0 . (cid:15) + 0 . (cid:15) , f (cid:39) . − . (cid:15) + 0 . (cid:15) + 0 . (cid:15) . (30)To determine which phase is thermodynamically preferred in a given ensemble, we must compare (30) with those forthe uniform black membrane. In the microcanonical ensemble, one keeps the energy density fixed and the dominantsolution is the one with higher entropy density. In the canonical ensemble, one keeps the dimensionless temperaturefixed and the dominant solution is the one with lower free energy density. For this analysis, it is thus important towrite the entropy density of the uniform black brane as a function of its energy density, and its free energy density asa function of its dimensionless temperature: σ H ( E ) (cid:12)(cid:12) unif = 4 π E , f ( τ H ) (cid:12)(cid:12) unif = 116 π τ H . (31)For the microcanonical analysis, we can now compute the entropy density difference ∆ σ H = σ H ( E ) − σ H ( E )between nonuniform and uniform black branes with the same energy density E to get the result ∆ σ H (cid:39) . (cid:15) FIG. 5: Phase diagram in the canonical ensemble: dimensionless free energy density difference ∆ f between the triangular (lowerblue curve) and hexagonal (upper dashed black curve) nonuniform black membranes and the uniform black membrane for agiven dimensionless temperature τ H . Solutions with the lower free energy are preferred. The green dot locates the zero-modeof the instability of uniform branes, with higher temperatures being unstable. presented in the main text. Similarly for the canonical ensemble, we can compute the free energy density difference∆ f = f ( τ H ) − f ( τ H ) at the same dimensionless temperature τ H to get the result ∆ f (cid:39) − . (cid:15) presented inthe main text.The microcanonical phase diagram was presented in the main text. For completeness, here we show the canonicalphase diagram ∆ f vs τ H in Fig. 5. Details of the perturbation theory for γ (cid:54) = 1 / Just as we did for γ = 1 /
2, in the γ (cid:54) = 1 / γ = 1 / x )and cos( y ). The equations of motion and boundary conditions at O ( (cid:15) ) imply k GL (cid:39) . q (1)1 ( z, x, y ) = (cid:32) z + 11 − z q (1)2 ( z ) + z (cid:0) − z (cid:1) − z q (1) (cid:48) ( z ) (cid:33) (cid:2) cos( x ) − cos( y ) (cid:3) ,q (1)2 ( z, x, y ) = q (1)2 ( z ) (cid:2) cos( x ) − cos( y ) (cid:3) ,q (1)3 ( z, x, y ) = − (cid:32) z + 11 − z q (1)2 ( z ) + z (cid:0) − z (cid:1) − z ) q (1) (cid:48) ( z ) (cid:33) (cid:2) cos( x ) − cos( y ) (cid:3) ,q (1) i ( z, x, y ) = 0 , for i = 4 , · · · , k GL and q (1)2 ( z ) are the (numerical) solutions of the eigenvalue equation (24). Note that with respect to the γ = 1 / q (1)2 ( z ) satisfies the same eigenvalue equation (24) and q (1) i ( z ) and k GL are the same for all γ ’s.The linear order solution of (24) and (32) now sources the solution at O ( (cid:15) ). The excited Fourier modes at thisorder come from squaring the sum of the two O ( (cid:15) ) modes. In total, they can be written as q (2) i ( z, x, y ) = f (2 , i ( z ) + f (2 , i ( z ) cos(2 x ) + f (2 , i ( z ) cos(2 y )+ f (2 , i ( z ) cos( x + y ) + f (2 , i ( z ) cos( x − y ) , for i = 1 , · · · , ,q (2) i ( z, x, y ) = 0 + f (2 , i ( z ) sin(2 x ) + f (2 , i ( z ) sin(2 y )+ f (2 , i ( z ) sin( x + y ) + f (2 , i ( z ) sin( x − y ) , for i = 7 , . (33)We need to determine the several f (2 ,α ) i ( z ), with α = 0 , · · · f (2 , , ( z ) = 0), for each γ (cid:54) = 1 /
2. Each Fourier mode α decouples from the others, so the equations reduce to an independent coupled ODE0system for each α (a second-order system of 6 ODEs for α = 0 and 8 for the others) and γ . These ODEs, as well asthose at higher orders in (cid:15) , take the form (26) where O ij is a second order differential operator and S ( n,α ) i , s ( n,α ) i aresourced by the lowest order solution (32).At the current order, n = 2, one has s ( n,α ) i = 0 for all the Fourier families ( i.e. for the five α ’s). In particular, thisimplies that k (1) = 0 (which was not the case for γ = 1 / x and cos y to (33). However, since these modes are not sourced by the two O ( (cid:15) ) modes, we find that theassociated ODE system is of the form (26) with S ( n,α ) i = 0. There is nevertheless a source term proportional to k (1) , i.e. s ( n,α ) i (cid:54) = 0, which results from the k -expansion in (22). However, the system only has the trivial solution k (1) = 0and vanishing eigenfunctions. The structure of the problem is such that k ( n ) = 0 for odd order n .The calculation at higher orders follows a similar process. At O ( (cid:15) ), we have a total of 8 Fourier modes, each ofthe form cos( A (3 ,α ) ), where A (3 ,α ) ∈ { x, y, x, y, x + y, x + 2 y, x − y, x − y } , (34)which are precisely those that follow from taking the cubic power of the sum of the two linear-order modes in (32).Namely, the perturbation functions take the form q (3) i ( z, x, y ) = f (3 , i ( z ) cos( x ) + f (3 , i ( z ) cos( y ) + f (3 , i ( z ) cos(3 x ) + f (3 , i ( z ) cos(3 y ) + f (3 , i ( z ) cos(2 x + y )+ f (3 , i ( z ) cos( x + 2 y ) + f (3 , i ( z ) cos(2 x − y ) + f (3 , i ( z ) cos( x − y ) , for i = 1 , · · · , ,q (3) i ( z, x, y ) = f (3 , i ( z ) sin( x ) + f (3 , i ( z ) sin( y ) + f (3 , i ( z ) sin(3 x ) + f (3 , i ( z ) sin(3 y ) + f (3 , i ( z ) sin(2 x + y )+ f (3 , i ( z ) sin( x + 2 y ) + f (3 , i ( z ) sin(2 x − y ) + f (3 , i ( z ) sin( x − y ) , for i = 7 , , (35)For each value of γ , we now have 64 functions f (3 ,α ) i ( z ) to solve. As before, each Fourier mode α decouples from allothers, so we have 8 independent systems of ODEs to solve (one for each α ), for each γ . These ODEs again takethe form (26). The systems corresponding to α = 1 ,
2, and only these, have s (3 ,α ) i (cid:54) = 0, and hence depend on k (2) .Therefore, for each of these two α ’s we solve the equations of motion to find f (3 ,α ) i ( z ) and k (2) . These two systems ofequations are independent but k (2) is unique, so we must get the same k (2) in each of them.There is an important difference between lattices with γ (cid:54) = 1 / γ = 1 /
2. The leading contributionto the entropy difference ∆ σ H occurs at order O ( (cid:15) ) for the γ (cid:54) = 1 / O ( (cid:15) ) as in the γ = 1 / O ( (cid:15) ). There are now 13 Fourier modes that are excited. Schematically,each of these is of form cos( A (4 ,α ) ), with A (4 ,α ) ∈ { , x, y, x, y, x + y, x − y, x + y ) , x − y ) , x + 3 y, x − y, x + y, x − y } , (36)which are precisely those that follow from taking the fourth power of the sum of the two linear-order modes. Solvingthe associated system of ODEs we find the corresponding Fourier coefficients f (4 ,α ) i ( z ) for α = 0 , · · · ,
12, and k (3) = 0(for reasons similar to those that give k (1) =0).To find the wavevector correction k (4) , which is required to compute the thermodynamic quantities up to O ( (cid:15) ),we still need to find the metric solutions at order O ( (cid:15) ) since k (4) only appears in the equations of motion at thisorder (in practice it is enough to analyse the cos x and cos y Fourier modes since k (4) only appears in the equationsassociated to this sector). At this order, a total of 18 Fourier modes are excited. Schematically, each of these is ofform cos( A (5 ,α ) ), with A (5 ,α ) ∈ { x, y, x, y, x, y, x + 2 y, x − y, x + y, x − y, x + 4 y, x − y, x + y, x − y, x + 3 y, x − y, x + 2 y, x − y } , (37)which are precisely those that follow from taking the fifth power of the sum of the two linear-order modes. Solving thecorresponding ODE system we determine the associated Fourier coefficients f (5 ,α ) i ( z ) for α = 1 , · · · ,
18 and k (4) (cid:54) = 0.Having reached the first corrections to the relevant thermodynamic quantities, we do not continue the calculationto higher orders.1 FIG. 6: Convergence of the O ( (cid:15) ) coefficient of the entropy difference for γ = 1 /
2. The convergence is power-law, which isconsistent with the non-smooth decay at infinity.
Numerical Convergence
Here, we show the convergence of our numerical method. We use a Chebyshev spectral methods which convergencesexponentially on smooth functions. However, the expansion at infinity is not smooth. We instead find a power-lawconvergence with grid-size 1 /N . , which is consistent with this decay. A plot of this convergence is shown in Fig. 6. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected][1] R. Gregory and R. Laflamme, Phys. Rev. Lett. , 2837 (1993), hep-th/9301052.[2] O. J. C. Dias, P. Figueras, R. Monteiro, J. E. Santos, and R. Emparan, Phys. Rev. D80 , 111701 (2009), 0907.2248.[3] O. J. C. Dias, P. Figueras, R. Monteiro, H. S. Reall, and J. E. Santos, JHEP , 076 (2010), 1001.4527.[4] G. S. Hartnett and J. E. Santos, Phys. Rev. D88 , 041505 (2013), 1306.4318.[5] O. J. C. Dias, G. S. Hartnett, and J. E. Santos, Class. Quant. Grav. , 245011 (2014), 1402.7047.[6] J. E. Santos and B. Way, Phys. Rev. Lett. , 221101 (2015), 1503.00721.[7] R. Emparan and H. S. Reall, Phys. Rev. Lett. , 101101 (2002), hep-th/0110260.[8] O. J. C. Dias, J. E. Santos, and B. Way, JHEP , 045 (2014), 1402.6345.[9] L. Lehner and F. Pretorius, Phys. Rev. Lett. , 101102 (2010), 1006.5960.[10] P. Figueras, M. Kunesch, and S. Tunyasuvunakool, Phys. Rev. Lett. , 071102 (2016), 1512.04532.[11] P. Figueras, M. Kunesch, L. Lehner, and S. Tunyasuvunakool (2017), 1702.01755.[12] S. S. Gubser, Class. Quant. Grav. , 4825 (2002), hep-th/0110193.[13] T. Wiseman, Class. Quant. Grav. , 1137 (2003), hep-th/0209051.[14] H. Kudoh and T. Wiseman, Phys. Rev. Lett. , 161102 (2005), hep-th/0409111.[15] J. M. Maldacena, Int.J.Theor.Phys. , 1113 (1999), hep-th/9711200.[16] S. S. Gubser, I. R. Klebanov, and A. M. Polyakov, Phys. Lett. B428 , 105 (1998), hep-th/9802109.[17] E. Witten, Adv. Theor. Math. Phys. , 253 (1998), hep-th/9802150.[18] O. Aharony, S. S. Gubser, J. M. Maldacena, H. Ooguri, and Y. Oz, Phys. Rept. , 183 (2000), hep-th/9905111.[19] E. Sorkin, Phys. Rev. Lett. , 031601 (2004), hep-th/0402216.[20] P. Figueras, K. Murata, and H. S. Reall, JHEP , 071 (2012), 1209.1981.[21] M. Rozali and A. Vincart-Emard, JHEP , 166 (2016), 1607.01747.[22] A. Donos and J. P. Gauntlett, JHEP , 148 (2016), 1512.06861.[23] M. Headrick, S. Kitchen, and T. Wiseman, Class.Quant.Grav. , 035002 (2010), 0905.1822.[24] P. Figueras, J. Lucietti, and T. Wiseman, Class.Quant.Grav. , 215018 (2011), 1104.4489.[25] T. Wiseman, Numerical construction of static and stationary black holes (in Black Holes in Higher Dimensions, CUP 2012,Ed. Gary T. Horowitz, 2011), 1107.5513, URL http://inspirehep.net/record/920553/files/arXiv:1107.5513.pdf .[26] O. J. C. Dias, J. E. Santos, and B. Way (2015), 1510.02804.[27] T. Harmark and N. A. Obers, JHEP , 043 (2004), hep-th/0403103.[28] M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, M. Triana, and M. Zilhao, JHEP , 129 (2017), 1703.02948.[29] R. A. Janik, J. Jankowski, and H. Soltanpanahi (2017), 1704.05387.[30] M. Hanada and T. Nishioka, JHEP09