BLACK HOLE QUASINORMAL MODE SPECTROSCOPY WITH LISA

An important source of gravitational waves in the frequency range 10⁻⁴–10⁻¹ Hz, which is the frequency band for the planned space-borne gravitational-wave detector Laser Interferometer Space Antenna (LISA), is the inspiral and merger of two supermassive black holes. The emitted gravitational waveform includes three major parts, corresponding to the inspiral, merger, and ringdown of the system. The last stage, the ringdown, is the settling down of the resulting black hole to quiescence as required by the "no-hair" theorem, and is characterized by (an incomplete) set of complex frequencies that depend only on the macroscopic parameters of the black hole, namely, its mass M and spin angular momentum J (Nollert 1999). The computation of the complex frequency is done using black hole perturbation theory (Chandrasekhar 1983), and its detection would be a smoking gun for the source being a black hole, and an accurate means for determining its parameters (Berti et al. 2006).

Naively, one would think that black hole quasinormal modes are not a strong source of gravitational waves for LISA. Indeed, the typical (real part of the) frequency for a Schwarzschild black hole (for quadrupole radiation ℓ = 2 and azimuthal mode m = 2 at the longest damping time n = 0) is f_{ℓ=2,m=2,n=0} ∼ 1.2 × 10⁴(M_☉/M) Hz whereas the e-folding time for the exponential decay of the amplitude (associated with the imaginary part of the complex frequency) is τ_{ℓ=2,m=2,n=0} ∼ 55(M/M_☉) μs. While for a stellar black hole this frequency is several orders of magnitude too high to be detected by LISA, the frequency for supermassive black holes falls inside LISA's good sensitivity frequency band. The quality factor Q associated with the oscillator is Q := πfτ, or for a Schwarzschild black hole Q ≈ 2, i.e., within two oscillations the amplitude of the waves is suppressed by a factor of e^2π ∼ 535. (While for rotating black holes the quality factor is higher, it is still low unless the black hole's spin approaches extremality. Indeed, the quality factor can be approximated by Q₂₂₀ = 0.7000 + 1.4187(1 − j)^−0.4990 ± 0.88% (Berti et al. 2006). A different parameterization of Q₂₂₀, which however agrees with this one to within 1%, was first given by Echeverria 1998.) Therefore, a black hole is expected to be a very poor resonator, and the waves difficult to detect. Regardless of these naive expectations, Flanagan & Hughes (1998) showed that the signal-to-noise ratio (henceforth S/N) from the ringdown stage could be comparable to that of the inspiral stage, so that waves emitted during the ringdown stage may be an important source for LISA.

The signal to be detected by LISA comes from sources that span a large parameter space, including both source parameters and detector parameters. To find the S/N from a large sample of ringdown sources as a function of the black hole's mass one needs to include our ignorance of the particular parameters. Specifically, we do not know the parameters θ, ϕ that are the spherical polar angles that refer to the direction to the source, and ψ which is the angle between the waves polarization plane and the plane ϕ = 0 that determine the detector pattern functions F^+,×(θ, ϕ, ψ) for the two gravitational-wave polarization states. Specifically,

$$\begin{eqnarray} F^{+}(\theta, \phi, \psi)&&=\case{1}{2}(1\,+\,\cos ^2\theta)\,\cos 2\phi \,\cos 2\psi \nonumber \\ &&-\,\cos \theta \,\sin 2\phi \,\sin 2\psi \nonumber \end{eqnarray}$$

$$\begin{eqnarray} F^{\times }(\theta, \phi, \psi)&=&\case{1}{2}(1\,+\,\cos ^2\theta)\,\cos 2\phi \,\cos 2\psi \nonumber \\ &&+\,\cos \theta \,\sin 2\phi \,\sin 2\psi \,.\nonumber \end{eqnarray}$$

We also are ignorant of the source's parameters ι, β that are the spherical polar angles in the frame of the source, defining the orientation of the source in its local frame (i.e., determining the orientation of the black hole's spin axis). Other unknown parameters are the black hole's dimensionless specific spin angular momentum j := cJ/(GM²), the phases of the two polarization states ϕ^+,×_ℓ,m,n, the total radiation efficiency _rd := M⁻¹c⁻²∫^∞₀( dE/ df) df (i.e., the fraction of Mc² being emitted in gravitational waves), and ⁺, the fraction thereof in the + polarization state, in addition to the distance to the source, which we take to be the luminosity distance D_L. Altogether, we have an 11-dimensional parameter space (three dimensions for the detector parameters, and eight for the source's; Thorne 1987).

The currently available method is to make the ignorance of the parameters manifest by averaging over some of the parameters (the detector parameters, i.e., θ, ϕ, and ψ, in addition to the two source parameters ι and β) and fixing other parameters (the source parameters j, ϕ^+,×_ℓ,m,n, _rd, ⁺, D_L). Consider the question of what is the S/N for ringdown sources of a particular mass. This question is important for detection purposes or for estimation of the population statistics. (This question is different from the question of what is the S/N for a ringdown source with particular parameters, which we address in what follows.) Specifically, as one black hole source may be in one direction in the sky with a certain value for the specific spin (and similarly for the other parameters), and another black hole source may be in a different direction with a different value of the specific spin (and the other parameters), one may assume that all black hole sources have average (or fixed) values for all parameters (note that Berti et al. 2006 did not average over the spin angular momentum j, but rather fixed it). Hereafter, we refer to this approach as the "all-sky average" (because it involves averaging over the detector parameters, although it also involves parameter fixing and not just averaging). In this approach, one would argue that a good estimate for the total S/N for ringdown sources in a given mass interval would be found by taking the spin of all such sources to be the average value, say j = 0.5 if the distribution function is taken to be uniform. The all-sky average approach was used recently by Berti et al. (2006) to find the S/N. The all-sky average approach clearly is useful as an approximation method to find the S/N, and as we show below does it to very high accuracy. Note that our ignorance of the parameters to be used is of two kinds: first, there are observational unknowns, which will become known after precise observations are made. Second, there are theoretical unknowns, which would become better known (including their distribution functions) after better theoretical modeling and understanding of the sources' physics is obtained. We do not distinguish here between the two types of unknowns.

An alternative approach is to integrate over a finite number of sources, and use Monte Carlo integration to reflect the ignorance of the parameters. That is, one uses random values of the parameters for each source, spanning the parameter space, and then averages over a large number of such sources. The range of the parameters and the distribution functions for the random sample then reflect the physics of the source and the gravitational-wave generation. In practice, we take here for simplicity all distribution functions to be uniform, but this assumption is easy to relax if evidence is presented for nonuniform distributions. Specifically, we take the sources to be homogeneous and isotropic (so that the uniformity of the distribution function is in D²_L and cos θ), and the sources' spin axes to be isotropic (i.e., the distribution function is homogeneous in cos ι). Recently, our understanding of gravitational-wave sources has improved dramatically. We now know that the radiation efficiency for binary inspiral, at least for the variables simulated using numerical relativity, is _rd ∼ 3% with circular polarization. However, it was argued that _rd ≳ 7% for some cases (Berti et al. 2007), and for head-on collision it is _rd ∼ 0.1% with + polarization (Berti et al. 2006). More recent numerical relativity simulations for the merger of comparable mass black holes have found a relationship of the gravitational-wave radiation efficiency to the spin of the initial black holes, and a relationship of the spin of the final black hole to the spin of the initial black holes (Campanelli et al. 2007). For the specific case of initially aligned or antialigned black hole spins, the relationship of the radiation efficiency in gravitational waves to the final black hole's spin j is given by

$$\begin{equation} \epsilon _{\rm rd}(j) = 0.994-0.438\,\sqrt{4.751-3.74\,j}-0.4545\,j\,, \end{equation} \tag{ 1 }$$

as obtained from a quadratic fit of the numerical relativity results of Campanelli et al. (2007). The final black hole's spin was found to be in the range 0.35 ≲ j ≲ 0.95, which implies radiation efficiency in the range 2.2% ≲ _rd ≲ 8.3%. The dependence _rd = _rd(j) reduces our parameter space from 11 to 10 independent dimensions. In what follows, we specialize our discussion to the case of quasinormal modes from black hole sources resulting from the merger of black hole binaries with initially aligned or antialigned spins. This specialization allows us to use actual numerical relativity values for the radiation efficiency in gravitational waves. Our method can be applied also for black holes resulting from less restrictive initial spin configurations when corresponding numerical relativity results become available.

One may naively expect that in the limit of very many sources, the Monte Carlo results reduce to the all-sky average results. However, they do not. The large number limit of the Monte Carlo integration approaches the all-sky average only if the parameters of the S/N are independent. When considering only Schwarzschild black holes, or even all black holes with a fixed value for the specific spin angular momentum, indeed all the parameters are independent, i.e., they live in different spaces. However, when one considers that for any given value for the black hole mass M the specific spin spans the range 0 ⩽ j ⩽ 1 (or, more specifically for the case considered here in detail, 0.35 ≲ j ≲ 0.95), these parameters are no longer independent.

Specifically, the S/N squared is given by (Berti et al. 2006)

$$\begin{eqnarray} \rho ^2 &=& 2\left(\frac{GM}{rc^2}\right)^2\int _0^\infty \frac{df}{S_h(f)}\Bigg \lbrace (b_+^2 + b_{-}^2)\Big[{A^+_{\ell mn}}^2{F^+}^2 \nonumber \\ &+& {A^\times _{\ell mn}}^2{F^\times }^2 - 2{A^+_{\ell mn}}{A^\times _{\ell mn}}{F^+}{F^\times }\sin \big(\phi ^+_{\ell mn}\nonumber \\ &-&\phi ^\times _{\ell mn}\big)\Big]\times \left| S_{\ell mn} \right|^2 \nonumber \\ &+& 2b_+b_{-}\left\lbrace \Re \left[\left({A^+_{\ell mn}}^2{F^+}^2e^{2i\phi ^+_{\ell mn}} - {A^\times _{\ell mn}}^2{F^\times }^2e^{2i\phi ^\times _{lmn}}\!\right){S_{\ell mn}}^2\!\right]\right.\!\!\!\! \nonumber \\ &+& \left. 2{A^+_{\ell mn}}{A^\times _{\ell mn}}{F^+}{F^\times }\Im \big[e^{i(\phi ^+_{\ell mn} + \phi ^\times _{\ell mn})}{S_{\ell mn}}^2\big] \right\rbrace \Bigg \rbrace. \end{eqnarray} \tag{ 2 }$$

Here, b_± = (1/τ_ℓmn)/[(1/τ_ℓmn)² + 4π²(f ± f_ℓmn)²] is the Fourier transform of the time domain single-frequency quasinormal signal that following Berti et al. (2006), we approximate using the delta-function (in frequency) approximation, and S_h(f) is the LISA noise spectral density. In practice, we calculate S_h(f) using the approximation given by Finn & Thorne (2000) and by Larson et al. (2000), which agrees well, except for oscillations at the high-frequency end, with the LISA Sensitivity Curve Generator results.³ The amplitudes in the two polarization states are denoted by A^±_ℓmn, and S_ℓmn(ι, β) are the spin-weighted spheroidal harmonics of spin weight 2, which depend on j through the overtone index n. Also the amplitudes A^±_ℓmn are functions of j, as

$$\begin{eqnarray} &&A^{+}_{\ell mn} = \sqrt{\frac{8c^3\,\epsilon _{\rm rd}^{+}(1 + 4Q_{\ell mn}^2)}{GMQ_{\ell mn}f_{\ell mn}(1 + 2Q_{\ell mn}^2)}},\nonumber \\ &&A^{\times }_{\ell mn} = \sqrt{\frac{4c^3\,\epsilon _{\rm rd}^{\times }(1 + 4Q_{\ell mn}^2)}{GMQ_{\ell mn}^3f_{\ell mn}}}, \end{eqnarray} \tag{ 3 }$$

and all three variables f_ℓmn, Q_ℓmn, and _rd are functions of j. Clearly, 〈A_ℓmn(j) S_ℓmn(j)〉 ≠ 〈A_ℓmn(j)〉 × 〈S_ℓmn(j)〉 when the averaging is done over j. We may not therefore expect the qualitative dependence of ρ² as a function of the mass M to remain unchanged when the all-sky averaging assumption is relaxed. In this paper, we study this question and find the S/N when the sources are obtained from a random sample. That is, we do not assume following Berti et al. (2006) that all black hole quasinormal mode sources have the average or fixed values for all parameters, and we let these parameters take random values, as one would expect from an actual sample of sources for LISA. The all-sky average approximation of Berti et al. (2006) assumed that ${F^{+,\times }}^2=\langle {F^{+,\times }}^2\rangle \equiv\frac{1}{5}$, F⁺F^× = 〈F⁺F^×〉 ≡ 0, and that $S_{\ell mn}^2=S_{\ell mn}^*\,S_{\ell mn}=\langle S_{\ell mn}^*\,S_{\ell mn} \rangle \equiv \frac{1}{4\pi }$, where an asterisk denotes complex conjugation. It was further assumed following Flanagan & Hughes (1998) that ϕ⁺_ℓmn = ϕ^×_ℓmn = 0 and that the amplitudes A⁺_ℓmn = A^×_ℓmn = A_ℓmn. The specific spin angular momentum j and the radiation efficiency _rd were fixed. To compare our Monte Carlo results with the all-sky average approach, we set the spin angular momentum in the all-sky average case to the value that gives the average radiation efficiency. Specifically, we find j that would solve $\epsilon _{\rm rd}(j)=\langle \epsilon _{\rm rd}(j)\rangle \equiv\,\int _{j_{\rm min}}^{j_{\rm max}}\epsilon _{\rm rd}(j^{\prime })\,dj^{\prime }/(j_{\rm max}-j_{\rm min})$. Over the range of j for which _rd(j) is valid, we find in practice that 〈_rd(j)〉 = 0.03824 and the corresponding spin angular momentum is j = 0.71882.

Our Monte Carlo approach is a simple application of Monte Carlo integration, as we are using uniform distribution functions for our random variables (see below for justification). The (pseudo-)random number generator is based on the Mersenne Twister algorithm (Matsumoto & Nishimura 1998). While not considered secure for cryptography applications, this generator is certainly reliable enough for our purposes. Specifically, we sample points uniformly from the integration region to estimate the integral and its error. Namely, we calculate the S/N by estimating the integral in Equation (2) by summing the integrand over the randomly chosen sample, dividing by the volume of the parameter space and by the number of iterations. The estimated error σ in the Monte Carlo integration is done by calculating the variance, according to σ² = (V/N)² ∑^N_i=1[ρ(x_i) − 〈ρ〉]², where V is the volume of the parameter space and N is the number of elements in the sample. Here, x_i denotes collectively the 11 parameters for the S/N for the ith member of the sample. The reason why this very simple approach to Monte Carlo integration is sufficient for our needs is that we choose the sources to be uniformly distributed across the parameter space. Specifically, we take the sources to be distributed homogeneously and isotropically, which seems to be a very reasonable assumption due to the cosmological nature of the sources. We also take the spin axis of the black hole members of the sample to be distributed isotropically. These assumptions mean that we take D²_L,  cos θ,  cos ι, ϕ, and β to be uniformly distributed. We have little motivation to prefer a nonuniform distribution function for the other parameters. Specifically, we currently do not have enough numerical relativity results to determine a realistic distribution function for the parameters ψ, ϕ^+,×, or ⁺ in addition to the population statistics of j. We therefore take these parameters to be distributed uniformly in their respective parameter spaces. The radiation efficiency _rd(j) is determined according to Equation (1).

Our approach also allows us to study sources with particular parameters instead of the all-sky averaged counterparts. Specifically, we consider the directional variations of fixed sources, and show in which directions in the sky LISA would be more sensitive, including finding LISA blind spots.

In all our simulations, we take a standard zero-curvature cosmological model with the 5-year WMAP values Ω_M = 0.279, Ω_Λ = 0.721, and H₀ = 70.1 km s⁻¹ Mpc⁻¹ (Hinshaw et al. 2009), and the luminosity distance–redshift relation to be given by (Carroll et al. 1992)

$$$ D_L(z)=(1+z)\frac{c}{H_0}\int _0^z\frac{\,dz^{\prime }}{\sqrt{(1+z^{\prime })^2(1+\Omega _M z^{\prime })-z^{\prime }(2+z^{\prime })\Omega _{\Lambda }}}\,. $$$

We first test our Monte Carlo integration for the simple case in which all black holes are assumed to be Schwarzschild black holes (j = 0). In particular, we confront our results with the all-sky average results, and study how well they agree. Figure 1(A) shows the S/N for the all-sky average and for the Monte Carlo simulation. We find that the two signals overlap, as indeed is expected. In particular, in the Schwarzschild case our Monte Carlo results agree to the expected accuracy level with the results of Berti et al. (2006).

Zoom In Zoom Out Reset image size

The Schwarzschild test does not test all sectors of our code, as it does not calculate the spin-weighted spheroidal harmonics. Fixing the value of j, we can test the remaining sectors of the code. Figure 1(B) displays the S/N for the all-sky average (for a fixed value of j) and the Monte Carlo results, which are in agreement to within the error bars.

We also test the convergence of our Monte Carlo code, letting all 10 independent variables be random. Figure 2 shows the S/N for different Monte Carlo runs, with different sample sizes. The convergence of the result is indeed as expected, namely, scaling with $\sqrt{N}$, with N being the number of sources in the sample (see a more detailed analysis of Figure 2 in Section 3).

Zoom In Zoom Out Reset image size

As we already pointed out, we do not expect the behavior of the S/N when j is let to vary to necessarily be identical to the case of j being averaged. In Figure 3, we fix all the variables, and let only j be a random variable. The distance of the sources is fixed at redshift of z = 1.509. Indeed, we find that the behavior is not identical. While for very low and very high masses the behavior is nearly indistinguishable, there is a range of masses, specifically M ∼ 4–7 × 10⁶ M_☉ for which the S/N in the variable case is significantly higher than in the all-sky average case. The ratio of the two S/Ns spikes at M ∼ 5.6 × 10⁶ M_☉, where ρ/ρ_all-sky ≈ 1.4 which would translate to an event rate greater by a factor of ∼2.75. This increase in the S/N is 14.2σ, so that it is not an artifact of the Monte Carlo error associated with the finite sample size.

Zoom In Zoom Out Reset image size

Next, we let all the variables—except for the distance, which we fix at redshift of z = 1.509—be random. Figure 4 shows the S/N, which preserves the behavior we found in Figure 3. Specifically, the sharp increase in S/N for masses in the range M ∼ 4–7 × 10⁶ M_☉ is retained. The S/N is increased by a factor of 1.6 that translates to an increased event rate by a factor of 4.1. This increase is 7.4σ. Indeed, this is the expected behavior as all other variables live in independent spaces. Note that as the number of dimensions here is higher than in Figure 3, the errors are commensurately larger.

Zoom In Zoom Out Reset image size

We present in Figure 2 the results for all variables being random, including the distance, from the distance to the Virgo cluster (18.0 ± 1.2 Mpc) up to redshift of z = 2. We find that the S/N is greatly increased compared with all sources being at the average distance (squared; redshift of z = 1.509) over a large range of masses, specifically for M ≳ 3.5 × 10⁶ M_☉. The increased S/N peaks at M ∼ 5 × 10⁶ M_☉, with an increase in the S/N by a factor of 2.4, which translates to an increase in the event rate by a factor of 13.8. This increase is 9.1σ. As LISA will detect sources at all distances, we believe that the increased S/N shown in Figure 2 is closer to the actual S/N for LISA.

Most importantly, our approach allows us to find the S/N in different directions in the sky. In Figure 5, we show the S/N as a function of the colatitude θ for various values of ϕ, along with the all-sky counterpart. We show that the same source (namely, same source's parameters) located in different directions have very different S/Ns projected on the LISA detector. In particular, there is a LISA blind spot at θ = π/2 and ϕ = π/4 (and other values of ϕ on the LISA equatorial plane given by the rotational symmetry, namely, at ϕ = 3π/4, 5π/4, 7π/4; our numerical results are within 0.1% of these values).

Zoom In Zoom Out Reset image size

This work can be extended in the following ways. First, one may relax the assumption of uniform distribution functions for the random variables, and base the variables on more realistic distribution functions. Second, we use here the delta-function approximation to the quasinormal mode frequency. This assumption may be relaxed (Berti et al. 2006), but the change in the S/N is not expected to be significant, and third, one may extend the single-mode analysis done here to multiple modes (Berti et al. 2006).

The authors are indebted to Vitor Cardoso, Leor Barack, and Emanuele Berti for discussions. This work was supported in part by grants NASA/GSFC NCC5-580 and NASA/SSC NNX07AL52A, by NSF grant PHY-0757344, and by a minigrant from the UAHuntsville Office of the Vice President for Research.