Quasi-normal Modes of Static Spherically Symmetric Black Holes in $f(R)$ Theory

We study the quasi-normal modes (QNMs) of static, spherically symmetric black holes in $f(R)$ theories. We show how these modes in theories with non-trivial $f(R)$ are fundamentally different from those in General Relativity. In the special case of $f(R) = \alpha R^2$ theories, it has been recently argued that iso-spectrality between scalar and vector modes breaks down. Here, we show that such a break down is quite general across all $f(R)$ theories, as long as they satisfy $f''(0) \neq 0$, where a prime denotes derivative of the function with respect to its argument. We specifically discuss the origin of the breaking of isospectrality. We also show that along with this breaking the QNMs receive a correction that arises when $f''(0) \neq 0$ owing to the inhomogeneous term that it introduces in the mode equation. We also discuss how these differences affect the"ringdown"phase of binary black hole mergers and the possibility of constraining $f(R)$ models with gravitational-wave observations.


I. INTRODUCTION
Recent observations of gravitational wave (GW) signals in LIGO and Virgo from compact object mergers have opened up a new chapter in the history of physics [1][2][3][4]. So far these detectors have already observed several binary black hole mergers and a binary neutron star merger [5]. The latter observation has put strong constraints on the equation of state of the matter inside neutron stars [6][7][8][9][10]. It has also provided a strong bound on the graviton mass [11] and the deviation of the velocity of GWs from that of the light [5,12]. As a result it has already helped us understand not only the mergers themselves but also fundamental aspects of the nature of gravity [12][13][14]. This scenario is expected to get even more interesting when LISA is launched in several years from now since the space detector will give us the opportunity to test gravity in a different frequency band [15][16][17][18][19].
These corrections come in various forms, such as functions of Ricci scalars and different combinations of the contractions of the Riemann tensor. These suggestions have gained importance due to other reasons as well. It has been shown that several actions with correction terms can drive inflation [31] and others have been successful in explaining the late time acceleration of the universe [32]. Such actions have also risen from the low-energy limit of quantum corrections or String Theory [33][34][35] and in Loop Quantum Gravity [36,37].
These results have initiated extensive research in such alternative theories of gravity in the last few decades. Now the availability of GW observations helps us examine if these theories can be subjected to useful tests.
As has been shown in Numerical Relativity, when the two black holes in a binary system merge, a distorted black hole is created. This remnant radiates GWs as it settles down into a Kerr black hole [38][39][40][41][42]. Asymptotically, these GWs can be expressed as superpositions of damped sinusoidal modes, termed as quasi normal modes (QNMs). These QNMs in GR depend only on a few parameters characterizing the black hole, namely, its mass and spin, for astrophysical black holes. This is a manifestation of the No-Hair theorem [43][44][45].
Therefore, observation of these QNMs in the future is an anticipated testing ground for the No-Hair theorem and alternative theories of gravity.
It is well known that in GR there are two modes contributing to GW observables. These are the odd (vector) and even (scalar) modes, and they share the net emitted gravitational energy equally. It has been argued in Ref. [46] that in the case of the (R+αR 2 ) action, where R is the Ricci scalar, QNMs emitted from black holes will not have the same spectrum as in GR. In this work, we show that such a break down in iso-spectrality is quite general across all f (R) theories, as long as they satisfy f ′′ (0) = 0, where a prime denotes derivative of the function with respect to its argument. We specifically discuss the origin of the breaking of isospectrality. We also show that along with this breaking the QNMs receive a correction that arises when f ′′ (0) = 0 owing to the inhomogeneous term that it introduces in the mode equation. We also discuss the resulting structural modification of the QNMs. Throughout the work we have adopted the {−, +, +, +} signature and set G = c = 1.
In Sec. II general results of an arbitrary f (R) theory and the modified Einstein equation in such theories are discussed. In Sec. III we present the perturbation equations in the context of GR around the Schwarzschild metric. In Sec. IV we discuss the perturbation of the modified Einstein equation in general f (R) theories. In Sec. IV A we deduce the different modes and the final perturbation equations in f (R) theories. Then in Sec. V we deduce the changes in the QNMs in comparison to GR. Finally, in Sec. VI we summarize our results.

II. GENERAL EQUATIONS FOR f (R) THEORIES
As mentioned above, here we focus our attention on f (R) theories of gravity. Our primary motivation for doing so is that they offer resolutions for certain limitations in the standard cosmological model and the fact that f (R) terms arise as correction terms from Loop Quantum Gravity. Another reason is that these theories are not affected by the Ostrogradsky instability [47]. Among all the f (R) theories we only focus on those that allow R = 0 solutions, for the obvious reason that such solutions exist in GR and are observationally relevant. We discuss this aspect in more detail later.
Let us begin by considering the gravity model that has a Lagrangian of the following form: where κ 2 = 8π. Varying this action with respect to the metric g µν yields the equation of motion [48], It is known that Eq. (2) has a constant curvature solution, R =R. In that case, the above equation takes the formR Taking its trace givesR Since our main goal here is to study perturbations of the Schwarzschild black hole, we focus only on such models of f (R) that haveR = 0 as a root of Eq. (4). It is straightforward to verify that GR is but one example of such theories.
Since our main goal is to look for deviations from GR, we will separately track the GR part of the f (R) action, namely, the R term. Therefore, we will express the Lagrangian as where φ(R) denotes terms in the gravity action beyond GR. Henceforth, we rename φ(R) as f (R). Thus, our working Lagrangian becomes R + f (R).

III. QNMS IN GENERAL RELATIVITY
For the purpose of our calculation we take the gauge invariant formalism discussed in [49][50][51][52]. Our focus is the Schwarzschild spacetime. In the case of a spherically symmetric spacetime, it is possible to split the manifold into an orbit space and unit sphere. We will be focussing on the four-dimensional black hole background. In that case the manifold splits into a two-dimensional orbit space and a two-sphere. We will use the co-ordinates (x a ) for the orbit space and (z A ) for the two-sphere. The covariant derivatives on the orbit space and the two-sphere are represented by D a and D A respectively. The d'Alembertian operator on the orbit space and the two-sphere are represented by˜ andˆ , respectively.
With proper choice of co-ordinates it is possible to cast the metric in the following form, where the overbar represents quantities of the unperturbed background, g(r) = 1 − 2M r , and dΩ 2 represents the metric on the unit two-sphere. The metric perturbation on the background can be expressed as follows [49][50][51][52], where ds 2 p represents the perturbed part of the line element. To simplify the problem, the metric functions are usually expanded in spherical harmonic basis. After linearizing the Einstein equations, they can be cast in terms of the co-ordinates of the orbit space alone. This problem has been well studied. We follow the formalism defined by Kodama, Ishibashi and Seto (KIS). With proper identification of the scalar and vector mode it is possible to show that for each multipole l ≥ 2 they satisfy the following equations [49][50][51]: where λ = (l−1)(l+2) 2 and the time dependence has been taken to be e −iωt . The subscripts S and V represent the scalar and the vector modes, respectively. Moreover, "GR" in the superscript distinguishes these modes from those in f (R).
These QNMs have been studied extensively in the literature. It has been shown that the transmission and the reflection coefficients of the V S and V V are equal. They make equal contributions to GWs asymptotically. It has also been shown that they share the same frequency spectrum. We will show below that this does not hold, in general, in f (R) theories.

IV. PERTURBATION IN GENERAL f (R)
Our main goal is to investigate the QNM structure of the Schwarzschild black hole.
Therefore, we take the background metric to be that given in Eq. (6). Now the perturbed metric can be written as whereḡ µν and h µν represent the background metric and its perturbation, respectively. Using the perturbed metric from Eq. (10) in the equation of motion Eq. (2), we find the equation for the perturbation as follows, where δR µν and δR are the perturbed Ricci tensor and Ricci scalar, respectively. Eq. (11) can be re-expressed as: where Taking the trace of the Eq. (11) it can be shown that where h is the trace of the perturbation. This is the massive scalar degree of freedom that is present in f (R) theories, unlike general relativity. Further details are available in Ref. [53][54][55]. This equation is equivalent to the equation of a massive scalar field with a source term that depends on the perturbation.

A. QNMs in f (R)
In the present work our main goal is to find the perturbed equation of Schwarzschild like solutions in f (R). For that reason we can fixR = 0 along side Eq. (4). This translates to This equation is a source-free massive scalar equation, with δR identified as the scalar field and the mass-squared being m 2 ≡ 1+f ′ 3f ′′ . For simplicity, from now on we will omit the 0 from the argument value of f (0) and all its derivatives. This massive scalar longitudinal mode [56,57] is absent in GR. Now the perturbed Einstein equation takes the following form: where where β = f ′′ 1+f ′ . Heretofore, we will use both expressions interchangeably. But important thing to notice is that there are possible f (R) theories, such as αR 3 , for which f ′′ = 0; in that case β becomes zero. Another key point is that, the information of a particular f (R) theory enters only through the β. So, apart from β, dynamically speaking all f (R) theories are same. Now we use the separation of variables to separate out the angular dependence. For this reason we write the perturbed Ricci scalar as, where S(z A ) are the scalar spherical harmonics. For radial coordinate we use the tortoise coordinate (r * ).
The time dependence is taken to be e −iωt . After simplification, Eq. (14) takes the following form: This result for f (R) = αR 2 theory was discussed in Ref. [46].
It is well known that each component of a tensor behaves differently under rotation group.
As result, the effective energy momentum tensor can be separated into scalar, vector and the tensor modes [49][50][51] as follows, where S, V A and T AB are the pure scalar, vector and tensor spherical harmonics, respectively. Rest of the tensors are defined from the pure spherical harmonic tensors. Further details, including the expression for Θ T can be found in Ref. [49][50][51] and in A.
We find the several components of the energy momentum tensor from Eq. (20) as listed below: where δR = Ω(x a )S(z A ).
The scalar and vector master equations get modified as follows: where and H ≡ k 2 − 2 + 6M r , It should be noted that as β → 0 these equations, appropriately, reduce to their GR counterparts.
We also find that the tensor mode gives rise to the following equation [51], where˜ represents the (r * − t) part of the d'Alembertian operator and Above, λ L is the eigenvalue of the Lichnerowicz operator [51]. Therefore we see that the tensor and the vector modes do not undergo any modification. This is understandable because the effective energy-momentum tensor arises from the different combinations of the derivatives of the Ricci scalar. It has been shown in Ref. [58] that the tensor harmonic functions are identically zero in transverse-traceless (TT) gauge. For this reason we will not pursue the study of the tensor mode any further.
Since the background in our studies is a R = 0 solution, its metric does not have any extra hairs. However, we notice from Eq. (22) that even though the background black hole metric has no additional hair and is similar to the black holes in GR, at the level of perturbations the Zerilli mode gets modified from GR due to an extra source term. This source term originates because at the level of perturbations the extra massive mode that was absent in the background now gets excited. The excited massive mode "generates" hair via β. This tells us that the perturbations will therefore get modified from GR. As result we can expect that the detection of QNMs through gravitational wave observations will, in principle, provide us with the opportunity to distinguish between GR and f (R) theories of gravity.

V. SOLUTION FOR QNM
A. Solution at radiation zone As we can see from the structure of the equations, the scalar equation gets modified by a source term. Therefore, it is understandable that the solution of that equation will have two parts: (a) One of these will be for the homogeneous part of the solution, which is identical to the GR solution; (b) The second part will be the solution of the inhomogeneous part.
Solution of the vector equation will be exactly equal to the GR solution. Therefore, we can This implies that from the observational perspective the nonzero inhomogeneous part will be the signature of any deviation from GR. If we take the r → ∞ limit then we find that, where, H = k 2 − 2. However, due to Eq. (18), the source term becomes Equation (30) has solution of the form e iχr * , with From Eq. (22) we find, This gives the form of the inhomogeneous term arising due to the f (R) theory: Therefore, one can see that there is a deviation from GR in this f (R) theory, howsoever tiny its effect might be on observables.
In Eq. (33) we derived the asymptotic deviation of the f (R) QNMs from GR. This provides the opportunity to constrain f (R) theories by using it. Theoretically, there exists the possibility of constraining β by tracking how much energy is lost by QNMs and if it is completely accounted for in GR. Solar system tests have already constrained f (R) theories (details can be found in appendix E). In comparison, the results found in the current work lend themselves to possible tests using GW observations. The inhomogeneous term arises due to the existence of the massive polarization mode, as seen from Eq. (15) and Eq. (17).
This occurs because the scalar mode and the massive mode are coupled, as was shown in Eq. (22). Through this inhomogeneous term, energy transfer can occur between the massive mode and the scalar mode. Therefore, tests for non-GR polarizations and energy content in each polarization may, in principle, constrain these f (R) theories.
In spite of the theoretical possibilities of constraining the f (R) theories discussed above, practically such tests may not be realizable. The solar system test suggests β to be very small (β ≤ 2 × 10 −9 m 2 ) and on top of that the inhomogeneous part gets highly damped for ω 2 ≪ 1 3β . Only when the mode frequency reaches of the order of m this damping will be lessened. Considering the current solar system tests, this range is ω ∼ 10 9 kHz. This frequency range is well beyond the observational bands of LIGO, VIRGO or LISA. Due to this damping it becomes practically impossible to observe these modes. If it were possible to compare the mass-energy of each constituent in the binary with that of the final remnant and the energy emitted in GWs then one might have been able to test for the existence of energy lost in other polarization modes. However, it is highly unlikely for these high frequency modes to get excited.
Another observational aspect of these findings is the breaking of the iso-spectrality. It has been shown in [42] that the scalar and vector mode have identical spectra in GR. But now due to the source term it becomes evident that in general f (R) theories this iso-spectrality will not be satisfied. The homogeneous solution will have an iso-spectral solution to the vector mode but the inhomogeneous part will not be iso-spectral. This point has been discussed in V B. They claimed this property will be there for all f (R) theories. We prove above that in general this iso-spectrality will be violated as long as the theory has nonzero β. The details of the breakdown of the iso-spectrality has been discussed in V B.

B. Broken Iso-spectrality
Breaking of iso-spectrality between the scalar and vector modes in the context of αR 2 theory has already been discussed in [46]. Even though they claim it to be a universal phenomenon for f (R) theories, we showed above that this is not true unless β = 0. For this reason we will investigate for the first time in detail exactly how the iso-spectrality breaks down when β = 0.
Iso-spectrality in GR was first discovered by Chandrasekhar [42]. We will follow his methods for this investigation. Say, Z 1,2 satisfy the following equations w.r.t. r * : where V 1 and V 2 are the corresponding potentials. Using the ansatz that Z 1 = pZ 2 + qZ ′ 2 and taking its derivative twice it is possible to show where the prime denotes derivative w.r.t. r * . Comparing it with Eq. (35) Chandrasekhar showed that Now, finding a solution for p and q establishes a relation between Z 1 and Z 2 . For the QNMs of the Schwarzschild solution in GR it is possible to identify, It is also possible to find p and q. As a result, it is one can show that where λ = (l−1)(l+2)

2
. Due to this relationship it becomes evident that if Φ S depends on the radial coordinate as e iωr * then so does Φ V . Therefore, they share the same spectrum [59].
By inspecting Eq. (35) it is clear that the reason for such a simplification is its homogeneous nature. If we modify Eq. (35) with a source term, as in then an ansatz of the previous form, Z 1 = pZ 2 + qZ ′ 2 , again leads to But comparison with Eq. (39) does not lead to equations of p, q that are independent of Z 1 .
So, in the sourced case it is not possible to write Z 1 = pZ 2 + qZ ′ 2 , resulting in the violation of iso-spectrality.
One key point is worth mentioning. It is always possible to separate Z 1 as and Therefore, it will be possible to write Z GR 1 = pZ 2 + qZ ′ 2 and, as a result, they will have an iso-spectral solution. Consequently, Φ S in f (R) theories will have a part that is isospectral to Φ V and another part that arises solely due to the source term that does not obey iso-spectrality.

C. Massive mode
As already shown above in Eq. (14) there exists a massive mode in f (R) theories. The equation for this scalar mode is exactly equivalent to the equation of a scalar field around a Schwarzschild black hole in GR. The QNM of a massive scalar field around a Schwarzschild black hole in GR has been studied extensively [62,63]. Because of the similarity between the two problems, the results for the QNM frequency of the massive mode relevant for our current work will be exactly the same as the ones found in those works. Naturally, the QNM frequency values found in those works depend on the mass of the scalar field. In our current work that mass depends on the specific f (R) theory chosen. Even though there are stringent restrictions on the value of this mass from solar system observations E, it is possible that further restriction may be found from GW observations. Current constraints on the mass imply that if an f (R) theory is the correct theory of gravity then, m ≥ 1.25 × 10 4 m −1 .
One key point is worth mentioning here. Usually ω is a complex number, which implies that χ should be such a number as well. Consequently, the QNMs are damped. On top of that, the large value of the mass of the massive mode will give rise to a strong damping pattern. As a result, the contribution of the inhomogeneous part in Eq. (33) will be very less. But massive scalar fields in Schwarzschild backgrounds can have purely real ω [64].
This originates mainly from the sub-dominant asymptotic contribution arising due to the irregular singularity at infinity [64]. These modes are called quasi-resonance frequency modes [63]. The requirement for the existence of such modes is ω QRM < m [64], where m is the mass of the massive mode. Further numerical study revealed that mass has a crucial role in the damping of that mode. They showed that the greater the mass the lesser is the damping rate. As result purely real modes originate corresponding to the non-damping oscillation.
They also showed that for a given mass of the field a certain number of lower overtones disappear. However, this disappearance happens only for the lower overtones, while the remaining overtones are still damped. As result due to the non-zero mass of the massive mode the inhomogeneous term will contain some quasi-resonance contribution that will not decay as fast as the other damped overtones.

VI. DISCUSSION
In this work we have focused on the QNMs of a Schwarzschild black hole in a general f (R) theory that allows the R = 0 solution. We came to the conclusion that in general f (R) theories the iso-spectrality between the scalar and the vector mode will not be valid.
Breaking of iso-spectrality was discussed in the work [46] for the specific case of αR 2 theory.
They claimed this result will be true for all other theories also. We found that this claim is not entirely true. In general, there will be breaking of isospectrality, but it is possible to have theories where it is not true, such as when f (R) = αR 3 . Indeed, it is possible to conclude that the isospectrality will be violated if f ′′ = 0.
Secondly, we discussed that the structure of QNMs gets modified at the asymptote. This modification arises as an inhomogeneous term alongside the GR contribution. Therefore, not only the iso-spectrality will be violated but also the QNM structure will get modified. But for f ′′ = 0, this contribution along with the breaking of the isospectrality will be absent. As a result one can in principle constrain β solely from GW observations. In practice, however, this may be difficult to realize.
One key-point is worth mentioning here. In our current work we have focused only on the perturbation around a vacuum R = 0 solution. We have used these similar assumptions while interpreting the solar system test. Owing to this choice we have not been able to pursue the study of Chameleon screening in f (R) theories [68]. This screening behaviour is the origin of the different behaviour of the theory at different scales. Due to the screening effect the mass of the massive scalar mode depends on the mass density of energy momentum tensor (T µν ) in the space-time. Effective mass of the massive scalar mode becomes heavier in the region with high mass density and lighter in the low density region. Therefore, in different environments the mass of the massive scalar mode will have different values depending on the density profile. This gives rise to the screening mechanism. Therefore, interpreting and using the solar system constraints should be done carefully. But unfortunately due to vacuum assumption (T µν = 0), we have not been able to explore this sector in our current work. Hence, a rigorous analysis with non-vacuum condition is an important endeavour that may be taken up in the future.
QNMs in scalar tensor theories have recently been calculated in Refs. [66,67]. It has been shown before that a class of scalar tensor theories are equivalent to f (R) theories.
Tattersall et al. have also found results similar the ones noted here in the context of scalar tensor theories. Due to the equivalence of scalar tensor theories and f (R) theories this is, of course, expected. The equivalence between Brans-Dicke with ω 0 = 0 and f (R) theory is valid under the condition f ′′ (R) = 0 (related to the sufficient condition for invertibility between the two theories). When f ′′ is not defined, or vanishes, the equality φ = f ′ (R) and the equivalence between the two theories cannot be guaranteed, although this is not a priori excluded by f ′′ = 0 [69]. Therefore, there are families of f (R) that in principle are not equivalent to scalar-tensor theories. Interestingly, in our work we see that this is precisely the sector where f (R) does not bring any modification w.r.t. GR, mainly due to the vanishing coupling between the GR mode and the extra massive scalar mode.

ACKNOWLEDGMENTS
We thank S. Shankaranarayanan for his useful inputs. This work is supported in part by the Navajbai Ratan Tata Trust. SD would like to thank University Grants Commission (UGC), India, for financial support as senior research fellow.

Appendix A: Spherical Harmonics
Any field on a spherically symmetric background can be classified as a scalar, vector and tensor depending on their transformation under the rotation group. For this reason three kinds of spherical harmonics are used to separate out the angular dependence.
One subset of these functions are the scalar spherical harmonics. Any scalar function on a spherically symmetric spacetime is expressible as a sum of the scalar spherical harmonics, Y lm (θ, φ). We call them S(Z A ).
A vector field U A can be expressed as, U A = V A + D A S, where S is the scalar spherical harmonic defined above. V A are called the vector spherical harmonics.
A tensor field X AB on spherically symmetric spacetime can be expressed as, where γ AB is the metric on the two sphere. V A and S are the vector and scalar spherical harmonics, respectively.L AB is the Lichnerowicz operator defined as,L AB = D A D B − 1 2 γ AB . S, V A and T AB are called the pure spherical harmonic functions.
They satisfy the following equations, where k 2 = l(l + 1), and l is a non-negative integer.
From the pure spherical harmonics several useful vector and tensor fields can be defined.
We only mention a few among them, namely, Appendix B: Separation of effective energy momentum tensor Let us assume that the metric of a manifold is as follows, In this case the Christoffel symbols take particular forms as discussed below. A hat represents connection on two sphere and 2 Γ a bc represents the connection in the orbit space. The nonzero components of the connection are as follows, Γ a bc = 2 Γ a bc , Γ a BC = − rD a rγ BC , (B2) Using these results for connection we can calculate various covariant derivatives. Various components of double co-variant derivatives of a scalar function F (x a , z A ) can be calculated to be as follows, It has been shown that the Einstein equation gets modified in our case. As result we get an effective energy momentum tensor of the following form, All other terms arising in Eq. [20] are zero.
The structure of the metric [54] gives rise to the potential of Yukawa type, as was already discussed in [54]. In fifth-force tests such potential has been studied extensively, where the considered form is, From the Eöt-Wash experiment [70,71] it is possible to put strong constraint, λ = m −1 ≤ 8 × 10 −5 m. Therefore, we are finding that for general f (R) this puts up the bound |β| = f ′′ 1+f ′ ≤ 2 × 10 −9 m 2 .