Quasinormal modes of black holes in Weyl gravity: Electromagnetic and gravitational perturbations

The recent reported gravitational wave detection motivates one to investigate the properties of different black hole models, especially their behavior under (axial) gravitational perturbation. Here, we study the quasinormal modes of black holes in Weyl gravity. We derive the master equation describing the quasinormal radiation by using a relation between the Schwarzschild-anti de Sitter black holes and Weyl solutions, and also the conformal invariance property of the Weyl action. It will be observed that the quasinormal mode spectra of the Weyl solutions deviate from those of the Schwarzschild black hole due to the presence of an additional linear $r$-term in the metric function. We also consider the evolution of the Maxwell field on the background spacetime and obtain the master equation of electromagnetic perturbations. Then, we use the WKB approximation and asymptotic iteration method to calculate the quasinormal frequencies. Finally, the time evolution of modes is studied through the time-domain integration of the master equation.

The recent reported gravitational wave detection motivates one to investigate the properties of different black hole models, especially their behavior under (axial) gravitational perturbation. Here, we study the quasinormal modes of black holes in Weyl gravity. We derive the master equation describing the quasinormal radiation by using a relation between the Schwarzschild-anti de Sitter black holes and Weyl solutions, and also the conformal invariance property of the Weyl action. It will be observed that the quasinormal mode spectra of the Weyl solutions deviate from those of the Schwarzschild black hole due to the presence of an additional linear r-term in the metric function. We also consider the evolution of the Maxwell field on the background spacetime and obtain the master equation of electromagnetic perturbations. Then, we use the WKB approximation and asymptotic iteration method to calculate the quasinormal frequencies. Finally, the time evolution of modes is studied through the time-domain integration of the master equation.

I. INTRODUCTION
From theoretical point of view, black hole is one of the most interesting and important solutions to gravity theories. The existence of such highly dense object in the cosmos has been also proved through the detection of gravitational waves (GWs) of black hole binary mergers [1] by the LIGO and VIRGO observatories and the captured image of the 'shadow' of a supermassive black hole by the Event Horizon Telescope collaboration [2,3]. In this regard, it is interesting to investigate the other kinds of black hole solutions such as ones which are constructed based on higher curvature theories of gravity.
Weyl gravity, which its Lagrangian is defined by the square of the Weyl tensor, is one of the successful and interesting theories in higher derivative gravity scenario [4][5][6][7][8]. This model of gravity is invariant under local scale transformation of the metric g µν (x) → Ω 2 (x)g µν (x), and thus, it is unique up to the choice of the matter source in order to keep the conformal invariance property. The Weyl gravity suffers the Weyl ghost that it might be possible to remove it under certain conditions [9][10][11][12][13][14][15][16]. In addition, one can consider this theory of gravity as a suitable model for constructing quantum gravity [17,18], since it is a higher-curvature theory of gravity which is power-counting renormalizable [19,20].
From the viewpoint of high energy physics, it is shown that the Weyl gravity arises from twister-string theory with both closed strings and gauge-singlet open strings [21], and it appears as a counterterm in adS/CFT calculations [22][23][24]. In addition, this theory can be examined as a possible UV completion of Einstein gravity [25,26]. It is worthwhile to mention that there is an equivalence between Einstein gravity and Weyl gravity by considering the Neumann boundary conditions [27,28].
The discovery of GWs produced by the merger of compact binary objects [1,29] added a totally different and new type of observations to the traditional astronomy based on electromagnetic waves. The emitted GWs contain a lot of information for fundamental physics and one can check the validity of the alternative theories of general relativity, such as massive gravity, scalar-tensor theory, and Weyl gravity. The signal of compact binary merger can be decomposed into three phases [30]. The first stage of a binary system is the inspiral phase which the frequency and amplitude of the signal chirp with time. During this phase, the signal is universal that just depends on the masses and the spins of compact objects and does not depend on the nature of source. The post-Newtonian approximation is a powerful tool for describing the inspiral phase [30]. The second stage is the merger phase and occurs after the inspiral phase with a rapid collapse of the two objects to form a black hole. The amplitude of GWs has a peak at this time and numerical relativity is used to compute the merger phase [31][32][33]. The ringdown phase is the final stage and describes the evolution of the new black hole. This new black hole is highly deformed due to the nonlinear dynamics of the collision. The perturbed black hole emits GWs in the form of quasinormal radiation and the perturbation theory can be used to calculate the quasinormal modes (QNMs) [34].
The ringdown phase of a compact binary merger might be better to study compared to the inspiral or merger phase depending on the physics being tested [35]. The nature of the binary components shows up only at high post-Newtonian order in the inspiral phase while investigating the merger phase requires time-consuming and theorydependent numerical simulations. On the other hand, the ringdown phase can be appropriately investigated by perturbation theory and it is relatively simple to model beyond Einstein gravity. The perturbation theory of black holes was started by the pioneering work of Regge and Wheeler [36] and was continued by Zerilli [37]. They have found the wave equations of axial and polar perturbations of the Schwarzschild black hole and examined the dynamical stability of the black hole under small perturbations. The frequencies of perturbations are called the QN frequencies (QNFs) and have been calculated by using analytic and semi-analytic approaches [38][39][40][41], and also, several numerical methods [42][43][44][45][46] (see [47][48][49] for reviews on QNMs). The electromagnetic and gravitational perturbations of the Schwarzschild black holes have been investigated in [50]. Besides, the QNMs for the gravitational and electromagnetic perturbations in modified gravity are calculated before [51]. In the case of conformal gravity, by imposing perturbations in Minkowskian spacetime, the GWs have been investigated and the effective energy-momentum tensor of the gravitational radiation is calculated [52]. The astrophysical GWs of inspiralling compact binaries have been investigated [53][54][55]. Moreover, the scalar, electromagnetic [56], and axial perturbations [57] of nonsingular black holes in conformal gravity have been studied. It was shown that it is possible to find the black hole solutions in Weyl gravity which are both thermally and dynamically stable under massive scalar perturbations, and also, the QNMs of this theory of gravity in asymptotically adS spacetime were obtained [58]. In addition, the nearly extreme black holes in Weyl gravity have been studied and an exact formula for the QNFs with an upper bound on them has been found [59].
In this paper, our main goal is studying the axial gravitational perturbations of singular black holes in Weyl gravity in order to investigate the QN radiation of the ringdown phase. We first obtain a master wave equation for an arbitrary metric which is conformally related to the Schwarzschild-(anti) de Sitter (Schwarzschild-(a)dS) spacetime. Then, by using the Weyl invariance property of the action, and also, the relation between the Schwarzschild-(a)dS black holes and Weyl solutions, we derive the wave equation of the axial gravitational perturbations of black holes in Weyl gravity. In addition, we consider both the scalar and electromagnetic perturbations in the background spacetime of these black holes. Then, we calculate the QNMs by employing the sixth order WKB approximation and the asymptotic iteration method (AIM). The time evolution of modes is also investigated by using the discretization scheme.

II. 4-DIMENSIONAL BLACK HOLES IN WEYL GRAVITY
Here, we give a brief review on the four-dimensional black holes in Weyl gravity and the connection between these solutions and the Schwarzschild-(a)dS black holes. The action of Weyl gravity is given by [60] where C µνρσ is the Weyl conformal tensor with the following explicit form and the field equations can be obtained by taking a variation with respect to the metric tensor g µν which W µν is the Bach tensor. It is straightforward to show that the following 4-dimensional line element satisfies all components of Eq. (3) where the metric function is as follows in which b, c, and d are three integration constants. It is notable that in contrast with the Einstein gravity that the cosmological constant should be considered in the action by hand, it is present as an integration constant in the Weyl gravity solutions. It is also worthwhile to mention that one can recover the Schwarzschild-(a)dS black hole by setting c = 1, d = −2M , and b = −Λ/3.
On the other hand, the line element describing the Schwarzschild-(a)dS spacetime in the radial coordinate ρ is given by in which the metric function is where M denotes the total mass of the black hole and Λ is the cosmological constant. One can show that there is a conformal relation between the black hole spacetimes in Weyl gravity (4) and Einstein theory (6). Indeed, these two spacetimes can be connected to each other by introducing a conformal factor S(ρ) so that ds 2 = S(ρ)ds 2 . Every spacetime like Schwarzschild-(a)dS case which is conformally related to (4) is also a solution of the field equations of Weyl gravity (3) since all the metrics that transform conformally are equivalent. By considering the conformal factor S(ρ) = (1 + qρ) −2 [60], one can find the relation ρ = r (1 − qr) −1 between the radial coordinates ρ and r. Multiplying (6) by the conformal factor S(ρ) and replacing ρ by r, we can obtain the following relations between the parameters where q is an arbitrary constant and we used ∂ ρ = (1 − qr) 2 ∂ r . Therefore, we have a spectrum of conformal solutions depending on the values of q. We shall use these relations to obtain the axial perturbation of Weyl gravity in the coming section.

III. GRAVITATIONAL PERTURBATIONS OF WEYL GRAVITY
Here, we are going to obtain a master equation for the axial gravitational perturbations of black holes described by the metric function (5). First, one may note that there is a conformal relation between the line element (4) and the Schwarzschild black holes in asymptotically adS spacetime. Indeed, the spacetime of Weyl solutions and the spacetime of Schwarzschild-(a)dS solutions are related to each other as ds 2 = S(ρ)ds 2 so that ds 2 is the line element of Weyl gravity (4), ds 2 is the line element of the Schwarzschild-(a)dS black holes (6), and S(ρ) is the conformal factor. Thus, if we obtain the master equation of black holes conformally related to the Schwarzschild-(a)dS black holes, it can also describe the gravitational perturbations of Weyl solutions by replacing the explicit form of S(ρ). In other words, if we construct the axial perturbations of S(ρ)ds 2 (for the general form of S(ρ)), it is equivalent to the master equation of Weyl solutions. It is worthwhile to mention that it is not possible to construct a second-order master wave equation by considering small perturbations in (4) and using the field equations (3) because there is fourth-order differential in the field equations and the linearized field equations are fourth-order. Here, since our goal is to find a second order wave equation for Weyl solutions, we shall use the conformal relation between Schwarzschild-(a)dS black holes and Weyl solutions.
The gravitational perturbations of black holes conformally related to the Schwarzschild-(a)dS spacetime, i.e. the perturbations of S(ρ)ds 2 , is derived in the appendix A. The master equation of the axial perturbations of the following spacetime (see APPENDIX A) is given by where ω is the QN frequency, ρ * = g −1 (ρ) dρ is the tortoise coordinate, and the effective potential V g (ρ * ) reads in which Z = ρ S (ρ). Note that the right-hand side of (9) is exactly equal to the Weyl gravity spacetime (4) under the conditions (8) whenever we consider the conformal factor S(ρ) = (1 + qρ) −2 . As a result, if we consider a coordinate transformation (in (9)-(11)) obeying this conformal factor, we can obtain the axial perturbations of singular black holes in Weyl gravity. Thus, we apply the coordinate transformation ρ to r so that S(ρ) = (1 + qρ) −2 and ρ = r (1 − qr) −1 . By considering ∂ ρ = (1 − qr) 2 ∂ r and ∂ r * = g (ρ) (1 − qr) 2 ∂ r , one can find the wave equation (10) converts into where r * is the new tortoise coordinate as dr/dr The effective potential (11) now is given by which we used Z = r and ∂ ρ = (1 − qr) 2 ∂ r . Therefore, Eq. (12) is the master equation of the axial perturbations of black holes in the Weyl gravity (5). In addition, Eq. (14) is the effective potential of perturbations and the free parameters b, c, and d of (5) are related to the parameters of (13) through the conditions (8). It is worthwhile to mention that for q = 0, the potential (14) reduces to which is the effective potential of the axial perturbations of the Schwarzschild-(a)dS black hole [61], as we expected.
In order to get rid of the free parameter q in (14), one can apply the conditions (8) in (14) to obtain the following effective potential which is a function of the free parameters of conformal gravity ( r * = f −1 (r) dr being the tortoise coordinate). Now, we can recover the axial perturbations of the Schwarzschild-(a)dS black hole (15) by setting c = 1, d = −2M , and b = −Λ/3 in (17).

IV. SCALAR PERTURBATIONS
Now, in order to ensure that our calculations of obtaining the axial perturbations of Weyl gravity are correct, we compare the effective potential of the massless scalar perturbations of black holes in Weyl gravity obtained by two methods; one is considering the evolution of a scalar field in the spacetime background (4) directly, and the other one is multiplying the Schwarzschild spacetime by a conformal factor (9) and obtaining the related effective potential.
The wave equation and the effective potential of a massless scalar perturbation in the spacetime background (4) are given by [59] where r * = f −1 (r) dr is the tortoise coordinate. On the other hand, the wave equation and the effective potential of the perturbative conformally related Schwarzschild-(a)dS black holes in ρ coordinate are as follows [56] where Z = ρ S (ρ) and ρ * = g −1 (ρ) dρ is the tortoise coordinate. Now, we apply the coordinate transformation ρ to r so that S(ρ) = (1 + qρ) −2 and ρ = r (1 − qr) −1 , as was described. Thus, the effective potential (21) reduces to It is straightforward to show that (1 − qr) 2f (r) is equal to the metric function of Weyl gravity (5) with the help of obtained conditions (8), and thus the effective potential (22) is equal to Eq. (19). Therefore, this comparison of scalar perturbations shows that our calculations of obtaining the axial perturbations of Weyl solutions given in the previous section are indeed correct.

V. ELECTROMAGNETIC PERTURBATIONS
The wave equation and effective potential of electromagnetic perturbation in the spacetime background (4) are given by (see APPENDIX B) where r * = f −1 (r) dr is the tortoise coordinate. On the other hand, the master equation of the perturbative conformally related Schwarzschild-(a)dS black holes in ρ coordinate is [56] with where ρ * = g −1 (ρ) dρ is the tortoise coordinate. By applying the coordinate transformation ρ to r such that S(ρ) = (1 + qρ) −2 and ρ = r (1 − qr) −1 , the effective potential (26) reduces to which (1 − qr) 2f (r) is equal to the metric function of Weyl gravity (5) by considering the conditions (8). Therefore, the effective potentials (24) and (27) are the same.
It is worthwhile to mention that scalar, electromagnetic, and gravitational perturbations of Weyl gravity can be collected and described by the following master equation with the potential where s = 0, 1, 2 is the spin of the perturbing field and we used the effective potentials given in Eqs. (17), (19), and (24) to obtain this equation. By inserting c = 1, d = −2M , and b = −Λ/3 into (29), one can obtain for the Schwarzschild-(a)dS black holes [48].

VI. QUASINORMAL MODES
Here, we consider the master equation (28) with the effective potential (29) for s = 0, 1, 2 as the results of perturbations of Weyl gravity. The spectrum of QNMs is the solution of the wave equation (28) and this spectrum becomes discrete when we impose some proper boundary conditions. The boundary conditions are applied to the modes Ψ (r * ) at r * → ±∞ which can be obtained by studying the flux of radiation detected by physical observers near the event horizon and cosmological horizon. The observers detect outgoing waves at the cosmological horizon and incoming radiation at the event horizon where r e is the event horizon and r c denotes the cosmological horizon. Now, we concentrate our attention on the conformal-dS solutions (b < 0) and calculate the QN frequencies by using the WKB formula as a semi-analytic approach and AIM as a numerical method. The WKB approximation is based on the matching of WKB expansion of the modes Ψ (r * ) at the event horizon and cosmological horizon with the Taylor expansion near the peak of the potential barrier. This method can be used for an effective potential that forms a potential barrier with a single peak. It was first applied to the problem of scattering around black holes [39], and then extended to the 3rd order [40], 6th order [41], and 13th order [63]. The WKB formula is given by where n is the overtone number, V 0 is the value of the effective potential at its local maximum, and Ω k 's denote the kth order of the approximation and depend on the value of the effective potential and its derivatives at the local maximum. It is notable that the explicit form of the WKB corrections is given in [40,41]. We shall use this formula up to the sixth order as a semi-analytical approach to obtain the QNFs of perturbations.
On the other hand, the AIM has been employed for solving the eigenvalue problems and second-order differential equations [64,65]. It was also shown that the improved AIM can be used as an accurate numerical method for calculating the QNMs [45,66]. Here, we will use this method up to 15 iterations as a numerical approach to obtain the QNFs of perturbations.
In addition, we can investigate the contribution of all modes by using the time-domain integration of the wavelike equation (28). The time-domain profile of modes shows the time evolution of modes at the ringdown stage and the behavior of the asymptotic tails at late times. In order to obtain the time evolution of modes, we follow the discrimination scheme given in [67]. The perturbation equation (28) takes the following form in terms of the light-cone coordinates u = t − x and v = t + x, and Ψ assumed to have time dependence e −iωt .
We can obtain the time-domain profile of modes by integrating this equation on the small grids from the two null surfaces u = u 0 and v = v 0 . One can obtain the evolution equation in the light-cone coordinates by applying the time evolution operator on Ψ (u, v) and expanding this operator for sufficiently small grids which ∆ is the step size of the grids. We shall obtain the time evolution of perturbations with a Gaussian wave packet on the surfaces u = u 0 and v = v 0 as initial data. The QNFs of scalar (s = 0), electromagnetic (s = 1), and gravitational (s = 2) perturbations are presented in tables I − IV . We calculated the lowest frequencies (tables I and II) and the second overtone (tables III and IV ) for some values of the free parameters b, c, d, and ℓ. The frequencies are written as ω = ω r − iω i where ω r (ω i ) is the real (imaginary) part of the frequencies. The obtained frequencies for gravitational perturbation show that the WKB formula and the AI method are in a good agreement. The modes of gravitational perturbation live longer with lower frequency compared to the scalar and electromagnetic perturbations. In addition, the all kinds of perturbations decay faster with more oscillations by increasing d and/or b. However, the effective potential is positive and all frequencies have a negative imaginary part which indicates that these kinds of perturbations will decay with time, and thus, the spacetime is stable.      The time-domain profile of modes for different perturbations is presented in Figs. 1 and 2 for some fixed values of the free parameters. According to the Fig. 1, we can observe that the QN oscillations of the wave function Ψ (t, r) at the ringdown phase of gravitational perturbations for early and intermediate times. This figure confirms that the wave function oscillates with a frequency that increases and decay faster with increasing b and/or d. In addition, the time evolution of modes for scalar and electromagnetic perturbations is illustrated in Fig. 2. This figure shows that the QN oscillations of gravitational perturbation live longer with lower frequency compared to the scalar and electromagnetic perturbations.
In order to obtain the deviations of Weyl solutions from the Schwarzschild-dS black holes, we compare the QNM spectra of the both cases. To do so, we first consider d = −2M and b = −Λ/3 in the metric function (5), and then calculate the gravitational QN modes. In this way, the metric function of conformal gravity takes the form and thus, the value of c characterizes the deviation. As c gets away from 1, both the real and imaginary parts of frequencies decrease which shows that the perturbations in the conformal black holes' background live longer compared to the Schwarzschild ones (see Fig. 3). It is notable that in order to have black hole solutions, the value of c cannot be chosen from c ≤ −0.5 and c ≥ 2, and also, the allowed range for nearly extreme black holes is −1 < c < 2 [59].
In addition, the exact relation for QNFs in the eikonal limit (ℓ → ∞) can be obtained by the first order WKB formula (32) and in this regime, the effective potential (29) is given by which is still a function of c due to the presence of f (r). Interestingly, we find that these black holes, unlike the nonsingular black holes in conformal gravity [57], can be distinguished from the Schwarzschild ones even in the eikonal limit.

A. Error estimation of QN frequencies
Although the WKB formula gives the best accuracy at ℓ > n and gives us an accurate and economic way to compute the QN frequencies, this method does not always give a reliable result and neither guarantees a good estimation for the error [41,68]. In addition, we cannot always increase the WKB order to obtain a better approximation for the frequency. In practice, there is some order that the WKB formula (32) provides the best approximation, and the error increases as the order of the formula increases. On the other hand, since the AIM relies upon a specialized barrier function, there is no systematic way to estimate the errors or to improve the accuracy [45]. In order to estimate the error of the WKB approximation, we can compare two sequential orders of the formula (32). However, since each WKB order correction affects either real or imaginary part of the squared frequencies, we should use the following quantity for the error estimation of ω k that is calculated with the WKB formula of the order k. This quantity allows determining the WKB order in which the error is minimal. In the case of the Schwarzschild black holes, it is shown that the error estimation (38) provides a very good estimation of the error order, satisfying [68] ∆ k |ω − ω k | , where ω is an accurate value of the QN frequency calculated by using numerical methods. Thus, we can use the error estimation (38) to find the order of the absolute error and determine the order which gives the most accurate approximation for the QN mode.
In the tables V and V I, we summarize the error estimation of the WKB formula for the fundamental and first overtone QNMs of fixed ℓ = 2, b = −0.08, c = 0.3, and d = −1 (related to the last line of tables I and III). From table V , we can see that the best order of the WKB formula for calculating the QN frequency of gravitational perturbations is 7th-order whereas the QN frequency of electromagnetic perturbations has the best accuracy with the help of 12th-order. In the case of scalar perturbations, the error of the WKB formula of 10th-order is minimal. Table V I shows the same results for the first overtone frequency, except for gravitational perturbations. The QN frequency of gravitational perturbations has the best accuracy with the help of 8th-order. If we assume that the band (39) is also correct for the Weyl solutions, since the maximum estimation of the error for the 6th-order WKB formula is of order 10 −6 and 10 −5 , up to 4 (3) digits of the frequency in the last line of table I (III) are reliable. Our calculations based on (39), not shown here due to economic reason, show that the frequencies given in tables I − IV are reliable up to a minimum of 3 digits.

VII. CONCLUSIONS
We have investigated the effects of both the axial gravitational and electromagnetic perturbations on a black hole system in Weyl gravity. We have derived the master equation, describing the QN radiation, from the conformal invariance property of the Weyl action, and also, a relation between the Schwarzschild-(a)dS black holes and Weyl solutions. We have found that the QNM spectra of the Weyl solutions deviate from those of the Schwarzschild black hole due to the presence of a linear r-term in the metric function. We have seen that, unlike the non-singular black holes in conformal gravity, this deviation was present even in the eikonal regime. Thus, it will be possible to test the Weyl solutions (or at least find a constraint on the free parameter c in order to recover the present universe after a phase transition where the conformal symmetry is broken) with the help of future gravitational wave detectors. Moreover, it was shown that the perturbations in the conformal black holes' background live longer compared to the Schwarzschild ones.
In addition, we have calculated the QN frequencies of scalar, electromagnetic, and gravitational perturbations through both the sixth order WKB approximation and the improved AIM after 15 iterations. For the obtained frequencies, the effective potential was positive and all the frequencies had a negative imaginary part. Therefore, one can obtain some stable black hole solutions in conformal gravity under these kinds of perturbations.
We have found that the QN modes of gravitational perturbations live longer with lower frequency compared to the scalar and electromagnetic perturbations. In addition, all kinds of perturbations decay faster with more oscillations by increasing the free parameters d and/or b. Furthermore, the time evolution of different perturbations for early and intermediate times is studied by using the time-domain integration of the master equation. The time-domain profile of modes confirmed the previous results mentioned above.

VIII. ACKNOWLEDGEMENTS
We wish to thank Shiraz University Research Council.
Appendix A: Gravitational perturbations of black holes conformally related to the Schwarzschild-(a)dS solutions Assume a black hole spacetime is conformally related to the Schwarzschild-(a)dS solutions so that Multiplying the Schwarzschild spacetime by a conformal factor can be described by an anisotropic fluid with the following effective energy-momentum tensor [57] where p 1 and p 2 are, respectively, the radial pressure and the tangential pressure, andρ is the energy density measured by a comoving observer. The explicit expressions ofρ, p 1 , and p 2 are functions of ρ and they can be derived by calculating the corresponding field equations constructed from the line element (A1). In addition, u µ is the timelike four-velocity and x µ is the spacelike unit vector (orthogonal to u µ and angular directions) and they satisfy in which the indices are raised and lowered by the metricĝ µν , and we also assumed u µ = [u t , 0, 0, 0] and x µ = [0, x ρ , 0, 0] in the comoving frame.
In order to obtain the master wave equation and the effective potential of the line element (A1), we follow Chandrasekhar and his notation is used [69]. The axial perturbations are characterized by introducing the non-vanishing parameters σ, q 2 , and q 3 in the unperturbed spacetime in the following form ds 2 = −e 2ν dt 2 + e 2ψ dϕ − σdt − q 2 dx 2 − q 3 dx 3 2 + e 2µ2 dx 2 2 + e 2µ3 dx 3 2 in which ν, ψ, µ 2 , µ 3 , σ, q 2 , and q 3 are functions of t, x 2 (radial coordinate ρ), and x 3 (polar angle θ) so that σ, q 2 , and q 3 are small quantities. The components of the unperturbed metric (σ, q 2 , q 3 = 0) are as follows e ν = S(ρ)g(ρ); e µ2 = S(ρ) g(ρ) ; e µ3 = ρ S(ρ); e ψ = ρ sin (θ) S(ρ), (A5) The proper field equations, describing this unperturbed metric, are where G µν is the Einstein tensor. Here, following [57], we use the tetrad formalism to derive the master equation because the axial components of the perturbed energy-momentum tensor in the tetrad frame also vanish in this case, and thus,ρ , p 1 , and p 2 have nothing to do with the master equation. The tetrad basis corresponding to the line element (A4) is defined as follows where the Greek letters (labelling the columns) are the tensor indices and Latin letters enclosed in parentheses (labelling the rows) are the tetrad indices run from 0 to 3. Based on this formalism, any vector or tensor field can be projected onto the tetrad frame to obtain its tetrad components in which η (a)(b) = e µ (a) e µ(b) is a constant symmetric matrix. In the tetrad frame, the perturbed field equations (A7) read where the perturbed energy-momentum tensor is as follows δT (a)(b) = (ρ + p 2 ) δ u (a) u (b) + (δρ + δp 2 ) u (a) u (b) By considering the constraints (A3) and u µ x µ = 0, we find that the axial components of the perturbed energymomentum tensor in the tetrad frame vanish We now consider the following decomposition where C −3/2 ℓ+2 (θ) is the Gegenbauer function governed by the following differential equation