Black ripples, flowers and dumbbells at large D

We explore the rich phase space of singly spinning (both neutral and charged) black hole solutions in the large D limit. We find several ‘bumpy’ branches which are connected to multiple (concentric) black rings, and black Saturns. Additionally, we obtain stationary solutions without axisymmetry that are only stationary at D → ∞, but correspond to long-lived black hole solutions at finite D. These multipolar solutions can appear as intermediate configurations in the decay of ultra-spinning Myers-Perry black holes into stable black holes. Finally, we also construct stationary solutions corresponding to the instability of such a multipolar solution.


Introduction
Black hole solutions in higher dimensional gravity show a far richer behavior than their counterparts in four spacetime dimensions. In higher dimensions, the rotation plays a significant role to fertilize a variety of new solutions. Since in D > 5, the (Newtonian) gravitational potential ∼ GM r D−3 falls off more rapidly than the centrifugal barrier ∼ J 2 M 2 r 2 , the horizon can be deformed to an extended shape at large angular momentum, and hence becomes vulnerable to a Gregory-Laflamme type instability [1,2]. This allows a family of non-uniform stationary solutions to branch off from the zero modes of the instabilities [3].
The increased number of degrees of freedom in a higher dimensional theory, however, complicate the construction of black hole solutions and analysis of their dynamics. To tackle this problem, several approximation techniques have been developed. One such approximation is the blackfold approach [4], which has been successful in elucidating the black hole phases in the ultra-spinning regime: for example for black (multi-)rings/Saturns in which the horizon has highly elongated shape , that allows to locally approximate them as loosely bent black strings/branes.
Another successful effective approach is the large spacetime dimension limit, or the large D limit [5,6], which has been proven to be useful in various problems involving higher dimensional black holes . This limit allows black holes to have a simple near horizon structure decoupled from the asymptotic region [31]. As a result, the Einstein's equation reduces to an effective theory on the horizon surface expanded in 1/D, namely the large D effective theory [13,14,32,33]. Different to the blackfold approach, the large D limit is naturally endowed with the separation of scales between gradients along and orthogonal to the horizon: the gradient orthogonal to the horizon becomes large compared to gradients along the horizon in the limit of large D as a result of the steepening of the gravitational potential. This enables us to formulate an effective theory without the requirement that the gradients along the horizon have to be infinitesimal, which makes the large D expansion a powerful tool to study the non-uniform 'bumpy' phases of black holes.
In this paper, we explore the phase space of compact stationary solutions with a single spin in the large D limit, specifically, we focus on the (non-)axisymmetric deformed families branching off from the Myers-Perry family. The instability of ultra-spinning MP black holes JHEP04(2020)108 In the ultra-spinning regime J /M > 2 the MP-BH develops instabilities and the corresponding zero modes appear at places marked with dots or crosses. For the analytically known black bar, we also study its non-uniform deformations ('dumbbells'), whose branches are shown in different shadings of a color to make them more distinguishable. and the existence of nearby 'rippled' solution was first conjectured in [34] and later, after the proof of existence of the zero modes and the instability [35][36][37][38][39], the rippled solutions were constructed numerically and identified as solutions that connect to black rings and black Saturns [4,[40][41][42].
Because of the strong suppression of gravitational radiation at large D [28], the effective large D description also admits stationary non-axisymmetric branches such as black bars [27] and other multipolar solutions. Here we apply the blob approximation developed in [27,29], where localized black hole solutions such as the Myers-Perry black hole are identified as stationary lumps ("blobs") on a membrane which share the same horizon topology as the black brane solution but nevertheless encode most of the physics pertaining to the localized solution. Figure 1 shows the full phase space plot of solutions we obtain. The solutions correspond to Myers-Perry solutions and their axissymmetric 'bumpy' deformations leading to black rings and black Saturns. We are also including stationary solutions without axisymmetry, which only can remain stationary at large D since gravitational radiation decouples. These solutions have been shown to play an important role in dynamical evolutions of the ultra-spinning instability [28,[43][44][45][46]. The first solution of this kind, a dipolar solution "black bar" was found analytically in [27]. Here we study its stationary deformations and also find its multipolar generalizations "black flowers". To illustrate features of the found solutions, we show plots of the mass density of four examples in figure 2. 1 JHEP04(2020)108 We observe that most of bumpy deformations remain tangential to their 'parent'branch until the deformation becomes comparable to the original solution and new blobs start to form. At some point, these blobs barely have any overlap and the branches enter a new asymptotic behavior for small Ω becoming completely separated. Some very short branches stick out non-tangentially above the parent-branch.
The paper is structured as follows: in section 2, we outline the derivation of our large D effective equations for black branes and describe how they also contain localized black hole solutions. In section 3, we construct perturbatively and numerically stationary 'bumpy' deformations of the MP black hole that lead to (multiple) black rings and Saturns. In section 4 and 5, we construct stationary non-axisymmetric solutions from multipolar deformations of MP black holes and deformations from black bars. Section 6 discusses effects of adding charge to obtain charged (but non-extremal) solutions. In the appendix we collect details of the perturbative calculation and describe our numerical procedure in greater detail.
2 Branes and localized black holes at large D

Large D effective equations
We study possibly charged black holes in Einstein-Maxwell theory in higher dimensions with n large and p a finite number. Ref. [18] developed an effective theory for fluctuations of p-branes along their extended directions σ i (i = 1, . . . , p) , where R = r n and The electric potential is The degrees of freedom of the effective theory are the mass density m(t, σ), the charge density q(t, σ) and the fields p i (t, σ). In the presence of charge it is convenient to introduce a new field v i (t, σ) defined by and the abbreviation The equations of motion of the effective theory are obtained by requiring that the Einstein-Maxwell equations are solve to leading order in a 1/D-expansion and take the form of conservation equations These equations simplify further if we consider only stationary configurations, that satisfy and v i is a time-independent killing vector i.e., Under these assumptions the equations of motion are reduced to a single master equation that is most elegantly formulated in terms of the area-radius R = ln m , (2.18) and is given by Using the scale invariance of the effective equations, which manifests itself in a shift symmetry of R, the above equation can be formally mapped to the uncharged equation by defining the charge rescaled velocity field 20) and shifting R to obtain the soap bubble equation [18] v 2 Which has the same form as the uncharged equation (i.e., eq. (2.19) with q = 0) but with the difference that the role of v 2 is now taken by the norm of the charged rescaled velocity field. Since the charged equation can be mapped to the uncharged one, solving eq. (2.21) for a given value of v q always gives a one parameter family of solutions, parameterized by the charge parameter q. In the case of non-vanishing charge, v q is not directly the physical velocity field and allows to study the effect of charging up the solution.

Black holes as Gaussian blobs on a membrane
Even though these equations were initially formulated to capture the dynamics of black branes. Ref. [27] found that this large D effective theory also contains localized black hole solutions when solved with different boundary conditions. We recapitulate here the findings of [27,29].
To capture effects of a single spin we consider the case of p = 2 and require the stationary solutions to have a purely rotational velocity field. Choosing angular coordinates for the spatial brane directions σ i = (r, φ), the only non-vanishing component of the (charge rescaled) velocity field can be set to v φ = Ω q and equation (2.21) becomes where Ω q is the charge rescaled angular velocity, according to eq. (2.20).

JHEP04(2020)108
The Myers-Perry (MP) black hole solution (and its charged Kerr-Newman counterpart described in [29]) corresponds to the axisymmetric solution (2.24) Since this corresponds to a Gaussian in the mass variable m = exp R, this solution is strongly localized in the directions σ i , but still shares the same horizon topology as the black brane (2.3). This feature of the solution is due to the fact that the rescaling of the spatial directions σ i → σ i / √ n assumed in eq. (2.3) leads for localized solutions to a magnification of the region around the center of one of its hemispheres. Since at large D most of the surface of the black hole is concentrated in this region, a description of it can capture most of the physics connected to the localized black hole.
The aforementioned localization of the mass density motivates the following definition of a localized black hole: we call a solution of eq. (2.22) a (stationary) localized black hole, if it has a finite mass M according to And it has an angular momentum given by where we used p φ = ∂ φ m + Ω r 2 m.

Axisymmetric sector: black ripples
First, we consider the axisymmetric deformation of the Myers-Perry, which leads to an infinite number of 'bumpy' solutions, or black ripples.

JHEP04(2020)108
Introducing a new radial variable z via the deformation equation becomes a Laguerre equation with a quadratic source term where we introduced the Laguerre operator L. We note that, in terms of the new variable, the MP-solution is now written as Perturbations of this solution should be normalizable in the sense of eq. (2.25), which means for the perturbed profile m = exp(R MP + δR) which is accomplished if the perturbation grows as a polynomial at each order, not showing exponential growth ∼ e z or any divergences. At leading order, the regular and normalizable perturbations are given by Laguerre polynomials [27], only if a 2 + 1 = 2N , for integer N . Non-trivial solutions have N ≥ 2. N has the interpretation of a 'radial overtone' number, i.e., it counts the number of oscillations along r. Since these zero modes correspond to 'bumpy black holes' [34,40,41], N can also be interpreted as the number of bumps in the cross-section of the corresponding solution.

Nonlinear perturbations
In the following, we study how to include higher order perturbations for these zero-modes obtaining better control over the phase space of stationary solutions and also to support the later numerical analysis. The general perturbative soution to eq. (3.4) is written as and for a leading order solution with a 2 + 1 = 2N , (N = 2, 3, 4, . . . ), the deformation equation (3.4) becomes at each perturbation order k. As usual, the source term S k (z) is expressed by the solution up to (k − 1)-th order. A similar higher order perturbation analysis has been performed in [12,21] for perturbations (non-uniformities) of black strings. It was found there, that the length of the JHEP04(2020)108 black string has to be renormalized to avoid secular terms that would break the periodic boundary condition. Here, for spinning localized solutions, it turns out that we have to renormalize the angular velocity Ω or the corresponding spin parameter a which changes the blob size, to avoid secular behavior that would break the normalization condition (3.6).

Resonance and secular perturbation
Secular behavior in perturbation theory is typically encountered when the dependence of some physical parameter on the perturbation parameter ε is ignored. A common example for this is the slightly anharmonic oscillator Note that if we assume x 1 the lowest order effect of the anharmonic term εx 3 is to modify the frequency: ω 0 → ω 0 + εω 1 . The appropriate ansatz accordingly should be x(t) = sin((ω 0 + εω 1 )t), but naive perturbation theory x(t) = x 0 (t) + εx 1 (t) leads to the solution (3.11) x 1 (t) = t · sin(ω 0 t) + . . . , (3.12) where the first correction grows unboundedly invalidating the perturbative ansatz and violating conservation of energy. Note here that the secular term (3.12) results from a resonance phenomenon between the zeroth order solution (3.11) acting as a resonant source for the first order correction. For our perturbative problem (3.9), a similar resonant behavior occurs. Assuming S k (z) can be decomposed into a linear combination of Laguerre polynomials L M (z), we have to distinguish two cases in (3.13) For M = N , the solution remains regular and normalizable, (3.14) However, for M = N , which we are going to call the resonant case, the solution is with B an integration constant and Ψ(N, 0, z) a Laguerre function of the second kind (see eq. (D.9)). Since Ψ(N, 0, z) has both a logarithmic divergence at z = 0 and exponential growth for z → ∞, the solution can never be regular and normalizable at the same time. This corresponds to secular behavior because the resonant term can be eliminated by a infinitesimal shift of a in eq. (3.4) since, JHEP04(2020)108

Recurrence formula
The perturbative solution can be obtained systematically by removing resonant terms in the sources order by order, which leads to an algebraic recurrence relation. For this, we assume both δR(z) and a are expanded in ε, where we set f 0 (z) = L N (z). (3.18) Plugging this into eq. (3.4) and expanding in ε, we obtain the perturbation equation for each order in ε, Assuming that f (z) are polynomials for < k, the source term also becomes a polynomial, and hence should be decomposed to the linear combination of the Laguerre polynomials, where M is a finite positive integer. After eliminating L N (z) from the source by using µ k , f k (z) can be expressed as a polynomial as well. And we can decompose the solution at each order into a finite linear combination of Laguerre polynomials The coefficients of the resonant term C k,N correspond to the reparametrizations of ε, and hence can be set to 0. So the problem reduces to determining the coefficients C k,I and µ k at each order. Substituting eq. (3.21) into the source term (3.19), we obtain (3.22) where X K I,J comes from the decomposition of the product of Laguerre polynomials [47],

JHEP04(2020)108
which is written as (3.24) The last line in eq. (3.22) is proportional to the resonant term, and hence should be removed by setting For non-resonant terms, the k-th order coefficients are determined by The coefficient of L N (z) is set to zero C k,N = 0 for k ≥ 1. With these recurrence equations, the perturbation equation can be solved algebraically.

Perturbation solution
To solve the recurrence equation (3.25), we first set Then, we have the solution for k = 1 Especially, the leading order shift in a is given by Here we note that µ 1 alternates in sign with N . For the first values of N , we obtain Using the relation to the Franel number (see appendix. C.1.1), one can show the amplitude of µ grows very rapidly with N , (3.32)

Phase diagram
Given the perturbative solution we can calculate the physical quantities M, J and the value at the origin R 0 = R(0) (which is used as an initial condition in the numerical analysis) perturbatively as follows.
Angular velocity and center thickness. By definition, the angular velocity has the expansion The center thickness is given by Which gives the gradient on the branching point is given by Since µ grows much faster than N , the gradient rapidly approaches to that of the Myers-Perry branch for the larger value of N . For the first few values of N , we obtain At higher order, the center thickness is given by where C k,I is the coefficients of the Laguerre expansion at each order in eq. (3.21). To compare with the numerical result (figure 3), we calculated the formula for (R 0 , Ω)-space up to ε 2 ,

JHEP04(2020)108
Mass and angular momentum. Provided that the perturbation is normalizable, the mass (2.25) and angular momentum (2.26) are easily obtained by where M MP is the mass of the Myers-Perry of the same a and z = L 0 (z) − L 1 (z) is used. Since these integrations take the form of the inner product of the Laguerre polynomials, it is convenient to use the expansion of the perturbative solution into the Laguerre polynomials, where C 0,M = δ M,N for the N -branch and M runs over some finite at each perturbative order k. Up to O ε 2 , one can expand as where we made use of the second order solution (3.27). Putting this into eqs. (3.42) and (3.43), we obtain in which a also should be expanded according to (3.17). We see that the ratio of angular momentum to mass only differs by O ε 2 from the Myers-Perry branch.

Numerical construction
To construct fully non-linear solutions we have to solve numerically the axisymmetric version of the soap bubble equation (2.22) which is a second order nonlinear differential equation for R(r). Since r is a radial coordinate, any physical solution of eq. (3.47) must satisfy the regularity condition R (0) = 0. This leaves the parameter R 0 ≡ R(0) as the initial condition that is needed to integrate the differential equation outwards radially. However, not all values of R 0 will result in physical solutions. In general, as a consequence of the nonlinearity of eq. (3.47), R(r) will become singular at a finite value of r = r s and only a discrete set of initial conditions will allow for solutions that that extend to r → ∞. To find these branches our numerical procedure consists in maximizing the value r s where the singularity appears. Solutions appear as singularities/ peaks of r s as a function of the initial conditions. See appendix A for a more detailed description of the numerical method. For fixed Ω ∈ [0, 1/2], the two branches (stable and unstable) of the MP black hole (2.23) correspond to two such solutions. In terms of the parameter a, the MP solutions describe an ellipse in the (R 0 , Ω) plane as

JHEP04(2020)108
(3.48) Apart from the MP black hole solutions, we find that multiple branches of bumpy solutions extend from every axisymmetric zero-mode. They can be represented in (R 0 , Ω) plane as curves that extend from the Myers-Perry ellipse, as shown in figure 3.
We observe that the bumpy branches fall in two distinct categories. Those branches that arise from even N zero modes, as defined in eq. (3.7), tend to R 0 → −∞ as Ω → 0 (asymptotically like R 0 ∝ − 1 Ω 2 ). This is equivalent to a rapidly decreasing mass density at the rotation axis as one moves along the branch. These bumpy branches connect the MP-branch to families of N − 1 concentric black rings. In figure 4, the mass density profiles m = e R are shown. On the other hand, for the zero modes with odd values of N , we have R 0 → 2, which means that the mass density at the center will closely approach that of a stable MP black hole. These branches will resemble black Saturns with N − 2 rings, as shown in figure 5.
As discussed in [40,41], every axisymmetric branch extends in both directions from the zero mode. This corresponds to the fact that linear perturbations of the Myers-Perry black hole can be added with either a positive or a negative amplitude. According to the convention in [40], branches adding the amplitude of the sign (−1) N +1 on the axis are called (+)-branches, which deform the MP-black hole towards the black rings or black Saturns, and the opposites, (−)-branches, which develop a singularity on the equator of the horizon. It is so far unclear if this (−)-branch connects to some singly spinning black hole solution.
Agreeing with this, we find that the negative side of the branches extends only for a very short interval, after which the allowed solutions cease to exist. This behavior is to JHEP04(2020)108 The (expected) pinching of the necks as we move away from the MP-branch follows a behavior described already in [40]: for multiple rings the pinching starts at the interior necks and later on the outer ones.  Figure 6. Phase diagram for axisymmetric solutions, we show the 10 first appearing branches: ring-branches are shown in purple, and Saturn branches in light-blue. The Myers-Perry and black bar solutions are also plotted by the black and red curves. We do not expect the Saturn branches to terminate, but they become harder to construct for low Ω. Figures 4,5,6 show that the bumpy branches for black rings and black Saturns seem to extend to arbitrary angular momentum 3 without encountering any conical singularities. For a sufficiently high angular momentum, the deformation ends up as multiple lumps/rings barely connected by exponentially thin necks. Figure 6 also shows this in a change of behavior of the curves: all branches show three phases of qualitative behaviors: in the first stage the branches are nearly tangential to the MP-branch. After that in an intermediate stage new (ringlike) blobs start to form until they reach a new asymptotic phase. In this final phase the blobs are practically separated and do barely deform further but the distance between the blobs keeps increasing, the angular momentum behaves asymptotically like For solutions with multiple ripples, we find that at low Ω the radii of ringlike blobs follow two different behaviors. The innermost ring has an approximate radius growing like Ω −1 , while the distance between the following outer rings increases slower than that and we estimate it to be ∼ | log Ω|. The Ω −1 -behavior agrees with the blackfold result for multi-rings if the separations of the rings are much shorter than the ring radius [48]. These observations on the far extended branches lead us to the expectation that our ring/Saturnlike bumpy solutions will be connected via a topology changing transition to the single bumpy rings/Saturns, not directly to multi-rings or higher Saturns. This picture is consistent with the numerical result in D = 6 bumpy Myers-Perrys [40]. JHEP04(2020)108

Multipolar zero modes
We study again perturbations of the MP-black hole, but this time allow for angular dependence of the perturbations Then, the deformation equation becomes where we defined It is convenient to expand the angular dependence as a Fourier series where each radial function is expanded in ε, With the Fourier decomposition, the linear part reduces to the generalized Laguerre equation which admits normalizable solutions for k = m when We also decompose the source terms into Fourier modes (4.10) Here the last line exists only for k > 0. JHEP04(2020)108

Nonlinear perturbation
For higher order perturbations, we proceed in almost the same manner as for the axisymmetric sector. The generalized Laguerre operators L (m) N also show resonant behavior if they are sourced by the corresponding resonant term L (m) N (z), provided N is a non-negative integer. Therefore, for the solution to be regular and normalizable, the resonant term has to be removed from the source for every mode by renormalizing the angular velocity as A new phenomenon we observe, is that some modes can not independently excited at linear order, otherwise the renormalization of the angular velocity becomes impossible. To show this, let us assume to the contrary that we start at linear order only with the zero mode corresponding to Then, this mode acts as a source for the neighboring perturbations f at next-to-leading order, If m is a even, eqs. (4.13) and (4.14) will contain resonant sources. 4 However, since we did not include the corresponding linear order term at leading order, the parameter renormalization cannot absorb the resonant terms. This implies that we are forced to include also the neighboring overtone modes at leading order Repeating the same argument for the new linear solution, one might be concerned that now we need an infinite tower of overtone modes to regularize the secular behavior. However, if N − (i − 1)m/2 < 0 for the i-th overtone, the equation ceases to produce secular behavior as long as the source term is a polynomial. Therefore, given m and N , the linear order solution should be a linear combination of its overtone modes whose overtone number does not exceed 2N/m + 1. 5 JHEP04(2020)108

Recurrence formula
Using the expansion of the spin parameter (4.11) we can derive a recurrence formula for all orders in perturbation theory. Eq. (4.2) can be rewritten as and S (k) (z) given through eq. (4.10). Under the perturbative expansion (4.6), we also expand the source term byS Using an inductive argument, the regular normalizable perturbations are shown to be polynomials to all orders of the perturbation. Therefore, we expand the radial functions at each order by the associated Laguerre polynomials, As discussed in the previous section, the linear order solution should include all the overtone modes with N − im/2 > 0, where η := 2N/m + 1. If m is odd, the even overtones are turned off. Using the reparametrization of ε, we set C (m) Substituting this expansion into eq. (4.18) , the source term can be decomposed into a resonant part and a normalizable part where T (k) p = 0 gives the normalization condition. 6 To extract the resonant term from the source, the following decomposition formula of the product of the associated Laguerre polynomials is used z where the coefficients are written by the integral of the triple product of the associated Laguerre polynomials becomes trivially zero.

JHEP04(2020)108
with This integration can be expressed through Lauricella's generalized hypergeometric functions (see appendix. C.2) [49]. 7 Since the LO-perturbation only contains the fundamental mode m and its overtones, also at NLO only m and its overtones are excited. To eliminate the resonant part in (4.23), we require for i = 0, . . . , η (again, only odd i if m is odd) where the last line only exists for i > 0. Other than the resonant terms, we also obtain the coefficients

Comparison to the numerical results
For later comparison with the numerical result, we derive an expression for the center value of each angular Fourier mode. As in the axisymmetric sector, the center thickness is defined by and for the multipoles, we define 8 (4.29b) 7 An English reference is found, for example, in [50]. 8 Which will serve as initial conditions in the numerical setup (4.70). JHEP04(2020)108

Even multipoles
The analysis for different fundamental modes (N, m) differs in important aspects, so we are going to distinguish several cases in the following. Let us begin with the case m even. As opposed to the axisymmetric modes, the normalization condition (4.27) already gives the coupled equation that determines the linear coefficients and the parameter renormalization, The nonlinear eq. (4.30) is hard to solve in general and we will further distinguish different cases.
Even multipoles with 2N < m. Here the leading order solution consists of only two modes f The normalization condition (4.30) becomes Setting α 1 = 0 immediately reproduces the axisymmetric result (3.30). Therefore assuming α 1 = 0, we obtain and Which has real solutions only if This leads to an upper bound for m (see figure 7). Since the sign of α 1 does not matter, we obtain  where we set α 0 = 1. The right hand side in eq. (4.39) monotonically grows in N , and for N ≥ 3, the bound (4.39) finally starts to exclude all of m > 2N . We will see that a similar bound appears also in the sector N < m ≤ 2N for N ≥ 3. This upper bound does not mean the absence of the higher multipole deformation, but rather implies such deformation should be coupled with the lower companions even in the linear order. For example, (N, m) = (0, 6) can be coupled with (N, m) = (2, 2) (together with (3, 0) and (1, 4)), which is in Lastly, we evaluate the center values and angular velocity in eq. (4.29) up to O (ε), and
Even multipoles with N < m ≤ 2N . Here three modes have to be excited at leading order The normalization condition (4.30) leads to a quadratic constraint for the relative amplitudes where the coefficients are given by and Setting α 1 = 0 immediately reproduces the previous analysis in which m is replaced by 2m. Therefore, we consider α 1 = 0 and (4.47b) reduces to Substituting this to the rest of eq. (4.47), we obtain two quadratic equations 4(I 1 − I 2 )α 0 α 2 + 2I 3 α 2 2 = I 3 α 2 1 .
(4.53) JHEP04(2020)108 I 1 and I 2 (and accordingly I 1 and I 2 ) have to have the same sign for fixed N and m. Thus eq. (4.52) and eq. (4.53) describe an ellipse and a hyperbola in the (α 1 /α 0 , α 2 /α 0 ) plane. The curves always have two (or no) intersections, which are shown to be identical by a constant shift in the angular coordinate φ → φ + π/m. Therefore, we have at most one branch for each (N, m) with N < m ≤ 2N . The radii of the ellipse from eq. (4.52) are proportional to (4.54) The positivity of this value is the necessary condition for the existence of the branch, which gives the upper bound for m ( figure 7). Since the last term in eq. (4.54) decays very quickly in N , the upper bound coincides with that from eq. where we also evaluated the amplitude of the overtone r 2 defined via

Odd multipoles with 2N < m
For odd m the leading order modes do not produce secular behavior at second order, but starting from third order it will also appear in this case. Here the LO-solution consists of a single mode, Iterating eq. (4.28) reveals that there are only even m modes for every odd order in ε, and vice versa. Which results in µ k = 0 for odd k. At third order, the normalization condition (4.27) becomes Different from the even cases, the normalization condition for the simplest odd multipoles does not lead to a bound for m. For the lower sector m ≤ 2N , we will have multiple overtones at linear order, which leads to coupled equations at third order as in the even modes. This may bound m as in the even modes.
In contrast to the case of m even, Ω and R 0 only have even powers of ε appearing in their expansion while R m is odd in ε, This means that odd branches go out from the Myers-Perry branch only in one direction. 9 The leading order corrections can be written as For N = 0 branches, eq. (4.62) gives µ 2 = 0 for any odd m,

Numerical construction
To obtain the fully non-linear multipole solutions numerically, we use a Fourier decomposition corresponding to overtones of a fundamental mode m Plugging this into the stationary master equation (2.22), we obtain a countable set of coupled equations for the fundamental Fourier mode R (m) (r) and its overtones R (n·m) (r) (n = 2, 3, . . . ). From the perturbative analysis, we know that close to the MP-branch higher overtones will only be weakly excited. So we truncate the Fourier series for some n max to obtain a finite dimensional problem. The resulting coupled ODEs can be now solved using the shooting method described in appendix A, where now the space of initial conditions is spanned by the amplitudes of the Fourier modes R (nm) (r) close to the origin, which we will denote as R 0 , R m , R 2m , . . . , R nmaxm . In figure 8, we show examples of branches extracted numerically with only the fundamental Fourier mode, i.e., n max = 1, and compare them to the perturbative result. We checked that the truncation n max = 1 is consistent for the beginning of the branch we show by comparing the results to a higher truncation with n max = 2 and finding good agreement of the results. To extend the branches further overtones should be included.
The inclusion of overtones however makes our numerical procedure much less efficient (see appendix A.3 for details), s.t. at this point we do not find conclusive results for odd multipole branches and even multipole branches corresponding to the opposite sign of the perturbation.
In figure 9, we show representative plots of mass densities for some of the branches. The profiles for even multipoles show a behavior that can be related to the perturbative result that modes of different N and m couple to each other: the black flower branches show mass profiles, which when averaged over the angular direction resemble the corresponding axisymmetric branch that starts at the same branching point, which results in a similar (J /M, Ω)-curve see figure 10.

Deformed black bars: dumbbells and spindles
As already studied in the previous section the large D effective equations allow for stationary solutions without axisymmetry around the rotation axis, the first (and so far only) analytically known solution is the dipolar black bar [27]. Like the other multipolar solutions, the black bar plays an important role in the decay of the ultra-spinning instability of MP-black holes [28,43,46]. At high enough angular momentum the bar gets very elon-JHEP04(2020)108 gated and develops a Gregory-Laflamme type instability. In this section, we are going to study the zero mode configurations corresponding to this instability.
The black bar is best studied in Cartesian coordinates in the co-rotating frame where it can be written as Note that for small Ω the bar becomes very elongated and in the limit Ω → 0 the solution connects to a non-rotating black string along the y-direction.

Co-rotating zero modes
We deform the bar perturbatively via R = R bar (x, y) + δR(x, y), where the deformation δR(x, y) satisfies At linear order, the regular solutions are given by Hermite polynomials where n x , n y are non-negative integers with Together with the constraint −2 ⊥ + −2 = 1, the regular and non-trivial perturbations are available only for n x = 0, n y = 2 ≥ 2. (5.7)

Nonlinear perturbations
Considering the linear result, we can assume only y-dependence even in the non-linear regime. Then, by rescaling  9) where H N is the Hermite operator defined by

JHEP04(2020)108
Given the value of , Ω and ⊥ is written by The corrections beyond the linear order can be derived in the same manner as the bumpy deformation of the Myers-Perry. First, we expand the deformation function by ε If we consider a branch bifurcating from the zero mode 2 = N on the black bar branch, one can set The length of the bar for the deformed branch should be expanded by ε, where the running coefficient µ k is determined so that f k (z) remains to be normalizable at each order. Expanding eq. (5.10) by ε, we obtain where the linear order solution is supposed to be C 0,M = δ M,N . Substituting this, the source term of each order becomes Using the properties of the Hermite polynomials, the source term can be decomposed to the resonant and non-resonant terms,

JHEP04(2020)108
where Q K I,J is given by eq. (C.4). Using C 0,M = δ M,N , the regularizing condition is given by (5.19a) and the non-resonant coefficients, For the resonant term, we simply set Using induction one can show for odd branches that f k (z) has only odd (even) power for the even (odd) order, and µ k vanishes for every odd order. Similarly, for even N , it can be shown that at each order only even powers appear.

Perturbation solution
By solving the recurrence equation with C 0,M = δ M,N , one can obtain the solution to arbitrary order. The result for O ε 2 is where Q N N,N = 0 for the odd N , giving µ 1 = 0 for the odd dumbbells.

Physical quantities
Once, given the deformation δR(z) as the physical quantities are calculated using properties of the Hermite polynomials.

JHEP04(2020)108
Value at the origin. Here we evaluate the center values R 0 = R(0) andR 0 = R (0), which are also used as the boundary condition for the numerical analysis. Due to the mirror symmetry in the even case,R 0 only exists for the odd branches. The center thickness R 0 of the deformed bar is given by With eq. (5.21), we obtain  This shows that one always need to spin up the black hole for the transition to an odd branch.

JHEP04(2020)108
Mass and angular momentum. The mass (2.25) and angular momentum (2.26) can be calculated by √ π e −z 2 exp(δR(z)), (5.36) and (5.37) where M bar = 2πe/Ω is the mass of the bar solution for the given Ω. Due to the orthogonal property of the Hermite polynomials, the integrals in M and J pick up H 0 (z) and H 2 (z) components in exp(δR(z)), respectively. Using the result in the previous section, the ratio of the angular momentum to the mass is given by where we note that Ω should also varies in ε. For the odd branch, both J /M and Ω become a function of ε 2 .

Numerical construction
In order to find fully nonlinear deformations of the black bar, we begin by considering equation (2.21) with the ansatz where we imply that R(y) ≡ R(0, y), and 2 ⊥ is defined by eq. (5.3). With this substitution, we are left with Since y is no longer a radial coordinate, the condition R (0) = 0 is no longer required. We can define R (0) ≡R 0 instead. Allowed solutions must extend regularly both to y → −∞ and y → ∞ simultaneously. If we start the integration from y = 0, the initial conditions are given by R 0 ≡ R(0) andR 0 ≡ R (0), which have to be tuned in order to get allowed solutions.
The branches arising from even N zero modes have a y → −y symmetry, soR 0 = 0. These bars only require R 0 to be tuned, so they can be found in the same way as the axisymmetric solutions. Nonzero values ofR 0 give rise to the branches originating in odd N zero modes. This requires a slightly more involved numerical algorithm, which is described in appendix A.
In figure 11, the first branches of deformed black bars are shown in the (R 0 , Ω) plane. In this case, there is a strong qualitative difference between even and odd N . Odd branches extend only in one direction. This is to be expected, since in this case, reversing the sign of linear perturbations is equivalent to the gauge change φ → φ + π. Surprisingly, for odd N branches, Ω increases as we move away from the zero modes, and these branches are also very short.  Figure 11. Branches of black bar deformations on the (R 0 , Ω) plane. The right plot is a close-up showing good agreement with the analytic expansions (orange) and also zooms in on the short branches. Different tones of green are being used for different branches for the sake of clarity. Even N branches result in the bar breaking apart in N/2 separated blobs. In (R 0 , Ω) plane, they behave in a way that is qualitatively similar to the axisymmetric case, and can therefore be classified in two types. If N is a multiple of 4, R 0 → 0 and the mass density approaches zero at the origin. If N is even but not a multiple of 4, then one of the blobs stays at the origin, with R 0 → 2. The profiles of the first two symmetric bars (N = 4, 6) are depicted in figure 12.

JHEP04(2020)108
Similar to the axisymmetric branches, even N branches can be extended far away from the black bar to the arbitrarily small Ω, in which the mass profile approaches to the multiple blobs located in the almost equal interval. Again, we observe these intervals grow very slowly at the same logarithmic rate as that of ring-like blobs in the axisymmetric branches. Therefore, one can expect these branches finally would pinch off to the array of binary black holes.

Effects of adding charge
Following the approach of [29] and as already described in section 2.1 we can easily construct the (non-extremal) charged solution corresponding to every uncharged solution. According to eq. (2.20) for a given charge parameter q = Q M and given Ω, the charged solution has the profile of an uncharged solution with rotation parameter The (J /M, Ω) phase diagrams for |Q| > 0 are thus the same diagrams as in the uncharged case with a rescaling of the Ω-axis by the factor 1 − 2q 2 −1/4 . Accordingly the bumpy branches will appear at the same J /M but at a lower Ω. As shown in the previous sections lower values of Ω correspond to more elongated/ further separated blobs, i.e., adding charge to the black holes leads to stronger deformations. This intuitively can be understood as charge repulsion deforming the horizon.

JHEP04(2020)108 7 Discussion
In this paper we have demonstrated that the hydro-elastic equations [18] contain a whole new class of 'rippled' stationary solutions, besides the already known black branes, their non-uniform deformations [21] and the non-deformed spinning localized black holes [27].
We have constructed solutions that branch off from the singly spinning Myers-Perry solution directly or indirectly via the black bar branch, which has been already identified in [27]. We found both axisymmetric and non-axisymmetric solutions, and only the former ones can remain stationary at finite D, since non-axisymmetric solutions will radiate gravitational waves. However, with increasing number of dimension the emission of gravitational waves becomes weaker, which will allow the non-axisymmetric solutions to be long-lived.
The axisymmetric solutions described in this paper, we have identified as ring-like and Saturn-like bumpy black holes, or black ripples in short. They bifurcate from the axisymmetric zero modes of Myers-Perry in the ultra-spinning regime. As in the numerical studies in finite dimensions [40,41], we found that all branches extend in two directions: either with a positive or a negative amplitude of the deformation. The direction that increases the angular velocity leads to a very short branch, the other direction extends indefinitely at large D. This implies that the former directions lead to singular solutions, as observed in previous numerical constructions [40,41].
Multipolar deformations can not be stationary in a fixed number of dimensions, but are indicative of ultraspinning instabilities of the Myers-Perry black hole. In high enough dimension they correspond to long-lived transient objects. We generically call them black flowers, the simplest case among them is the black bar and it has an analytic solution.
The black bar also has an infinite number of co-rotating zero modes, from which deformed branches develop: the dumbbells and the spindles. We classify the deformed bars by the parity of their zero mode as odd and even. Similarly to the ripples, the even branches go out in two directions. In the spin-down direction, the deformation grows a dumbbell-like profile with a distinct number of blobs for each branch, and hence we call them dumbbells. In the opposite direction, we could find only very short branches which we call spindles. Odd branches turned out very short as well. Odd branches and spindles correspond to solutions with increased angular velocity. One might expect that both the spindles and the odd branches end up forming a singularity.
It is very suggestive that the spindle branches correspond to the solutions that develop sharp pointy endings, as observed dynamically in [28,46]. These sharp endings of the deformed bar would be possibly affected by the Gregory-Laflamme instability, presenting a large number of zero modes close to the end of the short branch. The sharpened tips could, in principle, pinch off producing detached small black holes.
This process of a black hole developing long arms that end up pinching off has indeed been observed in [28,46], not only for the spindles but also for higher multipole deformations. We find it likely that these dynamical solutions would correspond to the short branches described above, i.e., those resulting from exciting the zero modes in the direction with increasing Ω. This would apply both to the spindle solutions and to multipolar deformations leading to multiple arms. This conjecture is supported by the fact that short JHEP04(2020)108 branches go in the direction of decreasing J /M, which should be favored in finite D simulations since gravitational radiation tends to decrease the angular momentum to mass ratio of the evolving object.
The method used to identify axisymmetric solutions should be exhaustive, and thus we do not expect the ripple branches to have their own secondary axisymmetric zero modes. We expect, on the other hand, that the axisymmetric solutions will become unstable to multipolar deformations. An indication of a ring-like ripple breaking apart into four black holes via an m = 4 deformation was already found at large D in [28]. Interestingly, black rings share the same kind of instabilities and subsequent pinch-offs [51][52][53][54]. Such instabilities would begin at zero modes along the branches of ripples. This fact leaves open the possibility of the 'long' multipolar branches actually merging with the ripple branches at these zero modes. No conclusive results have been obtained about this intriguing possibility in this paper.
We have found no evidence that the long multipolar branches have bifurcations. This possibility could be analyzed in future work, possibly with an improved numerical setup. The dumbbell branches end as an array of separated black holes and thus seem unlikely to have further zero modes.
The nature of the boundary conditions that are imposed in the blob formalism, together with the nonlinearity of the large D effective equations, leads to a remarkably challenging numerical problem. Ordinary relaxation and spectral techniques have not been shown to give reliable results so far. This fact is probably due to the requirement of imposing boundary conditions at spatial infinity, together with the equations of motion being numerically bad-behaved as r → ∞. Additionally, the equations are nonlinear, which rules out direct eigenvalue-finding standard algorithms. Fortunately, the shooting approach used in this paper, which consists in identifying sharp peaks in the radius where the numerical solution becomes singular, seems to be enough to find the right solutions. It is remarkable that this technique works even though the numerical method is usually able to integrate only to a finite value of r. Axisymmetric solutions are easily found this way. For the case of multipolar deformations, one encounters a multidimensional shooting problem with a scalar-valued output function (r s ), which becomes increasingly difficult as one increases the number of overtones. For this reason, an alternative method, possibly based on relaxation techniques, would be desirable in the future.
In the formalism employed here, the effect of the charge is simply incorporated in the effective angular velocity Ω q = Ω/(1 − 2q 2 ) 1/4 as in [29]. Therefore, with a given value of charge and Ω, the corresponding charged solution is immediately obtained from the uncharged one. Due to the factor 1 − 2q 2 −1/4 , the charged deformed branches will appear for the same J /M but for a lower Ω, which corresponds to more elongated/further separated blobs. This can be interpreted as the effect of the charge repulsion. Since all the analysis is written in terms of Ω q , one can take the extremal limit q 2 → 1/2 of all branches, keeping Ω q finite, resulting in a smooth limit, that leads to rather strange deformed 'extremal' branches, both with and without rotation. The proper large D limit of extremal horizons is however yet unclear, and a more careful analysis seems appropriate.

JHEP04(2020)108
Fate of far extended branches. All 'long' branches (corresponding to bulging deformations) extend far away from the original bifurcating points in the phase space, where they develop broad thin regions. Currently, very little is known about how to interpret these nearly zero thickness regions in the large D effective theory. In the case of spherical black holes the thickness falls off towards infinity as a Gaussian profile, which might be interpreted as the round tip of the black hole. Therefore, if the deformation develops a thin neck between blobs, and its size grows infinitely large, one can expect such deformation to end up as a pinch off of the horizon at finite D. This would correspond to a topology-changing transition.
We found that the ripple branches develop such long thin necks connecting Gaussianshaped ring blobs (with a central blob in the case of Saturns) at their final stages of deformation. Particularly, we observed that the separation process involves two distinct length scales. From the numerical solutions, we could easily estimate that the radii of ring blobs grow like Ω −1 as Ω → 0. The same behavior has been derived in the blackfold approach [4,48], which might imply that the blackfold approximation becomes already accurate in the pinch off phase, due to the localization of gravity at large D. Another scaling is that of the intervals between ring blobs, which are estimated as ∼ | log Ω|. Due to the hierarchy in these two scales, we expect the first pinch off to occur always on the axis, indicating a first topology change to a bumpy black ring/Saturn, before transitioning to the multi-rings/Saturns, as observed in the (+) 3 -branch of D = 6 bumpy black holes [40].
Dumbbell branches also extend far away from the black bar to arbitrarily small Ω, where the mass profile approaches that of multiple evenly spaced blobs. As opposed to the ripples, dumbbells show only a single scaling, which has the same logarithmic growth as the intervals between the ring blobs in the case of ripples. Therefore, one can expect that these branches would finally pinch off to multiple black holes. 10 Finite D effects. The blob coordinate is supposed to be identified as the small patch of the √ D-amplified entire coordinate. 11 Therefore, the blob approximation will break down if the length of the thin neck reaches ∼ √ D, when the 1/D corrections are included. This breakdown will give some information on the transition in phase space. For example, the pinch off from the ripples to black rings or Saturns will take place at Ω ∼ 1/ √ D. Actually, black rings are already constructed by using the large D effective theory approach in the same scaling [15,17]. This implies that one can use the effective theory result as the global setup to solve the local topology-change. For other logarithmic scalings ∼ | log Ω|, the break down will occur at much smaller spin Ω ∼ e −D . In the black string analysis, a similar type of breakdown is already seen after including 1/D corrections [21]. The black hole entropy is another important quantity to evaluate the stability of the solutions. Since the mass and entropy become degenerate at D → ∞, we would need to know the next-toleading order terms in 1/D expansion to calculate the entropy difference for a given mass.
Blob-blob interactions. For the ripples and dumbbells, we observed a universal scaling of the blob distance as | log Ω| at Ω → 0, implying an effective interaction between the JHEP04(2020)108 blobs (or ring-like blobs). This indicates the possibility to reconstruct the large D effective theory as a particle-like (or soliton-like) effective description of blobs weakly interacting via very thin necks. This possibility will be pursued elsewhere.
The origin of this logarithmic dependence, though very naively, might be understood as a force balance between the centrifugal force and the attraction between the blobs at large D. Assuming a black hole of radius r H and an orbiting particle, the gravitational force is approximated as (r H /r) D and the centrifugal force as Ω 2 r. The equilibrium is accomplished by r/r H ∼ 1 − 2D −1 log Ω. Therefore, the particle orbit exists very close to the horizon ∼ | log Ω|/D. This introduces the | log Ω| scaling in the near horizon region. Curiously, if we assume two adjacent black holes with the same mass, the equilibrium condition would be modified to r/r H ∼ 2 − 2D −1 log(e D/2 Ω) with e D/2 Ω = O (1) or | log Ω| ∼ D. This coincides with the value at which the neck length between blobs reaches | log Ω| ∼ √ D and the blob approximation breaks down.
Towards the topology change. The topology-changing transition at large D is described by the conifold metric which solves the Ricci flow equation [30]. Especially, the black string/black hole transition is completely solved by the King-Rosenau (KR) solution for the 2D Ricci flow. Some of the topology-changing transitions (Saturn-like ripples, dumbbells) can be reduced to the 2D Ricci flow problem in the co-rotating frame, since the transition occurs in a very narrow region. Hence, they should also be solved by the KR solution, due to the rigidity in 2D compact ancient flow [55]. 12 For the transition between ring-like ripples and black rings, we need a better understanding of the 3D Ricci flow.
Here we should note that, in the case of the black string/black hole transition, one just has to give the global configuration (such as the black hole (blob) radius and the compactification scale) as boundary conditions for the conifold metric, without considering the force balance condition. Now, for example, if we consider the transition between a dumbbell and binary black hole, we also have the rotation Ω, which will not appear in the large D conifold analysis after switching to the co-rotating frame. To relate Ω with the mass and separation, one needs to find the proper force balance condition at large D, as roughly estimated in the previous paragraph.
In the current formalism, we could only follow the (−)-ripple branches for a very short range. These (−)-branches are shown to develop a single-sided conical horizon on the equator when they approach the end of their branch [40]. Therefore, it should also be possible to study the ending phase of (−)-branches using the large D conifold metric and Ricci flow methods. Different from the usual pinch off problem, one may have to find the non-compact Ricci flow solution, in which only one side is the horizon.

A.1 Axisymmetric sector
Stationary axisymmetric black holes are regular solutions of eq. (3.47) that extend from 0 to r → ∞. Due to singular point at r = ∞ from the rotation term it is particularly difficult to use of spectral and relaxation methods. For this reason, the approach used in this paper is essentially a shooting method. By regularity at the origin the ODE can be generally integrated radially outwards with the initial conditions R(0) = R 0 and R (0) = 0. The numerical solution will generally become singular at some finite r = r s . In figure 14, the values of r s are shown as a function of the initial condition parameter R 0 , interestingly the appearance of singularities is (semi-) continuous in the space of initial conditions which makes it possible to look for singularities/ peaks where the solution extends to infinity. These peaks correspond to (approximate) locations of the allowed solutions. For this purpose, the (R 0 , Ω) plane is not a very suitable representation. This is so because the branches of solutions become very closely packed at low Ω, while the ring-like branches reach very large negative values. A numerical algorithm intended to find all these peaks with a high precision needs therefore an extremely high dynamic range of detection in R 0 , so it can both find widely spread peaks and resolve extremely close packed ones. This is solved by introducing the coordinates (α, β) as Ω = e β 2 sech α , R 0 = 2 − e α+β sech α . These coordinates both range from −∞ to ∞, and cover the (−∞, 2) × (0, ∞) region in (R 0 , Ω) plane. They are analytically invertible as In these new coordinates, the Myers-Perry black holes lay on the vertical axis (β = 0), with the Schwarzschild black hole corresponding to β = 0, α → −∞ (see figure 15). The ripple solutions become now much more suitable to be found numerically. In particular, the ring-like branches can be parametrized by β, and the Saturn-like by a polar angle θ such that α = ρ cos θ and β = −ρ sin θ.
When a branch ends, as for the negative amplitude modes, the peak that represents it becomes a local maximum, with no divergence whatsoever. This requires us to define a criterion for a local maximum to be considered a proper peak, or a vanishing peak. The criterion that has been taken for a peak to be valid is max {r s (α, β) − r s (α + δα, β), r s (α, β) − r s (α − δα, β)} > ∆ , where δα = 0.01 and ∆ = 3. When extracting the angular momenta of the solutions, it is also important to take into account that numerical error may result in extra (unphysical) oscillations of the R(r) profiles. These oscillations appear as additional bumps, or fake rings. These have to be properly removed before the angular momentum integration, since they could add an erroneous contribution to the integration result.

A.2 Black bar deformations
Deformations with even values of N are found in a way which is completely analogous to the axisymmetric case. In this case it is convenient to reparameterize the (R 0 , Ω) by the coordinates (γ, δ), Odd deformations of bars are described by solutions of eq. (5.40) that have a nonzero value ofR 0 = R (0). This increases the complexity of the problem, since it now requires to JHEP04(2020)108 tune both R 0 andR 0 in order to get a solution that extends to infinity both for the negative and positive sides of the y axis. This complication can be partially circumvented by noticing that, for the deformed black bars, the change y → −y is equivalent toR 0 → −R 0 . This means that, if (Ω, R 0 ,R 0 ) gives an allowed solution, then so does (Ω, R 0 , −R 0 ). This fact allows the right values of R 0 to be found by requiring the peaks in r s (Ω, R 0 ,R 0 ) to be located at opposite values ofR 0 . This is done by the secant root-finding method in a few iterations. Again, vanishing peaks and fake blobs are discarded in a similar way as in the axisymmetric case.

A.3 Multipole deformations
By using the ansatz (4.70) truncated at some Fourier mode cos(n max mφ), we obtain a set of n max + 1 coupled equations for the functions R (nm) (r). These equations, by imposing the regularity condition R (nm) (0) = 0 ∀n, can be solved by specifying the values of the radial functions at the origin. The problem reduces then to finding peaks in the singular radius r s (Ω, R 0 , R m , R 2m , . . . , R nmaxm ). Identifying peaks on a function with more than one variable is in general not an easy task, especially if there is no straightforward way of reducing the problem to one variable (as in the case of odd deformations of the black bar). For this reason, in this article we restrict ourselves to the fundamental Fourier mode, i.e., we maximize r s (Ω, R 0 , R m ). We use the Mathematica function NMaximize to identify the peak by incrementing Ω in small steps, and constraining the search in a small region around the result of the previous step.
Even with this method, the values of the R 0 , R m still are affected by small fluctuations (which are likely due to numerical error) around the branch. We correct this by subsampling the data points.

B Matching to the entire hemisphere
In general, blob solutions are thought to be identified as a polar cap of the compact black holes, in which the polar angle is stretched by √ D to match with the radial coordinate in the blob [27], r = √ Dθ. (B.1)

(B.2)
Here we show that this match is also consistent beyond the linear level, despite the increase in the degree of the polynomials in the higher perturbation order. The degree of each perturbation solution can be estimated from the recurrence formula ( where we rescaled the perturbation parameter byε = D 2N ε, so that the linear order remains finite at D → ∞. Therefore, the linear order match (B.2) turns out to be correct even up to the nonlinear order, and all the nonlinear perturbation will be matched with the subleading correction in 1/D, where the coefficients are given by It is worth noting that the coefficients in the above two formula become non-zero only if (I, J, K) satisfy the trigonometric inequality: any of the three cannot exceed the sum of the rest two.

C.1.1 Relation to Franel number
Interestingly, the renormalization coefficient µ in eq. (3.30), is related to the so called Franel number, which is known in combinatorics and number theory, Due to the identity, Using the large N approximation for the binomial coefficients, we can show the rapid growth in this number with respect to N , (C.8)

D Derivative of Laguerre functions with respect to the parameter
In this section, we study the infinitesimal parameter shift in the generalized Laguerre functions from the Laguerre polynomials.

D.1 Confluent hypergeometric equation
We start from reviewing the confluent hypergeometric equation, A solution is given by Kummer's confluent hypergeometric series where (a) k := Γ(a + k)/Γ(a) is the Pochhammer symbol. If b is not positive integer, the other solution is given by