Bumpy black holes

We study six-dimensional rotating black holes with bumpy horizons: these are topologically spherical, but the sizes of symmetric cycles on the horizon vary non-monotonically with the polar angle. We construct them numerically for the first three bumpy families, and follow them in solution space until they approach critical solutions with localized singularities on the horizon. We find strong evidence of the conical structures that have been conjectured to mediate the transitions to black rings, to black Saturns, and to a novel class of bumpy black rings. For a different, recently identified class of bumpy black holes, we find evidence that this family ends in solutions with a localized singularity that exhibits apparently universal properties, and which does not seem to allow for transitions to any known class of black holes.


Introduction and main results
In spite of the lack of effective solution-generating methods, the exploration of black hole solutions of the vacuum Einstein equations in D ≥ 6 has made significant strides through the complementary use of approximate analytical methods and numerical calculations. One line of study follows the observation that rapidly spinning Myers-Perry (MP) black holes in D ≥ 6 [1] approach black membranes, so they can be expected to admit, as black branes do, stationary deformations that ripple the horizon [2]. Such bumpy black hole solutions would naturally connect to black rings, black Saturns, and multi-ring solutions through topology-changing transitions in solution space [3,4]. Evidence for this picture has been provided in [5,6,7,8,9]. In this article we confirm, refine, and extend aspects of it through a detailed numerical investigation of bumpy black holes in D = 6.
Bumpy black holes, like MP black holes (in the singly-spinning case that will be the focus of this article), have horizon topology S D−2 with spatial symmetry group U (1) × SO(D − 3). What distinguishes them from the "smooth" MP black holes is that the size of the S D−4 symmetry orbits on their horizon varies in a non-monotonic fashion from the axis of rotation to the equator.
The different families of bumpy black holes are conveniently identified by the way they branch off the MP family. Refs. [5,6] identified linearized zero-mode perturbations of singly-rotating MP black holes, to which we can assign an 'overtone' number i = 1, 2, . . . 1 that, for fixed mass, grows with the spin. Let us conventionally fix the sign of the zeromode so that the i-th mode wavefunction at the axis of rotation has sign (−1) i . By adding or subtracting the zero-mode perturbation to the MP black hole, we obtain two different branches, (+) i and (−) i , of solutions emerging from the branching point. The evolution of the solutions along the (+)-branches was anticipated in [3]: the horizon develops bumps that grow until a S D−4 cycle pinches down to zero size, naturally suggesting a topologychanging transition to other black holes: black rings for i = 1, black Saturns for i = 2, and multi-ring configurations for i > 2. The (−)-branches of solutions were only identified recently in [8], and as we will see, they seem to terminate without any plausible connection to other black hole solutions.
Ref. [8] has studied the (+) 1 and (−) 1 branches in six and seven dimensions. Here we extend the analysis in the six-dimensional case to higher branches, i = 1, 2, 3, while pushing both the (+) and (−) branches closer to their ends in phase space at singular solutions. We also perform a detailed investigation of the geometrical properties of these bumpy black holes. Our main conclusions, partly illustrated in figs. 1 and 2, are:    for the largest deformations we have obtained. The value of R ⊥ gives the size of the spheres S 2 transverse to the rotation plane. The vertical axis u = 0 is the rotation axis, but u does not measure the radius of the rotation circles. We superimpose the embeddings of MP black holes with the same mass and spin (dashed black), and of the cones proposed for a local model of the critical singularity, eq. (3.9) (blue). The angular momentum j is normalized as in (2.19). In this and all subsequent plots, units are GM = 1. We give a simple argument to suggest that, as we move away of the MP solutions, the horizons in higher-i branches pinch-off in succession from the rotation axis to the equator. Let us also remark that the divergent length of the equatorial circle mentioned in point 3 is not visible in fig. 2, but will be made apparent in fig. 8 below.
These results are explained in detail in sec. 3, after having outlined in sec. 2 the construction of the solutions. In addition, in sec. 4 we compute the thermodynamic properties of these solutions and draw phase diagrams. We also analyze the spectrum of the Lichnerowicz operator, and relate the number of negative eigenvalues to the thermodynamic stability of the solutions. The details of our numerics are relegated to appendix A.
We remark that all these bumpy black holes are expected to be dynamically unstable; their importance lies in what they reveal about the possible geometries of black hole horizons in higher dimensions and the rich web of interconnections among them.

Construction of the Solutions
In order to construct deformed rotating black holes in six dimensions we solve the Einstein- Γ is the usual Levi-Civita connection compatible with the spacetime metric g andΓ is the Levi-Civita connection compatible with some reference metricḡ that satisfies the same boundary conditions as the spacetime metric g but needs not be a solution to Einstein's equations. This is a standard method used in numerical General Relativity to find static and stationary solutions [12,13,14]: the equations are manifestly elliptic and one can then use conventional numerical techniques for solving such partial differential equations. For asymptotically flat (AdS or Kaluza-Klein) static metrics [13] proved that the solutions to the Einstein-DeTurck equations must in fact be Einstein. This result is yet to be extended to the stationary case, and hence, since we are interested in Einstein metrics, we must check that the DeTurck vector ξ vanishes. For the solutions presented in this article we have checked that this is indeed the case, within our numerical accuracy. The solutions we study are stationary and with only one of the two possible rotations turned on. Thus the rotation group SO(5) is broken down to U (1) × SO(3), which act on the direction of rotation φ and on the spheres S 2 transverse to the rotation plane. The metric can then be written in the form and we denote the reference metric as The compact radial direction r ∈ [0, 1) covers the region from the horizon at r = 0, to infinity at r = 1. We seek solutions with horizons that are topologically S 4 , so we choose sections at constant t and r to also be topological S 4 's. The size of the φ-circles and of the symmetric S 2 's varies along the polar angular direction x ∈ [0, 1], with x = 0 corresponding to the rotation axis (where φ-circles shrink to zero) and x = 1 to the equatorial plane (where S 2 spheres shrink to zero). The "bumpiness" of the horizon corresponds to non-monotonicity (of, say, the size of the S 2 on the horizon) along this polar direction. We write the metric functions as and the reference metric is the MP metric with a small modification that enables us to control the temperature of the solutions. In order to obtain the functions (T 0 , W 0 , ... etc.) for the reference metric we begin with the single-spin MP metric in standard Boyer-Lindquist-like coordinates (herer,φ,t and θ) ds 2 = − dt 2 + r 3 0 rρ 2 dt + a sin 2 θdφ 2 + (r 2 + a 2 ) sin 2 θdφ 2 + ρ 2 ∆ dr 2 + ρ 2 dθ 2 +r 2 cos 2 θdΩ 2 (2) (2.5) where ρ 2 =r 2 + a 2 cos 2 θ, ∆ =r 2 + a 2 − r 3 0 r (2.6) and the dΩ 2 (2) is the line element of a 2-sphere. The horizon (r = r + ) is found by solving the mass and angular momentum are where Ω (4) is the area of a unit 4-sphere, and the temperature and angular velocity are (2.9) We perform the changes of coordinates (2.10) The first two are made so that the ranges of the coordinates are 0 < r, x < 1 and the third change, to co-rotating coordinates, is made because otherwise the W 0 function goes to zero too fast at infinity, which is inconvenient for numerical calculation. In co-rotating coordinates the function W 0 is 0 at the horizon and Ω H asymptotically. The MP metric then takes the form We will find bumpy black hole solutions with given values of the temperature and angular velocity. It is convenient to specify these in terms of parameters of the reference metric. In order to control the temperature, we introduce a parameter k in the reference metric so that the surface gravity κ of the reference metric is given by Obviously, whenever k = 1 the reference metric is not a solution of Einstein's equations but nonetheless it has a smooth horizon. However, k allows us to move along the branches of solutions by varying it as a parameter in the reference metric. We will choose boundary conditions on the Q's at the horizon in such a way that the surface gravity of the bumpy black holes is also given by (2.15). Note that by modifying k not only the temperature but also the mass and angular momentum of the solutions change. However, with the appropriate boundary conditions, Ω H remains unchanged.

Boundary conditions
The conditions we impose on the Q's at each of the boundaries of our domain in order to get regular solutions are Horizon (r = 0): The reference metric is already regular on the horizon. Since the spacetime metric is the reference metric multiplied by the Q's, we ensure regularity on the horizon by imposing Neumann boundary conditions ∂ r Q| r=0 = 0. In addition we impose , which fixes the surface gravity to the value (2.15).
Axis (x = 0): The reference metric is already regular on the axis of rotation so again we impose Neumann boundary conditions ∂ x Q| x=0 = 0. The φ circle goes to zero at this boundary and in order to avoid a conical singularity we impose Q 3 (r, 0) = Q 6 (r, 0).
Equator (x = 1): The boundary conditions are again Neumann ∂ x Q| x=1 = 0. Since here the radius of the S 2 shrinks to zero size, we impose Q 4 (r, 1) = Q 6 (r, 1) to avoid a conical singularity.
Infinity (r = 1): For asymptotically flat (AF) solutions, since the reference metric is already AF, we impose the Dirichlet boundary conditions so that the asymptotics are unchanged by the Q's. Since we are in co-rotating coordinates the horizon angular velocity relative to infinity is given by the asymptotic value of W (r, x). Then the condition Q 2 (1, x) = 1 ensures that Ω H is given by the same expression as in the MP black hole.

Physical magnitudes
Given our choice of boundary conditions, the temperature and the angular velocity at the horizon are easily extracted in terms of quantities present in the reference metric, namely, r 0 , k and r + , so that (2.17) Since we work with vacuum solutions, we can obtain the mass and angular momentum by evaluating their Komar integrals at the horizon, where χ is the 1-form dual to the asymptotic time-translation Killing vector ∂ t − Ω H ∂ φ , and ζ is dual to the axial Killing vector ∂ φ . In addition to these quantities we also compute the area of the horizon. In order to compare different solutions that have the same mass we use the dimensionless quantities with the numerical factors c a , c j , c ω , c t chosen as in [3].
Other geometric invariant quantities that are of interest for characterizing the solutions are the radii on the horizon (r = 0) of the circles parallel to the plane of rotation, R (x), and of the spheres S 2 orthogonal to it, R ⊥ (x). They are given by We render these dimensionless by dividing them by (GM ) 1/(D−3) without any additional factors.
We will often use j as the 'control parameter' that changes along a branch of solutions. The bumpy branches extend over rather narrow ranges of j. They originate at bifurcation points in the MP family given respectively by (±) 1,2,3 beginning : j = 1.20, 1.41, 1.57. (2.21) The (+)-branches initially extend towards larger values of j, but then bend backwards towards decreasing j, which we have followed down to

Geometry of bumpy black holes
In this section we explore the geometry of the solutions, in particular of their horizons. Since we have pushed the new branches close to their endpoints in solution space, one purpose is to examine whether the critical solutions of (+) branches have singularities modeled by Ricci-flat double-cone geometries that can mediate the transitions to black ring, black Saturn, and multi-ring solutions [10]. Another aim is to get a better understanding of the solutions in the (−) branches, in particular where and how these branches end.
The spatial horizon geometry, r = 0, t = constant, is In order to gain some intuitive understanding of these geometries, we perform two kinds of plots: embedding diagrams of sections of the horizon into Euclidean space, and plots of the invariant radii of the S 1 and S 2 symmetry cycles.
Embedding diagrams. Embeddings in Euclidean space provide useful and intuitive visualizations of the geometry. Here we use the same type of embeddings as ref. [15] presented for black rings. On the spatial horizon geometry we choose a section φ = const., and embed it in E 4 as a surface of the form so the induced geometry is The embedding is found by integrating which is possible since B(0, x) ≥ R ⊥ (x) 2 for all our solutions. In our plots we present The coordinate u does not have any invariant meaning as the radius of the rotational S 1 's, since this representation misses the information about R (x). For this, we employ a different type of plot.    Invariant-radii plots. These are plots of R ⊥ (x) versus R (x). Information about the length in the polar direction is lost now, which makes the horizon shapes in these plots look somewhat peculiar.

(+)-branch bumpy black holes
Representative solutions of these branches are depicted in embedding diagrams in fig. 3 and in invariant-radii diagrams in fig. 4. Observe that, contrary to what may seem from the embedding diagrams, the radius R of the S 1 near the equator is larger in MP black holes than in the bumpy solutions with the same mass and angular momentum.
Near the values (2.22) the solutions clearly approach configurations where a symmetric S 2 on the horizon pinches down to zero size, developing a singularity whose structure we analyze next.

Critical cone geometries
Depending on whether the singular pinch-off occurs along the rotation axis or on a circle away from the axis, the geometries are expected to be locally Lorentzian double-cones of with horizons at χ = π/2, and where L is the radius of the circle where the S 2 pinch to zero [10]. If we embed the section t = const, χ = π/2 of these geometries in Euclidean space as above, then it is easy to see that they are represented as the cones on-axis: We can superimpose these on the most deformed solutions we have obtained in these branches. Fig. 1 shows excellent agreement with the prediction of [10]. The invariant-radii plots probe complementary geometric aspects of the horizon. The geometries (3.7), (3.8) have slopes on-axis: which are also very well reproduced on-axis for (+) 1,3 , see fig. 5, but less well so off-axis for (+) 2 , reflecting (maybe unsurprisingly) a remaining small dependence of R on the polar angle that would become negligible only much closer to the critical singularity.    We also compare the Kretschmann scalar of both geometries 3 . For the cones, K depends only on the 'polar' coordinate z while for the black holes it depends not only on x but also on r. In order to make the comparison we must specify a way to map points between the two geometries, i.e., a function z(r, x). This involves a certain arbitrariness, which we fix by equating the radius of the 2-sphere in both geometries. Then (in six dimensions) on-axis: S(r, x) = z 2 4 , off-axis: S(r, x) = z 2 3 (3.12) and so K on-axis cone = 72 (3.13) These comparisons are dominated by how the size of the S 2 shrink close to the singularity, including away from the horizon, but they do not test the length of the equatorial S 1 , to which K is largely insensitive.
We have computed the discrepancy between the Kretschmann scalars of both geometries, K bh Kcone − 1 , for the three branches and it is less than 10% (often less than 5%) in the region near the singularity. Therefore, we conclude that the critical cones are locally a good description of the singular region.
Finally, we have also checked the appearance of a conical structure in the Euclidean time direction. Fig. 6 shows the rate at which the Euclidean time circle shrinks along the axis of rotation in our nearest-to-critical (+) 1 solution. 4 The slope in this curve fits well the slope of the conical solution over a range of distances close to the black hole. It departs from it very near the horizon, as it must since the cone is smoothed in our solution 5 .

(+) 3 : transition to bumpy black rings
It was naturally conjectured in ref. [3] that black holes along the (+) 1,2 branches would pinch to zero and transition to black ring and black Saturn phases, respectively. However, higher branches (+) i≥2 have multiple pinches and it was less clear what their fate could be. If pinch-down occurred first on a circle off-axis, then the branch (+) 3 would transition to a black Saturn configuration with a bumpy central black hole. However, the deformation of (+) 3 black holes is expected to be larger on-axis than off-axis. The reason is that in the black membrane limit of the MP black holes, and for small, linearized perturbations, the axisymmetric Gregory-Laflamme-type perturbation takes the form [2] δg µν ∼ J 0 (x)h µν (r), (3.14)  where x is the distance from the rotation axis in directions parallel to the horizon, and hence plays the role of the polar angle. The Bessel function J 0 (x) yields larger deformations close to the axis of rotation at x = 0, and decays away from it. Figs. 1 and 3 show that this behavior persists when the deformations are not small. This evolution of the (+) 3 branch has a natural end at a topology-changing transition to a branch of bumpy black rings, of horizon topology S 1 × S D−3 , with a deformed S D−3 . These have not been constructed yet, and our arguments are the first clear indication of their existence. It is also natural to expect that the bumpy black ring branch will connect, at its other end, to black di-rings. Indeed it seems implausible that they smooth out their deformations and connect to the known (smooth) black ring solutions, since these two branches are very far apart in solution space (see fig. 10 below).
The argument also suggests that the same behavior occurs in higher branches, with pinches being larger closer to the axis, and pinching-down sequentially at increasing values of the polar angle x. The conifold-type transition then connects them to new families of multiply-bumpy black rings, which eventually, through several transitions, connect to multi-ring configurations.

(−)-branch bumpy black holes
Figs. 7 and 8 show the previous two types of graphics for the horizon geometry of these black holes at their largest deformation, (2.23) (these were also shown in fig. 2).
From fig. 8 we see that these horizons spread in the rotation plane more than in the MP black holes of the same mass and spin. This could be anticipated near the bifurcation point, where the deformation is controlled by a zero-mode with i + 1 nodes: since the (−) i zero mode wavefunctions have sign (−1) i+1 at the rotation axis, then the wavefunction at the      fig. 7. R is the radius of circles parallel to the rotation plane and R ⊥ is the radius of the orthogonal S 2 . The black dashed curve shows a MP black hole of the same mass and angular momentum.  equator must always be positive, i.e., the bumpy black hole bulges out. 6 At least for the i = 1 solutions, we can also understand this in more physical terms: close to the branching point both solutions have the same mass, angular momentum and angular velocity. If the MP black hole is perturbed in such a way that some of its mass is concentrated closer to the axis of rotation, then in order to maintain the angular momentum constant (with the same angular velocity) some mass must also be moved farther along the rotation plane, preferrably around the equator.
Further along the branch the horizons stretch a lot on the rotation plane, see fig. 8, and get highly pancaked, R ⊥ R . Nevertheless, in contrast to ultraspinning MP black holes, they do not seem to approach black membranes in the limit, and in particular (as we will see in sec. 4) they do not develop the Euclidean negative modes that would signal the Gregory-Laflamme-type instability characteristic of black membranes. Fig. 9 strongly suggests that the length of the equatorial circle diverges in the limiting solutions -even though the radial distance to the equator remains finite. This behavior is known to occur for the extremal limit of the five-dimensional MP black hole, although in the latter case the extremal solution has zero temperature and area, whereas these remain finite in the critical (−) solutions.
The S 2 's on the equator shrink to zero in the limiting solutions in a singular way, causing the Kretschmann scalar to diverge. The effect seems to be the same in all three branches, being well reproduced on sections of constant t and φ on the horizon (such as are captured in fig. 7, and by the Kretschmann scalar) by the geometry which is also present in off-axis cones (3.8). This suggests that the local structure of the singularity at the equator in these solutions may be universal for all (−) branches: the S 2 shrink to zero along the horizon like in (3.15), while the length of the equatorial S 1 diverge.
Although we do not have a local model for the full singularity, it is not one of the conical geometries that effect a transition to another branch of black holes. In fact it seems unlikely that the singularity is a Ricci-flat scaling geometry. In view of this, and in the absence of a plausible candidate for a merger transition, we are led to conjecture that the (−) branches of black holes terminate in phase space without continuing into any other singly-spinning stationary black hole solutions.

Phase diagrams, thermodynamic stability, and negative modes
In fig. 10 we show the area, temperature and angular velocity as a function of the angular momentum for fixed mass. We can see the two different families of solutions branching off from each of the perturbative zero modes. Since in fig. 11 it is difficult to distinguish the two branches in the (j, a h ) plane, we also show plots of the area difference between the bumpy black holes and the MP solutions for values of j close to each branching point. The black ring phases obtained in [8] are also included in these plots, and it is apparent that the (+) 1 solutions tend to a merger point with the black rings. Although our results suggest that solution-trajectories inspiral close to this transition (which would lead to infinite discrete non-uniqueness of the kind found in [17,18,19]), our accuracy in this region is not enough to reach a definite conclusion. Thermodynamic stability of these black holes in the grand-canonical ensemble is obtained when the specific heat at constant angular momentum C j and the isothermal moment of inertia are both positive [20] Negative moments of inertia are possible for black holes since they are not rigid bodies. They can reduce their angular velocity while gaining angular momentum by spreading in the rotation plane. This is precisely what happens in ultra-spinning MP black holes. In this case it impossible for the black hole to remain in equilibrium with a co-rotating radiation reservoir. The specific heat and moment of inertia can be read off from the slopes of the solution curves in the (T , M ) and (Ω H , J) planes. The details of the plots for actual solutions are difficult to distinguish, so instead in fig. 12 we present sketches of them that capture their qualitative features.
In addition, we have also studied the spectrum of the Lichnerowicz operator, since its negative eigenvalues are directly related to the negative modes of the quasi-Euclidean     fig. 10). We only show the branches (±) 1 , but (±) 2,3 have the same behavior. The branch (+) 1 (brown) has two negative modes from O 1 to C 1 and three from C 1 onwards, while the branch (−) 1 (black) always has three negative modes. Color coding as in fig. 10. action. We have checked that the number of negative eigenvalues coincides with the expectations from thermodynamic stability. In particular, along the MP family of solutions in the direction of increasing j, initially the solutions have one negative mode that corresponds to negative specific heat (it is the MP extension of the Euclidean Schwarzschild negative mode), and acquire a second one at the cusp in the (Ω H , J) plane where the moment of inertia first becomes negative. This is also the minimum of the temperature (see fig. 10) which signals the entrance into the ultraspinning regime, and which coincides with the change of sign of . At higher j one encounters further zero modes that become negative ones. These are not associated to new thermodynamic instabilities, instead they are 'overtones' of Gregory-Laflamme-type negative modes.
The thermodynamic stability and negative modes along (+) branches are more complicated, as there are several points where the susceptibilities (4.1) change sign. Here we explain it for (+) 1 solutions (higher (+) branches exhibit the same qualitative behavior), referring to fig. 12: From O 1 to A 1 : Point O 1 is the bifurcation from the MP branch of solutions. The new branch bifurcates with higher area, hence the MP solution is expected to be less stable, and indeed it acquires an extra negative mode, while the bumpy solution keeps the number of negative modes present in the MP solutions just before O 1 . Both the specific heat and the moment of inertia are negative in this segment; accordingly, the Lichnerowicz operator on solutions from O 1 to A 1 has two negative eigenvalues.
From A 1 to B 1 : The point A 1 at which C j changes sign from negative to positive passing through zero corresponds to the cusp in the (j, a h ) plane in fig. 10, where the branch beyond A 1 has lower area and C j remains positive until B 1 . But there is no qualitative change in the spectrum of the Lichnerowicz operator at A 1 . We interpret the two negative eigenvalues present here as due to the negative and to the fact that there exists another solution with higher area for the same j. Observe that a given negative mode does not strictly correspond to just one instability.
From B 1 to C 1 : At B 1 the sign of C j changes from positive to negative and the sign of from negative to positive. Like before, the number of negative modes is preserved and the spectrum of the Lichnerowicz operator does not signal these changes in thermodynamic susceptibilities.
From C 1 : At C 1 the sign of changes from positive to negative going through infinity. The Lichnerowicz operator acquires a third negative mode.
We see that the (+)-branch solutions are always thermodynamically unstable, since either or C J or both are negative. The solutions are likely dynamically unstable to bar-mode perturbations, like MP black holes are at even lower values of j.
Regarding the (−) branches, they all have negative C j and . In addition they come out of the bifurcation with less area than the MP black holes. As expected from the arguments above, the Lichnerowicz operator on these solution has three negative eigenvalues. We also expect them to be dynamically unstable to bar-mode deformations.
work was done during the workshop "'Holographic vistas on Gravity and Strings" YITP-T-14-1 at the Yukawa Institute for Theoretical Physics, Kyoto University, whose kind hospitality we all acknowledge. MM is also grateful to the String Theory Group at NTU (Taiwan) for warm hospitality. RE

A Numerics
In this appendix we explain the details of our numerical construction of the bumpy black holes.
Plugging the ansatz (2.2) and the reference metric (2.3) into the Einstein-DeTurck equations gives a system of partial differential equations that we solve numerically.
First we discretize the system using Chebyshev points. We need more resolution in the angular x coordinate than in the radial r one, so we use conforming patches, see fig. 13 for an example. This is computationally cheaper than having one bigger grid and gives us the flexibility of increasing the resolution just where it is necessary. This type of patches coincide along one line of points (no overlapping regions), in the present situation they coincided along a line of constant x. We used 2 to 5 patches depending on various factors. Higher zero modes have more bumps (the Q's vary more along x) and we need more resolution. Close to transitions the functions become singular and therefore we need to concentrate more points in a specific part of the domain. We impose continuity of the functions and the first derivatives as boundary conditions between the patches.
Once discretized, we solve the system by an iterative Newton-Raphson method. Since this method needs a seed, we first solve the linearized problem, find the eigenvectors that correspond to the zero modes and use them (ref. metric perturbed with the eigenvector) as a seed for the first solution in each of the branches. Once we have solved the nonlinear problem, we move along the branch by using the previous solution as a seed and by changing the value of k in the background metric. We keep r 0 = 1 in all the solutions.
For the standard branches that connect the MP black hole with the black ring, black Saturn and black diring, we begin by increasing the temperature. At some point the branches reach a maximum of the temperature and in order to go past it we keep k fixed and vary r + instead. The solutions close to this maximum are tricky to obtain because the Lichnerowicz operator has a near zero mode, but once we pass it the following solutions are easily obtained by lowering the temperature (decreasing k with fixed r + ). The other type of branches do not have any extrema of the temperature and to obtain them we always decrease k. As for the resolution used, we began with two patches of 20 × 20 in the branches (+) 1,2 (heading towards the black ring and black saturn) and with four patches of 30 × 20 for the (+) 3 branch (heading towards the diring); we began with similar resolutions for the other branches. In order to know when to increase the resolution we estimated the numerical error in the physical quantities and if it was greater than a few percent we decided that more resolution was needed. We have also checked that our numerical solutions converge to the continuum limit according to our discretization scheme.