Non-uniform black strings and the critical dimension in the $1/D$ expansion

Non-uniform black strings (NUBS) are studied by the large $D$ effective theory approach. By solving the near-horizon geometry in the $1/D$ expansion, we obtain the effective equation for the deformed horizon up to the next-to-next-to-leading order (NNLO) in $1/D$. We also solve the far-zone geometry by the Newtonian approximation. Matching the near and far zones, the thermodynamic variables are computed in the $1/D$ expansion. As the result, the large $D$ analysis gives a critical dimension $D_*\simeq13.5$ at which the translation-symmetry-breaking phase transition changes between first and second order. This value of $D_*$ agrees perfectly, within the precision of the $1/D$ expansion, with the result previously obtained by E. Sorkin through the numerical resolution. We also compare our NNLO results for the thermodynamics of NUBS to earlier numerical calculations, and find good agreement within the expected precision.


Introduction
Black holes in the dimension D > 4 have totally different aspects than in D = 4. A remarkable property of higher dimensional black holes is the existence of dynamical instabilities in the long wavelength, which is originally discovered by Gregory and Laflamme for the uniform black string solution (UBS) whose horizon admits the isometry of S D−3 × S 1 [1]. The Gregory-Laflamme (GL) instability happens if the circumference of the compact dimension L is in the same order of r 0 , which is the radius of S D−3 . It is also known that the zero mode of the GL instability generates another static branch, called the non-uniform black string (NUBS), in which the translation symmetry in the compact dimension is broken. NUBS solutions are constructed perturbatively [2,3] and then numerically up to D = 15 [4][5][6][7][8].
The phase of NUBS for the fixed L is known to admit a critical behavior with respect to the dimension [3]. In D ≤ 13, NUBS becomes more massive than UBS at the GL critical point and has less entropy than UBS with the same mass. Therefore, the NUBS phase is unstable and cannot be the end point of the GL instability. On the other hand, for D > 13, NUBS becomes less massive and has more entropy, and then becomes stable.
Recently, the authors and other people have shown that, in the large number limit of the dimension D → ∞, the black hole dynamics can be solved analytically by the expanding in the inverse dimension 1/D [9][10][11][12][13][14]. These simplification can be explain by the fact that the large D limit of various black hole spacetimes universally admit the 2D string black hole structure, even in the existence of the cosmological constant and rotations [15].
The key feature of the large D limit is that the gravitational dynamics are divided into the two sectors each with the separated scales in terms of the black hole radius r 0 [12]: • Non-decoupling sector: the dynamics with the gradient large as ∼ D/r 0 , which describe the degrees of freedom of the gravitational wave propagating to the asymptotic region.
• Decoupling sector: the dynamics with the smaller gradient (≪ D/r 0 ), confined within the short distance (∼ r 0 /D) from the horizon, and decoupled from the asymptotic dynamics.
The dynamics reflecting the characteristics of each black hole, including its instabilities, is described by the decoupling sector. As for the GL instability of the black string, the 1/D expansion revealed that the threshold wavenumber approaches k GL ≃ √ D/r 0 at D → ∞, as suggested from the numerical result [3]. Its correction up to the fourth order in 1/D was also calculated and reproduced the numerical result, even for relatively lower dimensions D ∼ 10 [9,10,14].
From the decoupling nature, one can expect that the dynamics of the small gradient degrees of freedom reduces to an effective theory living in the near-horizon geometry. The authors and collaborators have shown that the Einstein equation for the decoupling sector is separated to the radial direction with the large gradient ∂ r ∼ D/r 0 and the dynamics along the horizon with the smaller gradient (≪ D/r 0 ) [16]. Under this separation, the radial dependence can be easily integrated out by solving the ordinary differential equation, and the remaining dynamics along the horizon is governed by a simple effective equation, which is expressed as the 'soap-flim' equation on the some constant surface of the radial coordinate placed at the short distance (∼ r 0 /D) from the horizon. This formalism enables one to solve the Einstein equation even in the non-perturbative regime, provided the gradient along the horizon remains not so large (≪ D/r 0 ). The authors of [17] also studied the large D effective theory including the slow time dependence.
In this paper, we study the effective equation for NUBS up to the next-to-next-toleading order (NNLO) in the 1/D expansion, extending the leading order result in [16]. We also demonstrate the matching to the far-zone geometry and evaluate the asymptotic charges. With these results, the thermodynamic variables are calculated up to NNLO. The phase diagram including the NNLO correction admits the critical behavior around D * ≃ 13.5. 1 Within the error of O(D −2 ), this value agrees with numerical observation [3,8].
The construction of this paper is as follows. In section 2, we explain the setup and outline of the analysis. In section 3, the far-zone geometry is studied. In section 4, the near zone is solved by the 1/D expansion and the effective equation up to NNLO is derived. In section 5, we study the effective equation. In section 6, the thermodynamics of NUBS is studied. Section 7 gives the summary of this paper. The appendices include the details of some lengthy equations used in the main part. We attached a Mathematica file containing the detail of the NNLO near solution which are too lengthy to show here. The file also contains the series solution for the effective equation.

Setup
In this paper, we consider the D = n + 4 static and spherically symmetric spacetime with one compact dimension. Such spacetimes are described by the following metric where we assume A, B, C → 1 at r → ∞ for the asymptotic flatness. In this conformal ansatz, we solve the near-horizon geometry by the large D effective theory approach [16]. Because the gravity leaking from the near zone falls off by 1/r n , we rather use 1/n as the small expansion parameter, instead of 1/D itself. For the near-zone analysis, we introduce the another radial coordinate suitable to describe the near-horizon structure with the large radial gradient (∼ n/r 0 ), We also introduce the new compact coordinate to capture the variation in the scale of In the literature, the critical dimension D * is usually defined as the highest dimension with the first order transition: D * = 13. In this paper, we adopt the definition of D * as the real number at which the order of the transition changes, as in the condensed matter physics. Because the large D analysis treats the number D just as a parameter, this is a rather natural definition.
This scaling is a natural choice if we want to keep the self-gravitating effect finite in the large n limit. 2 Assuming a bulge formed on the horizon with the wavelength λ and radius deformed by δr 0 ∼ r 0 , the mass excess from the original horizon is roughly estimated by Then, the extra gravitational potential made by this mass excess on its neighbor becomes which means that the gravitational interaction between bulges formed at intervals of λ has the finite limit at n → ∞ only if λ ∼ r 0 /n 1/2 . For more gentle deformations like λ ∼ r 0 , the deformed horizon needs other supporting forces in the limit. For example, the authors and collaborators have demonstrated the negative cosmological constant allows several deformed solutions: droplets and cavities [16], and the rotation also allows bumpy Myers-Perry solutions [19]. As for the Kaluza-Klein spacetime, we have localized black holes other than black strings [7,20], in which case the gravitational forces from distant black holes support the deformation. Though, such gravitational effect beyond the near zone will be rather contained in the non-decoupling sector with the large gradient. In our scaling (2.3), the static ansatz (2.1) allows the deformation only up to the magnitude of r 0 /n for the valid 1/n expansion [16]. However, this does not mean the analysis is perturbative. This magnitude of the deformation still produces the non-perturbative variation in the mass scale as (r 0 + δr 0 (z)/n) n ∼ r n 0 e δr 0 (z)/r 0 . We will see that this nonperturbative degree of freedom follows a nonlinear effective equation.
Matching Strategy For the evaluation of the asymptotic charges such as the mass and string tension, the author and collaborators proposed the general way of matching with the far zone through the quasi-local energy momentum tensor [16]. Instead, here we take the more primitive method. We will directly match the near and far solution in the overlap region 1 ≪ R ≪ e n . The strategy is as follows. The far solution is simply solved by the Newtonian approximation in the Minkowski background. The Newtonian potential Φ(r, z) is expanded in the overlap region as where the function φ(z) is identified with the mass deformation e δr 0 (z)/r 0 . Since δr 0 (z) is only the functional degree of freedom in the near region, terms higher in ln R/n will not provide any additional information. On the other hand, if we take r → ∞, the Newtonian potential Φ(r, z) admits the monopole behavior The important fact is that the monopole α is obtained by averaging φ(z) over z, Since the total mass (∼ αr n 0 L) will be the integration of the local mass density e δr 0 (z)/r 0 , this is a natural consequence. Therefore, once the near deformation δr 0 (z) is solved, the asymptotic charges are determined.

Far zone analysis
In this section, we examine the far-zone geometry, to which the near solution in the 1/n expansion will be connected. In the far zone r − r 0 ≫ r 0 /n, the gravitational potential ∼ r n 0 /r n becomes much smaller than any power of n. Then, the far-zone geometry should be solved in the expansion of the Newtonian order r n 0 /r n , instead of 1/n. We also study the boundary behavior of the far-zone geometry. The far zone includes two boundary regions: overlap and asymptotic region. In the overlap region, the far solution is further expanded by 1/n and matched with the near solution. In the asymptotic region, we merely take r → ∞ without taking the large n limit to obtain the asymptotic charges.

Newtonian approximation
Let us perturb the metric (2.1) from the Minkowski background. 3 where we write A = 1 + δA, B = 1 + δB, C = 1 + δC (3.2) and each correction is assumed to be in the Newtonian order O(r n 0 /r n ). Then, the equation reduces to the Laplace equation in the cylindrically symmetric space.
where Φ is introduced by Other components reduce to where β is an integration constant, and Ψ(r, Z) follows One can realize that Ψ(r, Z) is the degree of freedom of the conformal transformation. So we just set Ψ(r, Z) = 0. Assuming the periodicity Z → Z + L, the general solution with the regular asymptotic behavior is given by the Fourier decomposition Φ(r, Z) = αr n 0 r n + ∞ j=1 a j r n/2 0 r n/2 K n/2 (2πjr/L) cos(2πjZ/L) (3.6) where K n/2 (x) is the modified Bessel function of the second kind. The reflection symmetry Z → −Z is also imposed. In the limit r → ∞, only the first term can contribute to the asymptotic monopole, since K n/2 (x) falls off exponentially at x → ∞.

1/n expansion in the overlap region
In the overlap region, we introduce the rescaled wavenumber With r = r 0 R 1/n , the Debye expansion formula leads to the expansion in terms of 1/n, where f j (k) is defined up to NNLO in 1/n as This expansion will be rather formal which breaks down for the large j. But we believe such high frequency components are negligible in the small gradient assumption. Then, the 1/n expansion of eq. (3.6) becomes where we also inserted Z = r 0 z/ √ n. The leading behavior can be sumed to a functional degree of freedom f j (k)a j (k) cos(jkz). (3.13) With this function, eq. (3.12) up to the relevant order can be rewritten in this form where each Φ (i) is given by We replaced the factor j 2 k 2 in eq. (3.10) with −∂ 2 z . Therefore, the far solution (3.4) is determined only by φ(z) and β in the overlap region. In the following, we see that β is related to the tension.

Asymptotic charges
Once φ(z) and β are specified, the asymptotic charges are calculated as follows. In the asymptotic region r → ∞, the leading behavior of Φ in eq. (3.6) is given by α, which is equal to the average of φ(z) over z The far solution (3.4) gives the monopole at r → ∞, Therefore, the ADM formula determines the mass and tension as [22] M = Ω n+1 r n 0 L 16πG where Ω n+1 = 2π n/2+1 /Γ(n/2 + 1) is the volume of S n+1 .

Near zone analysis
In this section, we study the near-horizon geometry in the 1/n expansion, following the previous formulation [16] up to the one higher order.

Setup
The near-horizon geometry can be separated into the radial sector with the large gradient (∼ n/r 0 ) and the other sector with the smaller gradient (≪ n/r 0 ), Under this assumption, it is convenient to decompose the Einstein equation as where R µ ν is the intrinsic curvature for g µν and Our coordinate ansatz (2.1) leads to the intrinsic metric as Since the far solution (3.4) has only the normalizable corrections (∼ 1/R) in the overlap region, the boundary condition for R → ∞ is simply The large dimensionality of S n+1 and the rescaling (2.3) also requires for the valid 1/n expansion. 4 Hence, the metric components are expanded by 1/n as

Leading order
As was done in [16], the leading equation for the mean curvature K becomes If we expand K as the leading solution K (0) is given by where M (z) is an integration function denoting the horizon position. Solving eqs. (4.2) and (4.4), we obtain the leading solution and Adding to the condition (4.8), the regularity at R = M (z) is also imposed for B and C. The vector constraint (4.3) is trivially satisfied in this order. By setting M (z) = 1, this solution reproduces the uniform black string metric (A.10). Here we note that M (z) does not directly denote the deformation in the radius δr 0 (z), but the variation of the mass scale M (z) ∼ e δr 0 (z)/r 0 . Recalling the definition of R, the radius of the horizon embedding surface is given by Gauge choice As another gauge choice, one can set M (z) = 1 by the radial rescaling R → M (z)R and instead, use the integration function of C (0) as the degree of freedom [16], This choice makes the radial coordinate fit the equipotential surface, which leads to the simpler near horizon analysis. 5 In this 'equipotential' gauge, the near deformation function enters in the non-normalizable behavior in the overlap region, in exchange for the simplified near-zone analysis, the matching with the far-zone geometry requires a nontrivial coordinate transformation, involving both R and z. Instead, the current embedding condition (4.8) makes the near coordinate rather fit the asymptotic geometry. The deformation in the non-normalizable behavior C 0 (z) is now absorbed to the horizon shift M (z). This choice is convenient to keep contact with the asymptotic coordinate, and in which the matching between the near and far zone will not require any further coordinate transformation other than R = r n /r n 0 . 6

Next-to-leading order
The solution at NLO is similarly calculated. The horizon position is fixed to R = M (z) by the renormalization of M (z). In the appendix. B, we show the detail of A (1) , B (1) and C (1) . At this order, the vector constraint (4.3) admits the first nontrivial condition for M (z), This is equivalent to the differentiated version of the effective equation obtained in [16]. 7 As pointed out in [16], eq (4.17) is the system equivalent to the undamped Toda oscillator [24].
Effective equation from the matching Due to the assumption δr h (z) ∼ r 0 /n, the vector constraint (4.3) provides only the one lower order effective equation. This situation can be remedied by the use of the far-zone information. We can also derive the leading order effective equation only from the leading solution (4.14), just by using the far consistency relation (3.8). The leading solution (4.14) expanded up to O(R −1 ) is matched with the far Newtonian correction (3.4) as Substituting this into eq. (3.8), we obtain the leading order effective equation This is just the integrated form of eq. (4.17). We note that this derivation is due to the fact that the far-zone equation is solved explicitly, which depends on the asymptotic structure. More general asymptotics will not allow such simplification.

Next-to-next-to-leading order
Finally, we obtained the solution up to the next-to-next-to-leading order (NNLO). Because of its lengthy form, the detail of the NNLO solution is given in the attached Mathematica file. To this order, the vector constraint (4.3) gives 8 where the derivatives higher than M ′′ (z) in O(n −1 ) are eliminated by using this equation iteratively. This effective equation determines the horizon deformation only up to NLO.
We will see that the one higher order equation can be obtained from the matching.

Matching
Now, we identify the Newtonian correction in the far zone (3.2) with the near solution up to NNLO. Expanding the near solution up to O(R −1 ), we obtain where each coefficient consists of M (z) and its derivatives as in eq.
Therefore, we have completely determined the far-zone geometry up to NNLO in 1/n. 8 The relation between this equation and the soap-film equation in [16] is not trivial, because our radial coordinate is not normal to the horizon. We need the transformation to the normal coordinate to find the relation.
Consistency check Eq. (4.21) also contains the power of ln R/n. Since we have no more undetermined degree of freedom, the constraints (3.7) and (3.8) should be trivially satisfied for such terms. Actually, we found that our near solution satisfies up to the relevant order and A 1 − (n + 1)B 1 /n = 0, A 2 − (n + 1)B 2 /n = 0.

Analysis of the effective equation
In this section, we study the non-uniform black string solution by solving the NNLO effective equation (4.22).

Integrated form
After an integration, eq. (4.22) is arranged to the following form where β is the rescaled value of tension in eq. (3.18). Each coefficient is given by a 0 = 1 + a + 1 + 4a + 2a 2 + 4 ln 2 2n + 1 n 2 a + 2a 2 + 2 3 a 3 + 2(2a + 1) ln 2 , where a is an integration constant, whose definition is tuned for the later convenience. This is one of our main results. Here ln 2 in the coefficients is just originating from the current conformal coordinate, and does not affect the physical quantities. Eq. (5.1) is invariant under the rescaling, The definition of a is arranged so that this transformation works linearly on a. Then, the constant a is found to be the scaling degree of freedom of the spacetime, which can be set to an arbitrary value. The scaling at a = 0 is also fixed so that β = 1 just admits the uniform solution.

Leading order equation
Let us start with the leading order equation For the time being, we put aside the scaling a by setting a = 0. If we define ϕ(z) = M (z), then eq. (5.4) reduces to the potential problem of the classical dynamics, where the minus of the rescaled tension −β plays the roll of the energy. The minimum energy state with β = 1 is given by the uniform solution ϕ(z) = 1. For the positive tension 0 < β < 1, this potential also admits the oscillatory solution, which corresponds to the non-uniform black string (NUBS). Since the expansion breaks down at M (z) = 0, the solution for the negative tension β < 0 is not the physical branch. The solution also collapse to M (z) = 0 for the zero tension β = 0. However, the zero tension solution can be identified with localized black holes, which will be discussed in the later section. The maximum and minimum value for NUBS are given by where W {0,−1} (z) denote the upper and lower branches of the Lambert W function, defined as the inverse function for z(w) = we w . In the limit β → +0, we have Since the expansion is only valid for | ln M (z)| ≪ n, the range of the parameters is restricted to β ≫ e −n .
Non-uniformity parameter An useful dimensionless measure for the deformation is the non-uniformity parameter [2] λ = 1 2 where R max and R min is the maximum and minimum radius of S n+1 . This is expressed by the maximum and minimum value of M (z), Using eq. (5.6), the non-uniformity parameter is estimated for small β and large n as where we assumed β ≫ e −n . Thus, the 1/D expansion covers the deformation only with the small non-uniformity λ ≪ O(1).

Higher order corrections
Now, we study the effective equation including higher order corrections. For a = 0, eq. (5.1) is arranged to gives the UBS solution for β = 1 as in the leading analysis,  (2)) ln 2 + 720 ln 2 2 − 240ζ(2) − 84ζ(2) 2 + 180ζ(3))W {0,−1} (−b/e) + 2(5 + 120(−1 + 2ζ(2)) ln 2 + 240 ln 2  The inclusion of higher order corrections in eq. (5.1) is rather straightforward. We just expand M by 1/n as M = M 0 + M 1 /n + M 2 /n 2 , and the avoidance of the secular terms also adds the higher order corrections to the wavelength k(ǫ) = k 0 (ǫ) + k 1 (ǫ)/n + k 2 (ǫ)/n 2 . If we continue the analysis, the resulting solution can be written in the Fourier series M(z) = ∞ j=0 µ j (ǫ)ǫ j cos(jk(ǫ)z). (5.19) where the j-th harmonics begins from O(ǫ j ). We show only the first few terms in the appendix. D. Further terms up to O(ǫ 7 ) are given in the Mathematica file attached to this paper. The wavenumber k(ǫ) determined by the expansion up to O(ǫ 7 ) is given by (5.20) In the UBS limit ǫ → 0, k(ǫ) reproduces the higher order corrections to k GL shown by the linear analysis [9,10,14], except the factor 4 −1/n . Actually, we will see that the mass length of the corresponding UBS is given by ρ 0 = 4 1/n e a/n r 0 . Then, if we recover the scaling (5.12), the physical wavenumber at ǫ = 0 is given by √ nk(0) e a/n r 0 = k GL = √ n ρ 0 1 − 1 2n + 7 8n 2 .  The series solution in eq. (5.19) itself is convergent even around ǫ ≃ √ 2 where the rescaled tension b = 1 − ǫ 2 /2 reaches zero, regardless of the convergence in the 1/n expansion ( figure. 1). So, in the following analysis, we use this series solution, instead of solving eq. (5.1) numerically. Since the parameter a just involves the scaling, the shape of the deformation is governed only by ǫ. Figure 2 shows that as ǫ approaches to √ 2, the deformation tends to have the broader bulge and narrower neck as seen by the numerical calculations. Analytically, this can be shown by estimating the inflection point

Thermodynamics
In this section, we study the thermodynamic properties of NUBS in the 1/n expansion.

Variables
First, we compute the thermodynamic variables from the solution (5.19). The mass and tension are derived from eq. (3.18). Using eqs. (4.23) and (C.1), the average of φ(z) can where we eliminated ln M(z) and higher derivatives of M(z) by using eq. (5.1) and integrating by part. Using the series solution (5.19), the above averages can be computed as For the latter, we also used µ 1 (ǫ), µ 2 (ǫ), µ 3 (ǫ) and eq. (5.20). The form of µ 0 (ǫ), µ 1 (ǫ), µ 2 (ǫ) and µ 3 (ǫ) is given in the appendix. D. In the UBS limit ǫ → 0, we obtain the mass of UBS on the GL threshold point Comparing with eq. (A.2), the length scale of the UBS mass on the threshold becomes On the other hand, since the horizon position r h (z) = r 0 M (z) 1/n varies in z, the horizon area is given by where B H (z) and C H (z) are the metric components on the horizon. Using eq. (5.1) and integrating by part, the area just reduces to A = 4 1/n e a n Ω n+1 r n+1 where we used the form of eq. (6.1). The surface gravity is also computed by using eq. (5.1) The Smarr's formula just follows from eq. (3.18) and the above expressions for A and κ, In fact, the scaling factor can be factored out from these variables where we introduced the scale invariant quantities . (6.10) Substituting eq. (6.2) into eq. (6.1), φ 0 (ǫ) up to O(ǫ 6 ) is given by dA + T dL. (6.12) The first law for the scaling a is easily confirmed. For the variation with ǫ, the first law reduces to the condition From eqs. (6.11) and (5.20), this condition can be confirmed at least for up to O(ǫ 6 ) and O(n −2 ).

Circumference
So far, we have considered the scaling and deformation separately. However, we have to vary the scaling a together with ǫ to obtain the physical phase diagram. If we recover the scaling, the circumference of the extra dimension is given by L = 2πe a/n r 0 √ nk(ǫ) . (6.14) Therefore, the circumference also changes as the solution gets deformed. To obtain the sequence of NUBS with a fixed L, the scaling a should be tuned so that the running of ǫ is canceled. Since the running in k(ǫ) becomes comparable to e a/n only if ǫ 2 ∼ n −1 , we introduce a rescaled non-uniformity parameter where the scaling at s = 0 is absorbed in r 0 . Therefore, the physical phase of NUBS is expressed in terms of an O(1) non-uniformity parameter s.

Phase diagram
Now, we investigate the phase diagram for the constant L. With the assumption ǫ 2 = 2s/n and the scaling in eq. (6.16), L is now fixed to where k GL is given by eq. (5.21). From eqs. (6.11) and (6.16), eqs. (6.9) and (6.7) gives where the value at the GL point is given by In figure 3, we compare our result with and the numerical result in [8] at several dimensions. The corresponding non-uniformity parameter λ is determined by eq. (5.14) with b = 1 − s/n. The variables are plotted in the unit of L GL , where we writek GL = ρ 0 k GL and set G = 1. The formula in eq. (6.18) reproduces the numerical result within the error of O(n −3 ) for the small non-uniformity. For the large non-uniformity, the 1/n expansion breaks down and one can no longer expect the good accuracy.
In figure 3, the tension seems to have the relatively large error compared to other variables. Actually, the Smarr formula (6.8) implies that the tension has the error of one lower order than other variables. This is because the tension is sensitive to the inhomogeneity of the horizon, especially the neck part of the horizon which has relatively bad convergence (figure 1). On the other hand, the mass and area, which are the averaged quantities, will not be so sensitive to the neck geometry. The surface gravity, which is defined uniformly on the horizon, also will not have such sensitivity.

Localized black holes
We should note that our formalism also includes localized black holes as the zero tension solution. If the solution has the point M (z) = 0, the regularity in eq. (4.14b) requires From the effective equation (4.19), this is the case of the zero tension β = 0. For β = 0, eq. (5.1) admits an analytic solution where we omit the scaling for simplicity. This solution is not periodic nor convergent at the large z, so one can think this unphysical. But, this can be interpreted as the near-equatorial part of the spherical solution.
Let us consider D = n + 4 Schwarzshild spacetimes in the isotropic coordinate (dρ 2 + ρ 2 dθ 2 + ρ 2 sin 2 θdΩ 2 n+1 ). (6.23) We consider the transformation to the cylindrical coordinate r = ρ sin θ, r 0 z = √ nρ cos θ, (6.24) where z coordinate describes the region around the equatorial plane θ ≈ π/2. Then, we obtain where we set ρ 0 = r 0 and the following formula is used This is identical to what we have obtained in eqs. (4.14a) and (4.14b) with the leading order of eq. (6.22). As for the higher order, though it is difficult to find the coordinate transformation, Eq. (6.22) reproduces the mass spectrum of localized black holes, correctly up to NNLO. 9 Since we have β = 0, substituting eq. (6.1) into eq. (3.18) gives 32n 2 (6.27) 9 We thank Roberto Emparan for suggesting this possibility.
where we identified L M with the integration over Z = r 0 z/ √ n in the middle line. Using the following identity eq. (6.27) reproduces the mass of D = n + 4 Schwarzschild black holes in the isotropic coordinate (6.23) up to NNLO, despite the solution (6.22) itself breaks down at the large z. As discussed in [16], this is because the contribution from the equatorial region becomes dominant for the energy density in the large D limit. However, our construction fails to capture the thermodynamics and deformation in the localized black hole phase. Eq. (6.27) is only the mass of the spherical black hole without the deformation. If we compere the mass with the NUBS and UBS phase for the same L, any size of the localized black hole has the zero mass, because of the infinite period. Then, in the large D limit, the localized black hole phase is degenerate to a single point of the phase diagram: (M = 0, T = 0). This is because the deformation effect is nonperturbative in 1/n for localized black holes, at least, in the small gradient assumption. The deformation of the localized black hole with the radius ρ 0 and extra dimension size L is supported by the gravitational potential from itself at a distance of a half period L/2. This gravitational effect has the magnitude of (2ρ 0 /L) D , which is negligible in the 1/D expansion, unless the black hole is very large L/2 − ρ 0 ∼ ρ 0 /D.

Critical dimension
The phase diagram of NUBS admits a critical behavior between D = 13 and D = 14, across which the order of the transition changes [3]. This means that NUBS has the greater mass than UBS on the GL point for D ≤ 13, and the lower mass for D > 13. In eq. (6.18a), the gradient near the bifurcation point becomes 29) in which the gradient changes the sign at n 2 − 9n − 5 = 0. (6.30) This implies the existence of the critical dimension at n * ≃ 9.5 (D * ≃ 13.5), which agrees with the Sorkin's numerical result [3], within the error of O(n −2 ). Here we note that the n = 9 case admits a marginal behavior at NLO, in which case the NNLO correction plays a crucial role. 10 In figure 4, we show the critical behavior in the phase diagram, plotted in terms of the normalized mass and relative tension. The relative tension is a dimensionless quantity defined by τ = (n + 1)LT M = (n + 1) 2 (1 − ǫ 2 /2) (n + 2)nφ 0 (ǫ) + 1 − ǫ 2 /2 (6.31) which is normalized so that the UBS limit gives τ = 1. Substituting ǫ 2 = 2s/n, we obtain where the corrections can be determined up to O(n −3 ), because O(n −3 ) correction in eq. (6.31) should start with ǫ 2 for the normalization τ (ǫ = 0) = 1. Substituting this expression into eq. (6.29), we obtain the gradient in the phase diagram, The condition n 2 − 10n + 6 = 0 also gives n * ≃ 9.4. Since the gradient in eq. (6.32) does not affect the sign, this value is identical to n * ≃ 9.5 up to the error in O(n −2 ). , respectively. We plotted the results up to the NLO (gray) and NNLO (black) corrections, in comparison with the numerical data for n = 8 (circle), n = 9 (square) and n = 10 (triangle) in [8] (the data for the relative tension, not included in [8], is produced with the tension data provided by Pau Figueras). At NNLO, the mass seems to grow by the effect of the selfinteraction, which makes the n = 9 phase bellow the critical dimension.
The horizon area is also known to admit the critical behavior at the same dimension as the mass. Above the critical dimension, the horizon area of NUBS becomes greater than that of UBS with the same mass, which means that NUBS is more favored than UBS in the microcanonical ensemble. Below the critical dimension, NUBS becomes the unstable phase. However, we cannot observe the difference in the area between NUBS and UBS of the same mass to this order, where the horizon area of the UBS with the same mass is given by n a (n + 2)nφ 0 (ǫ) + 1 − ǫ 2 /2 (n + 1) 2 n+1 n .

(6.35)
This is expected since the first law guarantees the cancelation of the linear order dependence [2,3]. In this case, both A and A UBS are the function of ǫ 2 , and then the second order becomes ǫ 4 . In addition, the coefficient will be at most O(n −1 ), since there is no entropy difference in the large D limit [10]. Then, the leading behavior of the area difference should appear at much higher order than the mass The numerical calculation actually shows A/A UBS − 1 ≃ 10 −3 ∼ 10 −4 for the small nonuniformity in D = 11 to 14 (See figure 7 in [8]). As for the surface gravity in eq. (6.18d), though its thermodynamical interpretation is unclear, we see that the gradient at the GL point changes at which gives a smaller critical value n * ≃ 8.1 (D * ≃ 12.1). Then, the NUBS phase becomes hotter in D > 12 and cooler in D ≤ 12. This effect is also observed by the numerical calculation [3,8].

Summary
In the large D limit, the near-horizon dynamics with the small gradient (≪ D/r 0 ) along the horizon reduces to the nonlinear effective theory for some collective degrees of freedom of the horizon [16,17]. This effective approach will be a powerful tool for the various nonlinear problems in the higher dimensional gravity.
In this paper, we applied the large D effective theory approach to the study of the non-uniform black strings (NUBS) bifurcating from the uniform black strings (UBS) on the threshold of the Gregory-Laflamme (GL) instability. Extending the leading order analysis [16], we solved the near-horizon geometry up to the next-to-next-to-leading order (NNLO) in the 1/D expansion. We also analyzed the far-zone geometry by the Newtonian approximation and derived the matching condition between the two zones. As the result, we obtained the effective equation of the deformed horizon up to NNLO. We also obtained the thermodynamic variables up to NNLO from the matching result. The phase diagram produced by the large D analysis agrees with the numerical results for the small nonuniformity within the expected precision. We also show that our construction contains localized black holes which have the degenerate phase in the 1/D expansion. The phase diagram admits the critical behavior in the mass around D * ≃ 13.5 as shown numerically [3,8]. We cannot find the area difference between NUBS and UBS of the same mass in this order. It will be straightforward (though much harder) to obtain much higher order corrections, in which the critical behavior in the area difference will be observed.
In the higher dimension, various black objects admit the GL-like instabilities and their zero modes are expected to produce the new deformed branches. For example, Myers Perry black holes were conjectured to be unstable for the sufficiently large rotation [25], and it is numerically shown that each zero mode deformation leads to the new 'bumpy' solution [26,27]. AdS 5 × S 5 black holes also have the deformed blanches [28]. These deformed solutions will show the similar critical behavior. The large D effective theory approach will be a good analytical tool to observe such nonlinear phenomenon. The authors have already shown that the bumpy Myers-Perry solution can be obtained by the 1/D expansion [19].
As the further extension of this work, it is interesting to study the time evolution of the non-uniformity in the 1/D expansion. This will be partially realized by including the slow time dependence, which is already formulated in [17]. For D > 13, the NUBS phase becomes stable and the time evolution of the GL instability can be followed by the current small gradient effective approach. Since the numerical simulation becomes more difficult in the higher dimension by the resolution problem, the large D analysis will be a suitable approach to this problem. For D ≤ 13, the NUBS phase is no longer stable, and the GL instability develops to the self-similar structure which consists of the sequence of black holes and thin strings connected with each others [29][30][31]. Since the small gradient assumption will break down in such phase, we need to incorporate the large gradient (∼ D/r 0 ) dynamics which is not yet formulated.
The large gradient dynamics will be also required to understand static solutions supported by the interaction between the multiple horizons, because the small gradient degrees of freedom are confined in the near-horizon geometry. In section 6.3.1, we show that localized black holes have the degenerate phase in the large D limit. This is because the gravitational force becomes non-perturvative in 1/D, (2ρ 0 /L) D . To keep the deformation finite in the large D limit, one should assume very large black holes L/2 − ρ 0 ∼ ρ 0 /D, which will involve the large gradient dynamics. Near the merger point, the NUBS phase was conjectured to have the conical waist [32,33]. In such phase, the gradient along the horizon will be also large.
If the extra dimension is compactified with the circumference L, the mass and tension of the spacetime are given by where Ω n+1 = 2π n/2+1 /Γ(n/2 + 1) is the volume of S n+1 . This metric can be transformed to the conformal coordinate ds 2 = −F 2 (r)dt 2 + dr 2 + dZ 2 + ρ 2 (r)dΩ 2 n+1 . (A.3) by the transformation where r coincides with ρ at ρ → ∞. F (r) is implicitly given by In this coordinate, the horizon position F (r h ) = 0 becomes (A.6) Large D limit Though eq. (A.5) cannot be solved explicitly, the large n limit allows the expression in the 1/n expansion. Keeping R = 4r n /ρ n 0 finite, the limit n → ∞ in eq. (A.5) leads to where we introduced the horizon scale r 0 = 4 −1/n ρ 0 . Thus, in the conformal coordinate, the near-horizon geometry of the uniform black string in the leading order becomes ds 2 ≃ − R − 1 R + 1 2 dt 2 + dr 2 + dZ 2 + r 2 0 R 2/n 1 + We note that the horizon radius r h in eq. (A.6) slightly differs to r 0 in the higher order, (A.11)

B.2 Next-to-next-to-leading order
The NNLO solution is given in the Mathematica filed attached. We only note that B (2) includes the following integration which we could not express by known functions (C.1)