Numerical study of the quasinormal mode excitation of Kerr black holes

We present numerical results from three-dimensional evolutions of scalar perturbations of Kerr black holes. Our simulations make use of a high-order accurate multiblock code which naturally allows for adapted grids and smooth inner (excision) and outer boundaries. We focus on the quasinormal ringing phase, presenting a systematic method for extraction of the quasinormal mode frequencies and amplitudes and comparing our results against perturbation theory. The detection of a single mode in a ringdown waveform allows for a measurement of the mass and spin of a black hole; a multimode detection would allow a test of the Kerr nature of the source. Since the possibility of a multimode detection depends on the relative mode amplitude, we study this topic in some detail. The amplitude of each mode depends exponentially on the starting time of the quasinormal regime, which is not defined unambiguously. We show that this time-shift problem can be circumvented by looking at appropriately chosen relative mode amplitudes. From our simulations we extract the quasinormal frequencies and the relative and absolute amplitudes of corotating and counterrotating modes (including overtones in the corotating case). We study the dependence of these amplitudes on the shape of the initial perturbation, the angular dependence of the mode, and the black hole spin, comparing against results from perturbation theory in the so-called asymptotic approximation. We also compare the quasinormal frequencies from our numerical simulations with predictions from perturbation theory, finding excellent agreement. For rapidly rotating black holes (of spin $j=0.98$) we can extract the quasinormal frequencies of not only the fundamental mode, but also of the first two overtones. Finally we study under what conditions the relative amplitude between given pairs of modes gets maximally excited and present a quantitative analysis of rotational mode-mode coupling. The main conclusions and techniques of our analysis are quite general and, as such, should be of interest in the study of ringdown gravitational waves produced by astrophysical gravitational wave sources.

Figure 1

Illustration of the six-block grid structure and the cubed sphere coordinates that are used for the simulations in this paper. The left panel shows the distribution of grid points on a sphere of constant radius. The central panel shows a snapshot from a scalar wave evolution on an equatorial cut. The plot refers to an $\ell =m=2$ mode on the background of a Kerr black hole with spin $j=0.9$ at $t=92.2M$. Also shown are the locations of the interblock boundaries. The right panel magnifies the central region of the domain in the equatorial plane, showing the grid structure around the spherical excision boundary. The four dark lines mark the interfaces between blocks.

Figure 2

Event and Cauchy horizons for a Kerr black hole with spin $j\gtrsim 0.96$ in standard Kerr-Schild coordinates (as defined in the text), here shown in the $x$-$z$ plane. The horizons have an ellipsoidal shape; it is therefore not possible to fit a spherical excision region (dotted line) between the two horizons.

Figure 3

The left panel shows the $\ell =m=2$ component of the waveform extracted at radius $r=5M$ on a Kerr black hole background with a spin of $j=0.9$. The waveform is a superposition of the corotating and the counterrotating mode, and the beating of two different frequencies is clearly visible. The right panel shows the waveform for $t\ge 40M$ as well as a QNM fit with the fundamental $\ell =|m|=2$ modes. The interval used for the fit is $[74.5M,150M]$. The inlay shows the absolute value of the difference between the fit and the data. At times between the excitation of the QNM ($t\sim 25M$) and about $70M$ the differences are mainly due to the presence of the ($l=2$, $m=2$, $n=1$) mode, the exponentially decaying mode that can be seen in the inlay (a fit of this mode yields quasinormal frequencies in agreement with perturbation theory). At times $t\lesssim 25M$ the difference is due to the initial burst.

Figure 4

Critical uncertainty in the starting time, as defined by Eq. (21), assuming $ϵ={10}^{-2}$. In the left panel we give the critical ${\delta }_0$ for fundamental modes ($n=n^{\prime }=0$) with different angular dependence. For the first mode we assume $\ell =m=2$; the second mode has ${\ell }^{\prime }=2$ and different values of $m^{\prime }=1$, 0, $-1$, $-2$ (lines from top to bottom). In the right panel we show the critical uncertainty in the relative amplitude of the fundamental mode and first overtone, i.e., $n=0$ and $n^{\prime }=1$. Here we set $\ell ={\ell }^{\prime }=2$, consider all values of $m=m^{\prime }$ and once again we assume $ϵ={10}^{-2}$.

Figure 5

Residual in the fit, as defined in Eq. (22), as a function of the initial time for the fitting $T_i$. Looking at the minimum of the residual we can determine $T_i$ with high precision. This plot corresponds to a simulation with spin $j=0.5$, $\ell =m=2$ initial data with $A=0$, $B=1$ and a radial resolution $\Delta r=M/10$.

Figure 6

The left and right panels show the real and imaginary parts of the quasinormal frequencies extracted from the simulations of Fig. 5. From the optimal starting time determined by minimizing the residual, $T_i=(59.65\pm 0.025)M$ (see previous figure), we find ${\omega }_R=0.585887\pm 1\times {10}^{-6}$ and ${\omega }_I=0.0933851\pm 5\times {10}^{-7}$.

Figure 7

The left panel shows the extracted ($\ell =4$, $m=2$, $n=0$) waveform for three simulations with different spin parameters as seen by an observer at $r=5M$. The initial data are a pure ($\ell =m=2$) mode and are set up according to Eq. (10) with $A=0$, $B=1$, and $r_0=20M$. For zero spin the different multipole components of the solution should evolve independently and no modes besides the one in the initial data should be excited, while for nonzero spin modes with different $\ell$ but same $m$ do couple [26]. In the Schwarzschild case the ($\ell =4$, $m=2$, $n=0$) waveform differs from zero due to our grid structure and discretization errors, but it converges to zero with increasing resolution. This is illustrated by the right panel, which shows the extracted ($\ell =4$, $m=2$, $n=0$) amplitude for $j=0$ and $j=0.9$ from runs with two resolutions ($20\times 20\times 1000$ and $30\times 30\times 1500$ points per patch and outer boundaries at $100M$). Only for $j=0.0$ the mode converges to zero.

Figure 8

Results from a run with initial data parameters $\ell =m=2$ and spin $j=0.9$. The left panel shows the square of the amplitude of all modes up to $\ell =4$ which are not within the roundoff error. The right panel shows the sum over the square of those modes compared to the integral over a sphere of the full field squared: As expressed in Eq. (12) the two curves lie on top of each other, and there is no relevant contribution from higher modes.

Figure 9

Numerically obtained excitation amplitudes of the $\ell =|m|=2$ fundamental modes assuming an observer location $r_{\mathrm{obs}}=40M$ and a ringdown starting time $t_0=60M$. The left panel refers to a black hole with spin $j=0.5$. According to predictions from perturbation theory in the asymptotic approximation [cf. Eq. (24) and the following discussion] the maximum for $m=2$ should be located at ${\sigma }_{220}=2.445$, while the value that we obtain from our simulations is ${\sigma }_{220}=2.55\pm 0.05$ (the uncertainty describing the difference between consecutive values of $\sigma$ used in our simulations: $\Delta \sigma =0.05$). Similarly, for $m=-2$ the width at the maximum should be ${\sigma }_{2-20}=3.434$, while we obtain ${\sigma }_{2-20}=3.875\pm 0.075$. The right panel, in turn, refers to a black hole with spin $j=0.9$. In this case the theoretical (numerical) maxima are located at ${\sigma }_{220}=1.816$ (${\sigma }_{220}=1.85\pm 0.05$) and ${\sigma }_{2-20}=3.758$ (${\sigma }_{2-20}=3.85\pm 0.05$), respectively. The inset in the left panel is a zoom around the maximum for $j=0.5$ and $m=2$. As discussed in the text, an uncertainty in the excitation time of $0.09M$ would already explain the difference between the predicted location of the maxima and our numerical results.