Mode coupling of Schwarzschild perturbations: Ringdown frequencies

Mode coupling of Schwarzschild perturbations: Ringdown frequencies Enrique Pazos,1, 2 David Brizuela,3 José M. Martı́n-Garcı́a,4, 5 and Manuel Tiglio6 arXiv:1009.4665v1 [gr-qc] 23 Sep 2010 1 Center for Relativistic Astrophysics, School of Physics, Georgia Institute of Technology, 837 State Street, Atlanta, GA 30332-0430 2 Departamento de Matemática, Universidad de San Carlos de Guatemala, Edificio T4, Facultad de Ingenierı́a, Ciudad Universitaria z. 12, Guatemala 3 Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität, Max-Wien-Platz 1, 07743 Jena, Germany 4 Institut d’Astrophysique de Paris, CNRS, Univ. Pierre et Marie Curie, 98bis boulevard Arago, 75014 Paris, France 5 Laboratoire Univers et Théories, CNRS, Univ. Paris Diderot, 5 place Jules Janssen, 92190 Meudon, France 6 Department of Physics, Center for Fundamental Physics, Center for Scientific Computation and Mathematical Modeling, Joint Space Institute. University of Maryland, College Park, MD 20742, USA. Within linearized perturbation theory, black holes decay to their final stationary state through the well-known spectrum of quasinormal modes. Here we numerically study whether nonlinearities change this picture. For that purpose we study the ringdown frequencies of gauge-invariant second-order gravitational perturbations induced by self-coupling of linearized perturbations of Schwarzschild black holes. We do so through high-accuracy simulations in the time domain of first and second-order Regge-Wheeler-Zerilli type equations, for a variety of initial data sets. We consider first-order even-parity (ℓ = 2, m = ±2) perturbations and odd-parity (ℓ = 2, m = 0) ones, and all the multipoles that they generate through self-coupling. For all of them and all the initial data sets considered we find that —in contrast to previous predictions in the literature— the numerical decay frequencies of second-order perturbations are the same ones of linearized theory, and we explain the observed behavior. This would indicate, in particular, that when modeling or searching for ringdown gravitational waves, appropriately including the standard quasinormal modes already takes into account nonlinear effects. PACS numbers: 04.25.dk, 04.40.Dg, 04.30.Db, 95.30.Sf I. OUTLOOK AND MOTIVATION Black hole no-hair theorems [1, 2] state that within Einstein’s theory the end point of any system with enough gravitational energy to form a black hole is remarkably simple: it is uniquely characterized by one member of the Kerr family1 [3], which is described by only two parameters: the spin and mass of the final black hole. As a consequence, the details by which different systems decay to such endpoints have been of interest for many decades. Pioneering studies were done by a number of authors in the early seventies, starting with studies of linearized perturbations of non-rotating (Schwarzschild) black holes (e.g., [4–6]). Press realized that there is always an intermediate stage where the ringdown is dominated by a set of oscillating and exponentially decaying solutions, quasinormal modes (QNMs), whose spectrum depends only on the mass of the black hole and the multipole index ℓ of the initial perturbation [7]. This regime is followed by a power-law ‘tail’ decay due to backscattering [8]. In the case of gravitational perturbations of nonrotating black holes the relevant equations from which QNMs can be inferred are the Regge-Wheeler [9] and 1 Charge is expected not to play a significant role in most astrophysical scenarios. Zerilli [10, 11] ones. For rotating black holes the corresponding one (though based on a curvature formalism, as opposed to a metric one) is the Teukolsky equation [12]. Their QNMs were first studied by Teukolsky and Press [13]. See [14, 15] for comprehensive reviews on the rich area of QNMs. The QNM with the lowest frequency is called the fundamental one. Since the subsequent ones (overtones) decay much faster, the ringdown of Kerr black holes in linearized theory is in practice described by a few oscillating modes which decay exponentially in time, till they reach the tail regime. It is interesting to note that the tail decay problem for rotating black holes is still not completely understood [16–20]. From an observational point of view this universal ringdown spectrum is of great power: one can use a single QNM detection to infer the mass and spin of the black hole source, assuming General Relativity to be correct. Alternatively, through a two-mode detection one can test General Relativity and/or the assumption that a black hole is the source of the measured signal [21]. The main idea is that the QNM frequencies of both detections have to be consistent with respect to their inferred masses/spins. The LISA mission is expected to measure gravitational waves in the low-frequency spectrum: (10−5 − 10−1 ) Hz, such as those emitted in the collision of supermassive binary black holes (SMBBHs) [22]. Flanagan and Hughes [23] showed that, quite generically, the signal to noise ratio for these sources in the inspiral regime should 2 be comparable to that one in the ringdown. Therefore, detection of SMBBHs by LISA through the measurement of QNMs seems to be feasible. Assuming a lower cutoff of (10−4 − 10−5 ) Hz and requiring that the QNM signal lives long enough to travel once through LISA’s propagation arms places a constraint on the mass range of the SMBBH candidates: a few 105 M⊙ to (108 − 109 )M⊙ . A step beyond detection analysis is that one of parameter estimation. In Ref. [22, 24] it was found that through a single QNM detection LISA would be able to accurately infer the mass and spin of supermassive black holes: for black holes with mass M & 105 M⊙ the errors in mass and spin would be smaller than one part in 102 , and smaller than one part in 105 for the more optimistic case M & 5 × 105 M⊙ 2 . These predicted accuracies depend on the ringdown efficiency ǫrd , defined as the fraction of mass radiated in ringdown waves. In these references very conservative values were used: ǫrd ∼ 0.1%−3%. For example, it has been found in numerical simulations of two equal-mass, non-spinning black holes starting from quasi-circular motion that around ǫrd ∼ 2% − 3% of the total mass is radiated in the ringdown regime [25]. The inclusion of different masses and/or spin increases this value (see, for example, [26]). In [22, 24] it was also found that at least a second detection of either mass or spin should be possible for LISA. Resolving both spin and angular momentum (or, equivalently, both frequency and damping times associated with the QNM oscillation of this second mode) might require a very large critical signal to noise ratio, which might in turn need the second mode to radiate a significant portion of the emitted gravitational wave when compared to the first one. Whether this is feasible or not can only be established by giving precise predictions of the amplitudes for secondary candidates. Underlying in all these analyses is the implicit assumption that quasinormal modes and their spectrum of associated frequencies accurately describe the (intermediate) stage of the ringdown to a final Kerr stationary state. Which is certainly the case in linearized theory, but at the same time there is evidence that effects of self-interaction in gravitational waves might be observable with the expected sensitivity of LISA [27]. Similarly, quasinormal modes also play an important role in the semi-analytical modeling of intermediate mass black holes (IMBH), such as in the Effective One Body approach [28, 29], where the gravitational wave as modeled within this formalism is stitched to a ringdown one consisting of QNMs by enforcing continuity of the wave and its derivatives and using the values of quasinormal frequencies as expected from linearized theory [30–32]. In this context, it is worth recalling that in close-limit stud- 2 Only cases with M & 105 M⊙ are considered in [22, 22, 24] because otherwise the QNM signal would be short lived enough that special detection techniques might be needed. ies of binary black holes it was found that corrections from second-order perturbations were in some cases significant [33]. For IMBHs, which could have total masses in the range of ∼ 100M⊙ − 104 M⊙ , it is especially important to accurately model the merger and ringdown since they should fall in the frequency band of earth based gravitational wave detectors. Although the existence of IMBHs is still debatable, they could provide an interesting source for Advanced LIGO and VIRGO if they are present in dense globular clusters [34] (see also [35]). Recent observational evidence of an IMBH can be found in [36–39]. The previous discussions motivate us to study how nonlinear perturbations of black holes decay in time: do they do so just as in linearized theory or with a different spectrum of frequencies? We carry out our study through numerical simulations of first and second-order gaugeinvariant gravitational perturbations of Schwarzschild black holes. We find that for all practical purposes, and to high numerical accuracy, the complex decay frequencies of second-order perturbations are the standard quasinormal ones from linearized theory –in contrast to previous predictions in the literature [40–42]– and we explain why this appears to be the case. Essentially, in all our simulations we find that mode-mode couplings excite nonlinearities in the early stages of the perturbations before the quasinormal regime for the linearized perturbations kicks in. By the time the latter happens, those couplings have decreased to negligible values and the excited nonlinearities essentially propagate as in linearized theory; and, in particular, they decay with the standard QNM frequencies. The structure of the paper is as follows. Section II reviews the basics of Regge-Wheeler-Zerilli equations and quasinormal modes, and Sec. III the main features of the gauge-invariant approach here used for second-order perturbations. Section IV describes our numerical approach and setup for solving the first and second-order ReggeWheeler and Zerilli equations, and Sec. V presents and discusses our main results. II. FIRST-ORDER PERTURBATIONS OF SCHWARZSCHILD AND QUASINORMAL MODES Metric gravitational perturbations can be expanded in tensor spherical harmonics. At the linear level modes with different angular structure decouple from each other if the background spacetime is spherically symmetric, as is the case for the Schwarzschild metric. In addition, for each multipole, perturbations of Schwarzschild can be further split into two sectors with different parities, which also decouple from each other in linearized theory. For each multipole, each of these sectors is completely described by a master function which depends on time and radius. The Regge-Wheeler function contains all the relevant information of the axial or odd-parity sec- 3 tor, whereas the Zerilli one encodes the polar or even degrees of freedom. These functions satisfy master evolution equations, 102 (1) (2) {1} m Φℓ denote, respectively, the firstwhere {1}Ψm ℓ and order Zerilli and Regge-Wheeler functions for a given (ℓ, m) mode. The potentials that appear in these equations are given by zerilli function 10 {1} m  {1}Φm Φℓ = 0 , ℓ − VRW {1} m {1} m  Ψℓ − VZ Ψℓ = 0 , {1} Ψ l=2, m=2 0 10-2 10-4 10-6 10-8 2 VZ ≡ ℓ(ℓ + 1) 6M r λ(λ + 2) + 3M (r − M ) − 3 , r2 r (rλ + 3M )2 VRW ≡ ℓ(ℓ + 1) 6M − 3 r2 r (3) 50 100 150 t/M 200 250 300 Figure 1: Solution to the first-order Zerilli equation at observer location r = 51.8M . (4) and are, in particular, independent of the azimuthal index m. In these expressions λ ≡ 21 (ℓ − 1)(ℓ + 2), M is the mass of the Schwarzschild black hole background, r its areal radius, and  is the two-dimensional D’Alambertian operator corresponding to the time and radial sector of the background. The complete spectrum of QNMs can be numerically obtained by analyzing Eqs. (1,2) in the frequency domain. Computing in this way the amplitudes of QNMs for any given initial data is not so straightforward, though. Bearing in mind our motivation of studying the behavior of second-order perturbations, we instead solve Eqs. (1,2) in the time-domain. That is, we prescribe (a variety of) initial data, evolve them and analyze the solutions at different observer locations as functions of time. The early behavior of the solution depends on the type of initial data, followed by the QNM ringdown of the black hole. In Fig. 1 we show a typical solution of the Zerilli equation, for a fixed observer at r = 51.8M as a function of time. The initial data for this particular case consist of a Gaussian profile to the initial time derivative of the Zerilli function, centered at r = 20M with a width σ = 4M , and the initial value of the Zerilli function itself is set to zero. [This corresponds to what we call time-derivative initial data in the following sections, see Eqs. (12,13)]. To measure the complex QNM frequency we perform a numerical fit to a function of the form f (t) = A eat sin [b(t − t0 )] , 0 (5) where the parameters to be fitted are A, a, b and t0 and the choice of the starting time for the QNM regime is chosen as that one which optimizes the fit, as introduced and explained in Ref. [43]. The quasinormal frequency is therefore given as ω = a + bi. The expected value for the fundamental QNM for an ℓ = 2 perturbation is 0.37367 − 0.08896i [15], while our fit for this simulation yields 0.37077 − 0.08826i. III. SECOND-ORDER GAUGE-INVARIANT PERTURBATIONS OF SCHWARZSCHILD We study second-order gravitational perturbations of Schwarzschild black holes using a gauge-invariant formalism for arbitrary first and second-order perturbations [44]. The key feature of this formalism is being able to consider perturbations with arbitrary angular multipole structure, and has been possible mostly due to the development of a suitable theoretical framework [45, 46] and to the advance of very efficient symbolic algebra tools for tensor-type calculations [47, 48]. Next we very briefly summarize those results of the formalism presented in [44] which are relevant for the current work; see that reference for more details. Due to the intrinsic nonlinearities of General Relativity, any non-trivial solutions of Eqs. (1,2) generate second-order contributions which are solutions of Zerilli and Regge-Wheeler-type equations with source terms, {2} m  {2}Φm Φℓ = ℓ − VRW {2} m {2} m  Ψℓ − VZ Ψℓ = {2} SΦ , SΨ . {2} (6) (7) The sources {2}SΨ and {2}SΦ depend quadratically on the lower-order perturbations and their time and space derivatives from both first-order sectors. That is, in general the coupling of even (odd) parity modes generates odd (even) parity second-order modes; see Appendix A for a detailed description of the selection rules for this mode coupling. Second-order Regge-Wheeler-Zerilli (RWZ) type functions are not unique: one can add to them any quadratic combination of the first-order ones and they will still be gauge invariant and will still satisfy equations of the form (6,7) with different sources. For this reason, asking what are the ringdown frequencies of second-order RWZ functions is not an unambiguous question; as opposed to, say, asking what are the ringdown frequencies of the emitted gravitational power. It turns out, however, that, as dis- 4 cussed in [44] and made explicit below, since the secondorder quantities, like {2}Ψ, {2}Φ, {2}SΨ and {2}SΦ , that we use correspond to the ones labeled as regularized in Ref. [44], there is a simple relationship between emitted power and the RWZ functions. Under such choice of second-order RWZ functions it is therefore unambiguous and physically motivated to ask what are their ringdown frequencies. Reference [44] deals with the most general case for these sources and computation of the radiated energy. Here we quantitatively explore the ringdown frequencies First order Multipoles Parity (ℓ = 2 = |m|) even (ℓ = 2, m = 0) odd for some particular cases. Namely, we study first-order (ℓ = |m| = 2) even-parity and (ℓ = 2, m = 0) odd-parity perturbations and the modes that they generate through self-coupling. That is, for simplicity we do not consider the coupling between the mentioned even and odd-parity first-order modes. We do not explicitly list the sources {2} SΨ and {2}SΦ for the cases here considered because they are rather lengthy and complicated expressions, but they are available from the authors upon request. The generated modes and the radiated energy carried by them are described next and summarized in Table I. Second order Multipoles Parity (ℓ = 4 = m),(ℓ = 4, m = 0),(ℓ = 2, m = 0) even (ℓ = 4, m = 0),(ℓ = 2, m = 0) even Table I: First-order modes considered and second-order ones generated by self-coupling. A. CASE A: even-parity ℓ = |m| = 2 perturbations and generated modes As discussed in Appendix A, the self-coupling between these modes generates second-order (ℓ = 4, m = ±4) even-parity (polar) ones, whereas the coupling between them (m = 2 with m = −2) gives rise to second-order (ℓ = 4, m = 0), (ℓ = 2, m = 0) and (ℓ = 0, m = 0) even-parity (polar) modes, and (ℓ = 3, m = 0) and (ℓ = 1, m = 0) odd-parity (axial) ones. Since this paper only deals with different radiative aspects of the system, we can ignore modes with ℓ < 2. Furthermore, we assume that the Zerilli functions {1} ±2 Ψ2 that describe the first-order (l = 2, m = ±2) modes are real Ψ±2 2 ∈ R. {1} (8) This means that both modes are described by the Power(robs , t) = dE ǫ2 ∂t {1}Ψ22 = dt 12π 2 + ∗ same function, since generically the relation ( {1}Ψm l ) = m {1} −m (−1) Ψl holds, and, in essence, the system reduces to a problem of self-coupling. The second-order evenparity (polar) modes inherit this property, in such a way that we only need three functions to describe them: { {2}Ψ02 , Ψ44 } . {2} Due to the assumption (8), none of the second-order odd-parity (axial) modes are generated. This happens because the source for a m = 0 mode must be real. Schematically, the generic term of the source for this axial (l = 3, m = 0) mode can be written as i {1}Ψ22 {1}Ψ−2 2 and it is straightforward to see that its real part vanishes under the assumption (8). In this particular CASE A, the radiated power associated with the mentioned modes for a given observer located at robs as a function of time is given by 9ǫ4 n 2 ∂t {2}Ψ44 640π where ǫ is the perturbative parameter and all the expressions on the right-hand side are evaluated, naturally, at (robs , t). In principle this equation is valid only at null infinity but, as it is usually the case in computations, we evaluate it at a finite radius. Ψ04 , {2} B. 2 + ∂t {2}Ψ04 2 o + ǫ4 ∂t {2}Ψ02 96π 2 + O ǫ5
, (9) CASE B: odd-parity (ℓ = 2, m = 0) perturbations and generated modes In principle, following the selection rules, the odd parity (ℓ = 2, m = 0) mode would generate by self-coupling 5 Hence, none of the second-order odd-parity modes are excited and, up to this order, we are left with two radiative second-order modes, the second-order (ℓ = 0, m = 0), (ℓ = 2, m = 0) and (ℓ = 4, m = 0) even-parity modes, as well as the (ℓ = 1, m = 0) and (ℓ = 3, m = 0) odd-parity ones. However, in this axisymmetric case (m = 0 for all modes), the Clebsch-Gordan-like coefficients that appear in the sources {2}SΦ and {2}SΨ vanish for those second-order modes with an odd harmonic label ℓ (see Appendix A). 3ǫ2 ∂t {1}Φ02 Power(robs , t) = 32π IV. 2 3ǫ4 + 32π NUMERICAL APPROACH FOR SOLVING THE MASTER EQUATIONS We now describe in some detail our numerical approach for evolving the first and second-order RWZ equations, since in the past difficulties have been reported with the high-order derivatives in the sources of these equations. We numerically solve the first and second-order equations using a pseudo-spectral collocation method. The spatial derivatives are computed using Chebyshev polynomials and Gauss-Lobatto (GL) collocation points, and the system is evolved in time using a standard fourthorder Runge-Kutta scheme. We use a small enough timestep for the time integration so that the solution converges exponentially with the number of collocation points (see below). The accuracy of all the simulations presented in this paper are at the level of double precision round-off. GL collocation points are not equally spaced; instead, they cluster near the edges of the computational domain (equally spaced points would not give exponential convergence). For that reason it is standard to use a multidomain approach. Here we subdivide our radial domain in (non-overlapping) blocks of length 10M each, communicated through a penalty technique. At the interface each incoming characteristic mode u+ is penalized according to (see [49] and references therein) αN 2 δ + u˙+ = (. . .) − (u − v + ) rblock where v + is the value of the same mode at the interface point using the neighboring block, rblock is the size of the corresponding block (10M in these simulations), α is the associated characteristic speed, N the number of collocation points on that block and δ a penalty parameter chosen here to be δ = 0.6. At the outer boundary each characteristic incoming mode is similarly penalized to zero; though this is done simply to achieve stability, in our simulations the domain is large enough that our results are causally disconnected from the outer boundary. The singularity of the black hole is dealt with through excision. That is, by using Kerr-Schild coordinates for { {2}Ψ04 , Ψ02 } . {2} The radiated power is then given by  3 ∂t {2}Ψ04 20 2 1 + ∂t {2}Ψ02 9 2  + O(ǫ5 ). (10) the background spacetime and placing an inner boundary inside the event horizon. As an example, Fig. 2 shows a self-convergence test for the first-order (ℓ = 2 = m) Zerilli function, extracted at r = 51.8M 3 , both changing the number of collocation points as well as the timestep. The initial data used for such test correspond to the same one used in Fig. 1, Ψ22 = 0 , {1} Ψ̇22 (t = 0, r) = e−(r−r0 ) {1} 2 /σ2 , (11) with σ = 4M , r0 = 20M , and the spatial domain r ∈ [1.8M, 301.8M ]. In Fig. 3 we also show the result of a convergence test for the generated second order (ℓ = 4 = m) Zerilli mode. From these figures we see that using 30 collocation points per domain and a timestep ∆t = 0.001M gives a numerical error at the level of double precision round-off; from hereon we use such resolutions for all our simulations. In order to compare the magnitude of the errors with the solutions themselves, in Fig. 4 we show the absolute values of the first-order {1}Ψ22 and second-order { {2}Ψ02 , {2}Ψ04 , {2}Ψ44 } Zerilli solutions from the previous plots at their highest resolutions, all extracted at the same observer location. We note in passing that for most of the ringdown the order of magnitude of the second-order Zerilli functions appears to be comparable to (and in one case even larger than) the first-order one. There is no contradiction in this, since their contribution to the radiated energy is scaled by ǫ4 , while the contribution of the first order Zerilli function is scaled by ǫ2 , see Eq. (9). A. Setup of numerical simulations We could introduce non-vanishing second-order modes via initial second-order perturbations. However, we are interested in mode-mode coupling. Put differently, we 3 We place observers at the beginning/end of each domain: 1.8M, 11.8M, 21.8M , etc. 6 10-2 10-6 10-8 10-10 10-12 10-14 10-16 10-18 10-6 10-8 10-10 10-12 10-14 10-16 0 50 100 10-9 150 t/M 200 250 10-11 10-18 300 0 50 100 150 t/M 200 250 300 Figure 3: Numerical errors in the second-order Zerilli function {2} 4 Ψ4 for different spatial resolutions and a fixed timestep ∆t = 0.001M . The errors are to be interpreted as in the previous figures. For N = 30 (which are the ones used in the rest of this paper) and higher, they are at the level of double precision round-off. ∆t1=0.010 - ∆t4=0.0005 ∆t2=0.005 - ∆t4=0.0005 ∆t3=0.001 - ∆t4=0.0005 10-10 temporal resolution errors N15 - N60 N20 - N60 N30 - N60 N40 - N60 10-4 numerical error in {2}Ψ 10-4 spatial resolution errors 10-2 N10 - N60 N15 - N60 N20 - N60 N30 - N60 N40 - N60 10-12 10-13 10-14 10-15 10-17 l=2, m=2 l=2, m=0 l=4, m=0 Ψ l=4, m=4 10-2 0 50 100 150 t/M 200 250 300 Figure 2: Numerical errors for different spatial resolutions using a fixed timestep ∆t = 0.01M (top), and for different timesteps using a fixed spatial resolution of N = 60 points per domain (bottom). Both figures show the differences between several resolutions and the most accurate one, which is N = 60 for the top panel and ∆t4 = 0.0005M for the bottom one. In both cases the observer is located at r = 51.8M . We see exponential convergence and errors in the order of double precision round-off. zerilli function 10-18 {1} Ψ {2} Ψ {2} Ψ {2} 100 10-16 10-4 10-6 10-8 10-10 0 50 100 150 t/M 200 2. Time Symmetric (TS) In the following section we solve the first-order RWZ equations with four different types of initial data, varying both the location r0 and width σ of the initial data, 3. Approximately Outgoing (OUT) Ψ22 (t = 0, r) = M e−(r−r0 ) {1} Ψ̇22 (t {1} Ψ̇22 (t Ψ̇22 (t = 0, r) = e /σ2 , 2 /σ2 , = 0, r) = −(1 − 2M/r)∂r {1}Ψ22 (t = 0, r). 4. Approximately Ingoing (IN) Ψ22 (t = 0, r) = 0 , {1} 2 = 0, r) = 0 . Ψ22 (t = 0, r) = M e−(r−r0 ) {1} {1} {1} 2 −(r−r0 ) /σ 2 . 300 Figure 4: Absolute value of first and second-order Zerilli functions extracted at r = 51.8M . are interested in the particular solution of Eqs. (7,6) (that is, the one with vanishing initial data), since the homogeneous one will be exactly the same as at first order. Therefore, in this paper we always impose vanishing initial data for all the second-order modes and concentrate on those modes generated by first-order mode coupling. 1. Time Derivative (TD) 250 Ψ22 (t = 0, r) = M e−(r−r0 ) (12) {1} (13) {1} Ψ̇22 (t 2 /σ2 = 0, r) = (1 − 2M/r)∂r , Ψ22 (t = 0, r). {1} 7 V. DECAY FREQUENCIES OF SECOND-ORDER PERTURBATIONS {1} Ψ {2} Ψ source 10-5 {1} 10-10 10-15 ingoing ID 10-20 0 50 100 150 200 250 300 t/M {1} Ψ {2} Ψ source 10-5 10-10 {1} Ψ, {2}Ψ, source 100 10-15 outgoing ID 10-20 0 50 100 150 t/M 200 300 {1} Ψ {2} 100 Ψ source 10-5 10-10 {1} Ψ, {2}Ψ, source 250 10-15 time-derivative ID 10-20 0 50 100 150 200 250 300 t/M {1} Ψ {2} Ψ source Ψ, {2}Ψ, source 100 10-5 10-10 {1} Ioka and Nakano have put forward the suggestion that at second order new frequencies should appear in the ringdown spectrum, which would be given by the sum of different pairs of standard QNM frequencies. In particular, according to this prediction, the dominant frequency would correspond to the double of the standard fundamental one from linearized theory [40, 41]. This seems reasonable, since the sources for the second-order master equations are quadratic in the first-order modes, so one might well expect that frequencies get summed-up in Fourier space. The physical picture, however, appears to be at the same time more subtle and simpler: our numerical simulations indicate that in practice second-order perturbations decay with the standard QNM frequencies from linearized theory. Recall that the physical process here studied is the coupling of linear modes. That is, we initialize all secondorder perturbations to zero. The second-order master equations have sources which are, indeed, quadratic in the first-order modes. What we observe in all our simulations, though, is that those sources quickly excite the second-order perturbations and afterwards decay very fast in time. As a consequence, once the secondorder modes have been excited and reached the regime in which they oscillate with a constant complex frequency (i.e. what in linearized theory corresponds to the QNM regime), they essentially propagate with a vanishing source. In other words, as first-order perturbations do. And, in particular, they do oscillate and decay with the same, standard, QNM frequencies from linearized perturbation theory. We show this behavior in some detail for the four initial data types (Sec. IV A) of CASE A perturbations (Sec. III A), with r0 = 20M and σ = 4M . Figure 5 shows the first-order Zerilli function and the (ℓ = 2, m = 0) second-order one for different initial data profiles. In all cases the source decays much faster than the second-order solution itself and therefore its role in determining the behavior of the latter by the time it enters the QNM regime is negligible. From the same figure one can notice that the type of initial perturbation does determine the time by which the second-order Zerilli function enters the tail regime; but this is not surprising, since it already happens for the first-order one. In order to gain further insight into these observations we display the dynamics of first and second-order Zerilli functions and the source term {2}SΨ of Eq. (7), now for the (ℓ = 4, m = 0) second-order mode. Fig. 6 shows these three quantities as functions of radius at different times. The source term is dominant only during the first ∼ 20M , later decaying at a fast rate to several orders of magnitude below the second-order Zerilli function. Table II shows the fitted QNM frequencies from our numerical data, for the different initial-data profiles of CASE A perturbations, with r0 = 20M and σ = 4M . Ψ, {2}Ψ, source 100 10-15 time-symmetric ID 10-20 0 50 100 150 200 250 300 t/M Figure 5: CASE A simulations: first-order (ℓ = 2, m = 2) Zerilli function and second-order (ℓ = 2, m = 0) one, along with the source term for the second-order master equation [Eq. (7)], as functions of time for different initial data profiles. The source plays a role only at very early times in exciting the second-order modes, afterwards being much smaller than the first and second- order solutions. 8 102 {1} Ψ {2} Ψ source function absolute value 100 -2 10 10-4 10-6 ID TD TS IN OUT ℓ = 4, m = 0 0.80916 − 0.09418i 0.80920 − 0.09420i 0.80918 − 0.09416i 0.80931 − 0.09425i 2nd order ℓ = 2, m = 0 err % 0.37334 − 0.08883i 0.1 0.37335 − 0.08766i 0.3 0.37373 − 0.08945i 0.1 0.37074 − 0.08902i 0.8 10-10 t=0M -14 10 0 10 20 30 40 50 r/M 102 {1} Ψ {2} Ψ source 100 function absolute value 1st order ℓ = 2, m = 2 err % 0.37077 − 0.08826i 0.8 0.37353 − 0.08837i 0.2 0.37061 − 0.08887i 0.8 0.37107 − 0.08624i 1.0 10-8 10-12 10-2 10-4 10-6 10-8 10-10 10-12 10-14 t = 12 M 0 10 20 30 40 50 r/M 102 {1} Ψ {2} Ψ source 100 function absolute value ID TD TS IN OUT 10-2 10-4 10-6 10-8 10-10 10-12 10-14 t = 24 M 0 10 20 30 40 50 r/M Figure 6: CASE A simulations: snapshots of the first-order Zerilli function and second order (ℓ = 4, m = 0) one, along with the source term, for ingoing initial data. The generic behavior of the source is to rapidly decrease several orders of magnitude below the solutions themselves. [Note: In the first snapshot the second-order Zerilli function vanishes because it corresponds to t = 0 and our setup of the physical problem.] They agree quite well with those predicted by first-order theory for each of those modes. As described in Sec. IV, our numerical simulations are of very high accuracy (both at first and second-order): all the errors are at the level of double precision round-off, ∼ 10−14 −10−12 (see Figs. 2,3). Therefore, we do not consider lack of resolution as a possible reason for not finding traces of the predicted new second-order QNM frequen- 2nd order err % ℓ = 4, m = 4 0.003 0.80916 − 0.09418i 0.005 0.80920 − 0.09420i 0 0.80918 − 0.09416i 0.019 0.80931 − 0.09425i err % 0.003 0.005 0 0.019 Table II: Measured quasinormal frequencies from our numerical simulations (CASE A). They agree with those predicted by linearized theory, even for the second-order modes generated due to mode-mode coupling. The predicted QNM frequencies from standard linearized perturbation theory, as quoted in [15], are 0.37367 − 0.08896i for ℓ = 2 and 0.80918 − 0.09416i for ℓ = 4 (in linearized theory there is degeneracy with respect to the azimuthal index m. The relative errors for real and imaginary parts of the measured frequencies ω are computed as |ω − ωexact |/|ωexact |. cies in our simulations. Similarly, one might think that those predicted frequencies are in fact present, but with a very small amplitude. Figure 7 indicates that in practice that does not seem to be the case: the residual of the fit for the second-order Zerilli function [in the case shown it is the (ℓ = 2, m = 0) one, for TD initial data, first-order perturbations] —defined as the function minus its fit— does not appear at all to correspond to an oscillation and decay with twice the standard complex fundamental quasinormal frequency for that mode (or to an overtone). Instead, it appears to be the residual associated with the fact that QNMs are not complete. For completeness, in Appendix B we provide results of the fitted frequencies for the four initial data types of CASE A perturbations, now varying both the location and width of the initial data; all of them support the same conclusion. Finally, we briefly discuss the results of some CASE B [odd-parity (ℓ = 2, m = 0) first-order mode] perturbations, since the conclusions are identical. As discussed in Sec. III B and summarized in Table I, they generate both (ℓ = 4, m = 0) and (ℓ = 2, m = 0) even-parity second-order modes. The fitted frequencies from the numerical solutions (for a simulation of TD linear initial data with r0 = 20M and σ = 4M ) for these secondorder modes yield, respectively, 0.37441 − 0.08921i and 0.80932 − 0.09419i, to be compared with the expected values from perturbation theory: 0.37367 − 0.08896i for ℓ = 2 and 0.80918 − 0.09416i for ℓ = 4. 9 100 NSF PHYS 0925345, 0941417, 0903973, 0955825 to Georgia Tech, the Spanish MICINN Project FIS200806078-C03-03, the French A.N.R. Grant LISA Science BLAN07-1 201699, and by the Deutsche Forschungsgemeinschaft (DFG) through SFB/TR 7 “Gravitational Wave Astronomy”. residual solution -1 absolute value 10 10-2 10-3 10-4 10-5 10-6 120 130 140 150 160 We thank Emanuele Berti, Vitor Cardoso, Chad Galley and Bence Kocsis for helpful discussions and suggestions. We also thank Carlos Ann for inspiration, Tryst DC –where parts of this work were done– for hospitality, and Eric Cartman for proof-reading the manuscript and providing critical comments. t/M Figure 7: Residual from fitting a second-order Zerilli function to a complex frequency mode. The fitted frequency corresponds to the standard fundamental quasinormal one for that multipole index, and the residual does not appear to contain traces of twice that frequency (or any other). VI. FINAL REMARKS In this paper we have numerically evolved first and second-order self-generated gauge-invariant gravitational perturbations of Schwarzschild black holes with a variety of initial data sets, studying the oscillation and decay behavior of nonlinear modes and, more specifically, whether they correspond to the standard QNM frequencies or to a different spectrum. We have found, in all cases, that second-order modes decay through the standard QNM frequencies, and that the picture behind this is remarkably simple: first-order perturbations trigger high-order ones through source terms which afterwards rapidly decay in time. Besides, by the time the solutions reach the regime in which they oscillate and decay at a constant rate (the QNM regime in the case of linearized perturbations), the second-order modes for all practical purposes propagate as in linearized theory. Mode-mode coupling in the ringdown of black holes has been previously studied through numerical simulations of full Einstein equations [50, 51]; however, no conclusions seem to have reached or sought for in terms of the deviations in the ringdown spectrum from the linearized one (presumably due to lack of resolution). The fact that nonlinear aspects of Einstein’s equations in the ringdown of black holes appear to already be captured —at least in what oscillation and decay frequencies concerns— by linearized perturbation theory is somehow remarkable and should be of use both for modeling black holes in the ringdown regime as well as in data analysis searches of gravitational waves. Acknowledgments This research has been supported in part by NSF Grant PHYS 0801213 to the University of Maryland, Appendix A: Mode coupling: selection rules In this appendix we summarize, for completeness, the selection rules for mode coupling [45, 46]. A second-order (l, m)-mode gets a contribution from a pair of first-order modes (l̂, m̂) and (l̄, m̄) if three conditions are obeyed. First, the harmonic labels must be related by the usual composition formulas |l̂ − ¯l| ≤ l ≤ ˆl + ¯l, and m̂ + m̄ = m. (A1) Second, mode coupling must conserve parity. To any harmonic coefficient with label l, we associate a polarity sign σ such that, under parity, the harmonic changes by a sign σ(−1)l . Polar/even parity (axial/odd parity) harmonics have σ = +1 (σ = −1). Then, parity conservation implies the third condition: (−1)l̄+l̂−l = σσ̄σ̂, (A2) where σ̂ and σ̄ are the polarity signs corresponding to the modes (l̂, m̂) and (l̄, m̄) respectively. In particular, there is a special case in which the coupling of two modes satisfying Eqs. (A1) and (A2) does not contribute to a second-order mode, and the reason comes from the properties of the Clebsch-Gordan-like coefficients that appear in the product formula for the tensor harmonics [45]. In axisymmetry (m̄ = m̂ = 0) the mentioned ClebschGordan-like coefficients vanish if ¯l + ˆl + l is odd. This analysis can be extended to higher orders. In particular, the parity condition implies that a collection of k modes with harmonic labels {l1 , ..., lk } and polarities {σ1 , ..., σk } contributes to the mode (l, σ) only if (−1)l σ = Πki=1 (−1)li σi . 10 Appendix B: Numerical decay frequencies of second-order perturbations from time domain simulations r0 20 40 60 80 100 ℓ = 2, m = 0 0.37373 − 0.08949i 0.37147 − 0.08833i 0.37141 − 0.08834i 0.37102 − 0.08839i 0.37101 − 0.08840i Ingoing ℓ = 4, m = 0 0.80918 − 0.09416i 0.80916 − 0.09417i 0.80915 − 0.09413i 0.80917 − 0.09415i 0.80888 − 0.09414i ℓ = 4, m = 4 0.80918 − 0.09416i 0.80916 − 0.09417i 0.80915 − 0.09413i 0.80917 − 0.09415i 0.80872 − 0.09417i ℓ = 2, m = 0 0.37074 − 0.08902i 0.37079 − 0.08891i 0.37093 − 0.08844i 0.37140 − 0.08816i 0.37148 − 0.08782i Outgoing ℓ = 4, m = 0 0.80931 − 0.09425i 0.80911 − 0.09409i 0.80911 − 0.09417i 0.80913 − 0.09413i 0.80914 − 0.09426i ℓ = 4, m = 4 0.80931 − 0.09425i 0.80911 − 0.09409i 0.80911 − 0.09417i 0.80913 − 0.09413i 0.80914 − 0.09426i σ 2 0.37708 − 0.09116i 0.80918 − 0.09418i 0.80918 − 0.09418i 0.37365 − 0.08891i 0.80918 − 0.09416i 0.80918 − 0.09416i 4 0.37373 − 0.08949i 0.80919 − 0.09417i 0.80918 − 0.09417i 0.37079 − 0.08902i 0.80931 − 0.09425i 0.80911 − 0.09409i 8 0.37286 − 0.08893i 0.80875 − 0.09407i 0.80875 − 0.09407i 0.36887 − 0.08284i 0.80923 − 0.09599i 0.80921 − 0.09599i r0 20 40 60 80 100 ℓ = 2, m = 0 0.37334 − 0.08883i 0.37277 − 0.08893i 0.37378 − 0.08921i 0.37383 − 0.08995i 0.37387 − 0.08999i Time Derivative ℓ = 4, m = 0 0.80916 − 0.09418i 0.80919 − 0.09417i 0.80920 − 0.09418i 0.80915 − 0.09411i 0.80914 − 0.09409i ℓ = 4, m = 4 0.80916 − 0.09418i 0.80919 − 0.09417i 0.80920 − 0.09418i 0.80915 − 0.09411i 0.80914 − 0.09409i ℓ = 2, m = 0 0.37335 − 0.08766i 0.37376 − 0.08984i 0.37267 − 0.08948i 0.37344 − 0.08852i 0.37321 − 0.08823i Time symmetric ℓ = 4, m = 0 0.80920 − 0.09420i 0.80915 − 0.09413i 0.80917 − 0.09416i 0.80917 − 0.09413i 0.80923 − 0.09422i ℓ = 4, m = 4 0.80920 − 0.09420i 0.80915 − 0.09413i 0.80917 − 0.09416i 0.80917 − 0.09413i 0.80923 − 0.09422i σ 2 0.37360 − 0.08906i 0.80917 − 0.09417i 0.80917 − 0.09417i 0.37404 − 0.08972i 0.80917 − 0.09416i 0.80916 − 0.09416i 4 0.37334 − 0.08893i 0.80916 − 0.09417i 0.80916 − 0.09418i 0.37335 − 0.08766i 0.80920 − 0.09420i 0.80920 − 0.09420i 8 0.37315 − 0.08935i 0.80924 − 0.09342i 0.80924 − 0.09334i 0.37943 − 0.08760i 0.80867 − 0.09424i 0.80867 − 0.09423i Table III: For ℓ = 2 the expected fundamental frequency from linearized theory is 0.37367 − 0.08896i, while for ℓ = 4 frequency it is 0.80918 − 0.09416i [15]. Shown are the fitted values from our numerical solutions to the second order Zerilli equations, for four different sets of initial data. For each initial data configuration we vary the location r0 and width σ of the linearized perturbations, keeping one fixed (r0 = 20M and σ = 4M , respectively) while varying the other. The fitted frequencies appear to be insensitive to the choice of initial data and agree with the QNM frequencies from linearized perturbation theory. [1] M. Heusler, Living Reviews in Relativity 1 (1998), URL http://www.livingreviews.org/lrr-1998-6. [2] J. D. Bekenstein (1998), gr-qc/9808028. [3] R. P. Kerr, Phys. Rev. Lett. 11, 237 (1963). [4] C. V. Vishveshwara, Nature 227, 936 (1970). [5] C. V. Vishveshwara, Phys. Rev. D 1, 2870 (1970). [6] M. Davis, R. Ruffini, H. Press, and R. H. Price, Phys. Rev. Lett. 27, 1466 (1971). [7] W. H. Press, Astrophys. J. 170, L105 (1971). [8] R. Price, Phys. Rev. D 5, 2419 (1972). [9] T. Regge and J. Wheeler, Phys. Rev. 108, 1063 (1957). [10] F. J. Zerilli, Phys. Rev. Lett. 24, 737 (1970). [11] F. J. Zerilli, Phys. Rev. D. 2, 2141 (1970). [12] S. A. Teukolsky, Astrophys. J. 185, 635 (1973). [13] W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973). [14] E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), 0905.2975. [15] K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 1999 (1999), http://www.livingreviews.org/lrr-1999-2. [16] W. Krivan, Phys. Rev. D60, 101501 (1999), grqc/9907038. [17] E. Poisson, Phys. Rev. D66, 044008 (2002), grqc/0205018. [18] R. J. Gleiser, R. H. Price, and J. Pullin, Class. Quant. Grav. 25, 072001 (2008), 0710.4183. [19] M. Tiglio, L. E. Kidder, and S. A. Teukolsky, Class. 11 Quant. Grav. 25, 105022 (2008), 0712.2472. [20] L. M. Burko and G. Khanna, Class. Quant. Grav. 26, 015014 (2009), 0711.0960. [21] O. Dreyer et al., Class. Quant. Grav. 21, 787 (2004), gr-qc/0309007. [22] E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D73, 064030 (2006), gr-qc/0512160. [23] Éanna É. Flanagan and S. A. Hughes, Phys. Rev. D 57, 4535 (1998). [24] E. Berti, V. Cardoso, and C. M. Will, AIP Conf. Proc. 873, 82 (2006). [25] A. Buonanno, G. B. Cook, and F. Pretorius, PRD 75, 124018 (2007), arXiv:gr-qc/0610122. [26] E. Berti et al., Phys. Rev. D76, 064034 (2007), grqc/0703053. [27] B. Kocsis and A. Loeb, Phys. Rev. D 76, 084022 (2007). [28] A. Buonanno and T. Damour, Phys. Rev. D 62, 064015 (2000), [arXiv: gr-qc/0001013]. [29] A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 (1999), gr-qc/9811091. [30] T. Damour and A. Gopakumar, Phys. Rev. D 73, 124006 (2006), gr-qc/0602117. [31] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev. D75, 124018 (2007), gr-qc/0610122. [32] J. D. Schnittman et al., Phys. Rev. D77, 044031 (2008), 0707.0301. [33] R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, Phys. Rev. Lett. 77, 4483 (1996). [34] M. C. Miller and E. J. M. Colbert, Int. J. Mod. Phys. D13, 1 (2004), astro-ph/0308402. [35] J. Abadie et al. (LIGO Scientific), Class. Quant. Grav. 27, 173001 (2010), 1003.2480. [36] K. Wiersema et al. (2010), 1008.4125. [37] N. A. Webb et al., Astrophys. J. 712, L107 (2010), 1002.3625. [38] S. A. Farrell et al. (2010), 1002.3404. [39] S. A. Farrell, N. A. Webb, D. Barret, O. Godet, and J. M. Rodrigues, Nature 460, 73 (2010). [40] K. Ioka and H. Nakano, Phys. Rev. D76, 061503 (2007), 0704.3467. [41] H. Nakano and K. Ioka, Phys. Rev. D76, 084007 (2007), 0708.0450. [42] S. Okuzumi, K. Ioka, and M.-a. Sakagami, Phys. Rev. D77, 124018 (2008), 0803.0501. [43] E. N. Dorband, E. Berti, P. Diener, E. Schnetter, and M. Tiglio, Phys. Rev. D 74, 084028 (2006), grqc/0608091. [44] D. Brizuela, J. M. Martin-Garcia, and M. Tiglio, Phys. Rev. D80, 024021 (2009), 0903.1134. [45] D. Brizuela, J. M. Martin-Garcia, and G. A. Mena Marugan, Phys. Rev. D74, 044039 (2006), gr-qc/0607025. [46] D. Brizuela, J. M. Martı́n-Garcı́a, and G. A. M. Marugán, Phys. Rev. D76, 024004 (2007), gr-qc/0703069. [47] J. M. Martı́n-Garcı́a, Comp. Phys. Commun. 179, 597 (2008), URL http://metric.iem.csic.es/Martin-Garcia/xAct/. [48] D. Brizuela, J. M. Martı́n-Garcı́a, and G. A. Mena Marugán, Gen. Rel. Grav. (2009), 0807.0824. [49] J. Hesthaven, Appl. Numer. Math. 33 (1-4), 23 (2000). [50] G. Allen, K. Camarda, and E. Seidel (1998), grqc/9806036. [51] Y. Zlochower, R. Gómez, S. Husa, L. Lehner, and J. Winicour, Phys. Rev. D 68, 084014 (2003).