hep-th/0310096 PreprinttypesetinJHEPstyle-HYPERVERSION Caged Black Holes: Black Holes in Compactified Spacetimes II – 4 0 5d Numerical Implementation 0 2 n a J 9 1 Evgeny Sorkin, Barak Kol, Tsvi Piran 2 Racah Institute of Physics v Hebrew University 6 9 Jerusalem 91904, Israel 0 sorkin, barak kol, tsvi @phys.huji.ac.il 0 1 3 0 Abstract: We describe the first convergent numerical method to determine static black / h hole solutions (with S3 horizon) in 5d compactified spacetime. We obtain a family of t - solutions parametrized by the ratio of the black hole size and the size of the compact p e extra dimension. The solutions satisfy the demanding integrated first law. For small h black holes our solutions approach the 5d Schwarzschild solution and agree very well with : v new theoretical predictions for the small corrections to thermodynamics and geometry. i X The existence of such black holes is thus established. We report on thermodynamical r (temperature, entropy, mass and tension along the compact dimension) and geometrical a measurements. Most interestingly, for large masses (close to theGregory-Laflamme critical mass) the scheme destabilizes. We interpret this as evidence for an approach to a physical tachyonic instability. Using extrapolation we speculate that the system undergoes a first order phase transition. Contents 1. Introduction 1 2. Formulation 5 2.1 Choice of coordinates 5 2.2 Equations of motion 6 2.3 Boundary conditions and constraints 8 2.3.1 The z = 0 and z = L boundaries. 8 2.3.2 The r = 0 axis. 8 2.3.3 The horizon 9 2.3.4 The asymptotic boundary. 10 3. Thermodynamical and Geometrical Variables 11 4. Numerical Implementation 15 4.1 The scheme 15 4.1.1 Numerical lattice and discretization 15 4.1.2 Multigrid technique 18 4.1.3 Extracting measurables 19 4.1.4 Further developments 19 4.2 Testing the numerics 21 5. Results 29 6. Future directions 35 A. Equations of motion and boundary conditions on a d-cylinder 36 B. Asymptotic behavior. 39 1. Introduction Inbackgroundswithadditionalcompactdimensionstheremayexistseveralphasesofblack objects including black-holes and black-strings. The phase transition between these phases raises puzzles and touches fundamental issues such as topology change, uniqueness and Cosmic Censorship. Considerforconcreteness abackgroundwithasinglecompactdimension–Rd−2,1 S1. × We denote the coordinate along the compact dimension by z and the period by Lˆ. The – 1 – problemischaracterized byasingledimensionlessparameter1,e.g. thedimensionlessmass, µ =G M/Lˆd−3 whereG istheddimensionalNewtonconstantandM isthe(asymptotic) N N mass. Gregory and Laflamme (GL) [1, 2] discovered that a uniform black string – the d 1 − Schwarzschild solution times a line – becomes classically unstable below a certain critical value µ . They interpreted this instability as a decay of the string to a single localized GL black hole. Their discovery has initiated intensive research2 [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17] that attempted to trace out the fate of the unstable GL string: whether it settles at another intermediate stable phase as advocated in [3, 4], or whether it really decays to a single black hole. By now there is mounting direct evidence against the former possibility [8, 9], together with additional circumstantial evidence [5, 6]; and [15] which we also regard as evidence against stable non-uniform black string phase3. Here motivated by [8] we take another route, namely we address the question: what happens to a small localized black hole as its mass increases (by e.g. absorption of an interstellar dust)? Such a black hole grows and naively one expects that there is a moment when its “north” and “south” poles touch. Whether this is the case or not is yet to be established, but it is clear that some sort of instability will show up when the poles are getting closer. Put differently: is there a maximal mass, beyond which the black hole “does not fit into the circle” and there are no stable black holes. This maximal mass would be analogous to the GL critical mass, and would correspond to a perturbative, tachyonic instability . Yet another kind of instability may occur before that maximal mass is reached. Oncetheentropy of ablack hole equals theentropy of uniformblack stringwith the same mass, a transition between both phases will be allowed by quantum tunneling, or by thermal fluctuations. This first order phase transition is slower than the classical perturbative instability due to tunneling suppression. Noanalyticalsolutionforablack holeisknown. Moreover, eventhoughonecanexpect approximate analytic solutions to exist for very small black holes, the phase transition physics happens when the size of the black hole is comparable to the size of the compact dimension. Hence, in this work we take the numerical avenue. In our first article [23] we considered the theoretical background for the static d- dimensionalquasi-sphericalblackholes(BHs). Thereweoutlinedthegoalsofthenumerical study. Prime among these goals is to establish the very existence of the static black hole solutions. To our knowledge, there is no direct evidence in the literature that such BHs do exist4 though there are positive indications for that [14]. Among other goals is the study of such BH solutions in various regimes and dimensions. The ultimate and the most interesting aim is of course to determine the point of phase transition. 1Later we will use anotherparameter x definedin (2.4). 2Related research includes [18, 19, 20, 21, 22]. 3Notehoweverthattheauthorsof[15]didnotinterprettheirresultseitherassupportingorascountering the conjecture. 4Argumentslike: “inthelimitwhentheradiusofaBHissmallcomparedtothecompactificationradius theequivalenceprincipleimpliesthattheblackholemustbesimilartothe5dSchwarzschildsolution”,while intuitivearenotrigorously sufficient. Infact,thisargumentfails in4dwithoneofthespace-likedirections being curled to a circle, as there is no stable configuration of a periodic array of point-like sources. – 2 – Based on recent progress [9, 21] we develop a numerical scheme that allows us to find static axisymmetric BH solutions. Our scheme is dimension independent, provided that d > 4. As a first step, in this paper, we apply it to the 5d case, which is the example with the lowest dimension5 among spacetimes with extra dimensions. In 5d we construct numerically a family of static BHs, parametrized by x, which is the ratio of the size of the black hole to the size of the compact dimension, see (2.4). For small values of x the horizon region of our solutions approaches the 5d Schwarzschild solution which can be considered as the ’zeroth order in a perturbative expansion’ in powers of x. Moreover, in this limit our solutions satisfy the theoretical expectations for some next order corrections in this perturbative analysis [24] thereby allowing a confirmation of a new theoretical method through numerical “experiment”. This establishes the existence of static higher- dimensional BHs and shows that the Schwarzschild solution is the smooth limit of these solutions for x 0. → Wesucceedtocontroltheaccuracyofoursolutionsuptox . x 0.20(corresponding 1 ≃ to µ 0.047). Above this limit, up to the last value x 0.25 (corresponding to µ 1 2 2 ≃ ≃ ≃ 0.074), for which our solutions do not diverge the convergence rate was very slow and the numerical errors were not small. These values of µ should be compared with the critical GL mass µ 0.070. The slowdown of convergence and eventual divergence is mainly GL ≃ seen on one of our metric functions. By examining the equations of motion for our system we observe a ”wrong” sign in one of the equations (just like the plus sign in the following harmonic oscillator equation, ψ′′+ω2ψ =0), which is an indication of the presence of the o tachyon. Onecould expectthatthetachyonic behavior issuppressedforsmallxvalues and it is manifest for large x values, for which there are no static BH solutions6 7. However, the tachyonic behavior influences the numerics even before that critical x value and slows downtheconvergence. Webelievethattheproblematicvariableiscoupledtothetachyonic mode,andhencewhenthelatterdrivestheformertobehavepathologically, isanindication that the system is close to the phase transition point. In [23] we derived the d-dimensional Smarr formula, also known as the integrated first law, for the geometry under study (see also [16]). It is a relation between thermodynamic quantities at the horizon and those at infinity, relying on the generalized Stokes formula and the validity of the equations of motion in the interior. This naturally suggests to use this formula to estimate the “overall numerical error” in our numerical implementations. This method comes in addition to the standard numerical tests such as convergence rate, constraint violation etc. While it is possible that this is a lucky coincidence (though we believe it is not), for our solutions the Smarr formula is satisfied with 3 4% accuracy. − Moreover, it is intriguing enough that the Smarr formula is satisfied with the same 4% accuracy even for the problematic solutions for x& 0.20. This has to do with the fact that 5This is maybe the lowest dimensional example, but because of a very slow asymptotic decay, it is certainly not thesimplest to solve numerically [5]. We discuss this in detail later on. 6Consider a tachyon in a box: the mode can materialize only if its inverse mass is not less then the dimension of thebox otherwise the mode is suppressed. 7Thisisaclassical“revolutionarysituation”: a“poor”tachyonissuppresseduntiltheblackholebecomes toofat. Thenthetachyonrises,getsstronganddestroystheblackhole,headingtoanewfuture(toanother phase). – 3 – this formula relates only 3 variables of the 4 thermodynamical variables, characterizing the system. It turns out that the 4th variable is somewhat decoupled from the other three. However, as the inaccuracy in determining this 4th variable grows with x, this slows down and ultimately ruins the convergence. We believe that this is the variable which is coupled to the tachyonic mode. Even though the Smarr formula is satisfied to a good accuracy, we do have some larger inaccuracies in the solutions. One of the fields suffers from a certain convergence problems, andits asymptoticbehavior departsbysome30%fromsmallxpredictions. This is exactly the “fourth” asymptotic charge that does not appear in Smarr’s formula. Due to its approximate decoupling it is plausible that indeed we have good accuracy for all other measurements. Even better, we have indications that this field is reliable for a sub-range of x: 0.08 .x . 0.15. One can question what information can be extracted from knowledge of only three parameters for the entire sequence of BHs (x . 0.25), and knowledge of the fourth one for a smaller range (0.08 . x . 0.15). In particular, we show that the last black hole that we find (at x 0.25) deviates only slightly from being spherical, and moreover, its poles are ≃ quite distant from each other. In addition, one can ask whether there is a first order phase transition. We cannot establish this with certainty, since the entropy of our last black hole (at x 0.25) is still 2 ≃ larger than the entropy of the corresponding uniform black string. A naive extrapolation of our data to larger values of x indicates that the entropies will become equal just above the maximal BH that we find, namely at x 0.26 that corresponds to µ 0.082. It 3 3 ≃ ≃ is rather suggestive that µ ,µ and µ are all very close each to another. Since, all the GL 2 3 numbers in the system are expected to be of the same order, this fact may be regarded as an indication that we have found a real phase transition. Note that since µ µ we 2 GL ≃ come close to a first demonstration of a failure of higher dimensional uniqueness with two stable phases8. Finally we note, that generally in a first order phase transition one expects µ µ µ . This remains to be tested numerically. GL 3 2 ≤ ≤ While we expect that the instability we found corresponds to a physical one we stress that we cannot rule out the conservative possibility that it is a manifestation of imperfec- tions of the numerics. Since our numerical scheme is independent of the dimensionality of the problem provided d > 4, the immediate aim for the future work would be its applica- tion to higher dimensions, d 6, where the asymptotic fall off is faster and the solutions ≥ might be more stable9. In section 2 we describe our system. We employ the “conformal ansatz” and derive the equations of motion and the boundary conditions. A short excursion into theoretical background (summarized from [23]) is made in section 3. Our numerical implementation is described in detail in section 4, where we also describe various tests. The results are listed in section 5. We outline future directions in the final section 6. In appendix A, we derive 8Although we did not demonstrate that we assume that our BHs solutions are stable. 9In fact the preliminary results show that the picture that we find in 5d is qualitatively unchanged for d≥6. – 4 – the d-dimensional field equations and boundary conditions for the cylinder Rd−2,1 S1. In × appendix B we consider the asymptotic behavior of the equations. We also refer the reader to independent work by Kudoh and Wiseman who performed recently related calculations in 6d [25]. 2. Formulation In this section we focus on the five dimensional case – we derive the field equations and discuss the boundary conditions (b.c.). Equations and b.c. on a general d-cylinder are discussed in appendix A . The fifth spatial direction is denoted by z and it is compact with a period Lˆ, i.e. z and z+Lˆ are identified . We consider static localized BHs with an S3 horizon topology. We assume spherical symmetry (SO(3) isometry) of the 3 extended spatial dimensions and we denote the 4d radial coordinate by r. 2.1 Choice of coordinates We consider a static axisymmetric metric which is built out of three functions. We adopt a conformal (in the r,z plane) ansatz of the form { } ds2 = A2dt2+e2B(dr2+dz2)+e2Cr2dΩ2 , (2.1) − 2 where A,B and C are functions of r,z only and dΩ2 = dθ2+sin2θdφ2 2 To describe a BH it is convenient to transform to polar coordinates, defined by r = ρsinχ, z = ρcosχ, (2.2) since the BH horizon is represented by a closed curve in the r,z plane. The metric in { } these coordinates reads now ds2 = A2dt2+e2B(dρ2+ρ2dχ2)+e2Cρ2sin2χdΩ2 , (2.3) − 2 To simplify the numerical procedure it is desirable that the boundaries of the inte- gration domain10 lie along the coordinate lines. Note that by choosing the ansatz (2.1), or (2.3) we still did not fix the gauge completely. There is still a freedom to move the boundaries of the integration domain by a conformal transformation. It was shown in [21] that using this conformal freedom the horizon boundary could be set at a constant radius ρ , leaving the periodic boundaries along z = const lines. Thus the domain is h (r,z) : z L, r2+z2 ρ 2 , where for future use we define the half-period L = Lˆ/2 of { | | ≤ ≥ h } the compact circle, see Fig. (1). In addition, by fixing the ratio of the radius of the horizon to the period of the circle ρ h x:= , (2.4) L all residual gauge freedom is eliminated. In our implementation, we set ρ = 1, without a h loss of generality, and generate different solutions by varying L. 10Whatwecallhere“domainofintegration”couldbecalledalternatively“domainofdefinition”,“domain of relaxation” etc. By this term we refer to theregion of space-time where we solve our equations. – 5 – =0xis =0 xis ra χ a z=L z= L Conformal χ=π/2 Mapping z=−L z=−L black hole black hole (a) (b) Figure 1: A spacelike slice of the black-hole spacetime. (a) In the r,z plane the black hole’s { } horizon is a curve with a spherical S3 topology. (b) There is a conformal freedom to transform the domain to (r,z): z L, r2+z2 ρh2 . By fixing ρh/L the domain is uniquely specified [21]. { | |≤ ≥ } The fact that x cannot be changed freely for a given solution implies that x is a characteristic parameter analogous to the (normalized) total mass or the temperatureeven though it does not have a clear physical meaning. For example, if one enforces the horizon of a BH to be at a fixed radius set by some given x, it would be excessive to specify also the temperature. Conversely, specifying the temperature one does not have the freedom to constrain the location of the horizon [21]. Inpolarcoordinatesthereflectingboundaryofthecompactcircle,z = 0,isatχ = π/2, but the periodic boundary, z = L, does not lie along a coordinate line in the ρ,χ plane. { } Thetreatmentofthisirregularboundaryintroducesacertaincomplicationinthenumerical scheme as described in section (4.1). Nevertheless, we believe that it is preferable to work in polar coordinates (2.2) and to have an irregular boundaryat z = L, rather than work in rectangular coordinates r,z and have an irregular boundary at the horizon. Intuitively, { } this is because we expect that the region near the horizon would become the region of the ’activity’ as x increases. For numerical reasons it would be convenient to use another angular coordinate ξ =cos(χ) . (2.5) The benefit of using this coordinate is twofold. First, the irregular z = L boundary has a particularly simple representation, ρ= L/ξ. Second, as we explain shortly, the coordinate singularity at the axis, r = 0, becomes first order instead of second order. 2.2 Equations of motion Our basic equations are the five-dimensional time independentvacuum Einstein equations. There are five equations in 5d: two are equations of motion for A,C, while variation with respect to the metric in the (r,z) plane yields three additional equations. In the conformal ansatz oneofthemisanequation of motionforB whiletheother tworesultfromthegauge – 6 – fixing. The equations can be combined in a way that three of them will take the form of elliptic equations, which we call the interior equations. The other two combinations that contain ahyperbolicdifferential operator willbetermed ’theconstraints’. Theseconstraint equations are not independent as they are related to the interior equations via the Bianchi identities. In order to obtain the interior equations we can follow the general proceduredescribed in appendix A, or alternatively, use a suitable symbolic math application e.g. GRTensor [28] to evaluate the relevant quantities. In either route one obtains the interior equations which are the following combinations of the components of the Einstein tensor: θ + Gθ 1/2 χ+1/2 ρ 2 t ,2 θ 2 χ 2 ρ+ t and θ + χ+ ρ t. They can be written Gχ Gρ − Gt Gθ − Gχ − Gρ Gt Gθ Gχ Gρ −Gt respectively as 2∂ A 1 A + ξ ξ+(1 ξ2)∂ C +2∂ A +∂ C = 0, (2.6) △ ρ2 − − ξ ρ ρ ρ (cid:18) (cid:19) 2∂ A(cid:0)ξ 1 ξ2 ∂ C(cid:1) 2∂ C ξ 1 1 ξ2 ∂ C B + ξ − − ξ + ξ − 2 − ξ △ Aρ2 ρ2 − (cid:0) (cid:0) (cid:1) (cid:1) (cid:0) (cid:0) (cid:1) (cid:1) 2∂ C 2∂ A 1 1 e2B−2C ρ (∂ C)2 ρ +∂ C − = 0 , (2.7) − ρ − ρ − A ρ ρ − ρ2 (1 ξ2) (cid:18) (cid:19) − ∂ A ξ 1 ξ2 ∂ C 4∂ C ξ 1 1 ξ2 ∂ C C ξ − − ξ ξ − 2 − ξ + △ − Aρ2 − ρ2 (cid:0) (cid:0) (cid:1) (cid:1) (cid:0) (cid:0) (cid:1) (cid:1) 4∂ C ∂ A 1 1 e2(B−C) + ρ +2 (∂ C)2+ ρ +∂ C + − = 0 . (2.8) ρ ρ A ρ ρ ρ2 (1 ξ2) (cid:18) (cid:19) − Here we used the variable ξ instead of χ and the Laplacian becomes ∂2 +(1/ρ)∂ + △ ≡ ρ ρ (1/ρ2) 1 ξ2∂ 1 ξ2∂ . ξ ξ − − Thpe constrain(cid:16)tpequations(cid:17)expand to ∂ A 1 +∂ B ξ 1−ξ2 ξ ρ ρ 2ξ (∂ρB −∂ρC) + 2∂ξB +2∂ C∂ B Gρ ≡ ρ2 (cid:16) A (cid:17) − 1 ξ2 ρ ξ ρ p n − ∂ B∂ A ∂ A ξ ρ ρξ + 2∂ B∂ C 2∂ C∂ C + − 2∂ C = 0 (2.9) ξ ρ ξ ρ ρξ − A − ∂ A ξ 2 1 ξ2 ∂ B 4ξ (∂ B ∂ C) o ρ ξ ξ − − ξ + ξ − ξ Gρ −Gξ ≡ − Aρ2 ρ2 − (cid:0) (cid:0) (cid:1) (cid:1) 2 1 ξ2 (2∂ B ∂ C) ∂ C 1 ξ2 ∂2A − ξ − ξ ξ + − ξ + − ρ2 Aρ2 (cid:0) (cid:1) (cid:0) (cid:1) 2 ξ∂ C + 1 ξ2 ∂2C ∂ A 1 +2∂ B − ξ − ξ ρ ρ ρ + + + (cid:16) ρ(cid:0)2 (cid:1) (cid:17) (cid:16) A (cid:17) 1 ∂2A + 2 (2∂ B ∂ C) +∂ C ρ 2∂2C = 0 . (2.10) ρ − ρ ρ ρ − A − ρ (cid:18) (cid:19) β Assuming that the interior equations are satisfied, the Bianchi identities = 0, imply α;β G [9] the following relations between the constraint equations ∂ +∂ = 0 , ζ ξ U V – 7 – ∂ ∂ = 0 , (2.11) ζ ξ V − U ρ ξ where ζ = logρ., and we define the rescaled constraints = ρ√ g /2 , = U − Gρ −Gξ V ρ2√ g ξ with g det g . (cid:16) (cid:17) ρ αβ − G ≡ A nice feature follows[9]. The constraints and satisfy the Cauchy-Riemann (2.11) U V relations and hence each one of them is a solution of the Laplace equation. Hence, if one of the constraints is satisfied at all boundaries and the other at a single point along some boundary these constraints must be satisfied everywhere inside the domain. This fact will be referred hereafter as the “constraint rule”. In our implementation, following the choice in [9] we imposed = 0 along all boundaries and = 0 in the asymptotic region. It is V U important to check and confirm that the constraint and are satisfied everywhere for U V our numerical solutions, as we describe in section 4.2. 2.3 Boundary conditions and constraints The interior elliptic equations (2.6-2.8) are subject to boundary conditions. In this section we describe the boundary conditions that define the problem completely. The integration domain is defined by (r,z) : 0 z L, r2 +z2 ρ 2 , designated by the thick dashed { ≤ ≤ ≥ h } line in Fig. 1. The boundary conditions are specified on the axis, at the horizon, in the asymptotic region and at the reflecting and periodic boundaries z = L and z = 0. 2.3.1 The z = 0 and z =L boundaries. On the reflecting, z =0, and the periodic, z = L, boundaries we impose ∂ ψ = 0, ψ = A,B,C. (2.12) z While at the reflecting boundary this condition is simply ∂ = 0, at the periodic boundary ξ its implementation is not direct, see section 4.1.1. 2.3.2 The r = 0 axis. Regularity of the metric on the axis (absence of a conical singularity) requires B = C. (2.13) We use this equation as a Dirichlet condition for B. Equation (2.7) is not solved at the axis but it is only monitored there. For A and C the boundary conditions are automatic – on axis the (interior) equations for these functions become first order in derivatives normal to the boundary and have precisely the form of a b.c. Namely these equations are already incorporate b.c. and these do not need to be additionally specified. We term this an “automatic boundary condition”. This occurs because of our particular choice of the angular coordinate: we use ξ instead of χ. The axial symmetry of the problem dictates the ∂ = 0 condition for the metric functions, which translates to ∂ = 0 in r χ spherical coordinates and 1 ξ2∂ = 0 in our coordinates. But on axis ξ = 1 and hence ξ − this condition need not be imposed in our coordinates. While the coordinate singularity p at the axis is quadratic, sin(χ)−2 when using χ, it becomes linear (1 ξ)−1 when ∼ ∼ − – 8 – using ξ. With this advantage there is, however, a drawback: the metric functions are not differentiable at ξ = 1. This requires a modification of the numerical scheme there, replacing the second order normal derivative of the interior equations by a first order one due to the considerations above as described in subsection (4.1). 2.3.3 The horizon The horizon in our construction is located at ρ = 1. For static solutions various notions h of the horizon coincide – the event horizon (globally marginally trapped ), the apparent horizon (the outermost boundary of locally trapped surfaces) and the Killing horizon are all the same. The latter characterizes the horizon as a surface where A 0. (2.14) ≡ This implies that along the horizon ∂ A = ∂2A = 0. (2.15) ξ ξ Even though the horizon normally is not singular (in curvatures) our equations do become singular there as the function A vanishes. Now we describe how the physical regularity of the equations at the horizon gives boundary conditions for our functions11. Expanding Eqs. (2.7,2.8) at the horizon we obtain the condition ∂ C = 1 . (2.16) ρ − We still need a condition for B. We obtain this condition from the zeroth law of the black- hole mechanics (or thermodynamics), namely that for static solutions the surface gravity must be constant along the horizon (see for example [26]). The surface gravity along the horizon reads κ= e−B∂ A, (2.17) ρ and the derivative of κ along the horizon vanishes ∂ A ρξ ∂ κ ∂ B = 0, at ρ= ρ . (2.18) ξ ξ h ∼ − ∂ A ρ The upshot is that the boundary condition for B can be obtained in one of the forms: either ∂ A ρ B = C +log , (2.19) ξ=1 ∂ A |ρh (cid:18) ρ ξ=1(cid:19) from (2.17) and (2.13), or by integrating (2.18) outwards from the axis along the horizon. In our implementation we used the former form. However we have checked that a corre- sponding solution obtained by using the other option differs only slightly from our original one. Note that the condition (2.18) implies that eqn. (2.9) (or = 0) is guaranteed along V the horizon, and vice versa. 11Weassume hereafter that ∂ A| 6=0, i.e. the horizon is not degenerate. ρ ρh – 9 –