Renormalization group analysis of superradiant growth of self-interacting axion cloud

There are strong interests in considering ultra-light scalar fields (especially axion) around a rapidly rotating black hole because of the possibility of observing gravitational waves from axion condensate (axion cloud) around black holes. Motivated by this consideration, we propose a new method to study the dynamics of an ultra-light scalar field with self-interaction around a rapidly rotating black hole, which uses the dynamical renormalization group method. We find that for relativistic clouds, saturation of the superradiant instability by the scattering of the axion due to the self-interaction does not occur in the weakly non-linear regime when we consider the adiabatic growth of the cloud from a single superradiant mode. This may suggest that for relativistic axion clouds, an explosive phenomenon called the Bosenova may inevitably happen, at least once in its evolutionary history.


Introduction
String Theory provides us with an interesting scenario that our universe contains plenty of ultralight axion-like particles, which is called axiverse [1]. Detecting this type of particles will provide us with a lot of information about physics beyond the Standard Model, or even a clue to the string theory. In Ref. [1], the authors proposed that these axion-like particles may cause observable astrophysical or cosmological phenomena, such as the polarization of CMB or the step-like features in the matter power spectrum. In this paper, we will concentrate on the other possible phenomena involving a black hole.
We consider a massive scalar field around a Kerr black hole. As is well known, a rotating black hole possesses the ergoregion. If the negative energy field excitations allowed in the ergoregion fall into the black hole, energy conservation implies that the external field gains energy from the black hole. This is the superradiance, i.e., energy extraction from a black hole via waves, analogous to the famous Penrose process [2]. Now recall that the scalar field is massive. The mass makes the scalar field bounded around the black hole. This bound state continues to gain energy by the superradiance and keeps growing, which is called superradiant instability. Because we mainly concern QCD axion or axion-like particles from string theory in this paper, we refer to such an unstable bound state as an axion cloud [1,3].
The growth rate of an axion cloud was calculated approximately in Ref. [4,5] and numerically in Ref. [6]. These works show that the largest growth rate is ω I ∼ 10 −7 M −1 in units of G = c = = 1, where M is the mass of the black hole. ω I can be understood as the imaginary part of the scalar field oscillation frequency ω, with µ being the mass of the axion field. The largest growth rate is achieved when the effective gravitational coupling governed by α g = µM is ∼ 1. For example, if we consider a solar mass black hole, µ ∼ 10 −10 eV gives the largest growth rate. For the superradiant instability to be efficient within the age of the universe, the mass of the scalar field has to be in the range 10 −20 ∼ 10 −10 eV. Axion-like particles motivated by string theory naturally have a mass in this range [7].
If ultralight particles beyond Standard Model such as axion exist, an axion cloud forms around a black hole and interesting phenomena happen. Similar to the photon emission from a hydrogen atom, the axion cloud can emit gravitational waves by the transition between different energy levels. Also, axion pair annihilation can result in gravitational wave emission [3,8]. As the axion cloud grows, the black hole loses its spin and energy. As a result, a characteristic feature in black hole spin and mass distribution is expected [3,9].
Including the self-interaction of axion can cause dramatic effects on the axion cloud. The self-interaction of the axion is typically attractive, and thus when the axion cloud gets heavy enough, this attractive force would lead to a collapse of the axion cloud. This collapse is called bosenova and a burst of gravitational waves is expected during the collapse [3,10,11]. These phenomena have attracted great interest owing to the possibility of their detection [12,13]. To prepare for the future detection of axion clouds through observations of gravitational waves or astronomical observations of black holes, the precise analysis of the evolution of axion clouds is important.
To do so, we have to correctly take into account the self-interaction of the axion. Only a few works have discussed the self-interaction incisively. Numerical simulation [10,11] and analytic treatment within the non-relativistic approximation [3] were made to show that the self-interaction causes the bosenova. However, in the numerical simulations, there is a difficulty due to the huge hierarchy of time scales between the growth rate ω I and the dynamical time scale ω R . In Refs. [10,11], they used a large amplitude axion cloud as the initial condition to overcome this problem. However, in a realistic scenario, an axion cloud starts to evolve with a small amplitude seeded by fluctuation. To answer the question of whether bosenova happens or not in a realistic case, adopting a particular configuration with a large amplitude as the initial condition would be difficult to justify. Also, the most interesting case, in which the growth rate is maximized, is in the relativistic regime, so that non-relativistic approximation would not be satisfactory. Apart from bosenova, the self-interaction has one more important effect: the energy loss of the axion cloud due to the scattering of axion. Using the order estimate in Ref. [3], the energy loss due to self-interaction is typically more important than that due to the gravitational wave emission when the axion cloud grows by superradiant instability, but this effect is usually ignored in the literature. Our naive expectation is that when the cloud gets denser, this energy loss may balance the superradiant instability. This is because the time scale of the energy gain scales as ∝ |A| 2 , while the energy loss scales as ∝ |A| 6 , where |A| is the amplitude of the cloud. Thus, superradiant growth of the cloud should terminate in the weakly non-linear regime, and no explosive phenomena induced by the strong non-linear effects should happen.

2/31
Our main focus of this paper is to clarify the effects of self-interaction on the dynamics of an axion cloud avoiding the drawbacks in numerical simulation and non-relativistic approximation. We use perturbation theory to tackle the problem of axion self-interaction. When we apply the perturbation theory to a non-linear problem, we often encounter a secular term that destroys the validity of the perturbation theory. We encounter the same difficulty in the problem that we are concerned with. To overcome this difficulty, we use the renormalization group (RG) method for differential equations [14]. Using the RG method, we obtain the evolution equation that describes the long term behavior of the axion cloud. Our main conclusion is that the axion self-interaction would not terminate the superradiant instability by the energy loss within the validity range of our approximation, contrary to our naive initial expectation. What we will find is that the non-linearity accelerates the superradiant instability of the axion cloud because of the attractive nature of self-interaction.
This paper is organized as follows. In Sec. 2, we review superradiant instability and how a growing cloud around a rotating black hole is formed. In Sec. 3 we derive the evolution equation of the axion cloud by RG method. In Sec. 4 we analyze the evolution of the axion cloud by the equations derived in Sec. 3. Finally, we summarize our result and discuss the implication of our result in Sec. 5. Since the RG method is important for our analysis, we provide two simple examples of the RG method for differential equations in App.A that may help to understand the application of the RG method to the current problem. In this paper, we use G = c = = 1 unit, unless otherwise stated.

Superradiant instability
In this section, we review the superradiant instability of an axion field around a Kerr black hole. For a recent review on this topic, see Ref. [2].
We consider the following action for an axion field φ: Here, F a is the decay contant of axion and g µν is the metric of the Kerr spacetime, specified by with ∆ = r 2 − 2M r + a 2 , ρ 2 = r 2 + a 2 cos 2 θ .
There are two horizons in the Kerr spacetime specified by the solutions of ∆ = 0, which are As we are interested in axions, we adopt the cosine type potential induced by the quantum effect [15]. From the action (1), the equation of motion for the axion field is Here, g is the Laplacian on Kerr metric. 3/31 In this section, we solve the linearized version of Eq. (5), which is In Kerr spacetime, separation of variables is possible [16]. We take as an ansatz for the axion field, to get where Λ lm (ω) is the separation constant and has to be calculated numerically. Also, we normalize the angular solution as Following Refs. [4] and [6], we can show that this system develops instability, i.e., ω I = Im [ω] > 0. In Kerr spacetime, there is a timelike Killing vector ξ µ ≡ (∂ t ) µ . Using the conservation of T µ ν and the Killing equation, we obtain Here, we introduce the ingoing Kerr coordinates (t, r, θ,φ) defined by Integrating Eq. (12) overt =constant surface Σ, we obtain where n µ = −δ 0 µ and k µ = −δ 1 µ . By substituting the ansatz (7) into Eq. (14), we find where Ω H is the angular velocity of the outer horizon given by Ω H ≡ a/2M r + . Note that contribution to the right hand side from the boundary at r = ∞ disappears because we are considering a bound state, which vanishes exponentially there. We refer to the real part of ω as ω R . Equation (15) shows that the instability occurs when holds. Because the condition for the instability to occur is the same as the superradiance condition ω < mΩ H if ω I ω R , we will denote this unstable mode as the superradiant mode.

4/31
The superradiant mode grows exponentially around the black hole to develop a condensate of the axion field, which is called an axion cloud. There are several papers on the evaluation of ω I . In Ref. [4], the WKB method was used with the assumption µM 1 and in Ref. [5], the matched asymptotic expansion method was used with µM 1. These works show that ω I is the largest in the region µM ∼ 1. In this regime, ω I has to be calculated numerically [6,17]. For this purpose continued fraction method is used, which shows that ω I takes the maximum value ω I /M ∼ 1.5 × 10 −7 at l = m = 1, a/M ∼ 1, µM ∼ 0.42. We confirmed that our new code using the continued fraction methods to evaluate ω I consistently reproduces the results of Ref. [6].

Renormalization group analysis of an axion cloud
In this section, we formulate how to apply the RG method to the time evolution of an axion cloud. We are interested in whether the dissipative effect of self-interaction terminates the instability or not. Therefore, we will develop our formulation to second order perturbation theory, where the dissipative effect first appears. We assume that one can neglect the gravitational perturbation caused by the axion field, for simplicity.

Derivation of the Evolution equation
To consider the effect of the axion self-interaction, we solve Eq. (5) using perturbation theory. We assume that the amplitude of axion field is small enough, which allows us to approximate Eq. (5) as where λ = µ 2 /6F 2 a . This means that we have approximated the potential of the axion as The negative sign of the term proportional to φ 4 clearly indicates that the self-interaction of the axion is attractive. We solve Eq. (17) by expanding the solution in λ as Substituting Eq. (19) to Eq. (17), we obtain, Since we are interested in the dynamical evolution of an axion cloud, we take the fastest growing superradiant mode as the zeroth order solution: where S l0m0ω0 (θ) and R l0m0ω0 (r) are the solutions of Eqs. (8) and (9) with ω 0,I > 0. The label (l 0 , m 0 , ω 0 ) specifies the fastest growing superradiant mode that we are considering. A(t 0 ) is the amplitude of the cloud at t = t 0 , and c.c. denotes the complex conjugate. Here, 5/31 we introduce an arbitrary reference time t 0 , which is used in the RG method later. We first demonstrate the use of the RG method with the first order perturbation equation (21) in Sec. 3.1.1, and then proceed to the second order in Sec. 3.1.2.

First order perturbation
Using the retarded Green's function G ret , the first order equation (21) is formally solved as Here, the Green's function is defined to satisfy Separation of variables on Kerr space-time allows spectral decomposition of the Green's function as where and the function W lm (ω) is the Wronskian of R r+ and R ∞+ defined by The integration contour C is as shown in Fig. 1. As the integral over ω in Eq. (26) picks up poles of the integrand, which correspond to zeros of the Wronskian, we need to choose the integration contour C to pass above all the poles so that the Green's function satisfies the retarded boundary condition. We introduced R r+ and R ∞± as the solutions of Eq. (9) satisfying the boundary conditions where the symbol" * " denotes the complex conjugate of the variable and r * is the tortoise coordinate defined by dr * = (r 2 + a 2 )∆ −1 dr. The boundary condition for R r+ (see Eq. (29a)) implies that there exist only ingoing waves at the outer horizon and that for R ∞± (see Eq. (29b)) means purely outgoing/ingoing waves at infinity. The asymptotic forms at r → ∞ Fig. 1: Line C shows the contour of integration on the complex plane of the frequency ω. ω SR and ω QNM denote the frequency of the superradiant mode and that of the quasi-normal mode, respectively. The frequency Ω denotes the frequency of the source.
(see Eqs. (29a) and (29b)) determine the value of Wronskian as At the frequency of the superradiant mode ω = ω 0 , W lm (ω) vanishes. This is because superradiant modes satisfy ingoing boundary condition at the outer horizon and decaying boundary condition at infinity (we chose Im ω 2 − µ 2 > 0 branch). Substituting Eq. (26) into Eq. (24), we obtain Integration over ϕ gives δ mm . In the present case, the source φ 3 (0) (see Eq. (23) for the definition of φ (0) ) contains only m = m 0 and m = 3m 0 components. Therefore, the summation over m only picks up m = m 0 and m = 3m 0 . Notice that if we take the integration contours carefully so that both t and ω integrals converge, we can exchange their order. We first perform the integration with respect to t . When Im [Ω − ω] > 0, which is satisfied by taking the countour C as shown in Fig. 1, the integration over t gives Next, we integrate over ω. The form of this integration is At large |ω|, the integrand falls faster than |ω| −2 , if the mode functions R lmω and S lmω do not grow faster than |ω| 0 as |ω| → ∞. This |ω| −2 factor partly comes from 1/ω factor in the Wronskian. Now, we modify the integration contour C to a large semi-circle in the upper half plane and a small circle around the pole at ω = Ω. As we have stated above, the integrand falls 7/31 off as fast as |ω| −2 , the integral at a large semicircle vanishes. Therefore, all contributions to the integral come from the pole at ω = Ω.
In Eq. (31) the source contains frequencies with ω = ω 0 + 2iω 0,I and 3ω 0 , which are the only frequencies that appear in the first order particular solution: Here, f It is easy to see from Eq. (34) that the perturbative solution grows exponentially so that the solution would break the assumption that the amplitude of the perturbation is small. We use the RG method to avoid the breakdown of the perturbative expansion (see appendix A for review on the RG method).
To apply the RG method, we identify the term that diverges in ω 0,I → 0 limit. Note that we are considering the superradiant mode which satisfies ω 0,I /ω 0,R 1. Recall that the Wronskian W lm (ω) contained in the Green's function (26) has a zero at (ω, l, m) = (ω 0 , l 0 , m 0 ). From the expression (36), we see that the Green's function in f (1) l0m0ω0+2iω0,I (r) contains the factor 1/W l0m0 (ω 0 + 2iω 0,I ), which diverges for ω 0,I → 0. Taking the leading term in the Taylor expansion of the Wronskian around ω = ω 0 , we can extract the leading term in f (1) l0m0ω0+2iω0,I (r) in the limit ω 0,I → 0 as which we call a divergent term, although, strictly speaking, it is not divergent since ω 0,I is small but not vanishing. We can claim that the term of O ω −1 l0m0ω0 is purely imaginary. In the limit ω 0,I → 0, the equations for the radial and angular mode functions (9) and (8) are both real, and hence R l0m0ω0 and S l0m0ω0 can be chosen to be real. Also, α ω0 and A out , which are calculated through the mode functions, are real in this limit. Thus, the leading term 8/31 in C (1) l0m0ω0 is real except for the factor ω 2 0 − µ 2 , which is purely imaginary in the limit, ω 0,I → 0, because a supperradiant state is a bound state satisfying ω 0,R < µ. To summarize, we have Using the RG method, we cancel this divergence of Im C (1) l0m0ω0 by adding the homogeneous solution with an appropriate amplitude. This procedure is the same as inserting the counterterms in the renormalization in quantum field theory (QFT). After this procedure, the first order perturbative solution is given by We explicitly added δC (1) (= O(ω 0 0,I )) to express the arbitrariness in identifying the divergent part (see App.B for the relation between the choice of δC (1) and the definition of the amplitude A(t)). This is the same as choosing the scheme of renormalization in QFT. Now, we demand that the expression (41) satisfies the RG equation: We substitute Eq. (41) into the RG equation, to obtain the evolution equation for the amplitude as Redefining the amplitude to incorporate the exponentially growing factor contained in the zeroth order solution φ (0) will be more convenient. Namely, we define A(t) ≡ A(t)e ω0,I t , and then A satisfies the following equation, We decompose this equation into the amplitude and phase part by writing A in the form Substituting the solution of Eq. (45) to Eq. (41), and setting t 0 = t, we finally obtain the renormalized first order perturbative solution as where is the non-divergent part of f (1) l0m0ω0+2iω0,I (r). Two types of modes appear in Eq. (48). The first term in the parentheses has the frequency 3ω 0,R > µ, which means that the mode is unbounded and can dissipate to infinity. The second and third terms have the frequency ω 0,R < µ, and thus these are bounded and cannot dissipate energy to infinity. Also, an important property common to all modes is that they satisfy superradiant condition ω < mΩ H . Therefore, dissipating the energy of the cloud back to the black hole does not occur in our setup describing the adiabatic evolution starting with a single dominant superradiant mode.

Second order perturbation
Now, we proceed to the second order analysis to incorporate the dissipative effect from scattering via self-interaction. This calculation is almost in parallel to the first order one: we first solve Eq. (22) formally, then identify the divergent part in the formal solution, and finally apply the RG method to eliminate this divergence.
A formal solution of Eq. (22) is given by As we have seen before, the part that diverges in the limit ω 0,I → 0 originates from the Wronskian W l0m0 (ω 0 ). Using the spectral representation (26), we see that φ (2) contains the Wronskian in the following part of the solution: l m0ω0+2iω0,I (r ) 10/31 From this expression we identify the divergent part as where C (2) is defined by Same as first order solution (see Eq. (48) and discussion below), there are two contributions to C (2) . One is due to ω > µ modes, which dissipate energy to infinity and the other is due to ω < µ modes, which decay at infinity and conserve energy.
There is another type of source which produces divergence in the second order solution. This source comes from the homogeneous solution that we added to eliminate the divergent part in the first order solution. Using the Green's function, the contribution to the second order solution from this type of source is We identify the divergent part in this expression as − 9A|A| 4 e 2ω0,I t0 e 2ω0,I t 2C (1) l0m0ω0C We choose the initial condition for φ (2) at t = t 0 to eliminate the divergence derived above. This turns out to be adding the homogeneous solution to φ (2) . Here, δC l0m0ω0 represents the arbitrariness in the choice of the non-divergent part. Now, we can derive the amplitude equation by imposing the RG equation as in the first order case. Taking care of the time dependence of the amplitude, we obtain l0m0ω0 A|A| 2 e 2ω0,I t + 12λ 2 ω 0,IC l0m0ω0 A|A| 4 e 4ω0,I t , whereC (2) l0m0ω0 is defined bỹ 11/31 andĈ (2) l0m0ω0 is defined bŷ From expression (58), it seems that the divergence in the second order perturbation could be eliminated by adjusting the first order counterterm δC (1) l0m0ω0 . But we will see in appendix B that this is not true, at least for the real part of the divergence.
We see the cancellation of terms of O(ω −2 0,I ) inĈ (2) l0m0ω0 . In Eq. (53), the O(ω −2 0,I ) terms come from l = l 0 , m = m 0 mode, which are can be put together in the following form:

Behavior of the amplitude equation
Before calculating the numerical coefficients, we see some analytic feature of the evolution equation (57). In this section, we abbreviate the subscripts l 0 , m 0 , ω 0 , for brevity. First, we 12/31 comment on the validity of the perturbative renormalization. In the same way as Eqs. (46) and (47), we can rewrite Eq. (57) as the evolution equations of the amplitude and the phase by substituting A = |A|(t)e −ω0,I t−iΘ(t) , to obtain and neglecting the terms higher order in λ. Here, the frequency shift δω is defined by Explicit dependence on the ambiguous terms δC (1) and δC (2) From Eqs. (65) and (66), we find that the evolution of the amplitude and the frequency shift depends on the renormalization scheme, i.e., how to choose δC (1) and δC (2) . In App. B, we show these choices are related to the definition of the amplitude. δω is expected to have a scheme independent meaning. Therefore, the evolution equation for δω becomes scheme independent except for the terms of O(ω 0,I ). However, in order to obtain the evolution equation for δω up to O(λ 2 ), we need to proceed our calculation to one order higher, which turned out to be technically quite challenging. Below, instead of pursuing to obtain the scheme independent equation valid up to O(λ 2 ), we choose some schemes by adding conditions that specify δC (1) and δC (2) without the knowledge of higher order perturbation, and see the qualitative behavior of the cloud.

Minimal subtraction scheme
One scheme is to choose δC (1) = 0 and δC (2) = 0, which corresponds to minimally subtracting the divergent terms. This choice makes renormalization group equation quite simple 13/31 as one can see from Eqs. (65) and (66). After taking this scheme, evolution equations of amplitude and phase are given by δω A qualitative behavior of the axion cloud can be read off from Eq. (68). When the amplitude is Let us analyze the case when the saturation occurs in detail. As we shall confirm later, the attractive nature of the self-interaction implies Re C (1) < 0. This means that Re C (1) /λ < 0 irrespectively of the signature of λ, which is assumed in the following discussion. Then, we can assume Re Ĉ (2) < 0, since otherwise the right hand side of Eq. (70) is always negative. When Re Ĉ (2) − Re C (1) 2 , the saturation occurs at Here, we adopt the appropriate sign in the right hand side of Eq. (70) so as to make it positive.When − Re C (1) 2 Re Ĉ (2) < 0, the saturation of amplitude occurs at The saturation is likely to occur in the weakly non-linear regime, if there is some hierarchy in the magnitudes of coefficients, i.e., if Re Ĉ (2) Re C (1) 2 or Re Ĉ (2) Re C (1) 2 . However, we should recall that we assumed Re Ĉ (2) < 0 at the very beginning of this discussion, which we will find later not to be the case.

Dissipative scheme
Other than the minimal subtraction scheme, we can take the scheme in which the time evolution of the amplitude is totally governed by the dissipative energy loss to the infinity. Instead of putting δC (1,2) = 0, we demand that non-dissipative part in the right hand side of the amplitude equation (65) is set to 0. To do so, we must identify the dissipative part of 14/31 the C (2) . With the aid of the flux-balance equation, we identify as the dissipative part in the C (2) (see appendix C for the detail of the derivation). We observe that C (2)diss is real except for the imaginary part of ω 0 . Therefore, imaginary part is suppressed by ω 0,I and thus we treat C (2)diss as a purely real number. Now, the scheme which respects the dissipation is realized by taking Re Here, we defined the non-dissipative part ofĈ (2) aŝ Note that Re C (1) is also non-dissipative, and suppressed by ω 0,I compared with C (1) itself. This fact explains why we choose Re δC (1) so as to eliminate Re C (1) in Eq.(74). Because Im C (1) and Im Ĉ (2) start with O(ω −1 0,I ), we cannot eliminate these by δC (1) and δC (2) , which is O(ω 0 0,I ). Therefore, since the imaginary parts do not contribute to the evolution equation of the amplitude, we simply take Im δC (2) = 0 .
After taking this scheme, amplitude equations (65) and (66) are Eq. (79) clearly shows that the saturation of the amplitude occurs when which takes the same form as Eq. (71), butĈ (2) is replaced with the dissipative part C (2)diss . As seen trivially from Eq.(79), this saturation occurs irrespective of the sign of λ. However, we should be careful that obtained saturation is real or not. The validity of our approximation would be examined by checking if the non-linearity is already important or not at the saturation amplitude. For this purpose, we can compare the two terms in the right hand side of Eq. (80). If the higher order term dominates, i.e., if it would be a sign of breakdown of the perturbative expansion.

15/31
The above arguments tell us the qualitative behavior of the axion cloud and show the conditions for the saturation to occur in the minimal subtraction scheme and dissipative scheme. We will see in the next section that both schemes fail to show the existence of saturation. We understand that these arguments are still incomplete to conclude the absence of saturation, because there might be some other RG schemes that realize the saturation. However, we think that the absence of saturation in both distinctive schemes seems to be a strong evidence. One might be interested in up to which epoch our approximate treatment can be justified. To see it, we look at the size of the amplitude λ|A| 2 and the relative magnitude of the frequency shift δω/ω 0,R . As long as the perturbation theory is valid, these quantities must stay small.

Numerical results
In this section, we present the numerical calculation of the amplitude equation. In section 4.1, we concentrate on the attractive self-interaction and in section 4.2, we comment on the repulsive interaction. We used Mathematica for numerical calculation.

An example of time evolution of the amplitude equation
Here, we present the time evolution of the axion cloud with the fastest growing mode. First, we will solve the RG equation numerically with the minimal subtraction scheme and the dissipative scheme, then we compare the two schemes to see whether the results change. After examples of the time evolution, we search for the instability terminating region on the (µM, a/M ) plane.

Examples of time evolution
As an example, we present the time evolution of the frequency shift δω in the case when the mass of the axion and the spin parameter of the central black hole are µM = 0.42 and a/M = 0.99, respectively. We consider the cloud evolved from a single mode with l 0 = m 0 = 1. With lowest overtone number n 0 = 0, this parameter set gives the growth rate of the axion cloud as which is nearly the maximum and thus most interesting for the realistic observation.
l0m0ω0 and C  In Table 1, the calculated numerical coefficients that appear in the renormalization group equation is shown. We observe that the real parts of C (1) andĈ (2) are both suppressed compared to its imaginary part, as is expected from the order counting with respect to ω 0,I . 16  (84) The frequency shift blows up in the minimal subtraction scheme. This directly indicates the break down of the perturbation theory eventually occurs. On the other hand, the saturation occurs in the dissipative scheme, as expected by construction. However, the fractional frequency shift δω/ω 0,R is already larger than the unity when the saturation occurs. Therefore, although the qualitative behavior seems different, both two schemes indicate the breakdown of the perturbation. In the small amplitude regime, non-linearity is small and the frequency shift stays small. In this regime, the difference between the two schemes is negligibly small. However, eventually, non-linearity starts to affect the evolution. In the present case, non-linearity makes the cloud more unstable and finally, the frequency shift gets large. This is the consequence of Re Ĉ (2) l0m0ω0 > 0 in the minimal subtraction scheme and that of the condition (82) in the dissipative scheme, as noted in the previous section.
Positivity of Re Ĉ (2) l0m0ω0 could be understood as follows. Note that the m = m 0 mode behaves as an "almost" bound state as we have mentioned below Eq. (48) (we used "almost", because the boundary condition at the horizon is different from the bound state). Owing to the attractive nature of the axion potential (see Eq. (18)), the interaction energy of the axion 17/31 cloud is negative and quartic in amplitude A. This negative potential energy accelerates the increase of the amplitude, which results in the positive sign of Re Ĉ (2) l0m0ω0 . Since our calculation is based on the perturbation theory, time evolution in the large frequency shift regime is not trustable. One might be interested in to what extent our approximate solution can be used as the initial condition for the dynamical simulation. For this purpose, we estimate the time scale of the successive non-linear evolution beyond the weakly non-linear regime in the minimal subtraction scheme. As shown in Fig. 2 as a dotted vertical line, the time beyond which our renormalized solution cannot be trusted is around τ ∼ 12, where we adopt the criterion for the breakdown of perturbative solution that the difference in the frequency shifts in two schemes is greater than 10 %. Also, the time when the blow-up of the amplitude occurs is around τ ∼ 13.5 (the solid vertical line in Fig. 2). Therefore, the required time for the cloud to blow up after the perturbative approximation breaks down is estimated to be O(1) in the unit of the normalized time τ . This indicates that tracking the succeeding non-linear evolution by a numerical simulation with the initial data provided by the current method would be still computationally expensive.
On the other hand, the time scale of the evolutionδ ω/δω/ω 0,I is still O(ω −1 0,I ) during the above intermediately non-linear regime (see the right panel of Fig. 2). This is still much larger than the dynamical time scale ω −1 0,R . Hence, the adiabatic approximation is still valid at this epoch.
We changed the axion mass µM and the spin parameter a/M to see if the breakdown of the renormalized perturbation before the saturation of amplification is an universal feature of the relativistic cloud. Figures 3 shows the evolution of the frequency shift for l 0 = m 0 = 1, 2 modes adopting the minimal subtraction scheme and dissipative scheme. Table. 2 shows the frequency and growth rate of these modes as well as Re C (1) , Re Ĉ (2) and C (2)diss . We confirm that the perturbative expansion breaks down as the frequency shift gets large also in these cases. Table 2: Growth rate of the mode we have calculated. We also present Re C (1) l0m0ω0 , Re Ĉ (2) l0m0ω0 and C

Survey of the parameter space
From the result of the previous section, we find that the growth of relativistic axion clouds does not seem to saturate at least without getting into the strongly non-linear regime, which is beyond the scope of our current approximation. In this section, we investigate whether 18 this is really a generic feature of the self-interacting axion cloud around a rotating black hole or not. We present our results, focusing on the (l 0 , m 0 ) = (1, 1) modes here. We see that termination of the growth by saturation in the weakly non-linear regime is possible only when (1) Re Ĉ (2) is negative in the minimal subtraction scheme, or (2) the size of the left hand side of Eq.(82) is smaller than unity in the dissipative scheme. The second condition can be achieved when |C (2)diss | is sufficiently large. In other words, in both cases the dissipative effect must overwhelms the acceleration of the growth of cloud due to the attractive self-interaction. Thus, we explore the behavior of Re Ĉ (2) , the left hand side of Eq. (82), and |C (2)diss | on the (a/M, µM ) plane.
In Fig. 4, we show Re Ĉ (2) and the left hand side of Eq. (82) on the (a/M, µM ) plane. We see that Re Ĉ (2) is always positive and the left hand side of Eq.(82) is much larger than unity. Also, we show the behavior of |C (2)diss | and |C (2)diss |/Re Ĉ (2) on the (a/M, µM ) plane. We observe that |C (2)diss | gets smaller for the smaller a/M or µM , and 19/31 Re Ĉ (2) holds for any parameter set. We conclude that dissipation effect is relatively small. As a result, the inequality Eq.(82) holds.
In the non-relativistic region (M µ 1), although the cloud extends to a large radius, the dissipation due to the radiation to infinity gets less efficient. One may think that the multipole moments that become the source of radiation become large, and hence the radiative dissipation should be efficient. However, the wave length of the outgoing wave λ out = 1/(3ω 0 ) ∼ 1/(3µ) is much smaller than the size of the cloud r c ∼ M/(µM ) 2 , especially in the non-relativistic regime. Therefore, because of the phase cancellation, only the part of the cloud close to the horizon contributes to the flux integral j l3m03ω0 , which makes the dissipation to infinity inefficient. For this reason, C (2)diss gets small in the non-relativistic regime. These results suggest that the growth of the cloud to a strongly non-linear regime is a generic feature of self-interacting axion clouds.

Repulsive self-interaction
Finally, we will briefly discuss the case of repulsive self-interaction. For negative λ, the time evolution tracks of δω in the two schemes are shown in Fig. 6. Also in this case, the late time behaviors in the two schemes are different, which suggests the breakdown of the perturbative approach. In principle, there is a possibility that only the dissipative scheme maintains the validity of perturbative expansion. Even if we assume so, the significance of the non-linearity at the saturation amplitude can be investigated by examining the magnitude of δω/ω 0,R or the left hand side of (82). As seen in the previous subsection, the inequality (82) is always satisfied (see Fig. 5). This signifies that the breakdown of the perturbative expansion before the saturation occurs even in the repulsive self-interaction.

Summary and Discussion
Axion clouds around spinning black holes will induce interesting phenomena and allow us to test the existence of new particles via gravitational waves or other astrophysical phenomena. For further investigation of this possibility, understanding the precise dynamics of the axion cloud is important. In this paper, We examined the impact of the axion self-interaction on the dynamics of an axion cloud by applying the dynamical renormalization group method. Especially, we focused on the possibility that the energy loss through the emission balances with the energy input due to superradiance to realize a saturated quasi-stationary configuration. If the saturation occurs, the explosive phenomena important for the direct observation, such as bosenova, would not occur.
There are two types of effects caused by the self-interaction. One is a radiation of the axion field which dissipates the energy of the axion cloud. This effect is often ignored in the literature and it was one of our aim to clarify whether this effect saturates the growth of an axion cloud or not. Our calculation indicates that this dissipation does not work efficiently in any case, especially since the dissipation effect gets less important in the non-relativistic regime, µM → 0 and a/M → 0. The other type of effect is that the attractive self-interaction lowers the energy of the cloud to accelerate the instability.
Our results suggest that the saturation of the growth does not occur in the weakly non-linear regime, as long as we consider the φ 4 type self-interaction. More precisely, the 20  perturbative expansion breaks down before the dissipation starts to affect the cloud evolution. Inefficient dissipation due to the axion radiation to infinity is the main reason for this conclusion.
However, our analysis does not indicate that the explosive phenomena like a bosenova and hence bursts of gravitational waves always occur. Several aspects, which have not been taken into account in our current analysis, may terminate the instability. First of all, we cannot say anything about the evolution after the perturbative approach breaks down. In this paper, we truncated the cosine type potential at the φ 4 term. When the cloud gets denser, higher order terms in the potential come into play. The non-linear effects induced by the higher order terms in the cosine type potential may change the evolution.
The other possibility is the contamination of non-superradiant modes (such as m = −1 modes), which can dissipate energy efficiently by falling into the black hole. Within our current setup, there is no room for such modes to appear because the frequencies of the possibly excited modes are restricted to mω 0 , which all satisfy the condition for the superradiance. This situation is likely when we consider an adiabatic growth of the cloud from the fastest growing mode. However, once an explosive phenomenon occurs, the state with plural superradiant modes excited will be realized. In such a situation, the self-interaction can induce non-superradiant excitations. The impact of simultaneous excitations of plural superradiant modes on the evolution of an axion cloud has been already pointed out in Ref. [18,19]. Then, the cloud may develop into a quasi-stationary state, in which superradiant modes extract the mass and angular momentum of the black hole, while non-superradiant modes give them back to the black hole. If such a quasi-stationary state is realized with a relatively small amplitude, i.e., much smaller than the case without self-interaction, in which the gravitational wave radiation balances with the superradiant energy extraction, the spin parameter of the central black hole will change much more slowly. Then, the constraint from the population of black holes on the Regge plane [20][21][22], will be modified from the one without the self-interaction [23].
To tell whether bosenovae actually happen or not in the models with an axion field having an appropriate mass, we need an alternative numerical analysis method that can track the adiabatic evolution of the cloud in the strongly non-linear regime without the truncation of the potential. Full numerical simulation would be one option to treat such a problem [10,11,[24][25][26], although the computational cost is high. Our method will give a better initial data than just assuming a configuration in an ad-hoc manner, when we consider to start the simulation with a mildly large amplitude to reduce the computational cost.
anti-dissipation, which makes the oscillator unstable. On the other hand, the right hand side represents the dissipative effect, which corresponds to the radiative energy loss of the axion cloud.
We will solve Eq. (A13) using perturbation theory. We expand the solution x(t) as x = x 0 + x 1 + · · · . Substituting this expression to Eq. (A13), we get as the zeroth order equation. The solution to this equation is with ω = √ 1 − a 2 . We leave the integration constant C(t 0 ) arbitrary in order to apply RG method.
The first order equation is which has the following solution: +C * 3 a + iω 2(a + 2iω) e 3at e 4iωt e at e −iωt + c.c.
The terms with C 1 (t 0 ) represent the freedom of adding homogeneous solutions to the particular perturbative solution. We observe that some terms behaves as O(a −1 ) in the small a limit, and thus breaks the validity of the perturbation. The term diverging in a → 0 limit in Eq. (A17) becomes a secular term, which is proportional to t − t 0 , after subtraction of an appropriate counterterm as which corresponds to choosing C 1 as As in the case of the previous example, the problematic term can be renormalized. Now, we apply the RG equation to perturbative solution x = x 0 + x 1 + . . . . After some calculations, we can derive an amplitude equation or if we write C = |C|e iΘ . Solutions to an amplitude equations (A22) are Here, C 0 , Θ 0 are the amplitude and the phase determined by the initial conditions of x. After substituting Eqs. (A23) and (A24) to Eqs. (A15) and (A17), we obtain the renormalized first order solution to Eq. (A13) as follows: x = 2|C|e at cos(ωt − Θ) + 2(a 2 + 4ω 2 ) |C| 3 e 3at (−a cos(3(ωt − Θ)) + 2ω sin(3(ωt − Θ))) To see how good our RG method are, we compared numerical solution of Eq. (A13) and our renormalized solution Eq. (A25). Figure A1 shows the renormalized solution gives a good approximation. From this example, we believe that when the equation has a solution with a stable final state for which the perturbative correction remains to be small, we can obtain such a solution by the RG method, even though the zeroth order solution is unstable. 27/31 Let us summarize the recipes of the RG method to solve a differential equations. First, obtain a naive perturbative solution with leaving integration constants of the zeroth order solution arbitrary. The problematic terms, which have a secular growth or are parametrically enhanced to violate the perturbative expansion, will be absorbed by an appropriate choice of the integration constants. Then, choose initial conditions for higher order solutions to correctly cancel the problematic terms in the perturbative solution at some chosen time t 0 . Now, demand that the naive perturbative solution satisfies the RG equation. We can put t 0 = t because solution no longer depends on t 0 due to the RG equation. After all, we obtain the solution without divergence.

B. Scheme dependence of the amplitude equation
In this appendix, we show that the ambiguity of the amplitude equation (65) is related to the ambiguity in the definition of the cloud amplitude.
Now, assume that these two amplitudes are related by the following equation: A(τ ) = A(τ ) + a 1 λA(τ )|A(τ )| 2 + a 2 λ 2 A(τ )|A(τ )| 4 + . . . , Substituting Eq. (B3) into Eq. (B2) and using Eq. (B1), we obtaiñ This shows that the ambiguity in the choice of the homogeneous solution in RG method, which is shift in c i , is related to the definition of the amplitude. This result is general for any flow equation, such as rernomalization group flow in the QFT. In this case, amplitude corresponds to the coupling constant. Now we concentrate on the axion cloud case. In this case, we know that c 0 =c 0 = ω 0,I (see Eq. (65)). Therefore, as long as relation between two amplitudes (B3) do not contains O(ω −n 0,I ), n = 1, 2, . . . quantities, difference in the c 1 andc 1 is suppressed by ω 0,I . We can then identify a 1 = δC (1) . This shows the ambiguity in the counterterm is the ambiguity in the definition of the amplitude. Now we look at the second order equation. Eq. (B6) shows that difference in real part of c 2 is always suppressed by ω 0,I . Therefore, adjusting the amplitude by first order in λ won't eliminate the divergence in the second order solution.
C. Derivation of the dissipative effect C (2)diss l 0 m 0 ω 0 In this appendix, we derive Eq. (73). We start with the evolution equation of the particle number and convert this to that of the amplitude, then identify the coefficient which correspond to the C (2) in the RG equation. 28/31 The number density current J µ for the mode Ψ = Ae −i(ωt−mϕ) R lmω S lmω is given by and obeys conservation equation Following the derivation of Eq. (15), we obtain Here, Σ is constant t surface and ∂Σ is the boundary of Σ, which is r = r + and r = ∞. Therefore, N is the particle number on the constant t surface and the F is the flux through the boundary. When we take superradiant mode (l 0 , m 0 , ω 0 ) for Ψ and neglecting the self-interaction, then the the right hand side of Eq. (C3) is given by There is no contribution from the infinity due to the decaying boundary condition. Using the boundary condition (29a), we obtain ∆ R * l0m0ω0 ∂ r R l0m0ω0 − R l0m0ω0 ∂ r R * l0m0ω0 | r=r+ = 2i(ω − mΩ H )(r 2 + + a 2 ) .
Thus, F is given by On the other hand, left hand side of Eq. (C3) is given by because the frequency ω is complex. Substituting Eqs. (C9) and (C10) into Eq. (C3), we obtain the relation between the particle number and the amplitude as Now we include the effect of the self-interaction by evaluating the flux (C5) of the first order solution (48). The first order solution has two contributions to the flux. One is dissipation into infinity through (l, 3m 0 , 3ω 0 ) modes and other is accumulation via superradiant scattering by the rotating black hole. Because superradiant scattering is suppressed by tunneling through 29/31 angular momentum barrier, we only consider the disipation into infinity. The dissipative part in the first order solution is given by φ (1) (x) ⊃ e −3i(ω0t−m0ϕ) l S l3m03ω0 (θ) dr d cos θ (r 2 + a 2 cos 2 θ ) Using the boundary condition Eq. (29b), behavior around infinity is The coefficient of the |A| 5 term represent the dissipative effect due to the self-interaction on the evolution of the amplitude. Therefore, we define C (2)diss = 9ω 2 0 − µ 2 l |j l3m03ω0 | 2 2π(ω 0 − m 0 Ω H )(r 2 + + a 2 ) dθ|S l0m0ω0 | 2 .