Deci-Hz gravitational waves from the self-interacting axion cloud around the rotating stellar mass black hole
Hidetoshi Omiya Takuya Takahashi Takahiro Tanaka Hirotaka Yoshino






omiya@tap.scphys.kyoto-u.ac.jp
Department of Physics, Kyoto University, Kyoto 606-8502, Japan
ttakahashi1@rikkyo.ac.jp
Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan

t.tanaka@tap.scphys.kyoto-u.ac.jp
Department of Physics, Kyoto University, Kyoto 606-8502, Japan
Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

hyoshino@omu.ac.jp
Department of Physics, Osaka Metropolitan University, Osaka 558-8585, Japan
Nambu Yoichiro Institute of Theoretical and Experimental Physics (NITEP), Osaka Metropolitan University, Osaka 558-8585, Japan








Gravitational waves from condensates of ultra-light particles, such as axion, around rotating black holes are a promising probe to search for unknown physics. For this purpose, we need to characterize the signal to detect the gravitational waves, which requires tracking the evolution of the condensates, including various effects. 
The axion self-interaction causes the non-linear coupling between the superradiant modes, resulting in complicated branching of evolution. Most studies so far have considered evolution under the non-relativistic approximation or the two-mode approximation. In this paper, we numerically investigate the evolution of the axion condensate without these approximations, taking higher multipole modes into account. We also investigate the possible signature in gravitational waves from the condensate. We show that the higher multipole modes are excited, leading to the gravitational wave signal by the transition of the axion between different levels. The most prominent signal of gravitational waves arises from the transition between modes with their angular quantum numbers different by two. 
The gravitational wave signal is emitted in the deci-Hz band for stellar mass black holes, which might be observable with the future gravitational wave detectors.



Deci-Hz gravitational waves from the self-interacting axion cloud around the rotating stellar mass black hole
    Hirotaka Yoshino
    April 29, 2024
=============================================================================================================






§ INTRODUCTION


Black holes and gravitational waves offer a promising avenue to test physics beyond the standard model of particle physics, such as an axion. The interesting possibility suggested by string theory is that our universe may contain many species of axions, the so-called axiverse scenario <cit.>. In particular, axions in string theory typically have mass with their Compton wavelength around the astrophysical scale and the so-called axion decay constant, which determines the strength of the interaction involving the axion around the grand unification scale (∼ 10^16GeV) <cit.>. Besides the axiverse scenario, there are other motivations to assume the existence of axions. The QCD axion is needed to solve the strong CP problem <cit.>, while also axion is a dark matter candidate <cit.>.

Axions with the Compton wavelength about the size of astrophysical black holes (mass around 10^-20 - 10^-10eV) can be probed by observing the condensate of the axion around black holes <cit.>. The condensate of the axion is spontaneously formed by the superradiant instability <cit.>. Superradiant instability is a phenomenon such that the axion bounded by the gravitational potential exponentially grows by extracting the energy and angular momentum of the black hole. The rate of instability is fastest when the Compton wavelength of the axion is comparable to the size of the black hole. For the solar mass black hole, the rate is around one minute. As a result, the macroscopic number of axions will occupy the superradiant mode within the age of the universe.
In this paper, we call the condensate made solely by a single superradiant mode as an axion cloud and the superposition of axion clouds as an axion condensate.

One way to observe the axion condensate is to observe the gravitational waves emitted from it <cit.> (see for example <cit.> for other methods to observe axion condensates). The gravitational waves are emitted through two distinct processes: the pair annihilation of axions and the level transition between different superradiant modes. The characteristic feature of the emitted gravitational waves is that they are monochromatic. If a rotating black hole exists in our neighborhood, we have a chance to observe the gravitational waves from the associated axion condensate.

To better characterize the gravitational wave signal, a better understanding of the evolution of axion condensates under various effects, such as the axion self-interaction, is necessary (see  <cit.> for investigations on other effects such as tidal effects from the companion and the coupling to the photon). The self-interaction causes several effects  <cit.>, but the most significant one is the dissipation of the cloud by transferring the energy to another cloud  <cit.>. At a certain amplitude, the dissipation by the self-interaction is balanced with the growth due to the superradiant instability, resulting in a quasi-stationary configuration with various modes excited simultaneously.

So far, the evolution of the self-interacting condensate has been extensively investigated under the two-mode approximation with the non-relativistic approximation <cit.> or the three-mode approximation neglecting the spin down of the black hole <cit.>. These are not satisfactory to fully characterize the gravitational wave signals. The observationally attractive regime where the growth rate is the fastest, and the amplitude of the gravitational wave becomes large is in the relativistic regime. Furthermore, higher multipole modes are excited in the relativistic regime <cit.>. In addition, the inclusion of the spin down of the black hole is necessary since the spin of the black hole determines the lifetime of the condensate <cit.>.

In this paper, we numerically examine the co-evolution of the self-interacting axion condensate and the black hole, as well as the resulting gravitational wave radiation in the relativistic regime, including higher multipole modes. We use the adiabatic approximation to track the evolution, which is well-motivated in the current situation <cit.>. We find that the higher multipole modes are excited, but representative excited modes depend on the axion mass and the axion decay constant. In addition, we show that the excitation of the higher multipole modes will result in a characteristic gravitational wave signal different from the one predicted in Ref. <cit.>. More specifically, the transition, including the higher multipole modes, leads to a louder signal than the pair annihilation signal or the transition between the fastest and the second fastest growing modes for the GUT scale decay constant (10^17–10^16 GeV).

This paper is organized as follows.
In Sec. <ref>, we briefly review the axion cloud formed by the superradiant instability.
In Sec. <ref>, we formulate the evolution equations of the self-interacting axion condensate, including the multiple modes up to l=m=4 and the spin down of the central black hole.
In Sec. <ref>, we numerically solve the evolution equations formulated in Sec. <ref> and see how the excitation of the higher multipole modes proceeds.
In Sec. <ref>, we compute the gravitational and axion waves from the axion condensate.
Section <ref> summarizes our findings with some comments on the detectability of the gravitational waves.
In the rest of this paper, we take the unit with c=G=ħ = 1.



§ AXION CLOUD


We first review the formation of a condensate of an axion around a rotating black hole by the superradiant instability. For a detailed review of this topic, see <cit.>. In this paper, we employ the action

    S =   1/16 π∫ d^4 x √(-g) R       +  F_a^2 ∫ d^4 x √(-g){ - 1/2g^μν∂_μϕ∂_νϕ - V(ϕ)} ,

where we normalize the axion field  ϕ with the decay constant F_a. Note that F_a is measured in the unit of the Planck mass M_ pl = G^-1/2. Thus, we have F_a ∼ 10^-3 for the GUT scale decay constant. In this paper, we consider the leading order piece of the cosine type potential,

    V(ϕ) = μ^2 (1 - cosϕ) ∼μ^2/2ϕ^2 - μ^2/4!ϕ^4  .

The background metric is fixed to the Kerr metric given by

    ds^2    =  - (1 - 2 M r/ρ^2)dt^2 - 4 a M r sin^2θ/ρ^2 dt dφ        + [(r^2 + a^2)+ 2 M r/ρ^2 a^2 sin^2 θ]sin^2θ dφ^2            + ρ^2/Δ dr^2 + ρ^2 dθ^2 ,

with

    Δ   =(r-r_+)(r-r_-) ,    ρ^2    = r^2 + a^2cos^2θ ,

and

    r_±/M = 1  ±√(1 - χ^2) .

Here, M is the mass, and χ = a/M is the dimensionless spin parameter of the black hole. The radii r_+ and r_- correspond to the positions of the event horizon and the Cauchy horizon, respectively.

In the initial phase of the evolution, the amplitude of the axion is small, and we can neglect the quartic term in the potential (<ref>). In this case, the equation of motion for the axion is given by

    (□_g - μ^2)ϕ = 0 .

As shown in the literature <cit.>, we can solve the equation of motion by the separation of the variables,

    ϕ = e^-i(ω t - m φ)R_lmω(r) S_lmω(θ) +  c.c. .

Here, c.c. denotes the complex conjugate.
The radial part and the angular part of the equation become

    d/dr(ΔdR_lmω/dr)  + [K^2(ω)/Δ - μ^2 r^2 -λ_lm(ω) ]R_lmω = 0 ,


    1/sinθd/dθ(sinθdS_lmω/dθ)    + [c^2(ω) cos^2θ - m^2/sin^2θ]S_lmω          = - Λ_lm(ω) S_lmω ,

with

    c^2 (ω)    = a^2 (ω^2 - μ^2) ,      K(ω) = (r^2+a^2)ω - am ,λ_lm(ω)     = -2am ω +a^2ω^2 +Λ_lm(ω) ,

and Λ_lm(ω) is the separation constant that is reduced to l(l+1) in the limit a→ 0.  The normalization of the spheroidal harmonics S_lmω is taken to satisfy

    ∫ dcosθ S_lmω(θ)S_l'mω(θ) = δ_ll' .


Since we are interested in the bound states of the axion, we solve the radial equation (<ref>) with the boundary conditions

    R_lmω(r) →(r-r_+/M)^-i 2 M r_+/r_+ - r_-(ω - m Ω_H) ,     (r→ r_+) ,
    (r/M)^-1 -Mμ^2 - 2 ω^2/√(μ^2 - ω^2) e^-√(μ^2 - ω^2)r ,     (r→∞) ,

where  Ω_H = a/(r_+^2 + a^2) is the angular velocity of the event horizon. By the analytical <cit.> and numerical methods <cit.>, the solutions are found, and they are labeled by the three integers n,l, and m with n ≥ l+1, where n, l, and m is called the principal, azimuthal, and magnetic quantum numbers, respectively. We denote a solution as

    ϕ_nlm   = e^-i(ω_nlm t - m φ)S_lmω_nlm(θ)R_lmω_nlm(r) +  c.c. .

Here, ω_nlm is the angular frequency of the solution. We call the solutions with the smallest n  ( = l +1) as the fundamental modes and those with larger n as the overtone modes. Note that we normalize the radial mode function in such a way that the energy of the configuration calculated within the linearized approximation is equal to F_a^2 M.



The real and the imaginary parts of the frequency of the l=m=1 and l=m=2 fundamental modes are shown in Fig. <ref>. We observe that the frequencies possess the positive imaginary part, indicating the presence of instability. This instability is called superradiant instability. The superradiant instability occurs for <cit.>

    0<ω_nlm < m Ω_H .

 Because of the superradiant instability, the axion will occupy the superradiant states. As a result, the condensate of the axion will be spontaneously formed. In this paper, we denote the condensate made purely from a single superradiant mode as an axion cloud and the superposition of such clouds as an axion condensate.  



The radial configuration of the axion cloud is shown in Fig. <ref>. We observe that the configuration is quite similar to the radial wave function of the Hydrogen atom. In fact, in the small gravitational coupling regime, μ M ≪ 1,  the frequencies, and the mode functions coincide with those of the Hydrogen atom at a distant place <cit.>

    ω_nlm,R≡ω_nlm∼μ(1 - (μ M)^2/2n^2) ,


    S_lmω(θ)    ∼ Y_lm(θ) ≡√((l - m)! (2l+1)/2(l+m)!) P^m_l(cosθ) ,
    
    	R_lmω(r)    ∼ R_nl(r) ≡1/2√(π)μ M(2 (μ M)^2/n)^3/2√( (n - l - 1)!/ n^4
       (n + l)!)       ×(2 μ ^2 M  r/n)^l e^-μ ^2 M r/n L_n -l -1^2
       l+1(2μ ^2 M r/n) .

Here, L^k_n is the associated Laguerre polynomials.
By the asymptotic matching method, the approximate formula for the imaginary part of the frequency is found as

    
    ω_nlm,I≡ω_nlm∼   2 r_+/M(μ M)^4l+5(-μ+mΩ_H)2^4l+1(l+n)!/(n-l-1)! n^2l+4   ×(l!/(2l)!(2l+1)!)^2∏_j=1^l(j^2(1-a^2/M^2) + 4 r_+^2(μ-mΩ_H)^2) .


We refer to the formulas using the leading order approximations in μ M (Eqs. (<ref>) – (<ref>)) as the “non-relativistic approximation". The more accurate approximation can be found in Refs. <cit.>. As the figures show, the non-relativistic approximation breaks down in the μ M ∼ 1 regime. For example, the radial mode function becomes more compact than the one predicted by the non-relativistic approximation (see Fig. <ref>). 




§ EVOLUTION EQUATIONS OF THE AXION-BLACK HOLE SYSTEM


As the condensate grows, various effects would alter its evolution from the naive exponential growth predicted by the linear approximation in the previous section. The self-interaction <cit.> and the change of the spin of the central black hole <cit.> are particularly important. The self-interaction would change the evolution by causing the nonlinear mode coupling, which results in the saturation of the evolution <cit.>. 
In addition, the spin of the central black hole gradually slows down. 
Once the spin reaches a certain value, the growth of the condensate terminates because the superradiant condition (<ref>) is no longer met.

 Since the growth time scale of the cloud is much longer than its dynamical time scale and the light crossing time of the black hole, ω_nlm,I≪ω_nlm,R, M^-1_ BH, we can assume that the evolution is adiabatic. Under the adiabatic evolution, at each moment, the metric is approximated by the Kerr metric with a given black hole mass M_ BH(t) and an angular momentum J_ BH(t). The evolution of each axion cloud is captured by that of the normalized mass, M_ cl,i(t), where i specifies the mode. Since we have normalized the axion by the decay constant F_a, the actual mass of the cloud is expressed by the normalized mass as F_a^2 M_ cl,i(t).

Evolution equations for the clouds are derived from the local conservation laws of the energy and the angular momentum <cit.>. The results are summarized in the following form

 
    d M_ cl,i/dt =   
          + 2 ω_i,IM_ cl,i      - ∑_j,k(1 + δ_ij)ω_i,R/ω_i,R + ω_j,R -ω_k,R F^ℋ_ijk^*M_ cl,i M_ cl,jM_ cl,k + ∑_j,kω_i,R/ω_j,R + ω_k,R -ω_i,R F^ℋ_jki^*M_ cl,i M_ cl,jM_ cl,k     - ∑_j,k(1 + δ_ij)ω_i,R/ω_i,R + ω_j,R -ω_k,R F^ℐ_ijk^*M_ cl,i M_ cl,jM_ cl,k + ∑_j,kω_i,R/ω_j,R + ω_k,R -ω_i,R F^ℐ_jki^*M_ cl,i M_ cl,jM_ cl,k .


The first term represents the increase in the energy due to the superradiance.
The coefficients F^ℋ_ijk^* and F^ℐ_ijk^* are defined by 

    .dE_a/dt|_r=r_+ = F_a^2 F^ℋ_ijk^* M_ cl,iM_ cl,jM_ cl,k ,
    .dE_a/dt|_r→∞  = F_a^2 F^ℐ_ijk^* M_ cl,iM_ cl,jM_ cl,k ,

which represent the magnitude of the energy flux of the axion to the horizon and that to infinity, respectively, induced by the process involving modes i,j, and k. Note that F_ijk^*^ℋ,ℐ = F_jik^*^ℋ,ℐ, since the modes i and j are interchangeable in this process.
These processes are caused by the self-interaction, 
and see Fig. <ref> for the meaning of their subscript.
The terms with the Kronecker delta δ_ij are necessary to take into account that two axions are lost by a single transition when i = j.
See appendix A of <cit.> for details of the derivation. 
The factor F_a^2 comes from the normalization of the axion field.



The evolution of the black hole mass and the angular momentum is governed by the 
balance equations  <cit.>,


    d M_ BH/d t   = F_a^2(-∑_i 2 ω_i,IM_ cl,i +∑_(i,j,k) F^ℋ_ijk^*M_ cl,i M_ cl,jM_ cl,k)  ,
    d J_ BH/d t   = F_a^2( -∑_i 2 m_i ω_i,I/ω_i,R M_ cl,i+ ∑_(i,j,k)m_i + m_j - m_k/ω_i,R + ω_j,R - ω_k,RF^ℋ_ijk^*M_ cl,i M_ cl,jM_ cl,k) ,


where the right hand sides are the energy and angular momentum fluxes through the horizon, respectively. 
The first term in the parenthesis corresponds to the mass and the angular momentum loss (gain) due to the superradiance or absorption of axions. The second term corresponds to the energy and the angular momentum flux induced by the self-interaction (the left panel of Fig. <ref>). 



 §.§ Subleading effects


There are subleading effects that we do not take into account in the evolution. One is the other effect due to the self-interaction. As summarized in Ref. <cit.>, the effect of the self-interaction can be classified into three parts. One is the dissipation due to the low-frequency modes ω∼μ, which is considered in this paper. The second one is the dissipation due to the high-frequency modes ω≳ 2μ. These contributions are small when considering the evolution in the perturbative regime <cit.>. The third one is the deformation of the cloud. When the amplitude of the cloud is large, the condensate is deformed significantly <cit.>, and for a large amplitude, it leads to the collapse of the condensate <cit.>. However, it has been shown that the dissipative effects due to the self-interaction are large enough to prevent the condensate from entering the collapsing regime in the case the condensate starts with a small amplitude, e.g., that produced by the quantum fluctuation <cit.>. For these reasons, we only consider the dissipative effect due to the low-frequency radiation. 

 According to the cloud deformation, the eigenfrequencies of the cloud would change. It becomes essential when we consider the long-term observation of the continuous gravitational waves from the axion cloud <cit.>. When the amplitude of the cloud is not so large, the frequency shift can be calculated by the perturbative method both in the non-relativistic approximation <cit.> and in the relativistic approach <cit.>. When the amplitude of the cloud is large, it should be calculated numerically as in Ref. <cit.>. In this paper, we do not compute the frequency shift for ease of calculation, 
since the frequency shift does not change the evolution significantly. 

In addition, we neglect the gravitational wave radiation and self-gravity of the cloud. The energy loss rate of clouds due to the gravitational wave radiation is much smaller than the energy transfer rate due to the superradiant instability <cit.> and the scalar wave radiation induced by the self-interaction. For this reason, we neglect the effect of gravitational waves on the evolution.
The correction to the axion condensate due to the self-gravity is also small <cit.>.




§ TIME EVOLUTION WITH FOUR MODES


Now, we numerically solve the evolution equations (<ref>), (<ref>), and (<ref>). We analyze the situation in which the four fundamental modes l=m=1,l=m=2,l=m=3, and l=m=4 are excited simultaneously. 
In our previous paper <cit.>, it has been shown that the amplitudes of overtone modes are always subdominant. In the following we use i =1,2,3,4 to label the fundamental modes with (l,m) = (1,1),(2,2),(3,3),(4,4).  

Our four mode approximation would give a correct evolution unless the l=m=3 or the l=m=4 mode grows to a significantly large amplitude. In the case of the l=m=3 or 4 modes growing to a large amplitude, we need to include much higher multipoles. As we shall see below, the substantial growth of the l=m=3 or 4 mode only happens for a much later time when the l=m=1 and 2 modes have decayed. Therefore, when considering the gravitational wave signal involving the l=m=1 or l=m=2 mode, our four-mode approximation should give an accurate result.


When considering these modes, the fluxes of the processes that we need to take into account are the following eighteen (ten from the processes dissipating to the black hole and eight from the ones dissipating into infinity):


  Flux to the horizon


    F^ℋ_112^* ,  F^ℋ_113^* ,  F^ℋ_223^* ,   F^ℋ_114^* , F^ℋ_224^* ,    F^ℋ_334^* , F^ℋ_123^* , F^ℋ_124^* , F^ℋ_134^* ,  F^ℋ_234^* .



  Flux to infinity


    F^ℐ_221^* ,  F^ℐ_331^* ,  F^ℐ_441^* , F^ℐ_442^* ,  F^ℐ_231^* ,  F^ℐ_241^* ,  F^ℐ_341^* ,    F^ℐ_342^* .


Other processes become subdominant because of the following reasons: 1) the frequencies of the dissipative modes are much larger than the axion mass, or 2) for the modes dissipating to the black hole, the angular momentum barrier of the dissipative mode is large enough (l > 2)

 <cit.>. Note that when considering the flux to infinity, its frequency must be larger than the axion mass.

When solving the evolution equations, we need the numerical values of ω_i,  F^ℋ, I_ijk^* and so on for given (μ M,χ). We numerically calculate these quantities on the grid points in the parameter space of (μ M, χ) in the range between 0.02 ≤μ M ≤ 0.45 and 0.15 <χ < 0.998 with interval Δ (μ M) =0.01 and Δχ = 0.025. For the calculation of the time evolution, we use the value of ω_i and F^ℋ, I_ijk^* obtained by interpolating the values on the grid points. The whole calculation is done using Mathematica.



 §.§ Example of the time evolution


To demonstrate the behavior of solutions of Eqs. (<ref>), (<ref>), and (<ref>), we fix the initial mass of the black hole to M_ BH, ini = 10 M_⊙ in this section. Then, the superradiant instability of the l=m=1 mode occurs for the axion mass in the range μ≲ 6× 10^-12 eV. We take the initial spin of the black hole to be relatively large χ = 0.99. The initial mass of the cloud is taken to be M_ cl, i = 10^-10 M_ BH,ini, which is just set to a sufficiently small value. In reality, the cloud mass will start with a much smaller value, but in such a small mass (M_ cl, i∼ 10^-80 M_⊙ if the cloud starts with the quantum fluctuation with GUT scale decay constant), self-interaction is totally negligible, and hence our result should give the qualitatively correct results, even if we do not start with such an extremely small number. For the detail of the dependence of the evolution on the initial condition, see Ref. <cit.>.



We first show the evolution of the cloud for the axion mass with μ M_ BH, ini = 0.2, and F_a = 10^-3 in Fig. <ref>. This axion mass corresponds to μ∼ 2.6 × 10^-12 eV. The evolutionary track of the cloud is summarized as follows:


  * The l=m=1 cloud grows exponentially by the superradiant instability.

  * The energy of the l=m=1 cloud is transferred to the l=m=2 cloud. Then, the l=m=2 cloud grows exponentially. At a certain amplitude, the energy gained by the superradiance and the dissipation by the self-interaction balance result in a quasi-stationary configuration. The mass of each mode is roughly determined by solving the stationary conditions

    dM_ cl,i/dt   = 0 ,     (i   =1,2,3,4) .


Note that, around the maximum of the cloud mass, we cannot neglect the deformation of the mode function due to the self-interaction <cit.>. The clouds become compact at a large amplitude, and the frequencies of the clouds would be modified. However, the numerical analysis shows that the qualitative evolution does not change from the perturbative calculation.  The change in the cloud configuration might enhance the gravitational wave amplitudes, but we neglect this effect for simplicity. The modification to the frequencies is important for the actual observation, but we also neglect it.

  * Higher multipole modes are excited by the same mechanism. In the present example, the l=m=3 cloud is excited. Which modes are excited depends on the parameters of the axion, (μ, F_a), as we briefly discuss below.

  * The l=m=1 cloud is depleted because of the spin down of the black hole by the superradiance of the l=m=2 mode. As the obstruction due to the presence of the l=m=1 cloud disappears, the l=m=2 mode starts to grow by the superradiant instability. 

  * Once the energy of the l=m=2 cloud dominates, the energy transfer from  the l=m=2 cloud to the l=m=3 cloud becomes efficient. 
Since there is no process to transfer the energy of the l=m=3 mode to the l=m=2 mode, the l=m=3 mode continues to grow. 
For this case, we need to include much higher multipoles, such as l=m=6 mode.
Then, in the end, the energy is transferred to much higher multipoles and settles into the quasi-stationary state.
 



Now, let us vary the parameters of the axion, F_a and μ. In Fig. <ref>, we show the case with only the value of F_a changed from that of Fig. <ref>. Changing the value of F_a is identical to changing the relative size of the black hole spin down rate and the energy transfer rate of the cloud due to the self-interaction. Increasing F_a weakens the self-interaction and increases the relative spin down rate while decreasing F_a does the opposite. The top panel of Fig. <ref> is the case when F_a is increased. We observe that the saturation of the l=m=1 cloud occurs much before the excitation of the l=m=2 cloud. This saturation occurs because the black hole has spun down to the superradiance threshold of Eq. (<ref>). Although the saturation due to the spin down occurs, the energy transfer from the l=m=1 cloud to the l=m=2 cloud still occurs, leading to much earlier excitation of the l=m=2 cloud than the excitation by the superradiant instability alone. This regime corresponds to the “small self-coupling regime" in Ref. <cit.>. On the other hand, we observe a delay in the decay of the l=m=1 cloud when we decrease F_a (bottom panel of Fig. <ref>). This regime corresponds to the “large self-coupling regime" in Ref. <cit.>.





We now change the axion mass μ while the value of F_a is fixed as 10^-3. Figure <ref> shows the cases with μ M_ BH,ini = 0.1 and 0.4, which correspond to μ∼ 1.3 × 10^-12 and 5.2 × 10^-12eV, respectively. The behavior is similar to the case with μ M_ BH,ini = 0.2, except that the l=m=3 cloud is no longer excited before the depletion of the l=m=1 cloud. Since the simultaneous excitation of higher multipole modes will be important when we consider gravitational wave signals (Sec. <ref>), we discuss the parameter range in which the simultaneous excitation of higher multipole modes is possible. We define the effective growth rate of the l=m=3 mode, when the l=m=1 and the l=m=2 modes are in the quasi-stationary state as

    ω_3,I^ eff   ≡ 2ω_3,I + F^ℋ_113^*M_ sat,1^2 +  F^ℋ_223^*M_ sat,2^2  +    (F^ℋ_123^* - F^ℐ_231^*)M_ sat,1 M_ sat,2 .

Here, M_ sat,1 and M_ sat,2 are the masses of the l=m=1 and l=m=2 clouds when the growth is saturated by the self-interaction between the two clouds. We define ω_4,I^ eff in the same way. The region where ω_3,I^ eff and ω_4,I^ eff are positive is shown in Fig. <ref>.
The excitation of the l=m=3 mode is possible when 0.12 ≲μ M_ BH≲ 0.25, and the excitation of the l=m=4 mode is possible when 0.20 ≲μ M_ BH≲ 0.45. In the overlapping regime, the one with the larger effective growth rate will be excited, and the other mode is depleted by the processes F^ℐ_341^* and F^ℐ_342^*. In the most of the overlapping region, the effective growth rate for l=m=3 is larger than that for l=m=4.



Finally, let us comment on the evolution of the black hole spin. In Fig. <ref>, we show the evolution of the black hole spin with various μ and F_a. Since the superradiant rate depends on the axion mass μ, the time when the significant spin down occurs depends on μ. In addition, the spin down time scale becomes longer when one decreases F_a. This is evident from the equation governing the angular momentum evolution of the black hole, Eq. (<ref>). The spin down is initially driven by the l=m=1 mode and later by higher multipole modes because of the difference of the threshold spin and the hierarchy in the rate of the superradiant instability (see right panel of Fig. <ref>). However, we note that the spin down due to the higher multipole modes starts to function much earlier compared with the non-interacting case as long as the decay constant satisfies F_a ≲ 10^-1, since the growth of higher multipole modes is accelerated by the mode coupling due to the self-interaction. The spin up observed in the figure is due to the fallback of the lower multipole modes to the black hole.




§ OBSERVABLE SIGNALS






 §.§ Gravitational waves




  §.§.§ Calculation of the gravitational waves


Here, we consider gravitational waves emitted by the axion condensate.  
Gravitational waves are generated from the time-dependent part of the axion energy-momentum tensor T_μν(ϕ). When the multiple modes are excited simultaneously, the axion field is given in the form [Note that we normalized the axion field by F_a, so ϕ is non-dimension. In addition, M_ cl,i is normalized by F_a^2 M_ BH, so also nondimension.]

    ϕ = ∑_i√(M_ cl,i/M_ BH) e^-i(ω_i t - m_iφ)R_l_im_iω_iS_l_im_iω_i +  c.c. .

Since T_μν is quadratic in ϕ, the time dependent pieces behaves as ∝ e^-i(ω_i + ω_j) t or e^-i(ω_i - ω_j) t. 
Namely, the frequencies of gravitational waves are either the sum or the difference of the frequencies of two clouds. The former is called the “pair annihilation signal," while the latter is the “level transition signal."  The actual frequencies of the gravitational waves for the 10M_⊙ black hole mass with spin χ = 0.99 are shown in Fig. <ref>. The pair annihilation signal has a much higher frequency (∼ 1  kHz) than the level transition signal (∼ 1  Hz).

The energy flux and the amplitude of gravitational waves can be calculated by solving the Teukolsky equation. We briefly review the calculation of gravitational waves from the axion condensate. The details of the calculations can be found in Refs. <cit.>. 

The energy flux F_ gw and the amplitude h_ gw of a monochromatic gravitational wave 
with the angular frequency ω in a distant region from the black hole are described in terms of the Newmann-Penrose variable ψ_4 as <cit.>

    F_ gw   = lim_r→∞r^2/4πω^2∫ dΩ |ψ_4|^2 ,
    
    	h_ gw   = 2/ω^2|ψ_4| .

When the source of the gravitational wave has the form e^-i (ω t - mφ), the asymptotic form of the ψ_4 is shown to be
	
    ψ_4    →1/re^-iω t + i mφ e^i ω r∑_l=2^∞ Z_l^ out_-2S_lmω ,     (r →   ∞) ,

	where the amplitude is given by
	
    Z^ out_l   = 1/_-2W_lmω∫_r_+^∞ dr'  Δ^-2_-2R^ in_lmω(r') _-2T_lmω(r') .

	Here, _-2S_lmω is the spin-weighted spheroidal harmonics with the spin weight s = -2, _-2R_lmω^ in is the solution to the radial Teukolsky equation with the ingoing boundary condition at the event horizon, and _-2W_lmω is the Wronskian of the radial Teukolsky equation. The function _-2T_lmω is the source term of the radial Teukolsky equation, related to the energy-momentum tensor of the axion. The precise form of _-2T_lmω can be found in Ref. <cit.>. Note that _-2T_lmω depends on the mass of the cloud and the decay constant as ∝ F_a^2 √(M_ cl, iM_ cl,j), since the energy-momentum tensor is quadratic in ϕ. Therefore, the energy flux and the amplitude scale as

    F_ gw∝    F_a^4 M_ cl,iM_ cl,j ,
    
    	h_ gw∝    F_a^2 √(M_ cl,iM_ cl,j) .



	
In the left panel of Fig. <ref>, we show the energy flux of gravitational waves from a black hole with the mass 10M_⊙ and the spin χ = 0.99. In terms of the energy flux, the pair annihilation signal of the l=m=1 cloud and the level transition signal of the transition from the l=m=3 mode to the l=m=1 mode dominate if the masses of all clouds are assumed to be the same. The energy dissipation of the cloud due to the pair annihilation is much smaller than the energy gain by the superradiance <cit.>, which justifies neglecting the gravitational wave emission when we consider the evolution of the condensate.

In the right panel of Fig. <ref>, we show the amplitudes of respective gravitational waves. We observe the hierarchy among the amplitudes generated by different processes. First of all, the amplitude of the pair annihilation signal is the smallest. This is the case because the wavelength of the gravitational wave (∼ 1/μ) is much shorter than the size of the condensate (∼ 1/μ/(μ M)) in the situation we are considering. In contrast, the level transition signal is not suppressed since the wavelength is much longer (∼ 1/μ/(μ M)^2). The size of the level transition signal is the largest for the transition with |m_i - m_j| = 2, followed by |m_i - m_j| = 0, and then by |m_i - m_j| = 1. The transition with |m_i - m_j| = 1 is suppressed because the transition cannot occur through the mass quadrupole radiation because of the parity argument. The relative suppression of the signal with |m_i - m_j| = 0  to the one with |m_i - m_j| = 2 is due to the axisymmetry of the source for the case of |m_i - m_j| = 0. Note that the amplitude of the gravitational waves does not change significantly as one changes the black hole spin. This is because the configuration of the cloud does not depend on the black hole spin prominently.

 Other than the signals presented in Fig. <ref>, it is possible to consider pair annihilation between different modes. However, one can confirm that this “cross annihilation" process produces much smaller gravitational waves since it starts with the octopole radiation rather than the quadrupole one.



To summarize, we expect a much larger signal than the one previously considered, when
there simultaneously exist two modes whose magnetic quantum numbers are different by two, e.g., the l=m=1 and l=m=3 modes or the l=m=2 and l=m=4 modes. Let us estimate the amplitude of the gravitational wave emitted during the quasi-stationary state that is found in the previous section. In Fig. <ref>, we show the normalized masses of the clouds M_ cl, i when they settle into the quasi-stationary state. We observe that when the l=m=3 cloud is excited, the relation M_ cl,1∼ M_ cl,3∼𝒪(10^2)M_ BH holds regardless of the spin. The maximum of √(M_ cl,1 M_ cl,3) is ∼ 380M_ BH for χ = 0.99 and ∼ 290 M_ BH for χ = 0.7, around μ M_ BH∼ 0.15. From the numerical calculation, the value of  h_ gw around μ M_ BH∼ 0.15 normalized as in Fig <ref> is determined to be ∼ 3 × 10^-18 for both  χ = 0.99 and 0.7. Thus, the maximum amplitude of the level transition signal from the l=m=3 mode to the l=m=1 mode is given by 


    h_ gw∼
    	1.0 × 10^-21(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (√(M_ cl,1M_ cl,3)/380M_ BH) ,     (χ = 0.99, μ M_ BH∼ 0.15) ,
    
    	7.9 × 10^-22(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (√(M_ cl,1M_ cl,3)/290M_ BH) ,     (χ = 0.7, μ M_ BH∼ 0.15) .


For a black hole with mass 10M_⊙, the frequency of this gravitational wave is around 1 Hz, which is in the target frequencies of the future gravitational wave detectors such as DECIGO <cit.>, atomic interferometers <cit.>, TOBA <cit.>, TianGo <cit.>, DO <cit.>, AMIGO <cit.>, and LGWA <cit.>. The amplitude does not decrease significantly even if the spin of the black hole is reduced.



The amplitudes of gravitational waves from other processes can be estimated similarly. The maximum amplitude is given by


    h_ gw∼
    	1.0× 10^-24(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (M_ cl,1/36 M_ BH) ,     (χ = 0.99, μ M_ BH∼ 0.42) ,
    
    	2.1 × 10^-25(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (M_ cl,1/248 M_ BH) ,     (χ = 0.7, μ M_ BH∼ 0.2) ,

 for the l=m=1 pair annihilation, and 

    h_ gw∼
    	4.0 × 10^-23(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (√(M_ cl,1M_ cl,2)/54M_ BH) ,     (χ = 0.99, μ M_ BH∼ 0.42) ,
    
    	4.8 × 10^-24(1 kpc/r)(M_ BH/10M_⊙) (F_a/10^-3)^2 (√(M_ cl,1M_ cl,2)/64 M_ BH) ,     (χ = 0.7, μ M_ BH∼ 0.2) ,


 for the transition from the l=m=2 cloud to the l=m=1 cloud.
Thus, we find that the transition from the l=m=3 cloud to the l=m=1 cloud can provide the strongest signal for the GUT scale decay constant.

We show the typical time evolution of the gravitational wave amplitude for these processes in Fig. <ref>. After the excitation of the l=m=1 cloud, the gravitational waves from the pair annihilation in the l=m=1 cloud and the level transition between the l=m=1 and the l=m=2 clouds are radiated. This phase lasts for the time scale of superradiant instability, which is ≲ 1 yr in the current choice of the parameters. After a few years, the larger gravitational wave from the level transition between the l=m=3 and the l=m=1 clouds appears. The amplitude of this gravitational wave decays as the l=m=1 cloud depletes due to the spin down. 






  §.§.§ Dependence on the axion parameter


In the previous section, we estimated the amplitudes of gravitational waves for the GUT scale decay constant. However, the result does not scale simply as one changes the parameters, μ, F_a, the initial black hole mass M_ BH,ini and the spin χ_ BH, ini. 
Recall that the dependence of the evolution on these parameters is complicated  (see Sec. <ref>). Here, we investigate how the peak amplitude h_ peak and the duration of the signal τ_ gw depend on the parameters mentioned above by solving the evolution equation directly for various cases. We estimate the duration by the full width at the half-maximum (FWHM). 

As the initial parameters of the black hole, we choose the one similar to the Cygnus X-1 <cit.>, which has been proposed to be an interesting target for the axion search <cit.>, or to black holes produced from the compact binary coalescence such as GW170817 <cit.>. The corresponding parameters are M_ BH,ini = 14.8M_⊙, χ_ ini=0.99, and d = 1.82 kpc for Cygnus X-1-like black holes and  M_ BH,ini = 3 M_⊙, χ_ ini=0.7, and d = 40 Mpc for the black holes from the compact binary coalescence, respectively. Here, d is the distance to the black hole.

In Fig. <ref>, we show the peak amplitude and the duration of the gravitational wave from the pair annihilation of the l=m=1 cloud on the (μ,F_a) plane. The amplitude of the gravitational wave decreases as one decreases μ or F_a. The amplitude dependence on μ is evident from Fig. <ref>. 
The dependence on F_a is the consequence of Eq. (<ref>). The cut-off around μ M_ BH,ini∼ 0.2 observed in the lower left panel of Fig. <ref> is reflecting the superradiant condition for χ = 0.7 (see Eq. (<ref>)). Beyond this axion mass, the l=m=1 cloud will never be excited; thus, there is no emission of the pair annihilation gravitational wave.



Next, we consider the dependence of the duration on the parameters μ and F_a. For decay constant in the range F_a ≳ 10^-2, the duration becomes longer as one increases F_a or decreases μ. In this range, the duration is determined by the time when the l=m=1 cloud decays, which is determined by the spin down of the black hole due to the superradiance assisted by the presence of the l=m=2 cloud. Increasing the decay constant leads to relative suppression of the excitation of the l=m=2 cloud, which means that the spin down becomes slower. For a smaller axion mass, the superradiant instability rate is reduced, which leads to a slower spin down rate. The duration depends primarily on the axion mass for F_a ≲ 10^-2. In this region, the excitation of the l=m=2 cloud happens before the spin down to the threshold for the superradiant instability of the l=m=1 cloud completes. Duration determined by the FWHM is short since a sharp peak occurs for the amplitude of gravitational waves from the pair annihilation in the middle of the transition process (see Fig. <ref>).





Figure <ref> shows the case for the level transition between the l=m=2 and the l=m=1 clouds. The basic feature is similar to the pair annihilation signal from the l=m=1 cloud. The main difference is in the dependence of the amplitude on μ, which is much weaker than the pair annihilation signal. Another difference is in the size of the amplitude, which is about one order of magnitude larger for fixed (μ, F_a). The signal duration has the same order of magnitude as the l=m=1 pair annihilation process.

The case with the level transition process between the l=m=3 and the l=m=1 clouds is shown in Fig. <ref>. In contrast to the pair annihilation or the level transition between the l=m=2 and the l=m=1 clouds, the emission occurs for the narrow parameter range, i.e., F_a ≲ 10^-3 and 0.12 ≲μ M ≲ 0.24 (0.2) for χ = 0.99 (0.7). The region is determined by the condition that the l=m=3 mode is excited before the l=m=1 mode decays. The l=m=3 cloud is not simultaneously excited in the large parameter range of F_a since the spin down of the black hole is too fast and the l=m=1 cloud decays before the l=m=3 cloud appears (see Fig. <ref>). The region for the μ M agrees with the blue shaded region of Fig. <ref>.

Finally, we consider the case with the level transition process between the l=m=4 and the l=m=2 clouds (Fig. <ref>). We cannot expect the signal from this process for some parameter regions: small μ and F_a region, and the large μ region (only for χ_ ini = 0.7, see the region μ≳ 1.8× 10^-11 eV). The small μ and F_a region corresponds to the region where the l=m=3 cloud is excited. As pointed out in Sec. <ref>, when the l=m=3 cloud is excited, the l=m=4 cloud decays due to the mode coupling. The right edges in the lower two panels of Fig. <ref> correspond to the threshold of the superradiance condition of the l=m=2 mode.

The duration is determined in a complicated way (see right panels of Fig. <ref>). In the case of a relatively small initial spin (lower right panel of Fig. <ref>), the behavior is similar to the level transition between the l=m=2 and the l=m=1 clouds in the sense that the duration becomes longer as F_a is increased and μ is decreased (see right panels of Fig. <ref>).

Behavior is complicated for a larger initial spin (upper right panel of Fig. <ref>). We observe the transition around F_a ∼ 10^-3 and F_a ∼ 10^-4. This transition is determined by the competition between the decay of the l=m=1 cloud and the excitation of the l=m=4 cloud. In the large F_a regime, the black hole spin down is so fast that the l=m=1 cloud disappears before the excitation of the l=m=4 cloud. Therefore, in this case, the transition signal is emitted after the complete decay of the l=m=1 cloud and the black hole spin down.

As F_a is decreased, the spin down time scale due to the superradiance of the l=m=1 cloud becomes longer, and the transition signal between the l=m=2 and the l=m=4 clouds is emitted before the depletion of the l=m=1 cloud until F_a ∼ 10^-3. But still, the configuration is almost determined by the stationary conditions for the l=m=2 and the l=m=4 clouds, and the l=m=1 cloud is subdominant. The presence of the l=m=1 cloud reduces the peak amplitude by assisting the growth of the l=m=4 through the processes F_114^*^ℋ and F_124^*^ℋ. The earlier growth of the l=m=4 cloud results in the smoothing of the sharp peak.  Then, the peak amplitude becomes smaller than twice the amplitude at the plateau. 
In this case, the signal duration is determined by the spin down due to the l=m=2 superradiance, which greatly enhances the signal duration defined by the FWHM. 
For F_a ≲ 10^-4, the excitation of the l=m=4 cloud occurs before the depletion of the l=m=1 cloud. Then, the configuration is determined by the balance among the three modes, not between the two. Now, the signal duration is determined by the spin down due to the superradiance of the l=m=1 mode, which is much faster than the spin down due to the l=m=2 mode.
When we further decrease F_a, the spin down time scale becomes longer, and the duration increases.



 §.§ Axion wave


Besides gravitational waves, the self-interacting axion condensate will emit the axion to infinity. Since the mass of the axion is ultra-light and the frequency is close to the axion mass, the emitted axion has properties similar to the ultra-light dark matter. Here, we estimate the amplitude and the axion number flux from the saturated configuration. 

From the definition of F^ℐ_ijk^*, the energy flux of the axion is given by

    dE_a/dt = F_a^2 F^ℐ_ijk^* M_ cl,iM_ cl,jM_ cl,k ,

for each process. The axion field responsible for the energy flux (<ref>) behaves near infinity as 

    ϕ∼  A_∞ e^-i (ω t -m φ) S_mmω(θ) e^ i k r/r  +  c.c. ,

with

    A_∞   = √(F^ℐ_ijk^*M_cl,iM_cl,jM_cl,k/4 πω k) ,
    
    	k   =√(ω^2 - μ^2) .

 Here, the frequency and the magnetic quantum number of the axion are given by ω = ω_i + ω_j - ω_k and m = m_i + m_j - m_k, respectively, and used the fact that the contribution from the l=m modes dominates the condensate in the index of the spheroidal harmonics. The amplitude is converted to the axion number flux per unit area as

    Φ_a  ∼1/4 π r^2 ωdE_a/dt  = F_a^2 F^ℐ_ijk^* M_ cl,iM_ cl,jM_ cl,k/4 π r^2 ω .





When the l=m=1 cloud is present, the energy flux is dominated by the two processes, i.e., 2×|322⟩→|211⟩+() and |322⟩+ |433⟩→|211⟩ + (). In Fig. <ref>, we show the amplitude and the axion number flux for these processes, assuming the black hole with M_ BH = 10M_⊙ and χ=0.99 at 1kpc away from us. Note that if we assume the axion to be a dominant component of the dark matter, the amplitude (normalized by F_a) and the number flux in the solar system are estimated as

    |ϕ_ DM| ∼√(ρ_ DM/F_a^2 μ^2)   ∼ 2.6 × 10^-15(ρ_ DM/0.4  GeV/cm^3)^1/2       ×(10^-3/F_a)(10^-12 eV/μ) ,
    Φ_ DM∼v_ DMρ_ DM/μ   ∼ 1.2× 10^28[ cm^2/s](ρ_ DM/0.4  GeV/cm^3)       ×(v_ DM/10^-3)(10^-12 eV/μ) .

Therefore, the amplitude of the axion from the black hole that we consider here becomes smaller than that of the ultra-light dark matter, but the number flux can be comparable. This is because the velocity of the emitted axion (∼μ M_ BH) is one order of magnitude larger than that of the dark matter particle (v_ DM∼ 10^-3). Since the amplitude is not so small as the dark matter particle, we might have a chance to observe the axion emitted from the black hole <cit.>. We leave further investigation of the detectability of the axion wave from the black hole as a future work.



§ CONCLUSION


In this paper, we have numerically investigated the evolution of the self-interacting axion condensate and the observational signals of emitted gravitational waves in the relativistic regime (μ M ≳ 0.1), including the interaction between higher multipole modes up to l=m=4. We have found that the excitation of the higher multipole modes heavily depends on the axion mass μ and the axion decay constant F_a. These parameters control the relative size between the spin down rate and the energy transfer rate of clouds due to the self-interaction. As a result, the l=m=3 cloud can be excited only for the narrow parameter range: 0.12 ≲μ M ≲ 0.25 and F_a ≲ 10^-3. Except this parameter range, the l=m=4 cloud can be excited.

Owing to the presence of the higher multipole modes, we would expect gravitational wave signals in several frequency bands. The pair annihilation gives the signal in the high-frequency regime (kHz for a black hole with the mass 10M_⊙), and the level transition process gives the low-frequency signal (Hz for a black hole with the mass 10M_⊙). In this paper, we have shown that the gravitational waves from the level transition between the modes whose magnetic quantum numbers are different by two (i.e. Δ m = 2) can provide a larger amplitude than the pair annihilation signal or the level transition signal with Δ m = 1 (see Eqs. (<ref>)– (<ref>)). We have also estimated the flux of axion waves emitted from the condensate and have shown that it can be comparable to the flux of the dark matter.

Let us briefly consider the detection prospect of these level transition signals from black holes made by a binary neutron star merger. The frequency of the gravitational wave would be around Hz, which is the target of future detectors such as DECIGO <cit.>, atomic interferometers <cit.>, TOBA <cit.>, TianGo <cit.>, DO <cit.>, and AMIGO <cit.>. For simplicity, let us focus on the case of DECIGO. The current sensitivity target of DECIGO is ∼ 4× 10^-24 Hz^-1/2 at around 1 Hz <cit.>. If we could integrate the signal for 1 yr, we could observe the signal with an amplitude of h ∼ 7 × 10^-28. Therefore, the gravitational wave signal of the transition between the l=m=3 mode and the l=m=1 mode could be observed assuming the GUT scale decay constant and 40 Mpc as the distance to the source (see lower panels of Fig. <ref>). 
We can observe signals from black holes which are located much farther if we consider the heavier black holes. For example, we can observe it up to 800 Mpc for a black hole with the mass 60M_⊙ (such as the one produced in the GW150914 <cit.>). For the larger decay constant (F_a ≳ 10^-2), the dominant signal will be the gravitational wave from the level transition between the l=m=4 and l=m=2, which can be reached by the detectors with target sensitivity ∼ 10^-21 Hz^-1/2 around 1 Hz (see lower panels of Fig. <ref>).

However, to be more realistic, this estimation requires further study. We would need to calculate the frequency shift and fix the search strategy to estimate the detectability. The calculation of the frequency shift in the relativistic regime is already formulated in the literature <cit.> but requires some further computation. Hence, we leave the detailed investigation of the detectability of the level transition signals to future work.

In this paper, we only considered the stellar mass black hole and did not consider the supermassive black holes. A huge signal would be expected if we scale our result to the supermassive black holes due to their large mass. However, we cannot simply scale our result since the environmental effect must be addressed. We must especially consider the effect of the accretion. The accretion would spin up the black hole, which would affect the evolution of the cloud <cit.>. Furthermore, the tidal effect could potentially deplete the cloud <cit.>. We leave the investigation on the supermassive black hole for the future.


0
Note that <cit.> estimated that the advanced LIGO has a sensitivity to the continuous wave signal from the pair annihilation signal with h_ peak≳ 5 × 10^-26 and τ_ gw≳ 1 yr at the most sensitive frequency. 





0






H.O. thanks Jiro Soda and Hiroki Takeda for a valuable discussion.
H. O. is supported by JSPS KAKENHI Grant Numbers JP23H00110.  
H. Y is in part supported by JSPS KAKENHI Grant Numbers JP22H01220 and JP21H05189,
and is partly supported by Osaka Central Advanced Mathematical Institute
(MEXT Joint Usage/Research Center on Mathematics and Theoretical Physics JPMXP0619217849). T. Takahashi is supported by JSPS KAKENHI Grant Number JP23KJ1214. T.T. is supported by JSPS KAKENHI Grant Nos. JP24H00963, JP24H01809, JP23H00110 and JP20K03928.




apsrev4-1


