Quasinormal modes and their anomalous behavior for black holes in f(R) gravity

We study the propagation of scalar fields in the background of an asymptotically de Sitter black hole solution in f(R) gravity. The aim of this work is to analyze in modified theories of gravity the existence of an anomalous decay rate of the quasinormal modes (QNMs) of a massive scalar field which was recently reported in Schwarzschild black hole backgrounds, in which the longest-lived modes are the ones with higher angular number, for a scalar field mass smaller than a critical value, while that beyond this value the behavior is inverted. We study the QNMs for various overtone numbers and they depend on a parameter β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} which appears in the metric and characterizes the f(R) gravity. For small β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}, i.e. small deviations from the Schwarzschild–dS black hole the anomalous behavior in the QNMs is present for the photon sphere modes, and the critical value of the mass of the scalar field depends on the parameter β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} while for large β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}, i.e. large deviations, the anomalous behavior and the critical mass does not appear. Also, the critical mass of the scalar field increases when the overtone number increases until the f(R) gravity parameter β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} approaches the near extremal limit at which the critical mass of the scalar field does not depend anymore on the overtone number. The imaginary part of the quasinormal frequencies is always negative leading to a stable propagation of the scalar fields in this background.


Introduction
Modified theories of gravity in which the Einstein-Hilbert action is replaced with a generic form of f (R) gravity have been introduced in an attempt to describe the early and late cosmological evolution. Another motivation to study such theories is the understanding of the existence of dark energy and dark matter consistent with the recent observations, without introducing new material ingredients that have not yet been detected by experiments [1][2][3][4][5][6]. Modifying the action not only affects the dynamics of the universe, it can also alter the dynamics at the galactic or solar system scales. Therefore, modified theories of gravity with curvature corrections provide a deeper understanding of General Relativity (GR).
Assuming that the gravitational Lagrangian is not only a linear function of R, variable modified theories of gravity were considered describing the cosmic evolution at early times, in which the gravitational Lagrangian contains some of the four possible second-order curvature invariants. Also models in which higher order curvature invariants as functions of the Ricci scalar were introduced in the gravitational Lagrangian, resulting in various f (R) gravity models [7][8][9][10][11][12][13][14][15][16]. Although such theories exclude contributions from any curvature invariants other than R, they could also avoid the Ostrogradski instability [17] which proves to be problematic for general higher derivative theories [18].
One of the first modifications of the Einstein Lagrangian density was proposed in [19]. A more natural modification of the Einstein Lagrangian is to add terms R n , like the Starobinsky model f (R) = R + α R 2 [20][21][22]. For n < 0 such corrections become important in the late universe and can lead to self-accelerating vacuum solutions [23][24][25][26]. However, these models suffer from instabilities [27,28] and there are strong constraints from the solar system [29]. A wide range of phenomena can be explained by considering different f (R) functions. Some discussions have been performed, such as on gravitational wave detection [30,31], early-time inflation [32], cosmological phases [33][34][35], the singularity problem [36], the stability of the solutions [37][38][39], and various other branches have been studied [40].
Spherically symmetric static solutions have been studied in f (R) theories of gravity. In particular, it was shown that for a large class of models the Schwarzschild-de Sitter (SdS) metric is an exact solution of the field equations. However, in addition to the SdS metric, f (R) theories typically also have new different solutions [41]. Also, static spherically symmetric perfect fluid solutions have been studied by matching the outside SdS-metric to the metric inside the mass distribution that leads to additional constraints that limit the allowed fluid configurations [42], while in [43] black hole solution were found with and without electric charge, and exact spherically symmetric solutions were discussed in [44][45][46][47][48][49]. Also, exact charged black hole solutions in f (R) gravity theories with dynamic curvature in D-dimensions were discussed in [50].
A way to probe the behavior and the stability of black holes is the study of their quasinormal modes (QNMs) and quasinormal frequencies (QNFs) [51][52][53][54][55]. The QNMs depend on the black hole parameters and probe field parameters, and on the fundamental constants of the system and they are independent of the initial conditions of the perturbations. The QNMs are described by complex frequencies, ω = ω R +iω I , in which the real part ω R determines the oscillation timescale of the modes, while the complex part ω I determines their exponential decaying timescale (for a review on QNM modes see [53,56]).
The QNMs and QNFs have been calculated using various numerical and analytical techniques [57][58][59][60][61][62]. In the case of gravitational perturbations it was found that for the Schwarzschild and Kerr black hole background the longestlived modes are always the ones with lower angular number . In the case of a massive probe scalar field it was found [63][64][65][66], at least for the overtone n = 0, that, if the mass of the scalar field is small, then the longest-lived QNMs are those with a high angular number , whereas if the mass of the scalar is large the longest-lived modes are those with a low angular number . This inverted behavior appears when the mass of the scalar field exceeds a critical value, which corresponds to the value of the scalar field mass where the behavior of the decay rate of the QNMs is inverted and can be obtained from the condition I m(ω) = I m(ω) +1 in the eikonal limit, that is when → ∞. This anomalous decay rate for small mass scale of the scalar field was recently discussed in [67].
The study of the anomalous behavior of QNMs was extended to other asymptotic geometries, such as, Schwarzs -child-de Sitter and Schwarzschild-AdS black holes in [68]. It was found that the same behavior is present in the Schwarzschild-de Sitter background, i.e., the absolute values of the imaginary part of the QNFs decay when the angular harmonic numbers increase if the mass of the scalar field is smaller than the critical mass, and they grow when the angular harmonic numbers increase, if the mass of the scalar field is larger than the critical mass. Also it was found that the increase of the value of the cosmological constant results in the increase of the value of the critical mass of the scalar field. It was also found that the anomalous behavior is present in Schwarzschild-AdS black holes backgrounds; however, there is not a critical mass where the behavior in the decay rate is inverted. Also anomalous decay of QNMs in accelerating black holes was recently reported in [69], as well as for Reissner-Nordström black holes [70]. Furthermore, the anomalous behavior of QNMs was observed for fermionic field perturbations in Schwarzschild-de Sitter black holes [71].
Stability studies have been performed to f (R) gravity theories. However, the stability analysis appears to be very difficult to be applicable to f (R) gravity theories because it involves fourth-order derivative terms in the linearized equations [72]. One way to evade this problem is to transform f (R) gravity into the corresponding scalar-tensor theory to remove the fourth-order derivative terms [73]. However, it is well known that a non-minimally coupled scalar makes the linearized field equations very intricate when compared to a minimally coupled scalar in the context of General Relativity [74]. To avoid this problem one may use conformal transformations to find the corresponding theory in the Einstein frame where a minimally coupled scalar field appears. The resulting theory can be considered as a scalar-tensor theory with a geometric (gravitational) scalar field. Then it was shown in [75,76] that this geometric scalar field cannot dress a f (R) black hole with hair, therefore, the no-hair theorem is respected in these models. Matter was introduced in [77] and a R 2 gravity model was discussed having a conformally coupled scalar field with associated Higgs-like potential and a static spherically symmetric black hole was found. However, the conformally coupled scalar field could not yield black hole hair, because it is regulated by the same parameters as appear in the metric. A scalar field minimally coupled to gravity with a self-interacting potential in a f (R) theory was discussed in [78]. An exact hairy black hole solution was found similar to the BTZ black hole. Also the process of scalarization of f (R) gravity theories was investigated in [79].
As we already discussed, knowledge of the QNMs and QNFs can give us vital information of the behavior, properties and the stability of black holes. In this work we will consider a specific f (R) theory which accepts a black hole solution [80]. In this theory the Ricci scalar receives nonlinear correction terms which they introduce new parameters in the metric functions of the black hole solution, expressing the change of the curvature altering the strength of the gravitational force. In this theory we will consider a massless and a massive test scalar field scattered off the background black hole and we will calculate the QNMs and QNFs. The motivation of this work is to calculate the scalar field perturbations and study the behavior of QNMs and QNFs and see how they are effected by the presence of non-linear curvature terms which appear explicitly in the metric functions.
The manuscript is organized as follows: In Sect. 2 we give a brief review of the model considered for f (R) gravity and its black hole solution. In Sect. 3, we study the scalar field stability and we calculate the QNFs of scalar perturbations numerically by using the pseudospectral Chebyshev method for each family of modes. In Sect. 3.1.3, we perform an analysis using the WKB method to get some analytical insight. Finally, our conclusions are in Sect. 4.

f (R) modified gravity
We consider a generic action in f (R) gravity depending on the Ricci scalar R given by where S m is the matter content of the theory. In [80] a specific ansatz for the function f (R) was considered: where corresponds to the cosmological constant, R c is a constant of integration and R 0 = 6α 2 /d 2 , with α and d being free parameters. Then considering the spherically symmetric metric where A(r ) = B(r ) −1 , a vacuum spherically symmetric solution was found in [80] with B(r ) = 1 − 2M r + βr − r 2 3 , where β = α/d is a real constant. For the range R >> and R/R 0 >> 2/α the action reduces to and this is satisfied in the solar system range of r << d.
On the other hand, at large scales, for α << 1 and R ∼ R 0 ∼ , the action tends to f (R) ∼ R + and agrees with cosmological observations. For the metric (3) the curvature scalar is given by and the action tends to a constant value asymptotically, It is worth to mention that the solution in both regimes was obtained as perturbation around the Einstein-Hilbert action and it was shown that within this approximation the solutions are consistent with the modified gravity equations. Also, the parameters of the theory were constrained by using the pioneer data, the flat rotation curve of galaxies and the cosmic microwave background radiation and Supernova Type Ia gold sample data giving β ≈ 10 −26 m −1 , βd ≈ 10 −6 and βd << 10 −3 , respectively [80]. However, the spacecraft's heat loss is a more plausible explanation for the anomaly than the cosmological one [81]. Our analysis considers a range of values of β, such that the metric represents a black hole solution with an event horizon and a cosmological horizon as it shown in Fig. 1. In this figure for a fixed value of the cosmological constant, we show the transition among black holes, near extremal (r H ≈ r ) and naked singularity, when the β parameter goes from positive to negative values. Notice that, for β = 0, the spacetime is described by a Schwarzschild-dS black hole. So, for β positive (negative) the cosmological horizon radius is larger (smaller) than the Schwarzschild-de Sitter black hole, while for β positive (negative) the event horizon radius is smaller (larger) than the Schwarzschild-dS black hole, when the spacetime describes black hole solutions.

Scalar perturbations
In order to obtain the QNMs of scalar perturbations in the background of the metric (3) we consider the Klein-Gordon equation, with suitable boundary conditions for a black hole geometry.
In the above expression m is the mass of the scalar field ϕ. Now, by means of the ansatz the Klein-Gordon equation reduces to where = 0, 1, 2, . . . corresponds to the eigenvalue of the Laplacian on the two-sphere. Now, defining R(r ) = F(r ) r and by using the tortoise coordinate r * given by dr * = dr B(r ) , the Klein-Gordon equation can be written as which corresponds to a one-dimensional Schrödinger-like equation with an effective potential V eff (r ), parametrized as V eff (r * ) and given by and substituting the metric function B(r ) in the effective potential in Eq. (11) we obtain Now, by considering m = m c , as the value of the mass that cancels the divergence of order r 2 , i.e., m c = √ 2 /3, and = 0, the effective potential is Therefore we note that for β negative (positive) the potential diverges positive (negative) at infinity, and in principle the existence of a critical mass could not be present in this case, due to the arguments given in [68]. Now, in order to solve numerically the differential equation (9) by using the pseudospectral Chebyshev method [82], it is convenient to perform the change of variable y = (r − r H )/(r − r H ). Thus, the values of the radial coordinate are limited to the range [0, 1], and the radial equation (9) can be written as Also, in order to propose an ansatz for R(y) we consider its behavior in the vicinity of the horizon (y → 0), and at the cosmological horizon (y = 1), So, an ansatz that satisfies the requirements of only ingoing waves at the event horizon and at the cosmological horizon is Then, by inserting the above ansatz for R(y) in Eq. (14), it is possible to obtain an equation for the function F(y).
The solution for the function F(y) is assumed to be a finite linear combination of the Chebyshev polynomials, then it is inserted in the differential equation for F(y). Also, the interval [0, 1] is discretized at the Chebyshev collocation points. Then the differential equation is evaluated at each collocation point. Thus, a system of algebraic equations is obtained, and it corresponds to a generalized eigenvalue problem, which is solved numerically to obtain the QNFs. It is worth mentioning that, for β = 0, the spacetime is described by the Schwarzschild-de Sitter black hole. The complex family of QNMs for this geometry was determined in [83] by using the WKB and Pöschl-Teller method. Also, it was shown that the frequencies all have a negative imaginary part, which means that the propagation of scalar field is stable in this background, and the presence of the cosmologicalconstant leads to decrease of the real oscillation frequency and to a slower decay. High overtones were studied in [84]. Also, a family of purely imaginary modes were reported in [85] and an analysis of the photon sphere (PS) modes and de Sitter (dS) modes was recently performed in Ref. [68].

Photon sphere modes
The photon sphere QNMs are represented by damped oscillations whose decay rate (ω I ) is directly connected to the instability timescale of null geodesics at the photon sphere. The photon sphere is a spherical trapping region of space where gravity is strong enough so that photons are forced to travel in unstable circular orbits around a black hole. This region has a strong pull in the control of decay of perturbations and the QNMs with large frequencies. For instance, the decay timescale is related to the instability timescale of null geodesics near the photon sphere. For asymptotically dS black holes we find a family that can be traced back to the photon sphere and refer to them as PS modes. The dominant modes of this family are approached in the eikonal limit, where → ∞, for the mass of the scalar field lower than a critical value, and can be very well approximated with the WKB method. For this family of modes, higher angular momentum represents smaller decay rates, such that in the range of interest, the most representative modes are those for which → ∞ [86,87]. Now, in order to observe the behavior of photon sphere modes, we plot the behavior of the fundamental QNFs n P S = 0 1 of massless scalar fields as a function of Mβ with = 0; see Fig. 2. We can observe that, for small values of Mβ, the decay rate decreases when the cosmological constant increases. However, for large Mβ, the behavior is opposite, i.e., the decay rate increases. On the other hand, the frequency of the oscillation increases when Mβ increases and decreases when the cosmological constant increases. In Fig. 3 we consider = 1; note that the behavior of the QNFs is equivalent to the observed for = 0 and small Mβ; however, there is no inverted behavior, for the range of Mβ values considered. Thus, when the cosmological constant increases the decay rate decreases, for small Mβ and when Mβ increases the decay rate converges to the same value and it does not depend on the value of the cosmological constant. On the other hand, the frequency of the oscillation increases when Mβ increases and decreases when the cosmological constant increases, and the same behavior was observed for = 0.

Anomalous decay rate
Now, we will study if in f (R) modified gravity, in the presence of a cosmological constant, the photon sphere modes show an anomalous behavior, i.e. the longest-lived modes are the ones with higher angular number, and also the presence of a critical mass, where beyond this value the behavior mentioned is inverted. So, we plot the behavior of −I m(ω)M as a function of the scalar field mass for different values of the parameters Mβ, and n PS = 0 (left panels), and n PS = 1 with > n PS (right panels), 2 for a fixed values of the black hole mass and the cosmological constant; see  To explain this effect, we could say that the deviation of Schwarzschild-dS black hole, where the anomalous behavior and the critical mass have been observed, is given by the parameter Mβ which also appears in the metric function. So for small values of this parameter, i.e. small deviations Schwarzschild-dS black hole, the anomalous behavior occurs and the critical mass appears. However, for larger values of the parameter Mβ i.e large deviations from Schwarzschild-dS black hole, the QNMs show a different behavior, where the anomalous behavior and the critical mass are not present, this cut off point could numerically occur for Mβ ≈ 0.10; see Fig. 2, where the QNFs present an inverted behavior with respect to Mβ for massless scalar field.

Analysis using the WKB method
In this section, we use the method based on Wentzel-Kramers-Brillouin (WKB) in order to get some analytical insight of the behavior of the QNFs in the eikonal limit → ∞ [88][89][90][91][92][93]. It is worth mentioning that this method can be used for effective potentials which have the form of potential barriers that approach a constant value at the horizon and spatial infinity [55]. The method considers the behavior of the effective potential near its maximum value r * max . So, by using a Taylor series expansion, the potential around its maximum is given by where corresponds to the ith derivative of the potential with respect to r * evaluated at the maximum of the potential r * max . Using the WKB approximation up to sixth order the QNFs are given by the following expression [94]: where (2) ) 9/2 (77 + 188N 2 ) (2) ) 7/2 (51 + 100N 2 ) (2) ) 3/2 (5 + 4N 2 ) , and N = n PS + 1/2, with n P S = 0, 1, 2, . . . , is the overtone number. Now, by defining L 2 = ( + 1), we find that, for large values of , the maximum of the potential is approximately at The second derivative of the potential evaluated at r * max yields the expression where and the higher derivatives of the potential evaluated at r * max yield the expressions Then, using the expressions above with Eq. (20), we obtain the following analytical expression for the QNFs for large values of and small values of Mβ: where , The term proportional to 1/L 2 is zero at the value of the critical mass m c , which is given by (24) and it is valid for small values of β and n PS = 0. In all the above expressions we have performed a Taylor expansion around β = 0 in order to obtain expressions that will be easy to handle analytically. For β = 0 we recover the result of critical mass for Schwarzschild-dS black holes [68], and for M 2 = 0.04 and n PS = 0 we obtain from (24) Table 1, the values obtained with the two methods. Also, we show the relative error, which is defined by where ω 1 corresponds to the result from Eq. (23) and ω 0 denotes the result with the pseudospectral Chebyshev method. We can observe that the error does not exceed 0.0008 % in the imaginary part, and 0.0010 % in the real part. Also, as it was observed, the frequencies all have a negative imaginary part, which means that the propagation of massive scalar fields is stable in this background. As we mentioned, the Taylor series have been performed around β = 0, so the WKB analysis is true for small values of β. In Table 2, we show the error in the quasinormal frequency between the WKB and pseudospectral method, for different values of m M, and Mβ. Note that for Mβ = 0.20 the error is 18% in the imaginary part, and 31% in the real part, approximately, that mean a loss of robustness and diminished validity in the results. Also, we can observe that the error increases slightly, when m M increases.

The de Sitter family
The de Sitter family of modes are related to the accelerated expansion of the universe, which is related to the surface gravity of the cosmological horizon of pure dS space. They correspond to purely imaginary modes which can be very well approximated by the pure dS QNMs. The pure dS modes are determined by two branches [95,96], and the most important for our purposes corresponds to the lowest lying solution, given by where n is the overtone number. In order to visualize the behavior of the QNFs for the dS modes, in Fig. 7 we plot −I m(ω)M for = 1 and different values of the cosmological constant, note that there is a low decay rate when Mβ increases, and when the cosmological constant decreases. Also, as before, we will study the behavior of the dS QNFs as a function of the scalar field mass. So, we obtain some QNFs for different values of the parameters and n dS = 0, 3 for a fixed values of the black hole mass and    Table 3 for Mβ = 0.01. As we mentioned, this family is dominant for small masses and it can render a real part becoming complex, which depends on the scalar field mass and the angular number. Also, it is possible to observe in Table 3 that the decay rate always increases when the angular number increases, thereby the longest-lived modes are the ones with lower angular number and the anomalous behavior is not present in this family.

The dominance family
In order to observe the dominance between the photon and de Sitter family we plot in Fig. 8 the imaginary part of the QNFs The behavior is similar to that observed in Ref. [68]. We can recognize the two families for zero mass of the scalar field, a family of complex QNFs, and a purely imaginary family. As we have seen, the purely imaginary modes belong to the family of de Sitter modes, while the complex ones correspond to the photon sphere modes.
In order to interpret this figure, first observe the black points, and n dS = 0. This family is dominant for small masses and near m M = 0.17 it can render a real part becoming complex due to the continuity between the black and gray points, and also could be connected with the dS modes with −I m(ω)M ≈ 1. Also, for low values of the mass, it is possible to observe that, for n dS = 1 and n dS = 2, one can have a combination rendering a real part complex for m M ≈ 0.05; the same occurs for the overtone numbers n dS = 3 and n dS = 4, and for m M ≈ 0.025 becoming complex. Besides, for higher overtone numbers, there are dS modes purely imaginary for the whole range of mass considered. On the other hand, for the photon sphere modes, the gray points, we can see that for n PS = 0 and small mass this branch is subdominant; however, for m M ≈ 0.155 it begins to dominate. Now, we consider a bigger value of the parameter, Mβ = 0.20; see Fig. 9. As before, the black points correspond to the purely imaginary QNFs, and the gray points correspond to the complex QNFs. Note that for low values of the scalar field mass dS modes are dominant; however for larger mass the photon sphere modes are dominant. Also, it is possible to observe that dS modes with n dS = 0 could to connect with dS modes with n dS = 1 ( −I m(ω)M ≈ 1.10). Besides, for higher overtones numbers, there is a purely imaginary family of dS modes for the whole range of mass considered. Also, for the photon sphere modes, the gray points, we can see that for n PS = 0 and small mass this family is subdominant; however, for m M ≈ 0.13 it begins to dominate.

Conclusions
In this work, we considered an asymptotically de Sitter black hole solution with a specific f (R) modified gravity as background. In order to see if strong curvature effects, which arise from non-linear terms in R in the action, can modify the behavior of the QNMs, we studied the propagation of a probe scalar field and we analyzed the presence of anomalous behavior in the quasinormal modes spectrum. First, we showed that the QNM spectra consist of two families: the photon sphere modes and the de Sitter modes, and their   For the photon sphere modes we found that, for massless scalar field and small values of the parameter Mβ, the decay rate decreases when the cosmological constant increases; this behavior is similar to the observed for Schwarzschild-dS black holes [68]. On the other hand, for large values of the parameter Mβ the behavior is inverted, i.e. the decay rate increases when the cosmological constant increases. Also, we found that the photon sphere modes show an anomalous behavior that depends on the value of the parameter β and the mass of the scalar field. As in the case of the Schwarzschild-dS black holes, which corresponds to β = 0, there is a critical value of the mass of the scalar field, which depends on the cosmological constant, beyond which the QNMs show an anomalous behavior. However, in f (R) gravity this critical value depends not only on cosmological constant but also on the new parameter β, which expresses the presence of non-linear terms in the curvature. We found that the critical value of the scalar mass decreases when the parameter Mβ increases, and there is a critical value of Mβ where the critical mass is zero indicating that there is no anomalous behavior of the QNMs. Also, we saw that the critical value of the mass of the scalar field increases when the overtone number increases; however, when the parameter Mβ approaches a near extremal black hole, i.e. when the cosmological horizon approximates to the event horizon, the critical value of the mass of the scalar field does not depend of the overtone number.
On the other hand, for the dS modes we observed that for a massless scalar field they are purely imaginary, and for massive scalar field they can render an oscillatory part becoming complex for some value of the mass and higher. Also, we did not find an anomalous behavior for this family of modes.
Also, we found that the imaginary part of the quasinormal frequencies is always negative, leading to a stable propagation of the scalar fields in this background for both families of modes. The PS family or the dS family may dominate the evolution of the perturbation depending on the values of the parameters. We showed that for small values of M 2 the dS modes dominate for small values of the scalar field mass, while the photon sphere modes dominate for larger values of the scalar field mass. Also, the decay rate of the PS modes increases with Mβ, while the decay rate of the dS modes decreases, so it can occur that for low values of Mβ the PS modes dominate and for high values of Mβ the dS modes dominate, depending on the value of M 2 .
An interesting extension of this work is to calculate the QNMs and QNFs of the f (R) gravity theory if we allow the matter fields to backreact with the metric. Then we can study the properties of the resulting hairy black holes [78,79] and possibly see the interplay effects of the geometric scalar present in the f (R) theory and the matter scalar field minimally coupled to the modified gravity theory. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical paper without associated data.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .