Critical Kaluza-Klein black holes and black strings in D = 10

We construct static vacuum localized black holes and non-uniform black strings in ten spacetime dimensions, where one of the dimension is compactified on circle. We study the phase diagram of black objects with these boundary conditions, especially near the critical point where localized black holes and non-uniform black strings merge. Remarkably, we find that the merger happens at a cusp in the phase diagram. We verify that the critical geometry is controlled by a Ricci-flat double cone as previously predicted. However, unlike the lower dimensional cases, we find that physical quantities approach to their critical values according to a power law plus a logarithmic correction. We extract the critical exponents and find very good agreement with the predictions from the double cone geometry. According to holography, localized black holes and black strings are dual to thermal states of (1 + 1)-dimensional SU($N$) maximal Super-Yang Mills theory compactified on a circle; we recover and extend the details of the (recently found) 1st order phase transition in this system from the gravity side.


Introduction and results
In D = 4 spacetime dimensions, stationary asymptotically flat black hole solutions of the Einstein equation in vacuum of a given mass have spherical topology, are presumably unique and all the evidence suggests that they are dynamically stable. However, in D > 4 these properties change radically. The physics of black objects turns out to be much richer, allowing for nonspherical topologies, instabilities and non-uniqueness (and thus phase transitions). Black holes are fundamental objects in general relativity. In recent years the study of such objects in nonastrophysical settings has received much attention due to the intrinsic interest in understanding fundamental aspects of gravity as described by general relativity (see [1] for a review), and also because of the connections to string theory and the gauge/gravity duality [2][3][4][5][6]. In the latter context, it is natural to consider spacetimes that are asymptotic to M d × N n , where M d is ddimensional Minkwoski or Anti-de Sitter space and N n is an n-dimensional compact manifold so that the total number of spacetime dimensions is D = d + n > 4.
One of the most extensively studied models in this setting is that of M d = Mink d and N 1 = S 1 , a circle of length L. Since the compact dimension is flat, it is trivial to write down a black hole solution that is uniformly wrapped along the compact dimension: This is just given by the (D − 1)-dimensional Schwarzschild solution times a (compact) flat direction, Schw D−1 × S 1 . Such a higher dimensional black hole is known as the uniform black string (UBS). In [7], Gregory and Laflamme (GL) famously showed that thin enough black strings are unstable under linear gravitational perturbations with a non-trivial dependence along the S 1 -direction. 1 Determining the endpoint of such an instability has been the subject of intense studies during the past few years.
At the onset of the instability, the linear GL mode is time-independent (i.e., a zero mode) and can be continued to the non-linear regime. This indicates that there exists a new branch of black strings which are non-uniform in the compact direction and are thus known as non-uniform black strings (NUBS). NUBS were first constructed perturbatively in D = 5 by [8], subsequently constructed fully non-linearly in various spacetime dimensions using numerical methods [9][10][11][12][13][14][15] and, more recently, using the large-D expansion [16]. It turns out that in D < D * = 13(.5), NUBS have less entropy than UBS with the same mass and hence they cannot be the endpoint of the GL instability [17]. In fact, based on entropic arguments, [7] conjectured that unstable UBS would evolve into an array of localized black holes through a dynamical topology change transition; the latter can only happen through a singularity and hence the evolution of the GL instability of black strings could potentially constitute a counter-example of the weak cosmic censorship conjecture [18,19] around such spacetimes. This scenario was recently confirmed by [20,21], using numerical relativity techniques. 2 On the other hand, for D > D * , NUBS can be dynamically stable and hence be the endpoint of the GL instability, as [16,24] confirmed.
Apart from UBS and NUBS, spaces that are asymptotically Mink d × S 1 also admit static black hole solutions that are localized on the S 1 . These localized black holes (LOC) have been constructed numerically [12,[25][26][27] and perturbatively in the limit in which the black holes are 1 More precisely, the condition for the existence of a linear instability is that r 0 /L O(1), where r 0 is the mass parameter of the parent Schw D−1 solution and L is the asymptotic length of the compact circle. 2 Notice that this final fate is not exclusive of UBS and black holes with compact extra dimensions. Fully nonlinear time evolutions of analogous instabilities in asymptotically flat black rings or ultra-spinning Myers-Perry black holes spacetimes, similarly lead to violations of the weak cosmic censorship conjecture [22,23]. small compared to L [28][29][30][31][32][33]. Motivated by geometrical considerations, [34] conjectured that the NUBS and LOC branches should merge at a topology changing critical solution governed by a Ricci-flat double-cone. This conjecture was tested from the black string side in [35], and later in various dimensions in [13]. However, the most non-uniform black strings in these early constructions were still too far from the critical regime to provide conclusive results (see however [11]). Only recently, Kalisch et al. [15,27], in an impressive numerical construction, have managed to obtain NUBS and LOC in D = 5, 6 extremely close to the critical point, confirming the doublecone predictions to an unprecedented level of detail.
The goal of the present work is to construct NUBS and LOC in D = 10 very close to the critical point, where these branches of black holes merge. Critical solutions have only been previously constructed in D = 5, 6 [15,27]; for higher values of D, gravity becomes more localized near the horizon of the black object, which makes the numerical construction more challenging, especially very close to the critical point. Note that [36] previously constructed both NUBS and LOC in D = 10, but their solutions were very far from the critical regime since the aim of that paper was different (see below). Here D 0 = 10 is the critical dimension of the double-cone geometry [34] and D * = 13(.5) is the critical dimension in the microcanonical ensemble [17].
At the critical dimension D * the dynamical stability of weakly NUBS changes from being unstable for D < D * to being stable for D > D * . However, for D = 12, 13 [13] found that NUBS with a sufficiently large non-uniformity can also be dynamically stable. 3 This paper showed that there exists a turning point, i.e. a maximum of the mass/area, along the NUBS branch where the stability properties of the solutions change. On the other hand, in D = 5, 6 such a turning point is present along the LOC branch. In D = 10, as we move along both the LOC and the NUBS branches and approach the critical solution from both sides, we do not find any turning points on either of the branches. Therefore, the simplest picture that emerges from our work is that, sufficiently far from the critical solution, in D < 10 there should exist a turning point along the LOC branch, in D > 10 the turning point occurs along the NUBS, and in D = 10 there are no turning points at all. See Fig. 1.1. Recently, [16] confirmed the existence turning points in the phase diagram of NUBS in D < 14 using the large-D expansion, but the reliability of their approach breaks down at around D ≈ 9. Notice, however, that the methods of [16] did not allow them to study critical solutions in detail and therefore our results complement theirs.
Our numerical data suggests that in D = 10 the merger happens precisely at a cusp in the phase diagram. The study of the critical geometry in [34] showed that D = 10 is the critical dimension of the cone geometry that governs the topology change. For D < 10, the approach of physical quantities to their critical values is controlled by a (dimension-dependent) power law with infinitely many oscillations (i.e. turning points); this behavior has been beautifully confirmed in [27] for D = 5, 6. On the other hand, for D > 10 the approach to the critical point should be given by two independent power laws, with no oscillations. D = 10 is the marginal case and the approach to the critical point should be controlled by a power law with a logarithmic correction. In this paper we confirm this in D = 10.
In this paper we also compute the spectrum of negative modes of the Lichnerowicz operator, ∆ L , around the LOC and NUBS branches, restricted to modes that preserve the isometries of the background. Just as in [12,13], we find that NUBS posses two negative modes: one is continuously connected (as the non-uniformity parameter goes to zero) to the negative mode of the parent Schwarzschild black hole [37]. This mode diverges as the NUBS approach the critical solution. The other negative mode is the continuation of the GL zero mode to non-zero values of the non-uniformity parameter and our data suggests that it tends to a finite value at the critical solution. On the other hand, LOC have a single negative mode throughout the branch and it approaches the same finite value as the NUBS at the critical solution. See Fig. 3.9. Note that at an extremum of the temperature one can have zero modes corresponding to variations of the parameters of the solution that respect the boundary conditions (i.e. preserve the temperature). We do not find any evidence for new zero modes, which is consistent with the absence of extrema of the temperature along either branch of solutions.
Another motivation for the present work comes from the gauge/gravity duality [2][3][4][5][6]. The best well-understood example of this correspondence is between maximally supersymmetric Yang-Mills (SYM) theory in p + 1 dimensions and gauge group SU(N), and Type IIA (even p) or Type IIB (odd p) superstring theory containing N coincident Dp-branes in the decoupling limit. For p = 1 the duality is between 2-dimensional SU(N) SYM theory and type IIA or IIB string theory in the presence of D0-or D1-branes respectively [6]. At large N, strong coupling and finite temperature, the gauge theory is described by black hole solutions with D0-or D1-charge in the supergravity approximation, depending on the temperature (type IIA at low temperatures and type IIB at high temperatures respectively). In this paper we are interested in vacuum black hole solutions of the Einstein equation in 10 spacetime dimensions, one of which is compactified on a circle of length L. After a series of standard U-duality transformations [36,[38][39][40], these vacuum black holes can be given D0-or D1-brane charges.
According to the AdS/CFT correspondence, the black hole phase structure should be reproduced by the thermal phases of SYM on a circle at strong coupling and large N. Lattice simulations of SYM on a torus, with one of the circles being the thermal circle and with periodic boundary conditions for the fermions on the other spatial circle, have been performed. Most of the previous works in the past have focused on the p = 0 SYM quantum mechanics and agreement with the gravity predictions has been confirmed [41][42][43][44][45][46][47][48][49][50][51][52][53][54]. The case p = 1 has received less attention in the past [55], until the recent of work of [56]. This paper predicted the temperature at which a first order phase transition occurs from lattice simulations, in a regime where the latter should overlap with the supergravity calculations. The latter was only recently computed in [36] and found very good agreement with the lattice result. In this paper, as a by-product of our calculations, we recompute the value of the phase transition temperature (or energy, in the canonical ensemble); the values that we obtain are t crit = 1.09257 t GL and ε crit = 1.24181 ε GL for the temperature and energy at the phase transition measured with respect to the GL point. Our values differ with those found in [36] by less than a 0.25%. In addition, we have been able to locate the merger between the non-uniform and localized phases.
The rest of this paper is organized as follows. In §2.1 we start by reviewing some general aspects of black holes in Kaluza-Klein spaces. In subsections 2.2 and 2.3 we present our numerical construction of NUBS and LOC, respectively. In §3 we present our results. §3.1 contains the phase diagrams in the microcanonical and canonical ensembles, in §3.2 we consider the horizon geometry and in §3.3 we study in the detail the critical behavior of NUBS and LOC near the critical point and compare it with the predictions of the double-cone model. §3.4 is devoted to the computation of the spectrum of negative modes of the Lichnerowicz operator around the NUBS and LOC. §4 contains the results for the phase diagram of the supergravity solutions with D0-charge. We close the paper with a discussion in §5. Some technical details are relegated to the Appendices. In appendix A we give more details about the integration domain that we have used to construct the localized black holes and in appendix B we present some convergence tests. The mapping from neutral solutions to charged ones is presented in detail in appendix C.
Note added: while this paper was nearing completion, we became aware of [57], that has some overlap with ours and that has appeared on the arXiv on the same date.

Black objects in Kaluza-Klein theory
Consider vacuum Einstein's gravity in D = 10 spacetime dimensions with Kaluza-Klein (KK) asymptotic boundary conditions, i.e. Mink 9 × S 1 . For (ultra)static spacetimes, this theory contains three different families of static black holes, namely, UBS, NUBS and LOC. After fixing the overall scale by fixing the length of the asymptotic S 1 , these three different types of black holes can be parametrized by the temperature and one may distinguish them by the topology of the horizon and the isometries. Whilst UBS are translationally invariant along the S 1 and are known explicitly, for NUBS and LOC the translation invariance along the S 1 is broken and they have to be constructed numerically (or pertubatively). In this section we explain our numerical construction of such solutions. Since we are interested in studying the thermal phases, we will be working with the Euclidean form of the solutions where the Euclidean time τ is periodic, τ ∼ τ + β, with β being the inverse temperature. 4

Generic results and Uniform black strings
In this paper we are interested in Einstein metrics that asymptote to the flat Euclidean metric, where one of the directions corresponds the Euclidean time τ , times a KK circle of length L. As usual, the Euclidean time τ is compact and has period β. Ultimately we will consider ten dimensional spaces but for now we shall keep the total number of spacetime dimensions D general. Moreover, we will only consider spacetimes that preserve an SO(D − 2) subgroup of the full rotation group of the flat Euclidean metric in D − 1 dimensions. Therefore the asymptotic isometry group of the spaces that we shall consider is U(1) β ×SO(D − 2)×U(1) L , which is made explicit in the asymptotic form of the flat metric on the product space S 1 β × R D−2 × S 1 L , with τ ∼ τ + β and y ∼ y + L. For more general spaces, from the asymptotic behavior of the metric components one can extract two asymptotic charges, namely mass and tension, of the solution [58]: 5 From these quantities one can define the relative tension n = T L/M, which is bounded: n ≤ 0 ≤ D − 3. In addition to these charges, NUBS and LOC can be characterized using their own geometric quantities which are discussed in §3.
where (τ, r, Ω D−3 ) are the usual D − 1 Schwarzschild coordinates and y is the S 1 coordinate. The parameter r 0 labels each solution and it is directly related to the physical quantities: (2.5) (Notice that the uniform black string has C τ = r D−4 0 , C y = 0 and constant relative tension n = (D − 3) −1 .) Finally, recall that the topology of the horizon is S D−3 × S 1 .
Gregory and Laflamme [7] famously discovered that thin enough UBS, i.e. r 0 /L O(1), are dynamically unstable to clumping along the compact direction. More precisely, for fixed L there is a critical value r GL 0 below which there exist regular (linear) perturbations that grow exponentially with time and that break the translational invariance along the S 1 ; at precisely this critical value, the perturbations are time-independent thus signaling the existence of a linear solution of the Einstein equation which is not uniform along the S 1 . This linear solution can be continued into the fully non-linear regime, giving rise to the NUBS. For D = 10, the critical value of the horizon radius at the onset of the GL instability is: r GL 0 = 0.36671(3)L. Our aim in this work is to numerically construct vacuum NUBS and LOC solutions in D = 10. In practice, we will treat the different metrics as smooth Riemannian manifolds with a U(1) β Killing vector that vanishes at the horizon and solve the Einstein vacuum equation, R ab = 0, subject to certain regularity and asymptotic boundary conditions. As is well-know, due to the underlying gauge invariance of the theory, this equation does not yield a well-posed boundary value problem. Instead, we solve the Einstein-DeTurck equation, R H ab = 0, which is manifestly elliptic [12], where R ab is the Ricci tensor and ξ a is the so-called DeTurck vector. The last is formed from the usual Levi-Civita connection Γ compatible with the spacetime metric g, and a Levi-Civita connection Γ compatible with some reference metricḡ that we are free to prescribe. This has now become a standard approach in stationary numerical relativity and we refer the reader to the literature for more details [12,[59][60][61].
The equations are always discretized using pseudo-spectral methods on a Chebyshev grid, and we solve them by an iterative Newton-Raphson method; at each step of the iterative process the linear system of equations is solved using LU decomposition implemented by subroutine LinearSolve in Mathematica.

Numerical construction of Non-uniform black strings
NUBS wrap the KK circle, and, for regular solutions, the horizon S D−3 is finite everywhere. This implies that with our symmetry assumptions, the integration domain has the following effective boundaries: the horizon, asymptotic infinity and the periodic boundary. Due to the symmetry of first GL harmonic, NUBS have a Z 2 -symmetry and then one has an additional mirror boundary. Hence, a single coordinate patch is enough to cover the whole computational domain. In practice, to numerically construct highly non-uniform black strings near the critical point it is convenient to use more than one patch to get enough resolution in the regions of interest.
To find NUBS we consider the following ansatz for the metric: The UBS satisfies all the relevant boundary conditions that we will impose on our solutions (see below) and we shall use it as the reference metric in the Einstein-DeTurck equation. The compact radial coordinate x ∈ [0, 1) covers the region from the horizon (x = 0) to infinity (x = 1). Note that NUBS posses reflection symmetry along the S 1 direction. This allows us to consider only one half of the KK circle subject to mirror boundary conditions. Therefore, we take y ∈ [0, 1], where y = 0 corresponds to the reflection plane and y = 1 the periodic boundary. This implies that the asymptotic length of the KK circle is kept fixed to be L = 2.
The radius of the round S D−3 at the horizon is a good geometric invariant that can be used to describe NUBS; with our ansatz (2.7), this is given by whereas NUBS have λ > 0; the limit λ → ∞ corresponds to the merger point with the LOC branch, where R min → 0 while R max remains finite.
To obtain a well-posed boundary value problem that can be solved with elliptic methods we need to supplement the equations of motion with appropriate boundary conditions. These require regularity at the horizon, reflection symmetry, periodicity and KK asymptotics: • Horizon at x = 0: smoothness of the metric at the horizon implies that all Q's must be even in x and therefore we impose Neumann boundary conditions on all Q's, except the crossed term which must be Dirichlet. The condition Q 1 (0, y) = Q 2 (0, y) ensures that the geometry is free of conical singularities and fixes the surface gravity of the solution to be that of our reference metric.
• Reflection plane and periodic boundary at y = 0 and y = 1 respectively: all Q's must be even in the compact S 1 coordinate and thus we impose Neumann boundary conditions for all Q's, except for the crossed term which must be Dirichlet there.
To find NUBS, we start with the UBS close to the GL point and add a bit of the GL zero mode. This gives a good initial guess that allows us to find weakly non-uniform black strings. Once we have found a NUBS, we can move along the family varying the temperature; with our boundary conditions, the inverse temperature is given by We move along the branch of NUBS by using the previous solution as a seed and varying the value of the parameter r 0 ; we start at r 0 = 0.73450 which corresponds to λ = 0.04 (recall that r GL 0 = 0.73342(6) for L = 2) and, given our modest resources, we move up to a value of r 0 = 0.79184, corresponding to λ = 5.05.
For λ 1, the NUBS are relatively weakly non-uniform, not much resolution is required to construct the solutions accurately and one single patch is suffices. At this point, the solutions satisfy ξ a ξ a ≡ ξ 2 < 10 −10 , with estimated numerical error to be less than 0.01%. The Smarr's relation is satisfied up to the order 10 −7 . As we move along the branch of NUBS to greater values of λ, the function Q 4 develops very pronounced peaks near the origin, corresponding to the waist of the non-uniform black string, and some form of mesh-refinement there is needed to construct accurate solutions. We found that two conforming patches were enough to obtain good results, though the bound on the DeTurck vector goes up to ξ 2 < 10 −7 and the Smarr's relation is satisfied up to 10 −6 . Notice that our mesh-refinement introduces a new parameter x 0 , which is the coordinate location where the two patches meet. We also consideredỹ = mesh(y; 0, 1, χ), with the mesh-refinement function mesh(. . . ) given by (A.5); here χ is just a parameter that controls the density of the new grid points. Since the steep gradients move towards the origin as λ increases, we used two different setups with appropriate grid sizes x 0 (∼ 10 −1 , 10 −2 ) and values of χ (∼ 1, 10). It is possible that by choosing a different reference metric for highly NUBS one can achieve larger values of λ without losing accuracy.

Numerical construction of Localized black holes
To numerically construct LOC we follow the approach of Kalisch et al. [27] with minor modifications. Essentially, we considered a different compacitifaction of the radial coordinate so that we could extract the constants C τ and C y appearing in the conserved charges (2.3) by calculating 1st derivatives of our unknown functions. In this section we superficially discuss the actual numerical construction of LOC and refer the reader to [27] for further details.
We seek static axisymmetric black holes that are asymptotically KK and localized on a circle of (asymptotic) length L. We choose adapted coordinates so that symmetries of the spacetime become manifest. This implies that the actual boundaries of the computational domain are: the black hole horizon, the asymptotic infinity, the periodic boundary, the reflection plane and an axis of symmetry where the horizon S D−2 smoothly shrinks to zero size, which is exposed because the localization on the S 1 . From the point of view of finding these black holes numerically, since the integration domain has five boundaries, we naturally work with two coordinate patches: One patch adapted to a 'near' region (containing the horizon), and another one adapted to a 'far' region (containing the asymptotic infinity). The integration domain is schematically shown in Fig. 2.1.
As in [12], one can work with cartesian coordinates (x, y) in the far patch and polar coordinates (r, a) in the near patch, and the relation between them is simply given by the polar map: x = r cos a, y = r sin a. To transfer information between the two coordinate patches, one can use two overlapping domains and impose uniqueness of the solution [12]. This is simpler to implement if one uses finite differences. On the other hand, if one uses spectral methods, one can deform the two domains using some transfinite transformation and ensure that the two domains match along a curve; along this common boundary, one then imposes continuity of the functions and their normal derivatives. Alternatively, in the near region [27] introduce polar-like coordinates with a modified radial coordinate which naturally matches with the Cartesian coordinates sufficiently far from the black hole. This is the approach we follow in the near region. We recall the details of the integration domain and introduce the new compactification in appendix A. The ansatz for the metric in the far patch is: where the functions Q ≡ {Q 1 , Q 2 , Q 3 , Q 4 , Q 5 }(x, y) are our unknowns. As illustrated in the integration domain A.1, the coordinate x ranges from L/2, which is the boundary between the far and near regions, to infinity; on the other hand, y ∈ [0, L/2], where y = 0 is the reflection plane and L/2 is the periodic boundary. The boundary conditions we impose on the unknown functions Q in this patch are: • Asymptotic boundary at x = ∞: the metric must approach the KK space. This implies the Dirichlet boundary conditions, Q i (∞, y) = 1, ∀i = 5, and Q 5 (∞, y) = 0.
• Matching boundary at x = L/2: we impose continuity of the metric and its normal derivative.
• Reflection plane and periodic boundary at y = 0 and y = L/2 respectively: all Q's must be even in the compact coordinate y and thus we impose Neumann boundary conditions on all them except the crossed term Q 5 , which must be Dirichlet there.
The near horizon region ansatz covers the horizon and the symmetry axis; at the horizon, the Killing ∂ τ becomes null and at the symmetry axis the round S D−3 (and in fact the whole horizon S D−2 ) smoothly shrinks to zero. The ansatz we consider is: are the unknowns in this patch. This metric has a Killing horizon located at r = r 0 with surface gravity κ, and an axis at a = π/2; a = 0 is the reflection plane. Although the horizon is at r = r 0 , r 0 is simply a parameter in our ansatz and we keep it fixed throughout the calculation (we choose r 0 = 0.8 for convenience); the physical parameter labelling each solution is the surface gravity κ, and this the parameter that we vary to move along the branch of LOC. With the definitions given in Appendix A, the boundary conditions that we impose on the unknown functions Q ′ in this region are: • Horizon at r = r 0 : smoothness of the metric at the horizon implies that all Q ′ 's must be even in r and therefore we impose Neumann boundary conditions for Q ′ 1 , and Dirichlet for the crossed term Q ′ 5 . The condition Q ′ 1 (r 0 , a) = Q ′ 3 (r 0 , a) ensures that the geometry is free of conical singularities and fixes the surface gravity of the solution to be that of the reference metric (see below).
• Reflection plane at a = 0: all functions are Neumann except for Q ′ 5 , that vanishes there.
• Periodic boundary at r 3 (L/2, a) sin a = L/2: using the relation between the far and near coordinates and the relation between the far and near unknown functions, one can find the boundary conditions for the near horizon functions from the boundary conditions that the far region functions satisfy there.
• Matching boundary at r 2 (L/2, a) cos a = L/2: we impose continuity of the metric and its normal derivative.
In addition to the ansatz and the boundary conditions, the DeTurck scheme requires a global reference metric as part of the gauge fixing procedure. The reference metric must satisfy the same boundary conditions as the solution we seek. For the LOC, there is no known Einstein metric in closed analytic form that satisfies the required boundary conditions and hence one has to design it. In this paper we follow [12,27], and smoothly glue together two metrics, each of which satisfy the desired boundary conditions in each region, i.e. asymptotic Kaluza-Klein space in the far region and, for instance, the asymptotically flat Schwarzschild black hole in D dimensions: and where the function E(r) is given by ( 2.15) The reference metric depends on κ, r 0 , and r 1 , which is an additional parameter that can be adjusted (see appendix A). The function E(r) vanishes exponentially fast for r → r 1 , and the reference metric (2.12) tends to the KK space written in polar coordinates. For r → r 0 , E(r) ≃ 1 − κ 2 (r − r 0 ) 2 and (2.12) takes the form of the near horizon metric of the Schwarzschild black hole in D dimensions.
To find localized black holes we start with the reference metric as a seed with κ = 2.4 and L = 6 (we keep this value of L for all solutions). Recall that the convergence of Newton's method strongly depends on the choice of the initial seed and finding a first solution may be difficult. To stay within the basin of attraction at each iteration, we introduce a parameter α ∈ R + in the iteration loop so that the update is Q (n+1) ∼ Q (n) + αδQ, where α ∼ O(1/100) or O(1/10) during the first iterations and is O(1) towards the end. Once we have found the first solution, we use it as the initial guess to find the next solution with a slightly different κ while keeping α = 1. We kept the parameters of the integration domain and the coordinates fixed throughout the calculation and they are specified in Fig. A.1.
The most critical solution we found corresponds to κ = 1.262768. In this critical regime, the functions Q ′ 2 and Q ′ 4 develop steep gradients near the axis and the horizon; to resolve them, we redefine these two functions in the pure polar patch (blue and green dots in Fig. A.1): [27]. The boundary conditions for these redefined functions can be easily found from the original ones for Q ′ 2 and Q ′ 4 . All solutions we found satisfy ξ 2 < 10 −10 with numerical error less than 0.01% and the Smarr's relation is satisfied up to the order 10 −6 .

Results
In this section we present our results for both NUBS and LOC. We first consider the behavior of the various thermodynamic quantities along each branch of solutions and the phase diagrams, and then we focus on the horizon geometry. We then study the critical behavior near the merger point, and provide evidence that the double-cone geometry proposed by [34] does indeed control the merger. We finally compute their spectrum of negative modes of the Lichnerowicz operator.

Thermodynamics
The horizon temperature labels both NUBS and LOC and is given by (2.9) for non-uniform black strings and by κ/(2π) for localized solutions. The mass and the tension follow from (2.3). For the NUBS the asymptotic charges are computed by where in these expressions we first interpolate the numerical data and then perform the integration. For the LOC, these quantities are given in (A.3). The horizon area is found to be: In Fig. 3.1 we display the phase diagram in the microcanonical (top left) and canonical (top right) ensembles, and the behavior of the horizon area (middle) and tension (bottom) as a function of the inverse temperature (normalized by L). The behavior of the mass and the relative tension as a function of the inverse temperature is similar to that of the area and tension and we do not display the corresponding plots here. To make the microcanonical and canonical phase diagrams easier to visualize we plot the dimensionless differences ∆S/L 8 ≡ One of the remarkable features of the phase diagram in D = 10 is the lack of turning points away from the merger along any of the branches, either LOC or NUBS. This should be contrasted with the phase diagram in D = 5, 6, which exhibits a turning point along the LOC branch at some maximum mass and then there is a minimum of the temperature [12,27]. It is reasonable to expect that such a turning point (away from the merger) exists on the LOC branch for any dimension D < 10. This turning point switches to the NUBS branch in D = 12 (and presumably in D = 11), as shown in [13] and more recently in the large-D expansion in [16]. As we will argue below, the lack of turning points away from the merger in the phase diagram in D = 10 may be related to the nature of the merger in this specific number of spacetime dimensions. These plots reproduce and complete those shown in the appendix of [36]. Dimensionless horizon area A H /L 8 (middle) and tension T /L 6 (bottom) as a function of the dimensionless ratio β/L. The GL critical point is indicated with a solid black disc. The dimensionless mass and relative tension plots are very similar to the ones shown above.
In Fig. 3.2 we plot various physical quantities, normalized by their value at the GL point, against the normalized relative tension n/n GL . Close to the merger point, our results in D = 10 show that the physical quantities do not approach their critical values following a spiraling behavior, with presumably infinitely many turning points, as in D = 5, 6 [11,27]. Instead, the physical quantities of the NUBS and LOC branches merge at a cusp in the phase diagram, with no oscillations. As we discuss in §3.3, this behavior is precisely what the double-cone model of [34] for the merger predicts in D = 10. Notice that the physical quantities corresponding to both branches emerge from the cusp in the 'same direction'.

Horizon geometry
In this subsection we display the behavior of various geometric quantities defined on the horizon along the branches of solutions. Then we present the embeddings of the horizon geometry into flat space to help to visualize the geometry of NUBS and LOC. We characterize the NUBS using the non-uniformity parameter λ defined in [8], (2.8). In addition, we consider the proper length of the horizon along the S 1 : Following [12], for LOC, one can define R eq as the equatorial radius of the horizon round S D−3 , R eq = r 0 Q ′ 2 (r 0 , 0). Similarly, one defines L polar to be the proper distance from the 'south' pole to the 'north' pole along the horizon S D−2 , and L axis to be the proper distance between the poles along the axis: Recall that the asymptotically flat Schwarzschild solution in D dimensions is spherically symmetric and hence it enjoys the symmetry of the full rotation group SO(D − 1). On the other hand, LOC break the SO(D − 1) symmetry down to SO(D − 2), and only for very small localized black holes, i.e. high temperatures, the full SO(D − 1) is approximately recovered. We can characterize the deformation of the horizon geometry by comparing the area of the round equatorial horizon S D−3 , A eq ∝ R D−3 eq , and the area of the geodesic S D−3 on the horizon that contains both poles, A pol ∝ R D−3 pol with R pol = r 0 Q ′ 2 (r 0 , π/2). We compare these two areas defining the eccentricity parameter, A spherically symmetric black hole has zero eccentricity and ǫ diverges for the critical solution.
The behavior of these geometric quantities along each branch of solutions is displayed in Fig.  3.3. In the top row we display ǫ and λ as functions of β/L. At high temperatures, LOC are nearly spherically symmetric and the eccentricity is very small. In fact, ǫ remains quite small until pretty close to the merger with NUBS, where it diverges (notice that the vertical axis is in a log-scale). This explains why perturbation theory works so well for localized black holes in D = 10 [36], and it is another manifestation of the fact that gravity becomes more localized near the horizon as D grows. The behavior of the non-uniformity parameter λ for the NUBS is qualitatively similar. From these two plots it is clear that we managed to get closer to the merger from the LOC branch. From the behavior of ǫ and λ we can estimate that the merger occurs at β Merger ≃ 0.829L.
At the bottom of Fig. 3.3 we display the remaining geometric quantities as functions of the relative tension n normalized by its value at the GL point. We have added zooms of these plots to better appreciate the region where the various curves merge. In D = 10 the merger happens at a cusp, with the physical quantities of both the NUBS and the LOC coming out of the cusp in the same direction. This behavior should be contrasted with the D = 5, 6 case, in which a part from the shrinking spirals, the physical quantities for the NUBS and the LOC approach the merger from opposite sides. It would be nice to understand this behavior from the double-cone geometry. From the behavior of L axis /L and R min /L as they approach zero, we estimate the value of n/n GL at the merger to be n Merger ≃ 0.139n GL .
A useful way to visualize the geometry of τ = const. sections of the horizon is by embedding them into (D − 1)-dimensional Euclidean space E D−1 , with a flat metric For NUBS, the horizon geometry can be described as a surface X = X(y), Y (y) = R(y) in E D−1 , whilst for LOC one has X = X(a), Y (a) = r 0 cos a Q ′ 2 H . In each case, the embedding coordinate is given by NUBS: (3.10) In Fig. 3.4 we plot Y /L vs X/L for some representative solutions, including the most critical ones. We postpone the detailed comparison with the double-cone metric to the next subsection.

Critical behavior at the merger point
Kol argued that the merger between the NUBS and the LOC implies a topology change not only of the horizon geometry but in fact of the whole Euclidean manifold [34]. This is a much stronger statement than simply considering the change of the topology of the horizon. Moreover, . Note that the embeddings look 'rounder' or 'fatter' compared to the ones in lower dimensions; this is just a manifestation that gravity becomes more localized as D increases.
[34] conjectured that this topology change of the Euclidean manifold should locally be controlled by a Ricci-flat double-cone over S 2 × S D−3 : This double-cone arises as follows. Both the NUBS and the LOC possess an explicit SO(D − 2) spherical symmetry which must be inherited by the critical metric, i.e. it must contain a round S D−3 . The S 2 is less obvious and its origin is the following [34]: away from the waist, the Euclidean time, which is periodic to avoid a conical singularity at the horizon, is fibered over an interval whose endpoints are on the horizon, thus giving rise to a two-sphere. Such an S 2 is finite everywhere on the localized phase whilst it is contractible to zero size in the black string phase (see [62] for a nice depiction). On the localized phase, one can compute the radius of this sphere on the symmetry axis at the equidistant points from the poles of the horizon S D−2 . By symmetry, this corresponds to the equatorial radius of the S 2 and is given by One can compare it to the radius of this S 2 along the symmetry axis, R axis = L axis /(2π), along the branch of LOC. See Fig. 3.5. From this plot we see that R τ ∼ R axis as the solutions approach the merger and both radii tend to zero. This shows that the S 2 becomes round as it shrinks, just as the double-cone model of [34] predicts. Also shown in this plot is the minimum radius of the horizon S D−3 , R min , on the NUBS. This quantity also shrinks to zero at the merger. One can further test the double-cone model of the merger by considering the embedding of the τ = const. section of (3.11) into Euclidean E D−1 space. The embedding coordinates of the double-cone metric (3.11) are simply given by (3.13) In Fig. 3.6 we compare the embedding of the double-cone in D = 10 dimensions with the embeddings corresponding to the most critical LOC (red) and NUBS (blue) solutions that we have found. As this plot shows, the double-cone can be smoothed in two different ways, each one leading to one of the phases at each side of the transition. One can consider deformations of the double-cone metric of the form [34]: (3.14) The linearized perturbations satisfy the following equation of motion: For D < 10, the imaginary part of s ± causes oscillations in ǫ(ρ), while for D > 10 there are two independent (real) powers. Furthermore, [63] argued that the behavior of the deformations of the double-cone metric (3.16) should be reflected in the behavior of the physical quantities of NUBS and LOC sufficiently close to criticality. The argument goes as follows: if the zero mode ǫ(ρ) measures the deviation from the double-cone, then any physical quantity Q near the critical solution should behave as where δQ ≡ Q − Q c and ρ 0 is the typical length scale associated to the smooth cone. Recently, [27] has beautifully confirmed this prediction in D = 5, 6. The linearized solutions (3.16) degenerate in D = 10. Hence this is the critical dimension of the double-cone metric [34]. In this degenerate case, Frobenius' method gives two independent solutions of the form: ǫ D=10 (ρ) ∼ c 1 ρ 4 + c 2 ρ 4 ln ρ . (3.19) In the remaining of this subsection, we fit the different physical quantities of the near critical solutions that we have constructed according to the double-cone's prediction (3.19). Without loss of generality, for any physical quantity near the merger we have 20) where {a, b, c, d} are the fitting parameters and x measures the distance to the critical solution. We consider the following dimensionless quantities that tend to zero at the merger: where L is the length of the KK circle and r GL 0 is the horizon radius of the black string at the GL instability point given in §2.1. Any other definition of x should give equivalent results up to a rescaling. We use Mathematica's FindFit routine to carry out the fits.
In Fig. 3.7 we present the fits for the mass (normalized with respect to the values of a UBS at the marginal GL point) for the NUBS and LOC branches. The other physical quantities behave in a qualitatively similar way and we do not present the fits here. Note that in contrast to the D = 5, 6 cases, in D = 10 the physical quantities do not present any oscillations as they approach their critical values. In fact, the fits clearly show that the approach to the critical value is governed by a power law with a logarithmic correction, in very good agreement with the double-cone prediction (3.19).

Non-uniform black strings
Localized black holes In Table 1 we present the values of fitting parameters for the various physical quantities. To do the fits, we only have considered the solutions close enough to the merger, i.e. with small enough x; including more data points to perform the fit gives less accurate values of the critical thermodynamical values and exponent. For different physical quantities, the critical exponent coincides with the theoretical prediction of 4 with deviations of less than 0.05% in the worst case and the critical value of a given quantity coincide up to the 4 th or 5 th decimal number for both branches. We note that the critical values satisfy the Smarr's relation to the order 10 −6 and 10 −5 for NUBS and LOC respectively, which is consistent with the numerical error according to the values of ξ 2 we reached.  Table 1: Critical exponent and other parameters obtained from the fit of the nonuniform black strings (1st rows) and localized black holes (2nd rows) data points.
Only a couple of geometrical lengths do not follow the behavior (3.20), as it may be seen from Fig. 3.8. These are the horizon length L hor of the black string and the polar length L pol of the localized black holes. In lower dimensions this was also the case, and a linear term was introduced to get a proper fit [27]. In D = 10 the linear term appears naturally and the real critical exponent agrees to be one from both sides of the merger, just as in D = 5, 6. The equivalent plots to Fig. 3.7 for these lengths are shown in Fig. 3.8 and the extracted critical values and exponents are in Table 2. It would be interesting to better understand why these quantities do not follow the same critical behavior as the other physical quantities.  Table 2: Critical exponent and other parameters obtained from the fit of the NUBS's horizon length (1st row) and LOC's polar length (2nd row).

Spectrum of negative modes
In this subsection we present the spectrum of negative modes of the Lichnerowicz operator, ∆ L , around the NUBS and LOC solutions that we have constructed. The negative eigenvalues of ∆ L  are an invariant feature of the geometry and hence they provide another way to characterize the merger between NUBS and LOC. To compute the negative modes of ∆ L , we take advantage of the fact that, when using Newton's method to construct the solutions numerically, we have to linearize the Einstein-DeTurck operator as part of the iterative process. Around an Einstein metric, the linearized Einstein-DeTurck operator coincides with the Lichnerowicz operator [12]. It is then easy to readapt the code to find the low lying eigenvalues and eigenmodes of ∆ L , associated to (physical) metric fluctuations. Notice that with this approach we only find perturbations that are singlets under the action of U(1) β × SO(D − 3). We display the results in Fig. 3.9. We found that NUBS have two negative modes, as in lower D [12]. 6 The first one (green line in Fig. 3.9) corresponds to the continuation of the GL zero mode to a negative mode as one moves along the branch to larger non-uniformities. The other one (blue line in Fig. 3.9) is continuously connected to the negative mode of the UBS, which arises from the negative mode of Schwarzschild [37]. For the explored range of solutions in this work, no further negative modes appear on this branch. On the other hand, LOC have only one negative mode (red line in Fig. 3.9). For small black holes, this coincides with the negative mode of the asymptotically flat Schwarzschild solution in D = 10, as expected. No further negative modes appear or disappear along this branch. As we approach the critical solution from both sides, one of the modes of the NUBS diverges while the other appears to tend to a finite value; the latter seems to match the limiting value of the negative mode on the LOC branch. Notice that in D = 10, NUBS and LOC seem to have a different number of negative modes near the critical region. The reason is that there are no turning points along either branches, so modes cannot appear or disappear at a minimum of the temperature. However, modes can diverge as they approach the critical solution since it is singular. In D = 5, 6, there is a minimum of the temperature on the LOC branch, at which point ∆ L has a zero mode that continues to a (second) negative mode near the merger. At the merger, these two negative modes approach those of the NUBS and, in particular, a pair of them diverge at the critical solution [12].

Implications for Super-Yang Mills on T 2
In this section we (re)derive the thermodynamics of SYM on T 2 , see also [36,39,40,55,56], using the (neutral) KK black holes. We start with a lightening review of the different limits under which string theory can be described by its supergravity sector and then we use the solutions previously found to obtain the thermodynamics of those carrying D0-brane charge. Our results extend those in [36] and allow us to find the merger point.

Toroidal limits and type IIB/IIA supergravity duals
Consider (1+1)-dimensional SU(N) SYM at large N with 't Hooft coupling λ = Ng 2 YM . If the theory is at finite temperature T so that β = 1/T is the period of the thermal circle, and the spatial direction is also compactified on a circle of length L, then we can think of the theory as being defined on a 2-torus, T 2 = S 1 β × S 1 L . In these circumstances, we can define dimensionless quantities, t = T L, λ ′ = λL 2 to study different regimes of the theory. From the gauge theory perspective, phase transitions can be inferred by studying the expectation values of Wilson loops, P β and P L , wrapping the temporal and spatial circle respectively, which serve as order parameters. For SYM on T 2 the expectation value P L changes from zero (confined phase) to non-zero values (deconfined phase) upon heating the system, whereas P β is always non-zero at all temperatures (see [40,55,56] and references therein). Now we consider the dual gravity description. This is given by the near-horizon geometry of the spacetime sourced by a stack of N coincident D1-branes of type IIB string theory [6], with a periodic identification on one spatial coordinate. In particular, one is interested in nearextremal black D1-brane configurations of the gravitational theory in the decoupling limit, which were studied in [40]. The classical type IIB supergravity description is valid provided that N is large, to suppress string quantum corrections; α ′ -corrections are negligible when t ≪ √ λ ′ , while winding modes around the circle can be ignored when t ≫ 1/ √ λ ′ . The two conditions imply the window of validity of IIB supergravity description In this range one may use the type IIB supergravity solution to derive the thermodynamics of 2-dimensional SYM. In this regime, the thermal vacuum of type IIB supergravity is a black hole carrying D1-brane charge that uniformly wraps the compact circle. This solution is thought to be stable and corresponds to the uniform phase. At temperatures t ∼ 1/ √ λ ′ , stringy winding modes become unstable and the type IIB supergravity description is no longer valid. However, one can perform a T-duality transformation acting on the spatial circle, exchanging the theories IIB ↔ IIA and hence the charges D1 ↔ D0, so we can use type IIA black brane solutions to describe thermal states of SYM on T 2 on that range of temperatures. In this case, the requirement that the supergravity solution is valid gives Since D0-branes are point-like (instead of string-like, as the previous D1-branes), they can distribute the charge over the circle in various ways, either being uniformly or non-uniformly distributed, or localized on the compact circle. These three possibilities give rise to a non-trivial phase diagram, and it is then a dynamical question which case is preferred.
For higher temperatures (or lower energies) the charged UBS is thought to be dynamically stable. However, there exists a range of temperatures t GL < t < t PT for some t PT , where it is thought (and now known) that the uniform solution becomes globally thermodynamically less favored than the localized black hole solution. Then the temperature t PT represents a first order thermal phase transition between the uniform and localized phase. In the literature this has been termed the Gregory-Laflamme phase transition [55]. The natural interpretation of this picture on the dual gauge theory is a confinement/deconfinement phase transition.
In this section we find the temperature t PT and also t Merger , at which the non-uniform and the localized phase merge. Note, however, that the non-uniform phase never dominates any ensemble. So far, lattice simulations on the gauge side estimated the ratio t PT /t GL ∼ 1.5 [55].
To determine the precise ratio from the gravity dual theory one would need to construct the near-extremal charged solutions, take the near-horizon limit and extract their thermodynamic quantities. Clearly, solving the full system of supergravity equations is a formidable numerical task. However, it is possible to generate charged solutions from uncharged ones via a process of uplifting + boosting + KK reduction [36,[38][39][40]. Therefore, we can consider the neutral (vacuum) solutions we have previously found and from these obtain the thermodynamics that determine the phase structure of SYM theory under consideration.
Recent construction of KK black holes in D = 10 determined the energy or temperature to be ε PT = 97.067(N 2 /λ ′ 2 ) = 1.245ε GL or t PT = 2.451/ √ λ ′ = 1.093t GL [36]. Since the derivation of the mapping {Black hole thermodynamics on R 1,8 × S 1 } → {SYM thermodynamics on Ì 2 } was derived in there we do not include it here; for completeness, it is rederived in detail in the appendix C. Using (C.10), we have: Applying this map to the neutral UBS one gets the well-known results [40]: At r 0 = r GL 0 , these expressions correspond to the values ε GL and t GL previously discussed.

Thermodynamics
In this subsection we construct the phase diagrams in the microcanonical and the canonical ensemble describing the thermodynamics of SYM gauge theory on T 2 using the supergravity approximation. Both diagrams are shown below in Fig. 4.1 and reproduce and complete those in [36]. In addition, we determine for first time the merger point between charged NUBS and LOC. In Fig. 4.1 we plot dimensionless entropy or free energy difference between a given phase and the uniform one: Then the uniform phase is represented by a simple horizontal (black) line at the origin of the vertical axis. UBS are unstable for ε < ε GL in the microcanonical ensemble and for t < t GL in the canonical ensemble. The non-uniform phase, which exists beyond this point, is never dominant. On the other hand, the localized phase intercepts the uniform one at The ratios are: ε PT /ε GL = 1.24181 and t PT /t GL = 1.09257, and they are consistent with the predictions from the studies of SYM on the lattice. For energies or temperatures greater than this value, the uniform phase is dominant, and the LOC dominate the corresponding ensemble otherwise. The phase transition is first order, and it is thought to correspond to a confinement/deconfinement phase transition in the gauge side.
Our results allows us to determine, for first time, the merger between localized black holes and non-uniform black strings with D0-brane charge. From Table 1 (4.7) Note that since the non-uniform phase never dominates any of the ensembles, it may be difficult to test these numbers using lattice simulations.

Discussion and outlook
In this paper we have constructed NUBS and LOC in D = 10 and followed these two branches very close to the merger point. D = 10 is special from the point of view of this system since this is the critical dimension for the Ricci-flat double-cone metric that was conjectured to control the merger [34]. By fitting the physical quantities of both NUBS and LOC close to the merger point, we have shown that in D = 10, their approach to their critical values is governed by a power law plus a logarithmic correction, in accordance to the double-cone model. This result should be contrasted to the results in D = 5, 6 obtained in [27], which exhibit a spiraling behavior of the physical quantities as they approach their critical values. Moreover, we have found evidence that in D = 10, the merger happens at cusp in the phase diagram, and physical quantities belonging to the NUBS and the LOC emerge from the critical point in the 'same direction'. This feature should be related to the fact that D = 10 is the critical dimension for the double-cone metric and it would be very interesting to understand it in detail. To further confirm the double-cone model of the merger, one should construct NUBS and LOC very close to the critical point in D > 10 and verify that the physical quantities approach their critical values according to the predictions of the double-cone model. Work in this direction is underway.
In this paper we have not discussed the dynamical stability of NUBS and LOC. Ref. [24] considered the evolution of the GL instability of black strings in the large-D expansion and showed that they settle on a stable NUBS. More recently, [16] included corrections beyond the leading order term in the large-D expansion, and found that the endpoint depends on the thickness of the initial black string. For the cases where the black string is expected to pinch off, [24] could not follow the evolution all the way to the end. It would be very interesting to study the evolution of the instability of uniform black strings for large yet finite D. The techniques used in [23] seem appropriate and we are currently investigating this problem.
With the methods of [27] and those used in this paper, one can study the details of the mergers of other black hole systems of interest. In particular, for lumpy and localized black holes in AdS 5 × S 5 . Moreover, recently [56] obtained accurate results of the thermal phase diagram of 1 + 1 SYM on a twisted torus using lattice simulations. It would be very interesting to compare their results using the supergravity approximation, but to do so one needs to consider NUBS that are electrically charged with respect to a 2-form.

Acknowledgments
We would like to thank Toby Wiseman for enlightening discussions. PF would like to thank the Institute for Theoretical Physics (University of Amsterdam) for hospitality during the final stages of this work. BC and PF are supported by the European Research Council grant ERC-2014-StG 639022-NewNGR "New frontiers in numerical general relativity". PF is also supported by a Royal Society University Research Fellowship (Grant No. UF140319).

A Generic integration domain for localized solutions
In this appendix we describe the integration domain we have used to construct the localized black holes. Due to their nature, the numerical construction involves to work in two separate coordinates systems, one adapted to the asymptotic behavior and another one adapted to the near horizon behavior.
The construction considered in [27] takes the near chart (r, a) with five boundaries and divides it into three different subdomains. This encompasses the horizon, the axis, the boundaries of the internal space and a shared boundary with the far patch. The blue and green regions in Fig. A.1, say region 1 in the near patch, are covered by polar coordinates r ∈ [r 0 , r 1 ], a ∈ [0, π/2] whose relation with the far patch is simply given by x(r, a) = r cos a, y(r, a) = r sin a. The orange and yellow regions, say regions 2 and 3, are covered by polar-like coordinates with a modified radial coordinate which is parametrized in terms of v ∈ [r 1 , L/2]. In the orange region a ∈ [0, π/4], whereas in the yellow one a ∈ [π/4, π/2]. The precise relation with the far coordinates is: x(v, a) = r k (v, a) cos a, y(v, a) = r k (v, a) sin a, with Here r 0 and r 1 (< L/2) are parameters that we are free to specify, and L is the asymptotic length of the Kaluza-Klein circle. Notice that this construction assumes that the angular coordinate in the near patch 1 is further divided into two subregions: one patch where 0 ≤ a ≤ π/4, to match the density of grid points with that in the near patch 2, and another one where π/4 ≤ a ≤ π/2 to match the density of grid points with region 3. The far chart (x, y) covers the ranges L/2 ≤ x < ∞ and 0 ≤ y ≤ L/2. To deal with the infinity, the coordinate x is compactified introducing a new coordinate −1 ≤ ξ ≤ 1. Ref. [27] considers x(ξ) = L/(1 − ξ), such that ξ = −1 corresponds to the shared boundary x = L/2 with the near patch and ξ = 1 corresponds to asymptotic infinity. The problem with this is that to find the charges C τ and C y , i.e. the mass and the tension, one needs to consider the asymptotic expansion of the metric components up to (D − 4)th order, which implies to take D − 4 derivatives. Of course, this is problematic for D > 5, 6. We overcome this issue by considering the compactification where ∆ = (D − 4) −1 has been defined in §2.2. This way we still have x(ξ = −1) = L/2 and x(ξ = 1) = ∞, but it is sufficient to consider the metric expansion at infinity up to 1st order. In particular, the charges are given by 1st derivatives of the metric, Additionally, each patch has been further divided into other small subregions in order to be able to increase the grid resolution just where it is necessary. This is of particular interest since as we increase D gravity turns out to be more localized and the spacetime region close to the horizon needs special care. Moreover, close to the merger point with NUBS, some functions develop steep gradients. In practice, the radial coordinate in the near patch 1 is divided into two subdomains, and the compactified coordinate ξ in the far patch is divided into three subdomains. In the near patches containing the axis, the angular coordinate is also divided into two subregions. In total, this introduces four new parameters in the integration domain: r * , ξ * , ξ * * and a * , corresponding to the values where the different patches meet. At each shared boundary, either near-near, near-far or far-far patch, one must impose continuity of the functions and their first normal derivatives.
To impose these matching conditions one may consider the same grid point densities from both sides of a given shared boundary. Alternatively, one can still require continuity of the function and its normal derivative by performing the matching on an interpolation function. Unlike [27], we have opted to work with the same grid point densities. They are naturally always the same except at the shared boundary between near and far patches. We fix this by considering the coordinate y given in terms of a coordinate σ lying in the unit interval σ ∈ [0, 1]: Using Chebyshev grid points for σ, then the grid points along the y-direction are properly distributed.
Finally, we consider the same mesh-refinement as in [27], properly modified to take into account our redefinition of the angle a, near the axisã = mesh(a; π/2, a * , χ 1 ). We also use this type of mesh-refinement near the horizonr = mesh(r; r 0 , r * , χ 2 ), with the function mesh(. . . ) given by To check whether our code with the described modifications gives rise to reasonable solutions and, in particular, accurate values for the mass, we compare the obtained results, (i) for small localized black holes, with the mass of a Schwarzschild black hole in D = 10, or (ii) with the perturbative results. For small localized black holes one expects that the spacetime metric can be systematically expanded in a perturbation series with a small parameter ρ 0 /L, being ρ 0 the location of the horizon. The best available perturbative approximation for the thermodynamic quantities, with D arbitrary, are given in [28]. (We did not include these curves in our plots in §3.1 or §4.2 since they were not much clarifying.) For small enough black holes, i.e. with eccentricity ǫ < 10 −3 , our numerical values differ by less than a 0.05% when compared to those obtained by (i) or (ii). From the geometrical point of view, another check is to compare L polar defined in (3.7), with one half of the perimeter of a Schwarzschild black hole of the same temperature in D = 10. In this case the deviations are always less than 0.01%.

B Convergence tests
In this appendix we check that our numerical solutions converge to the continuum limit according to our discretization scheme. In the case of using pseudo-spectral methods, the error should be exponentially suppressed with increasing the grid size. To monitor it we use the squared norm of the DeTurck vector ξ 2 . We expect it to become zero in the continuum limit. Indeed, Fig. B.1 shows that our numerical implementation exhibits the expected behavior.
To produce this figure we picked up a reference solution of each branch, we interpolated it at different resolutions and then we filtered through the Newton-Raphson loop. For each output we computed the quantity of interest. For NUBS we just considered one single patch with resolution N = N x N y , being N x and N y the number of grid points in each direction. In the case of LOC, we considered the usual 12 patches and varied the mean resolutionN.

C D0-charge via Uplifting + Boosting + KK reduction
In this appendix we derive the mapping between the thermodynamics of neutral Kaluza-Klein black solutions and the thermodynamics of near-extremal D0-black branes on a circle of type IIA supergravity [36,[38][39][40]. This involves a M-theory lift-boost-reduce procedure.
Consider any static, axially symmetric metric solving R ab = 0 in D spacetime dimensions and approaching the direct product manifold R 1,D−2 ×S 1 asymptotically. The bar notation in S 1 (of lengthL) is to distinguish, in D = 10, the T-dual circle of type IIA supergravity from the original circle of length L of the type IIB theory. This solution can be written using isotropic coordinates: ds 2 = −f 2 dt 2 + g 2 dρ 2 + ρ 2 dΩ 2 D−3 + h 2 dy 2 , (C.1) with the generic functions f, g and h approaching 1 at ρ → ∞, and y ∼ y +L is the coordinate ofS 1 . If (C.1) is a black hole with a Killing horizon located at ρ = ρ 0 , then f (ρ 0 , y) = 0. We can construct the dimensionless quantity p 0 ≡ ρ 0 /L to label a given family of such metrics. Now uplift the solution adding a compact coordinate z, dŝ 2 D+1 = ds 2 + dz 2 . Boosting along z with rapidity parameter α yields a solution to vacuum general relativity in D + 1 dimensions. Upon dimensional reduction with respect to z we rewrite the metric as For D = 10, this procedure allows us to construct solutions of type IIA supergravity with only graviton, dilaton and 1-form field excitations turned on. The new metric, the non-trivial dilaton field and the 1-form gauge field are identified to be: f 2 H 2 dt 2 + g 2 dρ 2 + ρ 2 dΩ 2 D−3 + h 2 dy 2 , e φ = H 1/ζ , This works because momentum around the circle in the (D + 1)th dimension is reinterpreted as D0-brane charge from the lower-dimensional viewpoint [39]. We can start with (C.1) being a neutral non-uniform black string or localized black hole, and obtain the charged solution using (C.3) (which depends on the parameter α in addition to ρ 0 ). Since we are interested in their thermodynamics rather than the solutions themselves, we will proceed by expressing the quantities of interest of the new charged solutions in terms of the uncharged ones. For instance, it easy to see that the temperature and the entropy of the charged solution are simply shifted by a factor of cosh α with respect to the uncharged ones. To be precise, Now set D = 10. The derived quantities so far correspond to the thermodynamic quantities of the charged solution (C.3) of type IIA supergravity with a metric, a dilaton and a 1-form. At this point note that all physical quantities appearing in (C.5) and (C.7), reduce to those of the uncharged solution in the limit α → 0. The opposite limit α → ∞ corresponds to take the near-extremal limit. To obtain the desired mapping we have to take the decoupling limit of Ref. [6] of near-extremal configurations, which sends the string length ℓ s to zero while keeping g 2 YM fixed. According to holography, these are dual to the decoupled field theory at finite temperature. To do so, one needs the relation between the 10-dimensional Newton's constant and the string length and coupling constant, 16πG 10 = (2π) 7ḡ2 s ℓ 8 s , the T-dual relations [64]:L = (2πℓ s ) 2 /L, g s = (2πℓ s /L)g s , and the relation between the string coupling constant and SYM coupling. In the case of IIB string theory with D1-branes, this is g 2 YM = (2π) p−2 g s ℓ p−3 s with p = 1 [6]; to translate back the quantities above in terms of SYM variables recall that λ ′ = λL 2 = Ng 2 YM L 2 . The dimensionless energy above extremality, ε = LE = L(M −Q), in these limits corresponds to the energy density of the SYM theory. Since this becomes independent of ℓ s , the decoupling limit is trivial. The decoupling limit of the temperature and the entropy needs to be taken with more care. At the end, one finds that the dimensionless energy, temperature and entropy associated to a stack of N-coincident near-extremal D0-branes in the decoupling limit are: ε = 16 3 π 7 (4c t (p 0 ) − c y (p 0 )) N 2 λ ′ 2 , t = 2π 5/2 2c t (p 0 )t(p 0 ) 1 √ λ ′ , s = 16 3 π 11/2 s(p 0 ) 2c t (p 0 )