Probing Palatini-type gravity theories through gravitational wave detections via quasinormal modes

The possibility of testing gravity theories with the help of gravitational wave detections has become an interesting arena of recent research. In this paper, we follow this direction by investigating the quasinormal modes (QNMs) of the axial perturbations for charged black holes in the Palatini-type theories of gravity, specifically (1) the Palatini f(R) gravity coupled with Born-Infeld nonlinear electrodynamics and (2) the Eddington-inspired-Born-Infeld gravity (EiBI) coupled with Maxwell electromagnetic fields. The coupled master equations describing perturbations of charged black holes in these theories are obtained with the tetrad formalism. By using the Wentzel-Kramers-Brillouin (WKB) method up to 6th order, we calculate the QNM frequencies of the EiBI charged black holes, the Einstein-Born-Infeld black holes, and the Born-Infeld charged black holes within the Palatini R+αR2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R+\alpha R^2$$\end{document} gravity. The QNM spectra of these black holes would deviate from those of the Reissner-Nordström black hole. In addition, we study the QNMs in the eikonal limit and find that for the axial perturbations of the EiBI charged black holes, the link between the eikonal QNMs and the unstable null circular orbit around the black hole is violated.


Introduction
It is quite fair to say that one of the most enchanting events of recent discovery in modern physics is the direct detection of gravitational waves (GWs) from the coalescence of binary black holes [1,2]. The reason why the direct detections of GW signals are so appealing is that they not only confirm the predictions of Einstein's general relativity (GR) once a e-mail: b97202056@gmail.com b e-mail: mariam.bouhmadi@ehu.eus c e-mail: pisinchen@phys.ntu.edu.tw again, but also render GWs spectacularly a suitable tool for human being to hear deep into the sky far beyond the reach of electromagnetic signals. Not long after the first detection of GWs, the LIGO-VIRGO collaboration succeeded in detecting the GWs emitted from the merger of binary neutron stars, with an accurate localization of the source [3]. The prompt and accompanied electromagnetic signals emitted from the source were also detected. This outstanding accomplishment has initiated a new era of multi-messenger astronomy.
In addition, the direct detections of GWs could help us to test other gravitational theories, or even to falsify some extended theories of gravity [4]. In fact, one of the reasons to consider extended theories of gravity is that GR inevitability predicts the existence of spacetime singularities like the big bang singularity and the black hole singularity. To ameliorate these spacetime singularities, one may resort to some extended theories of gravity which modify Einstein equation at the large curvature limits, but reduce to GR when curvature becomes small. Within the new era of GW astronomy, one plausible way to test these extended theories of gravity, for instance, is via the speed of GWs, as was done in Refs. [5][6][7][8][9].
Another interesting aspect regarding GW detections could be the ringdown signals in the final stage of a merger event. Essentially, the final product of a merger event, no matter if seeded from binary black holes or from binary neutron stars, is usually a black hole. Before the final black hole settles itself, there is an intermediate stage where the distortion of the black hole is gradually relieved, with the emission of GWs. In practice, the ringdown stage can be described by the theory of black hole perturbations and the frequencies of the GWs are featured by quasinormal modes (QNMs). In this stage, the distorted black hole can be regarded as a dissipative system. The perturbations have a discrete spectrum and the QNM frequencies are complex numbers. The real part of the frequencies describes the oscillations of the perturbations and the imaginary part corresponds to the decay of the amplitude. Interestingly, these QNM frequencies merely depend on the parameters characterizing the black holes, such as the mass, the charge, and the spin. If there are additional parameters appearing in the underlying theory, these parameters should manifest themselves in the QNM spectra. Along this direction of research, the QNMs of black holes in several gravity theories have been investigated, such as in the Horndeski gravity [10][11][12][13][14], metric f (R) gravity [15][16][17], massive gravity [18,19], Einstein-dilaton-Gauss-Bonnet gravity [20][21][22][23], the Randall-Sundrum braneworld model [24], Hořava-Lifshitz gravity [25], higher dimensional black holes [26][27][28], and Einstein-aether theory [29], etc. Furthermore, the QNMs of some regular black holes [30,31] and the black holes with non-commutative geometry [32,33] have been analyzed in the literature. In addition, probing signatures of the black hole phase transitions in modified theories of gravity via QNMs has been shown to be possible [34]. See Refs. [35][36][37][38] for nice reviews on the latest progress of the field.
In this paper, as a further extension of our previous work [39] in which the QNMs of massless scalar field perturbations were studied, we consider the QNMs of the axial perturbations for the charged black holes in two Palatini-type gravity theories: (1) the Palatini f (R) gravity coupled with Born-Infeld nonlinear electrodynamics (NED) and (2) the Eddington-inspired-Born-Infeld (EiBI) gravity coupled with linear electromagnetic fields. To calculate the QNM frequencies, we use the WKB method up to 6th order [40][41][42][43]. We also calculate the QNMs in the eikonal limit in which the multipole number l → ∞. Furthermore, the QNM frequencies will be compared with those of the Reissner-Nordström (RN) black hole in GR. Note as well that for the merger events of binary neutron stars, the ringdown timescale is usually shorter than the timescale of charge neutralization of the black hole [44]. This justifies to some extent the validity of studying QNMs of charged black holes from the astrophysical point of view.
The charged black holes in the Palatini f (R) gravity coupled with Born-Infeld NED have been studied in detail in Ref. [45]. The black hole solutions are very close to the RN black hole at the exterior spacetime, while deviate from it inside the event horizon because of the Born-Infeld NED source and the nonlinear function f (R). It has been shown that there exist some regions in the parameter space where the singularity inside the event horizon is replaced with a finite size wormhole structure [46]. Moreover, one can construct the Einstein-Born-Infeld (EBI) black hole by choosing f (R) = R. The properties of this charged black hole have been widely studied in the literature [47][48][49][50]. Again, due to the Born-Infeld corrections from the NED source, the inte-rior structure of the black hole would change significantly as compared with that of the RN black hole in GR.
The EiBI gravity was formally proposed in Ref. [51] to resolve the initial big bang singularity [52]. This theory reduces to GR in vacuum but differs from it when matter is included. The exact expressions and some interesting properties of the charged black holes in the EiBI theory were studied in Refs. [53,54] (see Ref. [55] for a review on the EiBI gravity). Due to the Born-Infeld corrections from the gravity sector, the interior structure of the black hole could deviate from that of the RN black hole notably [56][57][58]. One can then compare the QNM frequencies of the EBI black holes and those of the EiBI charged black holes to see how the Born-Infeld structure from the matter and the gravitational sector affects the QNMs.
This paper is outlined as follows. In Sect. 2, we briefly review the tetrad formalism which will be used later to derive the master equations describing the axial perturbations of the black holes. In Sect. 3, the perturbed Maxwell equation for NED is obtained for the sake of later convenience. In Sect. 4, we derive sequentially the coupled master equations of the axial perturbations for charged black holes in the Palatini f (R) gravity coupled with Born-Infeld NED and in the EiBI gravity with linear electromagnetic fields. In Sect. 5, we calculate the QNM frequencies of these charged black holes by using the WKB semi-analytic method. The QNMs in the eikonal limit are discussed as well. We finally conclude in Sect. 6.

Tetrad formalism
To study the QNMs of the black holes of our interest, we shall consider the perturbations of a static and spherically symmetric spacetime. Without loss of generality, the perturbed spacetime can be described by a non-stationary and axisymmetric metric in which the symmetrical axis is turned in such a way that no φ dependence appears in the metric functions. In general, the metric can be written as follows [59]: where ν, ψ, μ 2 , μ 3 , σ , q 2 , and q 3 are functions of time t (t = x 0 ), radial coordinate r (r = x 2 ), and polar angle θ (θ = x 3 ). Because the system is axisymmetric, the metric functions are independent of the azimuthal angle φ (φ = x 1 ). In this work, the notation used in Ref. [59] is strictly followed. The only difference is that the metric function ω used in Ref. [59] is replaced with σ in Eq. (2.1) since we will use ω to denote the frequency of the perturbations later. Note that in the background spacetime which is static and spherically symmetric, we have σ = q 2 = q 3 = 0.
To study the perturbations of the spacetime metric (2.1), we will use the tetrad formalism in which one defines a basis associated with the metric (2.1) [59]: and where the tetrad indices are enclosed in parentheses to distinguish them from the tensor indices. The tetrad basis should satisfy Conceptually, in the tetrad formalism we project the relevant quantities defined on the coordinate basis of g μν onto a chosen basis of η (a)(b) by constructing the tetrad basis correspondingly. In practice, η (a)(b) is usually assumed to be the Minkowskian matrix In this regard, any vector or tensor field can be projected onto the tetrad frame in which the field can be expressed through its tetrad components: It has been shown in Ref. [59] that the master equations describing the gravitational perturbations of black holes (Schwarzschild, RN, etc) can be obtained by using the tetrad formalism in a straightforward and concise manner. One should notice that in the tetrad formalism, the covariant (partial) derivative in the original coordinate frame is replaced with the intrinsic (directional) derivative in the tetrad frame.
For instance, the derivatives of an arbitrary rank two object H μν in the two frames can be related as follows [59] where a vertical rule and a comma denote the intrinsic and directional derivative with respect to the tetrad indices, respectively. A semicolon denotes a covariant derivative with respect to the tensor indices. Furthermore, the Ricci rotation coefficients are defined by and their components corresponding to the metric (2.1) are given in Ref. [59].

Perturbed Maxwell equation for NED
In a gravity theory formulated upon the Palatini variational principle, the matter Lagrangian is assumed to be coupled with the physical metric g μν only. Therefore, the matter fields would follow the geodesics defined by this metric and the conservation equation of the energy momentum tensor follows the standard form with respect to g μν . In this section, we focus on the matter Lagrangian described by NED [45]: where we have set G = c = 1, and φ(X, Y ) is a function of gauge field invariants defined by [45]: where F * μν ≡ 1 2 μναβ F αβ is the dual of the field strength. The standard Maxwell electromagnetic fields are recovered when φ(X, Y ) = X ; For the sake of simplicity, we will assume a vanishing magnetic field, i.e., Y = 0, in the rest of this paper.
For a gravitational theory minimally coupled to NED with a purely radial electric field and no magnetic field, only the (t, r ) and (r, t) components, i.e. the (0, 2) and (2, 0) components of the field strength F μν appear at the background level. In the tetrad frame, the field strength F (a)(b) at the background level satisfies [45] F 02 = F (0) (2) where φ X = dφ/d X and Q * is an integration constant which can be regarded as the charge of the black hole. Note that the last equality in Eq. (3.3) can be obtained from the conservation equation of NED at the background level [45]. From Eq. (3.3), we get In the general case where the perturbations are taken into account, the metric functions and the field strength could depend on t, r , and θ . In this case, the Bianchi identity of the field strength F [(a)(b)|(c)] = 0 leads to Note that the comma here denotes the partial derivative with respect to the tensor indices. This derivative is related to the directional derivative shown in Eq. [59]. In addition, Eq. (3.5) is a redundant equation because it is just an integrability condition for Eqs. (3.6) and (3.7).
On the other hand, the conservation equation for NED η (b)(c) (F (a)(b) φ X ) |(c) = 0 can be written explicitly as follows Again, it can be shown that Eq. (3.9) is a redundant equation since it is an integrability condition for Eqs. (3.10) and (3.11).
To linearize the equations above, the scalar field φ X and the gauge field invariant X should be decomposed into the background and the 1st order parts: where δφ X = φ X X δ X = 2φ X X F (0) (2) δ F (0) (2) . (3.14) In the above expressions, we have decomposed F (0) (2) as (2) in which the background F (0) (2) is given in 1 Eq. (3.4). Note that the quantity X at the background level is given by X = F 2 (0) (2) . After the decomposition, the linearized Eqs. (3.6), (3.7), and (3.8) are On the other hand, the linearized Eqs. (3.10), (3.11) and (3.12) can be written as where δν, δμ 2 , δμ 3 , and δψ are the 1st order parts of the metric functions ν, μ 2 , μ 3 , and ψ, respectively. For the sake of abbreviation, we define the field perturbation Recall that in the derivation of Eq. At the end of this section, we would like to write down the perturbed energy momentum tensor of NED, which will appear on the right hand side of the perturbed gravitational equations later. The energy momentum tensor of NED in the tetrad frame is Furthermore, the perturbed energy momentum tensor reads Finally, we can write down its components explicitly as follows: (3.28)

The master equations
As we have mentioned previously, the master equations describing the gravitational perturbations of a charged black hole are two coupled equations. This is because of the coupling between the gravitational field and the electromagnetic field in the system. So far we have derived one of the coupled equations, i.e., Eq. (3.22), from the Maxwell equation of the NED source. In this section, we will carry out the derivation of the other coupled equation from the gravitational equation of the theory. We will first consider the Palatini f (R) gravity coupled with NED and obtain the master equations in this theory. After that we will turn to deduce the master equations of the EiBI gravity coupled with linear electromagnetic fields.

Axial perturbations of charged black holes in Palatini f (R) with NED
In this subsection, we will consider the Palatini f (R) theory coupled with NED. The action of the theory reads [45] where the matter Lagrangian S m is given in Eq. (3.1). We would like to emphasize again that the theory is formulated within the Palatini variational principle in which the metric g μν and the affine connection are independent variables. For a nonlinear function f of the Ricci scalar R ≡ g μν R μν ( ), the equations of motion would be different from those in the metric f (R) theory.
In addition, after deriving the master equations we will consider a particular NED model, that is, the Born-Infeld NED: The background solutions of the charged black holes in the Palatini f (R) gravity coupled with the Born-Infeld NED Lagrangian (4.2) have been studied in Ref. [45]. The QNMs of a massless scalar field of such black holes have been discussed in Ref. [39]. The most general form of the metric functions of these black holes have been derived in Ref. [45] as well (see also Eqs. (3.24) and (3.25) in Ref. [39]). Specifically, if we focus on Einstein gravity ( f (R) = R) coupled with the Born-Infeld NED, we would get the Einstein-Born-Infeld (EBI) black hole whose deviations from the RN black hole result purely from the matter sector. Its exact metric functions at the background level read [39,[47][48][49][50] where F(.., ..; ..; ..) is the hypergeometric function [60] and r m is defined as r m ≡ √ Q * /β m . Note that we have used the following dimensionless rescalings: where r s /2 ≡ M 0 denotes the mass of the black hole seen by an observer infinitely far away. For the sake of convenience, we will use these rescalings in the rest of this paper.
As mentioned in Ref. [39], it is interesting to compare the QNMs of the EBI charged black hole with those of the charged black holes within the EiBI gravity coupled with Maxwell electromagnetic fields. One can then compare directly the QNMs of charged black holes within two theories, one with the Born-Infeld correction from the matter sector and the other one with this kind of modification from the gravitational sector.

Perturbed field equations
In a gravitational theory constructed on the Palatini variational principle, the variation of the gravitational action with respect to the affine connection determines the auxiliary metric q μν which is compatible with the affine connection: In the Palatini f (R) gravity, the auxiliary metric and the physical metric are conformally related Therefore, their metric functions are related as follows On the other hand, the variation of the action with respect to g μν reads It should be emphasized that the symmetric part of the Ricci tensor R (μν) in Eq. (4.7) is defined solely by the affine connection . Because the field equation ensures the compatibility between the auxiliary metric q μν and the affine connection, it is fair to say that the Ricci tensor is defined by the auxiliary metric, that is, In order to recast Eq. (4.7) into the tetrad frame, we rewrite Eq. (4.7) as follows μ is a tetrad basis mapping the auxiliary metric q μν onto the tetrad frame. More explicitly, it satisfies This additional tetrad basis is related to the physical tetrad basis according to the conformal relation between the two metrics: In this regard, Eq. (4.7) can be rewritten as where the expressions ofR (a)(b) are given in Ref. [59], in which all quantities should be replaced with their tilde counterparts. Furthermore, the scalar curvature R can be written as Since we are focusing on the axial perturbations (odd parity perturbations) which change sign when φ → −φ, we only consider the (1, 3), (1, 2) and (0, 1) components of the linearized Eq. (4.11): Then, we define with which Eqs. (4.13) and (4.14) can be rewritten as  3 . (4.20) Equations (4.19) and (4.20) form the coupled master equations describing the evolutions of the perturbation fields B and Q, which correspond to the perturbations of the matter field and the gravitational field, respectively.

Effective potentials
For the sake of later convenience, we would like to recast the coupled master equations (4.19) and (4.20) into a Schrödinger-like form, which in practice is more suitable for the calculations of QNMs with the WKB method. We consider the ansatz [59] Q(r, θ) = Q(r )Y (θ ) , B(r, θ) = B(r )Y ,θ / sin θ, (4.21) where Y (θ ) is the Gegenbauer function satisfying [60] d dθ where μ 2 = (l − 1)(l + 2) and l is the multipole number. From Eq. (4.22), it can be proven that where we have used Eq. (3.4) for the background field F (0) (2) and the Fourier decomposition ∂ t → −iω.
We introduce the following definitions where Z ≡ r f 1 .

(4.29)
It can be seen that the coupled master equations have been recast into a Schrödinger-like form and they can be written in a matrix expression as follows where V i j is given in Eqs. (4.28) and (4.29). According to the coupled master equations (4.28) and (4.29), one can see that: 1. When f R = 1 and the NED model is assumed to be the Born-Infeld NED given by Eq. (4.2), the master equations reduce to those of the EBI black hole given in Ref. [ i , the master equations reduce to those of the RN-dS(AdS) spacetime given in Refs. [62][63][64][65][66]. 3. If φ = X , and f R = 1, we have The master equations turn out to be those of the RN black hole [59]. (iv) If φ = X , Q * = 0, and f R = 1, we have V 12 = V 21 = 0, and Therefore, the potential for pure electromagnetic perturbations and for pure axial gravitational perturbations of the Schwarzschild black hole (the Regge-Wheeler equation [67]) are recovered, respectively.
As we show in the appendix A, the coupled master equation (4.30) can be decoupled within a WKB approximation in the cases that we are studying in this work. Therefore, to proceed we will decouple the master equation by diagonalizing the matrix V i j to obtain its eigenvalues V 1 and V 2 : If φ = X , Q * = 0, and f R = 1, it can be seen that V 1 and V 2 reduce to Eqs. (4.34) and (4.35), respectively. From the discussion above, we have proven that the potential terms in the master equation (4.30) reduce to their RN counterpart in the correct limits. In the presence of nonlinearity of NED and the gravitational function f (R), the potentials would change significantly and, consequently, alter the QNM frequencies.
We will discuss this issue later in Sect. 5.

Axial perturbations for EiBI charged black holes
In this subsection, we will consider the EiBI theory coupled with linear Maxwell fields. The total action is given by [51] where the matter Lagrangian is described by the linear Maxwell fields: φ X = 1. In the above action, = ±1 indicates that one can freely choose the Born-Infeld coupling constant to be either positive or negative. The dimensionless constant λ is related to an effective cosmological constant via = β 2 g (λ − 1). In the rest of this paper, we will assume a zero effective cosmological constant (λ = 1) and focus on black hole solutions which are asymptotically flat. It should be stressed that only the symmetric part of the Ricci tensor R (μν) ( ) is considered to respect the projective symmetry of the theory.
Since only the linear Maxwell fields are considered, the perturbed Maxwell equation can be obtained by simply rewriting the NED equation On the other hand, one needs to take the perturbed gravitational equation into account to complete the derivation of the coupled master equations of the EiBI theory. The perturbed gravitational equation contains the perturbed energy momentum tensor, which can be rewritten from that of the NED with φ = X as: δT (2)

Perturbed field equations
Because the EiBI gravity is also formulated within the Palatini variational principle, there is an auxiliary metric q μν which is compatible with the affine connection as before. The line element of q μν can be similarly expressed as in Eq. (4.5).
According to the variation of the action, the auxiliary metric satisfies the following Eq. [51]: At the beginning, the symmetric part of the Ricci tensor is assumed to be constructed from the affine connection. However, we can recast it as a function of the auxiliary metric, that is, R (μν) (q), because the compatibility between the affine connection and the auxiliary metric.
In the tetrad frame, we can rewrite Eq. (4.43) as Using a similar procedure to the one we have carried previously in the Palatini f (R) gravity, we have constructed another tetrad basisẽ In the tetrad frame, this equation reads where e = √ −g andẽ = √ −q. Even though the two tetrad bases are not related explicitly as in the Palatini f (R) gravity, that is Eq. (4.10), in the EiBI theory these two bases are still related implicitly via Eqs. (4.45) and (4.47). We will immediately show how to derive the master equations of the axial perturbations by using these two equations. The We then have where we have defined We then have (2) , (4.58) which leads to e 3ψ−ν−μ 2 +μ 3 σ ,2 −q 2,0 ,2 + e 3ψ−ν−μ 3 +μ 2 σ ,3 −q 3,0 ,3 = 4e μ 2 r Q * sin 2 θ F (1) (2) . (4.59) To simplify the equations, we definẽ  3 . (4.63) We have obtained one of the coupled master equations (4.63). Now we need to consider the other equation which comes from the perturbed Maxwell equation, that is, Eq. (4.38). It is necessary to replace the right hand side of Eq. (4.38) with its tilde counterpart because of the definition (4.60). If we express Eqs. (4.52) and (4.57) more explicitly by writing down the metric functions, we obtain (4.67)

Effective potentials
Similar to what we have done in the previous subsection, we substitute (4.21) into the master equations. Consequently, Eqs. (4.63) and (4.67) can be rewritten as follows With the further definitions where W ≡ r √ σ + and S ≡ √ σ − /σ + , and after introducing the tortoise radius we can obtain These coupled equations can be written in a matrix form similar to the one given in Eq. (4.30). Furthermore, since the eigenvectors of the matrix V i j are approximately constant as those for the EBI black hole shown in the appendix A, we can diagonalize the matrix as we did in Eq. (4.36). Because the expressions are so complicated, we did not write down the explicit form of V i in this paper. It is still important to check whether the potentials reduce to those of the RN black hole in the proper limits. Specifically, if β 2 g → ∞ we have σ + ∼ σ − ∼ 1 and the master equations reduce to those of the RN black hole. Moreover, if Q * = 0 the master equations of the Schwarzschild black hole are recovered.
In addition, we should highlight a crucial result following from the master equations (4.72) and (4.73). It can be seen that in the second term on the right hand side of Eq. (4.72), i.e., the term containing μ 2 , there is a factor σ + /σ − . On the other hand, in the second term on the right hand side of Eq. (4.73), there is a factor σ − /σ + . Actually, these factors play a crucial role when the QNMs in the eikonal limit (l → ∞) are considered. In that limit, we will show later in Sect. 5 that, because of these factors, the QNMs cannot be calculated directly from the associated quantities of the unstable photon sphere of the black hole and the correspondence proposed in Ref. [68] is not satisfied for the EiBI charged black holes (for more fundamental illustration on the photon sphere, see Ref. [69]).
The derivation of the exact metric functions given in Eqs. (4.74) and (4.75) were first obtained in Refs. [53,54]. One can also refer to Ref. [39] in which we recast the metric functions in a simpler form for calculating the QNM frequencies. Note that in the EiBI gravity, there are some regions of the parameter space where no black hole solution exists [53,54]. In this paper, we will only focus on the cases where the black holes exist and calculate their QNM frequencies.

QNM frequencies: the 6th order WKB method
The evaluation of the QNM frequencies is essentially based on treating the master equations of the perturbations as an eigenvalue problem with proper boundary conditions. In the literature, there have been several methods to calculate the QNMs, ranging from numerical approaches [70,71] to semianalytic methods (see Refs. [35][36][37][38] and references therein).
In this paper, we will use a semi-analytical approach firstly formulated in the seminal paper [40]. This approach is based on the WKB approximation and the QNMs can be calculated by just using a simple formula once the potential terms in the master equations are given. In Refs. [41,42], the 1st order WKB method was extended to the 3rd and 6th orders WKB approximation, respectively. Recently, a further extension of the WKB method up to the 13th order has been proposed with the help of the Padé transforms [43]. The WKB method is expected to be accurate as long as the multipole number l is larger than the overtone n [36]. On the other hand, for astrophysical black holes, the fundamental mode n = 0 has the longest decay time and therefore dominates the late time signal of the ringdown stage. At this regard, we will mainly focus on the fundamental mode. The formulation of the WKB method to calculate the QNMs is essentially based on the fact that the master equations can be written like a Schrödinger wave equation in quantum mechanics. The potential term, in most cases (including ours), has a finite value when r * → ∞ (spatial infinity) and r * → −∞ (at the event horizon). Furthermore, the potential has a peak at some finite r * . One can then treat the problem as a quantum scattering process through a potential barrier after suitable boundary conditions for the problems are imposed. At spatial infinity, only outgoing waves moving away from the black hole exist. On the other hand, there can only exist ingoing waves moving toward the black hole at the event horizon.
The idea of the WKB method to encompass the aforementioned boundary conditions is to consider a quantum scattering process without incident waves, while the reflected and the transmitted waves have comparable amounts of amplitudes. The peak value of the effective potential V eff (r * ) ≡ −ω 2 +V is required to be slightly larger than zero in the sense that there are two classical turning points near the peak. The solutions far away from the turning points (r * → ±∞) are solved by using the WKB approximation up to a desired order and the boundary conditions should be taken into account. At the vicinity of the peak, the potential is expanded into a Taylor series up to a given order, and one uses a series expansion method to solve the differential equation. Finally, the numerical values of the QNM frequencies ω can be obtained by matching the solution near the peak with the solutions deduced from the WKB approximation simultaneously at the two turning points.
In the 6th order WKB method, the WKB formula for calculating QNMs is [40][41][42] i where the index m denotes the quantities evaluated at the peak of the potential. V m is the second order derivative of the potential with respect to r * , calculated at the peak. i are constant coefficients resulting from higher order WKB corrections. These coefficients contain the value and derivatives (up to the 12th order) of the potential at the peak. 2

Fundamental QNMs
In Fig. 1, we consider the EBI black hole and show its QNM frequencies calculated from V 1 and V 2 , respectively [see the expressions of these potentials in Eq. (4.36)]. We consider the multiple l = 2 and the fundamental mode n = 0. For the sake of convenience to highlight the deviations due to the Born-Infeld corrections, we present the QNMs ratio of the EBI black hole and the RN black hole. In Fig. 2, we consider the EiBI charged black holes and exhibit their QNM frequencies calculated from V 1 and V 2 , respectively. The solid curves and the dashed curves cor- 2 The explicit expressions of i are given in Refs. [41,42] (see Eqs. (1.5a) and (1.5b) in Ref. [41], and the appendix in Ref. [42]). respond to the results when the EiBI coupling constant is positive ( = +1) and negative ( = −1), respectively.
For a more general case, we shall consider the Palatini R + α R 2 gravity coupled with Born-Infeld NED. In this case, there is no exact expression of the metric functions. The metric functions can only be written in an integral form [45]. In Fig. 3, we rescale α as α/r 2 s → α, fix the value of the charge Q * = 0.2, and exhibit the QNM frequencies with respect to α (we shall mention that the qualitative behaviors of our results remain unchanged when we alter the values of Q * as long as the charge is smaller than its extremal value). It can be seen that when β m gets large, the frequencies remain almost constant when changing α. This is expectable because in this case, the NED reduces to linear Maxwell fields and the Palatini R +α R 2 reduces to GR in absence of a cosmological constant. Therefore, the QNMs reduce to those of the RN black hole.

Eikonal QNMs
In Ref. [68], it has been shown that in GR the QNMs in the eikonal limit (l ≈ μ → ∞) of any stationary, spherically symmetric, and asymptotically flat black hole can be deduced from the properties of the unstable null circular orbit around the black hole. More precisely, the QNM frequency in the eikonal limit can be expressed as [68] ω ≈ c l − i(n + 1/2)|λ c |, (5.2) where c is interpreted as the angular velocity of the null circular orbit and the parameter λ c is the Lyapunov exponent quantifying the instability of the orbit. The derivation of Eq. (5.2) is related to the fact that for these black holes, the potentials in the master equations within the eikonal limit can be approximated as It can be proven that the peak of this potential (5.3) coincides with the radius of the null circular orbit. After inserting the potential (5.3) into the the 1st order WKB formula, we can derive Eq. (5.2). It can be shown that this equation is valid as well in some modified theories of gravity. In fact, it can be seen from Ref. [39], and from the master equations (4.28) and (4.29) that the massless scalar field perturbations and the axial perturbations of a charged black hole in the Palatini f (R) gravity coupled with NED satisfy the approximation (5.2). The same is also valid for the massless scalar field perturbations of an EiBI charged black hole [39]. However, Eq. (5.2) may not be valid for the axial perturbations of the EiBI charged black holes. According to the master equations (4.72) and (4.73), the potentials in the eikonal limit Note that these approximated expressions are valid for both = ±1. Because of the factors β 2 g r 4 ± Q 2 , the relation between the eikonal QNMs and the properties of the null cir-cular orbit around the black hole would be violated. Instead, the QNMs of the axial perturbations described by the potentials (5.4) and (5.5) can be expressed as where i = 1, 2 and the index p denotes the quantities calculated at the peak of the potentials. Note that at the location In Fig. 4, we exhibit the eikonal QNM frequencies of the EiBI charged black holes in terms of 1/β g . The charge is fixed to Q * = 0.4 (we choose this value to amplify the effects of the charge on the QNMs). The blue (red) curves correspond to a positive (negative) EiBI coupling constant. The dashed and the dotted curves are the eikonal QNMs described by the potential V 1 and V 2 , respectively. We also present the eikonal QNMs for the massless scalar field perturbations which can be described by the potential (5.3) in the sense that V s = V (see Ref. [39]). It can be seen that the eikonal QNMs of the axial perturbations for the EiBI charged black holes (dashed and dotted curves) deviate from those corresponding to the unstable null circular orbit (solid curves).
Before closing this subsection, we would like to mention that the violation of Eq. (5.2) for the axial perturbations of the EiBI charged black holes could be due to the non-trivial coupling between the electromagnetic and the gravitational fields in this theory. On the other hand, if we assume that the electromagnetic perturbations do not alter the spacetime geometry, the electromagnetic perturbations will be described by the master equation (4.38) without the metric perturbation terms Fig. 4 The real part (upper) and imaginary part (lower) of the QNMs of the EiBI charged black holes in the eikonal limit (l → ∞) are shown with respect to 1/β g . The solid curves are included as well, indicating the QNMs of the massless scalar field perturbations which are described by the potential V s = V given in Eq. (5.3) on the right hand side. In this regard, the potential describing the electromagnetic perturbations in the eikonal limit can be approximated as Eq. (5.3), and the correspondence proposed in Ref. [68], i.e., Eq. (5.2) is satisfied.

Conclusions
In this paper, we consider specifically two gravitational theories within the Palatini formulation and study the QNMs of the axial perturbations for the charged black holes in these theories. These theories of gravity are, respectively, the Palatini f (R) gravity coupled with Born-Infeld NED and the EiBI gravity coupled with linear electromagnetic fields. One of our goals is to see how the Born-Infeld structures from the gravitational sector and from the matter sector change differently the QNM frequencies. Therefore, we pay special attention to the comparison between the QNMs of the EBI black holes and the EiBI charged black holes. The QNMs of the Born-Infeld charged black holes in the Palatini R + α R 2 gravity are also discussed. In fact, our paper can be regarded as a further extension of our previous work [39] in which we studied the QNMs of the massless scalar field perturbations to these different charged black holes.
By using the tetrad formalism, we have derived the coupled master equations describing the axial perturbations of the charged black holes. In the two theories that we are considering, the coupled equations reduce to those of the RN black holes when the ratio of the charge and the Born-Infeld coupling constant Q * /β m (or Q * /β g ) is small. The QNM frequencies of the charged black holes are evaluated by using the WKB method up to the 6th order, which is accurate for modes whose multiple is larger than the overtone l > n. In this paper, we mainly focus on the QNMs of the fundamental modes (n = 0), since these modes have the longest decay time and would dominate the late time ringdown signals from an astrophysical perspective. Our results indicate that the charged black holes are all stable against the axial perturbations. Besides, the QNM frequencies would deviate from those of the RN black hole when nonlinearity of matter fields (Born-Infeld NED) or modification of the gravitational theory (EiBI or f (R)) are considered. For instance, both the absolute value of the real part and the imaginary part of the QNM frequencies for the EBI charged black holes would increase with the value of 1/β m . On the other hand, we show that by increasing the value of 1/β g , the real part of the QNM frequencies and the decay time (∝ 1/|Im ω|) would increase (decrease) for the EiBI charged black holes with = +1 ( = −1).
Furthermore, we study the QNMs of these black holes in the eikonal limit (l → ∞). Interestingly, we find that the QNM frequencies in this limit for the EiBI charged black holes cannot be described by the properties of the unstable null circular orbit around the black hole. In other words, the QNM formula (5.2) proposed in Ref. [68] is not valid for the EiBI charged black holes. This violation could be an extra possible imprint from the EiBI corrections on the QNMs, aside from the QNM spectra, that may be detectable in the future.
In addition to the axial perturbations, it is necessary to investigate the QNMs of the polar perturbations (even parity perturbations) for the charged black holes considered in this work. For the Schwarzschild [72] and the RN charged black holes [59] in GR, it is well-known that their axial and polar perturbations are isospectral. This means that the potential terms in their master equations satisfy a certain relation in such a way that the QNMs of the axial and polar perturbations have identical spectra. The isospectrality could be violated in the presence of, for instance, nonlinearity in the matter source [73,74], or modifications of the Einstein-Hilbert action [16], and so on. The violation/fulfillment of the isospectrality for the charged black holes in the Palatini-type gravity theories could be an additional tool to test the underlying theories and we shall leave this interesting issue for a coming work. The ratio |∂ r * P c |/|P RN | for the two eigenvectors of the matrix V for the EBI charged black holes is presented with respect to r and β. We assume here Q * = 0.3 and l = 2. Those values are chosen to estimate the largest errors in our calculations. Note that the domain of r we are considering here is r ≥ r H , where r H is the horizon satisfying e 2ν = 0 can be decoupled by diagonalizing the matrix V to obtain its eigenvalues V 1 and V 2 : This is because the eigenvectors P = P RN of the matrix V are non-vanishing constant vectors for the RN black hole. Therefore, one can use a similarity transformation H (−) = PH (−) to rewrite the matrix equation (A1) as follows P ∂ 2 r * + ω 2 H (−) = VPH (−) .
By multiplying the above equation by P −1 , the equation can be decoupled. On the other hand, for the extended charged black holes considered in this paper, the eigenvectors of the matrix V are not constant vectors anymore. Instead, they would depend on r and, strictly speaking, one is not able to diagonalize V to decouple the master equations. It can be proven that for these extended charged black holes, the eigenvectors can be expressed as where P c (r ) stands for the correction term. If we use the same similarity transformation H (−) = PH (−) , the matrix equation (A1) can be rewritten as P ∂ 2 r * + ω 2 H (−) + 2 ∂ r * P c ∂ r * H (−) + ∂ 2 r * P c H (−) = VPH (−) . (A3) In general, Eq. (A3) cannot be decoupled due to the last two terms on the left hand side. However, we will argue that for our present work and for the parameter space of interest, the last two terms on the left hand side of Eq. (A3) are actually very small as compared with the other terms. The arguments are the following: (i) In this work, we calculate the QNM frequencies with a semi-analytical approach, which is formulated within the WKB approximation. For the cases where this approach is valid, the wave functionsH (−) can be solved with the WKB approximation. It can then be proven that H (−) and its derivatives (∂ r * H (−) and ∂ 2 r * H (−) ) have the same order of magnitudes (note that we have normalized all relevant quantities with respect to r s , so the magnitude of the frequencies would be of order one, which is also consistent with our results shown in Sect. 5). (ii) It can be shown that the magnitude of ∂ r * P c is very small as compared with the magnitude of P RN outside the event horizon and in the parameter space of our interest. In Fig. 5, we assume Q * = 0.3 and l = 2, and exhibit the smallness of the ratio |∂ r * P c |/|P RN | for the EBI black hole, with respect to r and the value of β. For the Born-Infeld black holes with f (R) being a quadratic function and for the EiBI charged black holes, this ratio is also very tiny.
According to the arguments above, the last two terms on the left hand side of Eq. (A3) are very small as compared to the other terms. Therefore, when studying the QNMs of the extended charged black holes, we shall omit these two terms and decouple the master equation (4.30) by diagonalizing the potential matrix V as given in Eq. (A2).