Gravitational radiation from an axion cloud around a black hole: Superradiant phase

Motivated by possible existence of string axions with ultralight masses, we study gravitational radiation from an axion cloud around a rotating black hole (BH). The axion cloud extracts the rotation energy of the BH by superradiant instability, while it loses energy through the emission of gravitational waves (GWs). In this paper, GWs are treated as perturbations on a fixed background spacetime to derive the energy emission rate. We give an analytic approximate formula for the case where axion's Compton wavelength is much larger than the BH radius, and then, present numerical results without approximation. The energy loss rate of the axion cloud through the GW emission turns out to be smaller than the energy gain rate of the axion cloud by superradiant instability until nonlinear self-interactions of axions become important. In particular, an axion bosenova must happen at the last stage of superradiant instability.

Motivated by possible existence of string axions with ultralight masses, we study gravitational radiation from an axion cloud around a rotating black hole (BH). The axion cloud extracts the rotation energy of the BH by superradiant instability, while it loses energy through the emission of gravitational waves (GWs). In this paper, GWs are treated as perturbations on a fixed background spacetime to derive the energy emission rate. We give an analytic approximate formula for the case where axion's Compton wavelength is much larger than the BH radius, and then, present numerical results without approximation. The energy loss rate of the axion cloud through the GW emission turns out to be smaller than the energy gain rate of the axion cloud by superradiant instability until nonlinear self-interactions of axions become important. In particular, an axion bosenova must happen at the last stage of superradiant instability.

Introduction
In a few years, the ground-based detectors, Advanced LIGO, Advanced VIRGO, and KAGRA, are expected to begin operation with sufficiently high sensitivity to detect gravitational wave (GW) signals from binary mergers of black holes (BHs) or neutron stars. We will have a very exciting era of general relativity, and many interesting sciences will be done. One interesting possibility is to find signals from unknown sources or effects of unestablished hypothetical theories. Do we have a possibility to detect signals from string theory through GW detectors? Naïvely it seems difficult because the effect of string theory is expected to appear at very high energy scales (e.g, the string scale). However, the authors of Refs. [1,2] answered "Maybe yes," if there are string axions with ultralight masses. Since string theories (or M-theory) require our spacetime to have 10 or 11 dimensions, the extra dimensions other than four have to be compactified. Then, the extra dimensions have dynamical degrees of freedom of changing their shape and size which are called moduli, and they effectively behave as fields in our four-dimensional spacetime. One of the moduli will be the QCD axion. We can also expect the existence of other pseudoscalar fields with ultralight masses, called string axions, whose expected number is typically from 10 to 100. The allowed mass range is from 10 −10 eV to 10 −33 eV and further below. Depending on their masses, they cause a variety of phenomena that could be observed in the context of cosmophysics (=cosmology and astrophysics). This is called the axiverse scenario (See also Ref. [3] for an overview).
Kerr parameter a is introduced by J = M a. In this paper, we often use the nondimensional rotation parameter a * = a/M . In the Boyer-Lindquist coordinates (t, r, θ, φ), the metric is given by with Σ = r 2 + a 2 cos 2 θ, ∆ = r 2 + a 2 − 2M r.
The horizon is located at r + = M + √ M 2 − a 2 and it rotates with the angular velocity Ω H = a/(r 2 + + a 2 ) as seen by observers at spatial infinity. We consider a real scalar field Φ in the Kerr spacetime. Here, Φ is treated as a test field and its gravitational backreaction is neglected. Because we study the situation much before the bosenova, the effect of the nonlinear self-interaction is less important. For this reason, we ignore the nonlinear term and assume Φ to satisfy the massive Klein-Gordon equation In the Kerr spacetime, the separation of variables of the field Φ is possible as Throughout the paper,f indicates the complexified quantity of a real function f , i.e., f = ℜf . In the following, we denote R = R ωµ ℓm (r) and S = S ωµ ℓm (θ) for simplicity. See Appendix A for the equations that S(θ) and R(r) satisfy.
The superradiant instability occurs only for quasibound states. Such a quasibound state is obtained by assuming the field to decay at infinity and to be ingoing at the horizon. The procedure is similar to the calculation of the bound states of a hydrogen atom in quantum mechanics. In the case of the hydrogen atom, one imposes the wave function to decay at infinity and to be regular at the origin. Then, the energy levels are discretized. In the same manner, the frequency ω is discretized in the case of the BH-axion system. Since there is an ingoing flux at the horizon, the frequency ω becomes complex, ω = ω R + iω I . The imaginary part gives the exponential behavior, and the field decays for ω I < 0 and grows for ω I > 0. In the case ω I > 0, the bound state is unstable, and this happens when the superradiance condition ω R < mΩ H is satisfied. Physically, the instability happens because the energy density in the neighborhood of the BH horizon becomes negative for ω R < mΩ H , and negative energy falls into the horizon. Then, due to the energy conservation, the amplitude of the outside field becomes larger. This is the superradiant instability.
Approximate solutions for superradiant bound states were constructed in the parameter region M µ ≪ 1 by the matched asymptotic expansion (MAE) method [4] and in the parameter region M µ ≫ 1 by the WKB method [5]. The solution in the MAE method is closely related to our analysis in Sec. 3. In this method, two approximate solutions are constructed in a near-horizon region and in a distant region, and they are matched in the overlapping region where both approximations are valid. The equation in the distant region is identical to the Schrödinger equation for a hydrogen atom except that the electric potential −e 2 /r is  replaced by the Newton potential −M µ/r. Therefore, the BH-axion system can be regarded as a huge gravitational atom, which is characterized by the angular quantum numbers (ℓ, m) and the radial quantum number n r = 0, 1, 2, .... In the parameter range M µ ∼ 1, numerical calculations are required to determine the bound-state levels. Such analyses were done by several authors [6][7][8]. In Ref. [8], a detailed analysis was reported using the continued fraction method developed by Leaver [13]. We also reproduced the consistent result using our independent code. See Appendix A for technical details. Figure 1 shows the imaginary part M ω I of the eigenfrequency as functions of M µ for the BH rotation parameter a * = 0.90 and 0.99. The results for the modes ℓ = m = 1, 2, and 3 and n r = 0 are shown. The unstable range of M µ becomes larger as m is increased because of the superradiance condition µ ≈ ω R < mΩ H . The peak value of M ω I becomes smaller as ℓ = m is increased and as a * is decreased. Figure 2 shows an example of the radial profile of R for the bound state solution with ℓ = m = 1 (the left panel) and the corresponding configuration in the equatorial plane (the right 5/31 panel). The system parameters are adopted to be M µ = 0.40 and a * = 0.99. The function R is shown as a function of the tortoise coordinate r * defined by Figure 3 shows the same information for the mode with ℓ = m = 2 and the system parameter M µ = 0.89. In these two figures, we chose the modes which are reduced to gravitational atoms with the radial quantum number n r = 0 for M µ ≪ 1. For these modes, the field Φ has one maximum and one minimum for ℓ = m = 1, and two maxima and two minima for ℓ = m = 2 in the equatorial planes (i.e., right panels) outside of the ergoregion. In the left panels, ingoing waves which behave as R ∼ e −i(ω−mΩH )r * can be seen. These waves carry negative energy density into the BH. In Sec. 4, we use these numerical solutions for mode functions of the quasibound states as the source term of GWs.

GWs from an axion cloud
The energy-momentum tensor of the axion field is where T µν (A, B) is defined as We study GWs generated by this energy-momentum tensor T µν . For this purpose, we consider the perturbationg whereg µν is the spacetime metric, g µν is the background Kerr metric, and h µν is a small perturbation. We introduce ψ µν by Then, in the de Donder gauge 6/31 the Einstein equations for ψ µν become where ∆ L is the Lichnerowicz operator that reduces in the vacuum case to Equation (11) can be rewritten as the equation for h µν , Here, we introduce one approximation. The field Φ is proportional to e ωI t by the superradiant instability, and the energy-momentum tensor gradually grows larger. However, the time scale of the instability is very long, 10 7 M , and its effect on the GW flux must be small. For this reason, we ignore the imaginary part ω I and calculate GWs from an oscillating source without exponential growth. This also means that we ignore the ingoing flux of the field Φ to the BH horizon.
In addition to this, there is another subtlety. Originally, the superradiant unstable mode for Φ is obtained as a complex functionΦ proportional to e −iωt+imφ . Then, the energymomentum tensor should be estimated by using its real part Φ = ℜΦ as where z * denotes the complex conjugate of z. Clearly, the left-hand side of Eq. (14) is not the real part of T (Φ,Φ)/2 whenΦ are sum of several oscillating modes proportional to e −i(2ω)t+i(2m)φ because of the presence of the term T (Φ,Φ * ). In fact, this term can be interpreted as the source for the GW emission by the level transition [2]. However, for a monochromatic eigenmode of Φ, the stationary terms corresponding to T (Φ,Φ * ) generate only stationary perturbations h µν , and we are not interested in such perturbations. Hence, we neglect this contribution and focus on GWs generated by the oscillating terms of Eq. (14). For this purpose, it is convenient to introduce the complexified energy momentum tensorT defined by SinceT is proportional to e −2iωt+2imφ , the frequencyω and the azimuthal quantum number m of GWs satisfyω = 2ω,m = 2m.
Hereafter, tilder indicates a quantity for GWs. Once a solution for the mode function is given, the GW radiation rate toward null infinity can be calculated by dE whereḣ TT µν denotes the time derivative of a metric perturbation in the transverse-traceless (TT) gauge, and dS is the area element of the 2-surface in a distant region specified by t and r = r out . Here, the conditions for the TT gauge are ∇ µ h TT µν = 0 and g µν h TT µν = 0, (18) and the existence of the TT gauge has been shown for a vacuum perturbation of a vacuum spacetime (e.g., p. 186 of Ref. [14]). On the other hand, the formula for the energy flux 7/31 dE (bh) GW /dt absorbed by the horizon of a Kerr BH is derived in Ref. [15] in terms of the Newman-Penrose quantity ψ 0 or ψ 4 . A consistent formula is found using the Isaacson effective energy-momentum tensor for high-frequency GWs in Ref. [16]. In this paper, we do not explicitly calculate the GW energy flux across the horizon with the following reason: In Refs. [15,16], dE (bh) GW /dt has been shown to be proportional toω −mΩ H . As stated in Eqs. (16), the frequency and the azimuthal quantum number of GWs satisfyω = 2ω and m = 2m. Then, if the axion cloud satisfies the superradiance condition ω < mΩ H , emitted GWs also satisfy the superradiance conditionω <mΩ H . Therefore, the GW energy flux to the horizon is negative, dE (bh) GW /dt < 0, and hence, the total radiation rate is This means that dE (out) GW /dt gives the upper bound of the total energy loss rate dE GW /dt of the axion cloud due to the GW emission, and if dE (out) GW /dt is smaller than the energy extraction rate dE a /dt of the axion cloud, we can safely declare that the GW emission does not stop the growth of the axion cloud by consuming all the energy extracted from the BH by the superradiant instability.
In order to estimate the GW emission rate toward null infinity with the aid of the formula (17), in principle, we have to solve Eq. (11) first and then find a gauge transformation that puts the solution into the TT gauge. However, this is a rather intricate task. Therefore, we adopt the following different method.
First, we pay attention to the fact that the solution h µν of Eq. (13) satisfies the vacuum Einstein equations outside the finite source region. Therefore, it is related to the homogeneous out-mode solutions 2û (j) µν at future null infinity I + aŝ Here, we introduce a complex metric perturbationĥ µν ∝ e −iωt+imφ with h µν = ℜ[ĥ µν ] as before, andj indicates the collection of indexes to specify each mode, i.e. the angular quantum numbersl andm and the "polarization state" P = ±1 (in the definition of Chrzanowski [16]). The out-modeû µν is a solution to the homogeneous equation satisfying the boundary condition of vanishing flux across the future horizon H + (i.e., the absence of ingoing waves at r * → −∞). Since we consider monochromatic waves, the outmode solutionsû (j) µν are also in proportion to e −iωt+imφ . By a gauge transformation, the perturbation can be expressed in the TT gauge, in whichĥ TT µν can be expanded in terms of u (j)TT µν with the same expansion coefficients Cj. By inserting this expansion in the TT gauge into (17), we obtain where X implies the time average over a time sufficiently larger than 1/ω. Next, we determine Cj. First, we apply Green's theorem to the identity This relation holds for an arbitrary domain D. Let us adopt as D the domain surrounded by t = t 0 and t = t 0 + ∆t, r = r out , and r = r in . Here, r = r out is located at a sufficiently far position from the axion cloud, and r = r in is located close to the horizon.
In the left-hand side of the relation (25), the surface integrals at t = t 0 and t = t 0 + ∆t remain finite and we can safely neglect their contribution if we take a sufficiently large ∆t. Then, the both sides of this relation are proportional to ∆t. Besides, because of the boundary condition for u * (j)TT µν imposed at the future horizon, the surface integral at r = r in vanishes.
where the "inner product" of u µν and T µν is defined as Note that because u TT tµ ≈ u TT rµ ≈ 0, only the angular component h IJ (I, J = θ, φ) appear in the right-hand side of Eq. (26). From this, it follows that the right-hand side is invariant under the gauge transformation h IJ → h IJ + δh IJ with δh IJ = ∇ I ξ J + ∇ J ξ I . This implies that h µν can be replaced by its counterpart in the TT gauge, h TT µν . Hence, substituting the expansion formula (20) in the TT gauge into right-hand side of Eq. (26), we obtain Here, we used the fact thatû (j)TT µν behaves as ∼ e −iω(t−r) /r for large r. Furthermore, in spite of the fact that the relation (26) was derived assuming the TT gauge conditions, the inner product u (j)TT , T is invariant under the gauge transformation u TT µν →û µν =û TT µν + δû µν with δû µν = ∇ µ ξ ν + ∇ ν ξ µ because δu, T = 0 holds. This can be shown by rewriting the integral of (∇ µ ξ ν )T µν = ∇ µ (T µν ξ ν ) in the domain D into the surface terms by the Gauss law and using the absence of the energy flux at the boundaries r = r in and r out . Therefore, the inner product û (j) , T can be calculated without taking a special care to the gauge choice, and it is sufficient to find a solution u (j) µν satisfying the TT gauge 9/31 condition only at future null infinity. In the Kerr background, we can explicitly construct such out-mode solutions as we will see later.
Substituting Eq. (28) in an arbitrary gauge into Eq. (22), the energy emission rate is given by with Jj defined in Eq. (23).

Methods for calculating the energy emission rate
In order to evaluate the radiation rate using Eq. (29), we have to calculate the homogeneous solutionû µν for a gravitational perturbation and perform the integration of the inner product u, T . Also, the values of Jj have to be calculated. We explain the methods for calculating these quantities.

Homogeneous GWs.
We use the Teukolsky formalism [18] in order to calculate the homogeneous solution of a gravitational perturbation of the Kerr spacetime. By setting M = a = 0, this formalism also can be applied to a perturbation for a flat spacetime studied in the next section. The Teukolsky formalism realizes the separation of the variables (t, r, θ, φ) of the perturbative equations using the Newman-Penrose formalism [19]. In the Teukolsky formalism, the Kinnersley null tetrads are adopted: In terms of these tetrads, the metric (2a) of the Kerr spacetime is The equations for the Newman-Penrose quantities ψ 0 and ρ −4 ψ 4 are given by with s = ±2, respectively. Here, the energy-momentum tensor is set to be zero because we need to generate homogeneous solutions. By writing ψ = e −iωt e imφ s Rω ℓm (r) s Sω ℓm (θ), the 10/31 Teukolsky equation (32) is separated as The function s Sω ℓm (θ)e imφ is called a spin-weighted spheroidal harmonics, and s Alm is its eigenvalue. For a spherically symmetric spacetime, the spin-weighted spheroidal harmonics reduce to the spin-weighted spherical harmonics s Ylm(θ, φ) with s Alm =l(l + 1) − s(s + 1). We adopt the standard normalization condition s Sω ℓm (θ) GW /dt that gives the upper bound on the total radiation rate. In order to evaluate this quantity, we need to generate the out-mode solutionû µν for a gravitational perturbation of a Kerr spacetime. The asymptotic behavior of R for the out-mode is where r * is the tortoise coordinate defined in Eq. (5). In this paper, we choose s = +2 for a technical reason that is explained in Sec. 4.1.
Once the functions +2 R(r) and ±2 S(θ) are obtained, the next step is to reconstruct the metric perturbationû µν from the Teukolsky functions. The function ψ 0 or ψ 4 is known to completely determine the metric perturbation [20]. Intuitively, this is because ψ 0 or ψ 4 is a gauge-invariant quantity, and its real and imaginary parts correspond to the two degrees of freedom of GWs. The explicit formula was derived by Cohen and Kegeles [21,22] and Chrzanowski [16] under some assumptions, and its complete proof was given by Wald [23]. Here, we adopt the following metric formula given by Chrzanowski in the outgoing radiation gaugeû µν n ν =û ν ν = 0: 3 where A, B, and C are operators Here, α, β, γ, µ, and π are the Newman-Penrose variables [19] (see [18] for expressions in the Kinnersley tetrad), and ∆ = n µ ∇ µ , δ = m µ ∇ µ . The polarization-state parameter P takes a value +1 or −1, and the perturbations with P = ±1 are reduced to the even-type and odd-type perturbations of a Schwarzschild spacetime in the nonrotating case, a * = 0, respectively [25,26] (or the scalar and vector modes in [27]). For outgoing GWs, the outgoing radiation gauge agrees with the TT gauge at the distant region r ≫ M : Here, we wrote only the I, J = θ, φ components because the other components are subdominant. The explicit formulas for S IJ and V IJ are Note that in the case of a spherically symmetric spacetime, these tensors S IJ and V IJ are identical to those given in Ref. [27] for the transverse-traceless part of gravitational perturbation except for overall factors. For this asymptotic form ofû µν , the value of J defined in Eq. (23) is calculated as for both P = ±1. Substituting Eq. (41) into the formula for the radiation rate (29), we have 2.3.2. Energy-momentum tensor. Once the out-mode GW solution u µν and the value of Z out are obtained, we can calculate the inner product using the energy-momentum tensor (6). Because the homogeneous solutionû µν of Eq. (37) is spanned only by n µ n ν , n (µ m ν) , n (µ m * ν) , m µ m ν , and m * µ m * ν , the terms proportional to g µν do not contribute to generation of GWs. For this reason, we only considerT ′ µν = (1/2)∇ µΦ ∇ νΦ , and the components necessary for 12/31 our calculation arê T ′µν m µ m ν = e −2iωt e 2imφ 4(r + ia cos θ) 2 In the following two sections, we apply this formalism to the analytic approximation in the flat background spacetime and fully numerical calculations in the Kerr background spacetime, respectively.

Gravitational radiation in the flat approximation
In this section, we derive approximate formulas for the GW radiation rate in the case M µ ≪ 1. In this case, most of the axion cloud distributes in a region far from the BH. Detweiler [4] constructed an approximate solution by the MAE method for this case. As stated earlier, the outer solution agrees with the wave function of a hydrogen atom. The solution iŝ Here, L 2ℓ+1 n−ℓ−1 (x) represents the Laguerre polynomial, and k is defined by The value of k is discretized as where n = 1, 2, ... is the principal quantum number defined by with the radial quantum number n r = 0, 1, 2, .... Note that the gravitational analogue of the Bohr radius is given by a 0 = 1/M µ 2 = n/k. In Eq. (44), we added the extra factor √ 2E a /ω compared to the standard wavefunction of a hydrogen atom in order to normalize the axion field so that is satisfied, where dΩ = sin θdθdφ. We ignore the near-horizon solution because it is small and gives a minor contribution to the generation of GWs. In this approximation, the axion field is bounded by the Newton potential −M µ/r, and the spacetime is approximately flat. For this reason, we approximate the background metric g µν as the flat spacetime metric η µν in calculating the homogeneous GW solution and the inner product u, T . We call 13/31 this approximation the flat approximation. In this flat approximation, we replace g µν of the energy-momentum tensor (6) by η µν , and raise and lower the tensor indices µ, ν by η µν .
Here, we point out a potential problem in this approximate method. The approximate solution (44) for the axion field satisfies the equation The right-hand side comes from the Newton potential, or in other words, the contribution from the static perturbation generated by the BH mass M . On the other hand, in calculating the homogeneous GWs and the inner product, we completely ignore the background curvature and use the flat metric η µν . As a result, we have and the conservation of the energy and momentum is weakly violated. Because of this, if we adopt the gauge transformation δu µν = ∂ µ ξ ν + ∂ ν ξ µ , we have and therefore, the gauge invariance of the inner product is not guaranteed. For this reason, we have to check carefully to what extent the "error" in Eq. (51) can be large. We will come back to this point at the last part of this section.

Solution to Teukolsky equation
In the flat spacetime, the radial Teukolsky equation becomes Here, the rescaled radial coordinate x =ωr is introduced. This equation can be solved analytically. Choosing s = +2, the general solution is given as where F and U denote the confluent hypergeometric functions of the first and second kinds, respectively. We choose the integral constants to be A = 0, B = 1; so that the radial function becomes regular at the origin, x = 0. The asymptotic behavior at large x is Here, the first term represents outgoing waves. From this asymptotic behavior, the coefficient Z out for the asymptotic outgoing waves in Eq. (36) is read to be 14/31

Energy emission rate
The remaining task is to obtain the homogeneous solutionû (j) µν and calculate the integral u (j) , T , where the index (j) is a shorthand for (l,m, P ) to label each GW mode as noted in Sec. 2.2. Since the calculation of the inner product is rather tedious, we sketch the calculations in Appendix B, and just present the results here. It is convenient to discuss the results for the odd-type GWs (P = −1) and for the even-type GWs (P = +1) separately.

Odd-type perturbation.
In the case of odd-type perturbation, P = −1, the inner product u, T vanishes after integration with respect to the angular coordinates (θ, φ). Therefore, we find that GWs of the odd-type modes are not radiated in our BH-axion system in the flat approximation. This is a natural result because the oscillating part of the energymomentum tensor of the axion cloud has just the even-type mode. The stationary part of the energy-momentum tensor generates the odd-type perturbation that corresponds to the gravitational angular momentum, and this part has been ignored in our analysis.

Comparison with the superradiant growth rate
Here, we compare the GW radiation rate with the energy extraction rate of the axion cloud by the superradiant instability. For the situation M µ ≪ 1, an approximate formula for the growth rate by the instability has been derived by Detweiler [4] with the MAE method. The approximate growth rate is given by Note that the definition of n in Ref. [4] is different from ours: it corresponds to our radial quantum number n r . With this ω I , the energy extraction rate by the superradiant instability is given by Eq. (1). Because the energy extraction rate and the GW radiation rate are proportional to (M µ) 4ℓ+5 and (M µ) Qℓ with Q ℓ = 4ℓ + 10, respectively, the GW radiation rate has a higher power of (M µ). Since our approximation holds for M µ ≪ 1, the GW radiation rate is expected to be much smaller than the energy extraction. Figure 4 shows the GW radiation rate dE GW /dt and the energy extraction rate dE a /dt by the superradiant instability normalized by (E a /M ) 2 as functions of M µ for ℓ = m = 1, 2, and 3. The BH rotation parameter is fixed to be a * = 0.99, and we show the two cases where the energy of the axion cloud is E a /M ∼ 10 −3 and 10 −1 . These two values correspond to the energies when the bosenova happens for the choice of the decay constant f a = 10 16 and 10 17 GeV, respectively [9]. The figure shows that the energy extraction rate is much larger. For other values of the BH rotation parameter a * , we have found that the result does not change in the region where the superradiant instability works effectively. Therefore, in the region M µ ≪ 1, the GW radiation rate is smaller than the energy extraction rate, and the GW emission does not hinder the growth of the axion cloud by the superradiant instability for a wide range of system parameters.

Reliability of the flat approximation
Now, we discuss the problem of the gauge dependence of the inner product in the flat approximation, Eq. (51). Here, we consider the gauge transformation δu of u generated by a vector field ξ µ ∝ e −iωt . Since δu ∼ ωξ and ∂ t Φ 2 ∼ ωΦ 2 , we have Here, we have introduced the small parameter β := k/ω ≈ µM/n. On the other hand, since the order of the tt component of the energy-momentum tensor is T tt ∼ ω 2 Φ 2 , readers may expect that u, T ∼ (ω 2 /∆t) uΦ 2 √ −gd 4 x, and thus, δu, T ∼ β u, T ≪ u, T . However, this is incorrect because the leading order terms with respect to β cancel out in the calculation of u, T as explained just before Eq. (B21) in Appendix B. The fact is that because of this In all cases, the energy extraction rate is larger than the GW radiation rate. cancellation, and therefore, u, T ∼ δu, T . The change in the inner product caused by the gauge transformation contribute to the leading order of u, T . Although this problem of the gauge dependence would be solved by introducing the static perturbation from a flat spacetime that corresponds to the Newton potential, the analytic calculation will become much more difficult in such an analysis. From the observations above, the value of the coefficients C nℓ of Eq. (57c) is not reliable and may be changed by a factor. In fact, we also calculated the inner product u, T using the gauge-invariant formalism for perturbations in spherically symmetric spacetimes in general dimensions [27] in the radiation gauge u 0µ = u µ µ = 0 and found that the prefactor C nℓ for the radiation formula (57a) for the even-type mode differs from Eq. (57c) by a factor of (ℓ − 1) 2 /4 (although the result for the odd-type mode is unchanged). However, we would like to stress that the power dependence on (µM ) must not be changed by the gauge choice: We can trust the power Q ℓ of (µM ) given by Eq. (57b). We will confirm this statement in the fully numerical calculation in the Kerr background spacetime in the next section.
Note that this gauge dependence just appears in the case of the flat approximation. In the full numerical calculations in the Kerr background below, the gauge invariance of the inner product is guaranteed because the energy conservation ∇ µ T µν = 0 is fully satisfied.

Gravitational radiation in Kerr background
In this section, we explain the numerical method of computing the GW radiation rate from an axion cloud in the Kerr background spacetime and present the numerical results. This analysis can be applied for the range of the system parameter M µ ∼ 1.

Numerical method
The numerical method of calculating the axion bound stateΦ = e −iωt e imφ R(r)S(θ) was already explained in Sec. 2.1. The bound-state frequency and the radial function R(r) are obtained by the continued fraction method. The angular function S(θ) is approximately calculated with an expansion formula S(θ) = i S (i) c i with c 2 = −a 2 k 2 up to the sixth order for each (ℓ, m). See Appendix A for details. Once these functions are obtained, the necessary components of the energy-momentum tensor are calculated with Eqs. (43a)-(43e).
In order to generate the homogeneous solutions u (j) µν , we first calculate the functions and quantities that appear in the Teukolsky function: The spin-weighted spheroidal harmonics s Sω ℓm (θ)e imφ and its eigenvalue s Aω ℓm , and the radial function s Rω ℓm (r).
For the eigenvalue s Aω ℓm , we adopt the approximate formula for small c := aω given in Refs. [28,29] that is applicable for arbitrary values of (s,l,m). This approximate formula is given in the form of the series expansion with respect to c = aω up to the sixth order. Since approximate formulas for the functions s Sω ℓm have not been given in an existing literature, we have derived the series expansion s S(θ) = i s S (i) c i up to the sixth order for each value of (s,l,m). The reason why we use the approximate formula is that regularization at poles is necessary when we calculate the metric from the Teukolsky functions. This regularization procedure is much more difficult for numerically generated data. The radial function +2 Rω ℓm (r * ) was numerically generated using the fourth-order Runge-Kutta method starting from r * = −200 to outward up to r * /M = 1000 or 4000 depending on the situation. Here, we choose s = +2 because the ingoing mode behaves as s Rω ℓm ∼ ∆ −s e −i(ω−mΩH )r near the horizon and it is easy to kill this mode for a positive s. In order to calculate the radiation rate, we have to determine the value of Z out from the numerical data. This is done by fitting the numerical data with the asymptotic expansion formula of the form Here, the coefficients A 1 , ..., A 5 and B 1 , B 2 are determined from the asymptotic expansion of the radial Teukolsky equation (33), and the fitting parameters are Z in and Z out . Although the outgoing mode decays more rapidly compared to the ingoing mode at large r, the value of Z out can be evaluated fairly accurately by increasing the numerical accuracy of +2 Rω ℓm (r). The fit was done in the range 100 ≤ r * /M ≤ 300, because the numerical noise increases as the value of r * /M is further increased.

General properties
The modes of radiated GWs can be restricted by the symmetry properties of the (spinweighted) spheroidal harmonics. Assuming that the axion cloud is in the mode ℓ = m, we have S ℓm (θ) = S ℓm (π − θ).
In the flat approximation in the previous section, all vector modes (P = −1) vanished. But in the case of the rotating Kerr background, our numerical data show that the P = −1 modes are nonzero as we see below. The reason can be understood by considering the slow rotation case. The effect of the BH rotation is given by the odd-type perturbation, and the coupling between the even-type axion distribution and the odd-type rotation becomes the source for the odd-type GWs at the second order.

Numerical results
Now we show the numerical results. for M µ ≪ 1. Figure 5 shows the GW radiation rates dE In each panel of a * = 0.50, 0.90, and 0.99, the vertical line indicates the threshold of the superradiant instability: In the left-hand side of the line, the axion cloud grows by superradiant instability at least if the GW emission is neglected. In the right-hand side, no superradiant instability occurs, and the axion cloud, even if it were produced by some mechanism, would simply shrink gradually due to infall into the BH. We find the general tendency such that the radiation rate decreases for very large M µ. The reason is as follows. In the superradiant regime, there is always a potential minimum of the axion field and the quasibound state is formed. Although this potential minimum also exists for M µ that is not much larger than the threshold, it disappears at some point as M µ is further increased. In such a situation, the field cannot form a bound state and it just falls into the horizon. Then, the field is concentrated near the horizon and the GW emission becomes inefficient, because the redshift effect becomes more and more significant as M µ increases further. The GW radiation rate decreases rapidly with M µ.
In the panel of a * = 0.00, the result of the flat approximation is also shown. For a small value of M µ, our numerical data obey the same power dependence on M µ as the approximate formula, dE GW /dt ∝ (M µ) 14 . There is a shift by a constant factor from the curve of the approximate formula, because the value of C nℓ includes the error by a factor as discussed in Sec. 3.4. Figure 6 compares the total GW radiation rate dE GW /dt and the energy extraction rate dE a /dt normalized with (E a /M ) 2 as functions of M µ for a * = 0.90 and a * = 0.99. Here, the total GW emission rate is approximated by sum of the GW radiation rates for the first four 20/31 Since the power dependence on (E a /M ) of the two rates are different from each other, we have to specify this value to compare them. As we have done in the flat approximation, we consider the two cases, E a /M = 10 −3 and 10 −1 . Each curve of the energy extraction rate crosses the curve of the GW radiation rate at the point where the former curve suddenly drops near the threshold of the superradiant instability. Therefore, the GW emission scarcely affects the growth of the axion cloud due to the superradiant instability.
4.3.2. ℓ = m = 2. We turn our attention to the results for the axion cloud in the ℓ = m = 2 mode with the radial quantum number n r = 0. Figure 7 shows the radiation rates dE   of the GW radiation rates for the first four modes with respect tol, and the curves of the energy extraction rate for the cases E a /M = 10 −3 and 10 −1 are shown. Similarly to the case ℓ = m = 1, the GW radiation rate is much smaller than the energy extraction rate except for a very small region near the threshold of the superradiant instability. Therefore, the GW emission scarcely affects the growth of the axion cloud due to the superradiant instability also in the case ℓ = m = 2.
4.3.3. ℓ = m = 3. Finally, we show the results for the axion cloud in the ℓ = m = 3 mode with the radial quantum number n r = 0. Figure 9 shows the radiation rates dE    Fig. 10: The same as Fig. 6 but for the axion cloud in the ℓ = m = 3 mode. The total GW energy emission rate is evaluated by summing the first four GW modes ℓ = 6, 7, 8, and 9.
and 0.99. The general features are same as the previous two cases. Although only the mode (l,m) = (6, 6) is radiated in the nonrotating case, the contribution of other modes become important as the value of a * is increased. The GW radiation is suppressed for a very large value of M µ. In this case, the numerical data coincide well with the flat approximation for M µ ≪ 1 as shown in the panel of a * = 0.00. Figure 10 compares the total GW radiation rate (approximated by sum of the radiation rates for the first four modes with respect tol) and the energy extraction rates for E a /M = 10 −3 and 10 −1 as functions of M µ for a * = 0.90 and 0.99. Again, the GW radiation rate is much smaller than the energy extraction rate except for a very small region near the threshold of the superradiant instability.

Summary of the numerical result
We have calculated the GW radiation rates to infinity, dE (out) GW /dt, as functions of M µ from the axion clouds in the modes ℓ = m = 1, 2, and 3. Our numerical data have the same power dependence as the flat approximation for small values of M µ. We compared the numerical results of dE GW /dt with the energy extraction rates dE a /dt of the axion cloud. The value of dE (out) GW /dt is much smaller than dE a /dt for the both cases E a /M = 10 −3 and 10 −1 . Therefore, similarly to the flat approximation, the axion cloud grows by the superradiant instability until the effect of the nonlinear self-interaction becomes important also when the background spacetime is adopted to be a rotating Kerr BH spacetime.

Conclusion
In this paper, we have studied GW emissions from an axion cloud in the superradiant phase around a rotating BH by combining analytic methods and numerical calculations. First, for M µ ≪ 1 for which the Newtonian approximation is good, we have derived the analytic formula (57a)-(57c) in the flat spacetime approximation. This formula can be translated into the expression for the observed non-dimensional GW amplitude h as 23/31 with C nℓ given in the Table 1, where d is the distance to the source from us, α g is defined by and we have recovered G and c. This formula cannot be used when α g becomes of the order of unity, because the axion cloud approaches the BH horizon and relativistic effects become significant. Therefore, for the parameter range α g ∼ 1, we have computed numerically the GW emission rate from the axion cloud corresponding to the modes ℓ = m = 1, 2, and 3 by solving the perturbation of the Kerr background spacetime exactly. As shown in Figs. 5, 7, and 9, the results of the numerical calculations indicate that the GW emission rate for each GW mode deviates from the above analytic formula when α g approaches unity, and rapidly decreases beyond some critical value of α g due to relativistic effects. This critical value of α g depends on the angular momentum parameter a * = J/M 2 and the quantum number ℓ characterising the superradiant mode of the cloud and increases with a * and ℓ. Interestingly, up to this critical value, the above analytic formula reproduces the numerical results rather well when a * is large if the contributions of all emission modes are summed up, although the deviation becomes significant at smaller α g for a * ≪ 1.
These GW emission rates were compared with the energy extraction rate by the superradiant instability. When the axion decay constant is in the range f a = 10 16 -10 17 GeV [9] and the angular momentum is in the range a * = 0.90-0.99, we have found that the GW emission rate is much smaller than the energy extraction rate except for a small region near the threshold of superradiant instability, as illustrated in Figs. 4, 6, 8, and 10. Therefore, the GW emission does not hinder the growth of the axion cloud through the superradiant instability for a wide range of the system parameters.
In our previous paper [9], we simulated the time evolution of an axion field around a Kerr BH taking account of nonlinear self-interaction, and showed that "axion bosenova" happens as the result of superradiant instability. Our results in this paper indicate that we do not have to change this picture even if we take into account the backreaction of the GW emission from the axion cloud. Since the GW emission rate is sufficiently low, the superradiant instability continues until the bosenova happens by nonlinear self-interaction.
Because we can safely declare that bosenova happens, our next task is to calculate the GW emission during bosenova. In our previous paper [9], we gave a preliminary estimation of the emission rate by the quadrupole formula. The result can be expressed as where ∆E a is the change of the axion cloud energy during the bosenova collapse. ∆E a includes the energy that falls into BH, which is around (0.05-0.2)M . Thus, the main difference between the estimation based on the flat approximation and that obtained by the quadrupole formula comes from the numerical coefficient and the power of α g . Here, for small ℓ, we have found that C N is of a similar order to or less than C Q . The power of α g for the flat approximation is always larger than that for the quadrupole formula. Hence, the quadrupole formula always gives a larger emission rate than that during the superradiant phase, which is detectable by the advanced LIGO, advanced VIRGO and KAGRA if the source is inside our galaxy. This also implies that bosenova can produce strong GW emissions for small α g if the instability growth time is shorter than the age of the central BH.
In the present paper, we have studied the case in which only a single mode grows by superradiance. In realistic situations, however, there may exist several unstable modes. Although the instability growth rates for them are largely different, the amplitudes of all of them can become large simultaneously if the instability growth times are all shorter than the age of the BH. Then, GW emissions in the superradiant phase will be more complicated because the emitted GWs will be the superposition of those from the two-axion annihilation processes in several bound-state levels and from the process of axion level transitions. This multi-level occupation will also make the bosenova collapse and the associated GW emissions much more complicated. Further, we should also carefully estimate the detectability of GW emissions in the superradiant phase because it lasts for much longer time than the bosenova collapse. Clearly, GW emissions in such realistic situations can be correctly evaluated only by combining the method developed in the present paper and numerical simulations of bosenova in our previous paper [9]. This program is now in progress.
In Sec. 3, we presented the formula for the gravitational radiation rate. The radiation rate vanishes for the odd-type modes, and is given by Eqs. (57a)-(57c) for the even-type modes. These results are derived from Eqs. (42) and (56) after calculating u, T . In this section, we sketch the calculation of this inner product u, T .
For simplicity, we denote the radial function and spin-weighted spherical harmonics for the Teukolsky function asR := +2 Rlm(r) and sỸ := s Ylm(θ, φ), respectively. On the other hand, R = R ℓm (r) and Y := Y ℓm (θ, φ) represent the radial function and spherical harmonics for the axion field.
Then, the functions appearing in the Chrzanowski formula (37) are calculated as On the other hand, the necessary components of the energy-momentum tensor, Eqs. (43a)-(43e), becomeT