Superradiant (In)stability, Greybody Radiation, and Quasinormal Modes of Rotating Black Holes in non-linear Maxwell f (R) Gravity

The research of superradiant instability in the realm of quantum gravity is a well-known topic, with many physicists and astronomers studying the potential impact it can have on gravitational waves, the structure of the universe, and spacetime itself. In this work, we investigate the superradiant (in)stability of a rotating black hole obtained from the nonlinear Maxwell $f(R)$ gravity theory. In this study, the evaluation of stability/instability is going to be based on non-existence and existence of magnetic field, when the magnetic field constant becomes $c_{4}=0$ and $c_{4}\neq 0$, respectively. The analyzes of greybody factor (GF) and quasinormal modes (QNMs) are investigated in the stationary black hole spacetime both in the absence and presence of the magnetic field parameter. To this end, we first consider the Klein-Gordon equation for the complex scalar field in the geometry of that rotating black hole. In the sequel, the obtained radial equation is reduced to a one-dimensional Schr\"{o}dinger-like wave equation with an effective potential energy. The effects of the nonlinear Maxwell $f(R)$ gravity theory parameters ($q$, $c$, and $c_{4}$) on the effective potential, GFs, and QNMs are thoroughly investigated. The obtained results show that even though the factors $q$, $c$, and $c_{4}$ all affect the effective potential, this phenomena, surprisingly, is not valid for the GFs and QNMs. With the proper graphics and tables, all outputs are depicted, tabulated, and interpreted.


I. INTRODUCTION
Superradiance is a term used to describe a radiation amplification system that includes a scattering mechanism. Superradiance plays a noteworthy role in the studies of relativity, astrophysics, quantum mechanics, and optics [1]. This phenomenon can be considered as a quantum aspect of black hole (BH) physics. Namely, in the context of quantum gravity phenomenology, superradiance can play a crucial role in the study of BHs. For example, the emission of Hawking radiation (or greybody radiation) [2,3] from a BH can be enhanced by the presence of surrounding matter, leading to a process known as BH superradiance.
This phenomenon has been proposed as a possible explanation for the observed properties of BH systems, such as the large amounts of energy emitted from the centers of galaxies.
Because it is influenced by the BH's area theorem, tidal forces, the Penrose process, and Hawking radiation [4]. The concept of superradiance was first introduced by Dicke in 1954 [5]. Subsequently, in 1971, this phenomenon was further understood and developed by Zel'dovich [6,7] during his investigation of reflected wave amplification from a rotating BH (Kerr metric). Zel'dovich proved that if the frequency ω of the ingoing wave having the plane wave structure with e −iΩt+imφ (Ω represents the angular velocity, φ is the cyclic coordinate, and m denotes the magnetic quantum number) convinces ω < mΩ, the scattered waves amplify in such a way that the waves coming out of the BH to become more than the ones entering to it. Namely, the dispersed wave is amplified.
Combination of superradiance with a confining procedure to force the wave to consistently interact with the BH to undergo exponential increase known as "BH bomb" [8], which has been an active case since 1970s [9]. This phenomenon can be viable by considering a mirror around the rotating body that could make the system unstable [7,10]. In [11], it was mentioned that such a mirror can be recognized by applying a charged massive scalar field breeding in the Kerr-Newman spacetime. Moreover, BHs can be transformed into efficient particle detectors by applying strong constraints to ultralight bosons via the superradiant instabilities of spinning BHs. On the other hand, instability formation and whether or not its nonlinear time evolution follows the linear intuition are, nevertheless, topics that are not well-understood, yet [12].
f (R) theories of gravity are straightforwardly generated by replacing Ricci scalar R in the Einstein-Hilbert action. Namely, we have a generic action for f (R) as follows where k = 8π denotes the gravitational constant and g is the determinant of metric. Throughout the paper, unless otherwise noted, we shall work in natural units with c s = G = = 1. Since general relativity (GR) has had many unresolved problems, includ-ing the existence of dark energy and dark matter, deflection from Einstein's theory allows us to estimate the fundamental matters and extension of GR (modified gravity). Based on different formalisms, there are three types of f (R) gravity models: 1) metric, 2) Palatini, and 3) metric-affine, f (R) gravity [39][40][41]. The use of f (R) gravity in many contexts is significant; for example, see its astrophysical perspective in [42][43][44][45], the cosmological models with f (R) in [46,47], and the derivations of novel BH solutions in [48][49][50][51]. Moreover, we refer the reader to the monographs [52,53] for some good reviews about f (R). In particular, spherically symmetric BH solutions in f (R) gravity have been receiving special attentions (see for instance [54], in which the exact static spherically symmetric solutions in f (R) gravity coupled with nonlinear electrodynamics derived by Hollestein and Lobo [55]).
Searching for alternative gravitational theories to conventional Einstein's general relativity (GR) is supported by difficult challenges ranging from quantum gravity to dark energy and dark matter. Indeed, there are a lot of unresolved problems with GR, such as singularities, the nature of dark energy and dark matter. All of these problems motivate researchers to improve or modify GR in order to address the challenges at the UV and IR scales [56].
However, the obtained feasible modified/extended theories should be consistent with the present observational/experimental restrictions. With the awareness of this issue, an ambitious study has recently been carried out on the nonlinear Maxwell f (R) gravity [54].
By using dynamical Ricci scalars that asymptotically converge to flat or (A)dS spacetimes, Nashed and Saridakis [57] have derived a new charged rotating BH solutions, which will be our main reference metrics in this study.
Due to the quantum effects, a BH can act as a blackbody object, which emits thermal waves [2,3]. The mass of a BH decreases during its HR, which can lead to complete BH evaporation. As a matter of course, the emitted particles are affected by an effective potential originating from the curvature of the spacetime. As a result, although some waves penetrate the potential barrier and extend to infinity, the remainder is reflected back to the BH. Due to the structure of the effective potential, the radiation spectra are altered and different from that near the event horizon. As a result, the term GF [58] refers to a quantity that measures the deviation of the radiation spectrum from the blackbody radiation. At the event horizon, the BH emission rate [59] is defined as follows by which ω represents the wave frequency, T H and k denote the Hawking temperature and surface gravity, respectively. The relation between emission rate and GF is given by [60] γ where |A l,m | 2 represents the GF. There are various methods to compute the GF, such as matching technique [60][61][62], WKB approximation [63,64], finding Bogoliubov coefficients method [65][66][67][68], the Miller-Good transformation method [58,69], and the rigorous bounds [70]. effective potential. We also study the behaviors of the obtained effective potential under the influence of charge q and magnetic field constant c 4 . As being two subsections, in Sec.
(IV), we examine the superradiant instability for zero and non-zero c 4 values. Sections (V) and (VI) are reserved for the analysis of the greybody radiation and QNMs, respectively.

II. ROTATING BHs IN NONLINEAR MAXWELL f (R) GRAVITY
In this section, we briefly review both new static and rotating BH solutions obtained in nonlinear Maxwell f (R) gravity [57], whose the total action is given by where k is a gravitational constant, which can be considered as k = 1, without loss of generality, in this study. Moreover, L(F) indicates a general gauge-invariant electromagnetic Lagrangian where the usual antisymmetric Faraday tensor is F = 1 4 F αβ F αβ .
√ −g represents the determinant of the metric g µν . The corresponding gravitational field equations of the action (4) can be derived as follows by which nlem µυ denotes the energy momentum tensor and f ≡ df ( ) d . Using the following ansatz for a spherically symmetric line element: in Eq. (5), after making some tedious calculations, the following metric function was finally obtained by Nashed and Saridakis [57] where c is a positive constant, which can take limited values: 0 < c ≤ 2, q and M stand for the charge and mass, respectively (see Ref. [57] for the details). In Fig. 1, we show the behavior of the metric function H (r) under the influence of varying parameters q and c. It is clearly seen that the spacetimes exhibit flatness at the asymptotic distances, independent of the values of q and c. In that case, the horizon radius is obtained as, Hence, it will be Schwarzschild black hole if we set c = 2 and q = 0. Also, extremal limit is given by the following condition, to the static and spherically symmetric metric (6). In Eq. (10), a is the rotating parameter and Ξ = √ 1 + a 2 . Thus, we have (11) in which H(r) is provided by the static solution (7) previously derived. On the other hand, the general gauge potential is defined by [57] where q(r), s(r), and n(φ) are 3 free functions generating the electric and magnetic charges in the vector potential as follows in which c 4 represents the magnetic field constant. In the following sections, our investigation will consider cases in which the magnetic field constant does not exist (c 4 = 0) and exists (c 4 = 0).

III. SCALAR PERTURBATION
In recent decades, perturbations of BHs and stars have arisen as one of the main topics of relativistic astrophysics. Furthermore, perturbations are hot subjects right now because of their functions in gravitational waves. In this part, we employ the charged KGE to arrive at the Schrödinger wave equation in one dimension. The effective potential to be obtained in this section is crucial for studying superradiance, greybody radiation, and QNMs.
Let us consider the charged and massive KGE: where Q and m are the charge and mass of the scalar field (spin-0), respectively. Moreover, √ −g represents the determinant of the metric. Here, for metric (11), we consider the following ansatz for the spinor field: where ω represents the frequency of the wave and k is azimuthal quantum number. During the derivation of the scalar wave equation, we will consider the dyonic case. Plus, we set s(r) = n(φ) = 0 in Eq. (12). Thus, the components of the vector potential read Throughout the paper for more convenience in our calculation, without loss of generality, we shall consider qQ → q 2 . After substituting Eq. (16) and ansatz (15) into the massive charged KGE (14), one can obtain where λ is the eigenvalue whose value can be found with the help of the angular part: In the mean time, throughout the paper, a prime (dash) symbol is used to denote the derivative of a function with respect to its argument. By considering the definition of the tortoise coordinate: and in the sequel applying the transformation R(r) = U (r) r to Eq. (18), one can acquire one-dimensional Schrödinger like wave equation as follows: in which 2 = ω 2 (1 + a 2 ) = ω 2 Ξ 2 and the effective potential is given by The behaviors of the effective potential (21) when the charge parameter q is changed for various values of c are depicted in Fig. 2.

IV. SUPERRADIANCE PHENOMENON
A. For c 4 = 0 Here, we investigate the stability of the rotating BH obtained from the non-linear Maxwell f (R) gravity. To this end, we consider the method prescribed in Ref. [87,88]. After applying the transformation U (r) = e −i r * ψ (r) to the Schrödinger equation seen in Eq. (20), one Now, let us replace the tortoise coordinate (19) with the naive radial coordinate r and multiply Eq. (23) by ψ * . By imposing H (r + ) = 0 and ψ (∞) = 0, one can solve the final differential equation by performing the well-known integration by parts method. Thus, where the second term in the integrated can be expanded to Therefore, the integration (24) recasts in It is also possible to write Eq. (26) as In Eq. (27), where the potential is approximated to and correspondingly Thus, in Region I, the solution of the radial equation (20) is obtained as which slightly away from the event horizon yields [88] where the near horizon tortoise coordinate is defined as r * = dr H(r) ≈ 1 H (r h ) ln (r − r h ). Between the event horizon and the distant regions, Region II acts as an intermediate zone and is defined as: Therefore, the radial equation (20) reads By comparing the solutions in the first and second regions, we define the constants as as A = B, and C = −Ai (ωΞ + (q 2 + ak)).
To find the solution in Region II, we take into consideration an asymptotic series for tortoise coordinate i.e., r r h . Hence, we get The third region (Region III) is the asymptotic zone (r r h ,), where the conducting terms of the effective potential become Thus, one can get the Region III solution as Then, after matching the solution of Region II with the solution of Region III, we get where η = m 2 c 4(ωΣ + (q 2 + ak)) .
To obtain the reflection coefficient and the GFs and to complete our assessment of superradiant stability/instability, we employ the flux expression as follows Therefore, one can obtain the near-horizon flux as Similarly, the asymptotic flux at spatial infinity becomes where ξ = (ωΞ + (q 2 + ak)) 2 + m 2 c 2 . By considering D 1 =D 1 +D 2 and D 2 = i(D 2 −D 1 ), the asymptotic incoming and outgoing fluxes (42) can be written as follows and By employing the definition of the reflection coefficient and GF, we get and Now, based on Eqs. (45) and (46), we are able to determine the superradiant condition.
To this end, either the reflection coefficient should be greater than 1 or the GF should be negative: By taking into account superradiant instability conditions, in the relevant sub-cases, we now review the behavior of the effective potential in Fig. 2. It is obvious from Fig. 2 whereω = H(ω + qc 4 φ) and the effective potential is determined as a complex expression.
In Region 3, its real part can be expressed as and its imaginary part reads In Eqs. (49) and (50), the terms including c 2 4 and c 4 r 2 are ignored. In addition, our analysis has shown us that it is reasonable to consider only the real component of the effective potential.
The wave solutions of regions 1 and 2 for the existing magnetic fields of the BH are the same as the non-existing ones, Eqs. (30) and (34), but the wave solution of Region 3 (r >> r h ) with c 4 = 0 is different than the c 4 = 0 one: where in which To determine the flux expressions at the horizon and spatial infinity, we apply the same method followed in the previous sub-section. Thus, we have and where By substituting D 1 =D 1 +D 2 and D 2 = i(D 2 −D 1 ) in Eq. (56), one can obtain Therefore, the reflection coefficient of the rotating BH with small c 4 reads and the corresponding GF becomes γ = F hor F ∞in = (ωΞ + (q 2 + ak))|A| 2 Since the result is almost the same as Eq. (46), Moreover, no explicit well in the effective potential behavior in Fig.3, the interpretation for the stability in the presence of a magnetic field will be the same as the nonexistent one.

V. SEMI-ANALYTICAL GREYBODY RADIATION
In this section, we shall follow the method, which was reviewed in [58] (and references therein) to analyze the greybody radiation for both cases of c s = 0 and c 4 = 0.
The general semi-analytic bounds for GFs are given by [89] where σ l represents the GF and ℘ is formulated as follows by which h is a positive function that satisfies the following condition: h (−∞) = h (+∞) = ω. Normally, one follows the method of replacing the V parameter with the potential obtained in Eq. (36) and then employs the tortoise coordinate to evaluate the GF (62). On the other hand, that method is not always the best course of action to take.
In fact, this method also fails in our situation. So, to overcome this discrepancy, we set . This allows us to rewrite the expression for GF (62) as which corresponds to One can immediately observe that Eq. (65) is valid for ω 2 > V peak , the peak value of the potential [58]. Besides, Eq. (65) can be rewritten as To find the maximum or the peak of the potential, first we derive the r peak from the effective potential equations by taking derivative with respect to r , as is well-known, one should find where the graph shifts from increasing to decreasing. To find out the rate at which the graph shifts from increasing to decreasing, we look at the second derivative and see when the value changes from positive to negative. Depending on the values of V peak , the GFs are obtained. For instance, by setting M = m = 1 and λ = 2, in relation to the variables q and c, the V peak expression is given by infinitesimal c 4 case, V peak is found to be

VI. QNMs
QNMs are important in the study of BH perturbation because they provide a way to understand the behavior of a BH in the presence of external perturbations such as a scalar, electromagnetic, gravitational, etc. When a BH is perturbed, it will respond by emitting gravitational waves that have a characteristic frequency. This frequency, known as the QNM, is determined by the properties of the BH, such as its mass, charge, and spin. QNMs of a BH is characterized by complex numbers. The complex frequency of a QNM is given by a real part and an imaginary part. The real part represents the oscillatory frequency of the mode, while the imaginary part represents the rate of decay or growth of the mode. By observing the QNMs in gravitational waves, astronomers can learn about the physical characteristics of the object that produced the waves.
In this section, for the scalar perturbations, we consider a semi-analytical approach to derive the frequencies of the QNMs of the charged rotating BHs in nonlinear Maxwell f (R) gravity. To this end, we employ the WKB (Wentzel-Kramers-Brillouin) approximation, which is a mathematical method used to solve differential equations (DEs) with a largescale parameter. This approximation allows for a simplified solution to the differential equation and provides an estimate of the energy levels (frequency) within a certain accuracy.
Conventionally, the WKB approach is based on the assumption that the solutions can be expressed as an exponential power series, where the coefficients of the series are determined by solving a set of recursive equations: The DE has a general form as where Q (r) = V (r) − ω 2 , which contains two turning points. The boundary condition of waves is chosen to be Ψ (r) = Z out Ψ (r) out for outgoing waves as r * → +∞ and Ψ (r) = Z in Ψ (r) in for incoming waves while r * → −∞. First, Mashhoon [90] invented this approach and applied it to the BHs in 1983. Then it was developed by [91,92]. The WKB approximation can be extended from the third to sixth order. The sixth-order WKB approximation, also known as the Konoplya approximation, is a method used to approximate the solutions of differential equations with complex potentials. The sixth-order WKB approximation by Konoplya can be found in his seminal papers [84,93]. The Konoplya approximation uses a series expansion of the solution to the DE in powers of the small parameter , which is the wavelength of the solution. Konoplya approximation for obtaining the complex frequencies of the QNMs is given by the following expression [84]: where and In Eqs. (71)- (73), the primes and superscripts (4, 5, 6; for the higher order derivatives) denote the differentiation with respect to r * and α = n + 1 2 , where n denotes the tone number. By considering the effective potentials of both solutions, the results are tabulated in Tables (I) and (II) for the zero and non-zero magnetic field constants, respectively.
The behaviors of the QNMs for c 4 = 0 are illustrated in Figs. 6. One can observe that both parts (real and imaginary) of the QNMs decrease by increasing the charge parameter, q.
Moreover, for n = 0, QNMs increase by growing the c parameter and then start to decrease.
In addition, when n = 1 and q increases, the oscillation frequencies rise and the damping mode steadily declines.   Tables I and II under different physical parameter changes and shown that the results support Figs. 7 and 8.
In our findings, we have discovered that all unstable modes exhibit superradiance and all stable modes do not, consistent with the superradiant condition (47). This means that scalar waves can experience superradiant amplification by extracting charge from the BH, indicating that the BH geometry is unstable. Additionally, superradiance can also be used to probe the fundamental principles of quantum gravity, such as the behavior of quantum