1 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_{\mu \nu }(x)\rightarrow \Omega ^{2}(x)g_{\mu \nu }(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 theory-dependent 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.

2 Four-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]

$$\begin{aligned} I=\frac{1}{16\pi }\int d^{4}x\sqrt{-g}C^{\mu \nu \rho \sigma }C_{\mu \nu \rho \sigma }, \end{aligned}$$
(1)

where \(C_{\mu \nu \rho \sigma }\) is the Weyl conformal tensor with the following explicit form

$$\begin{aligned} C_{\lambda \mu \nu \kappa }= & {} R_{\lambda \mu \nu \kappa }+{\frac{1}{6}}R_{\alpha }^{\alpha }\left[ g_{\lambda \nu }g_{\mu \kappa }-g_{\lambda \kappa }g_{\mu \nu }\right] \nonumber \\&-{\frac{1}{2}}\left[ g_{\lambda \nu }R_{\mu \kappa }-g_{\lambda \kappa }R_{\mu \nu }\right. \nonumber \\&\left. -g_{\mu \nu }R_{\lambda \kappa }+g_{\mu \kappa }R_{\lambda \nu }\right] , \end{aligned}$$
(2)

and the field equations can be obtained by taking a variation with respect to the metric tensor \(g_{\mu \nu }\)

$$\begin{aligned} W_{\rho \sigma }=\left( \nabla ^{\mu }\nabla ^{\nu }+\frac{1}{2}R^{\mu \nu }\right) C_{\rho \mu \nu \sigma }=0, \end{aligned}$$
(3)

which \(W_{\mu \nu }\) is the Bach tensor. It is straightforward to show that the following 4-dimensional line element satisfies all components of Eq. (3)

$$\begin{aligned} ds^{2}=-f(r)dt^{2}+f^{-1}(r)dr^{2}+r^{2}\left( d\theta ^{2}+\sin ^{2}\theta d\varphi ^{2}\right) , \end{aligned}$$
(4)

where the metric function is as follows

$$\begin{aligned} f\left( r\right) =c+\frac{d}{r}+\frac{c^{2}-1}{3d}r+br^{2}, \end{aligned}$$
(5)

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=-\Lambda /3\).

On the other hand, the line element describing the Schwarzschild-(a)dS spacetime in the radial coordinate \(\rho \) is given by

$$\begin{aligned} d{\tilde{s}}^{2}=-g(\rho )dt^{2}+g^{-1}(\rho )d\rho ^{2}+\rho ^{2}\left( d\theta ^{2}+\sin ^{2}\theta d\varphi ^{2}\right) , \end{aligned}$$
(6)

in which the metric function is

$$\begin{aligned} g\left( \rho \right) =1-\frac{2M}{\rho }-\frac{\Lambda }{3}\rho ^{2}, \end{aligned}$$
(7)

where M denotes the total mass of the black hole and \(\Lambda \) 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(\rho )\) so that \(ds^{2}=S(\rho )d {\tilde{s}}^{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(\rho )=\left( 1+q\rho \right) ^{-2}\) [60], one can find the relation \(\rho =r\left( 1-qr\right) ^{-1}\) between the radial coordinates \(\rho \) and r. Multiplying (6) by the conformal factor \(S(\rho )\) and replacing \( \rho \) by r, we can obtain the following relations between the parameters [60]

$$\begin{aligned} b=q^{2}\left( 1+2Mq\right) -\frac{\Lambda }{3};\ \ \ c=1+6Mq;\ \ \ d=-2M,\ \end{aligned}$$
(8)

where q is an arbitrary constant and we used \(\partial _{\rho }=\left( 1-qr\right) ^{2}\partial _{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.

3 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(\rho )d {\tilde{s}}^{2}\) so that \(ds^{2}\) is the line element of Weyl gravity (4 ), \(d{\tilde{s}}^{2}\) is the line element of the Schwarzschild-(a)dS black holes (6), and \(S(\rho )\) 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(\rho )\) . In other words, if we construct the axial perturbations of \(S(\rho )d {\tilde{s}}^{2}\) (for the general form of \(S(\rho )\)), 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(\rho )d{\tilde{s}} ^{2}\), is derived in the Appendix 1. The master equation of the axial perturbations of the following spacetime (see Appendix 1)

$$\begin{aligned} ds^{2}= & {} S(\rho )d{\tilde{s}}^{2}=S( \rho ) \left[ -g (\rho ) dt^{2}+\frac{d\rho ^{2}}{g (\rho ) }+\rho ^{2}\right. \nonumber \\&\quad \times \left. \left( d\theta ^{2}+\sin ^{2}\theta d\varphi ^{2}\right) \right] , \end{aligned}$$
(9)

is given by

$$\begin{aligned} \frac{d^{2}\Psi ^{\left( -\right) }\left( \rho _{*}\right) }{d\rho _{*}^{2}}+\left[ \omega ^{2}-V_{g}\left( \rho _{*}\right) \right] \Psi ^{\left( -\right) }\left( \rho _{*}\right) =0, \end{aligned}$$
(10)

where \(\omega \) is the QN frequency, \(\rho _{*}=\int \ g^{-1}\left( \rho \right) d\rho \) is the tortoise coordinate, and the effective potential \( V_{g}\left( \rho _{*}\right) \) reads

$$\begin{aligned} V_{g}\left( \rho _{*}\right) =\ g\left( \rho \right) \left[ \frac{\ell \left( \ell +1\right) }{\rho ^{2}}-\frac{2}{\rho ^{2}}-Z\frac{d}{d\rho } \left( \frac{\ g\left( \rho \right) }{Z^{2}}\frac{dZ}{d\rho }\right) \right] , \end{aligned}$$
(11)

in which \(Z=\rho \sqrt{S\left( \rho \right) }\). 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(\rho )=\left( 1+q\rho \right) ^{-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 \(\rho \) to r so that \(S(\rho )=\left( 1+q\rho \right) ^{-2}\) and \(\rho =r\left( 1-qr\right) ^{-1}\). By considering \(\partial _{\rho }=\left( 1-qr\right) ^{2}\partial _{r}\) and \(\partial _{r_{*}}=g\left( \rho \right) \left( 1-qr\right) ^{2}\partial _{r}\), one can find the wave equation (10) converts into

$$\begin{aligned} \frac{d^{2}\Psi ^{\left( -\right) }\left( r_{*}\right) }{dr_{*}^{2}}+ \left[ \omega ^{2}-V_{g}\left( r_{*}\right) \right] \Psi ^{\left( -\right) }\left( r_{*}\right) =0, \end{aligned}$$
(12)

where \(r_{*}\) is the new tortoise coordinate as \(dr/dr_{*}=\left( 1-qr\right) ^{2}\ {\bar{f}}\left( r\right) \) with

$$\begin{aligned} {\bar{f}}\left( r\right) =1-\frac{2M\left( 1-qr\right) }{r}-\frac{\Lambda r^{2} }{3\left( 1-qr\right) ^{2}}. \end{aligned}$$
(13)

The effective potential (11) now is given by

$$\begin{aligned}&V_{g}\left( r_{*}\right) =\ \left( 1-qr\right) ^{2}{\bar{f}}\left( r\right) \nonumber \\&\quad \times \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}-\frac{2}{r^{2}}-r \frac{d}{dr}\left( \frac{\ {\bar{f}}\left( r\right) }{r^{2}}\left( 1-qr\right) ^{2}\right) \right] , \end{aligned}$$
(14)

which we used \(Z=r\) and \(\partial _{\rho }=\left( 1-qr\right) ^{2}\partial _{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

$$\begin{aligned} V_{g}\left( r_{*}\right) =\left( 1-\frac{2M}{r}-\frac{\Lambda r^{2}}{3} \right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}-\frac{6M}{r^{3}} \right] , \end{aligned}$$
(15)

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

$$\begin{aligned} V_{g}\left( r_{*}\right)= & {} \ f\left( r\right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}-\frac{2}{r^{2}}-r\frac{d}{dr}\left( \frac{\ f\left( r\right) }{r^{2}}\right) \right] \end{aligned}$$
(16)
$$\begin{aligned}= & {} f\left( r\right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}+\frac{ 2\left( c-1\right) }{r^{2}}+\frac{c^{2}-1}{3rd}+\frac{3d}{r^{3}}\right] , \end{aligned}$$
(17)

which is a function of the free parameters of conformal gravity ( \(r_{*}=\int \ f^{-1}\left( r\right) 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=-\Lambda /3\) in (17).

4 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]

$$\begin{aligned}&\frac{d^{2}\Psi \left( r_{*}\right) }{dr_{*}^{2}}+\left[ \omega ^{2}-V_{s}\left( r_{*}\right) \right] \Psi \left( r_{*}\right) =0, \end{aligned}$$
(18)
$$\begin{aligned}&V_{s}\left( r_{*}\right) =f\left( r\right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}+\frac{f^{\prime }\left( r\right) }{r}\right] , \end{aligned}$$
(19)

where \(r_{*}=\int \ f^{-1}\left( r\right) 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 \( \rho \) coordinate are as follows [56]

$$\begin{aligned}&\frac{d^{2}\Psi \left( \rho _{*}\right) }{d\rho _{*}^{2}}+\left[ \omega ^{2}-V_{s}\left( \rho _{*}\right) \right] \Psi \left( \rho _{*}\right) =0, \end{aligned}$$
(20)
$$\begin{aligned}&V_{s}\left( \rho _{*}\right) =g\left( \rho \right) \left[ \frac{\ell \left( \ell +1\right) }{\rho ^{2}}+\frac{1}{Z}\frac{d}{d\rho }\left( g\left( \rho \right) \frac{dZ}{d\rho }\right) \right] , \end{aligned}$$
(21)

where \(Z=\rho \sqrt{S\left( \rho \right) }\) and \(\rho _{*}=\int \ g^{-1}\left( \rho \right) d\rho \) is the tortoise coordinate. Now, we apply the coordinate transformation \(\rho \) to r so that \(S(\rho )=\left( 1+q\rho \right) ^{-2}\) and \(\rho =r\left( 1-qr\right) ^{-1}\), as was described. Thus, the effective potential (21) reduces to

$$\begin{aligned}&V_{s}\left( r_{*}\right) =\ \left( 1-qr\right) ^{2}{\bar{f}}\left( r\right) \nonumber \\&\quad \times \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}+\frac{1}{r}\frac{d }{dr}\left( {\bar{f}}\left( r\right) \left( 1-qr\right) ^{2}\right) \right] . \end{aligned}$$
(22)

It is straightforward to show that \(\left( 1-qr\right) ^{2}{\bar{f}}\left( r\right) \) 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.

5 Electromagnetic perturbations

The wave equation and effective potential of electromagnetic perturbation in the spacetime background (4) are given by (see Appendix 2)

$$\begin{aligned}&\frac{d^{2}\Psi \left( r_{*}\right) }{dr_{*}^{2}}+\left[ \omega ^{2}-V_{e}\left( r_{*}\right) \right] \Psi \left( r_{*}\right) =0, \end{aligned}$$
(23)
$$\begin{aligned}&V_{e}\left( r_{*}\right) =f\left( r\right) \frac{\ell \left( \ell +1\right) }{r^{2}}, \end{aligned}$$
(24)

where \(r_{*}=\int \ f^{-1}\left( r\right) dr\) is the tortoise coordinate. On the other hand, the master equation of the perturbative conformally related Schwarzschild-(a)dS black holes in \(\rho \) coordinate is [56]

$$\begin{aligned} \frac{d^{2}\Psi \left( \rho _{*}\right) }{d\rho _{*}^{2}}+\left[ \omega ^{2}-V_{e}\left( \rho _{*}\right) \right] \Psi \left( \rho _{*}\right) =0, \end{aligned}$$
(25)

with

$$\begin{aligned} V_{e}\left( \rho _{*}\right) =g\left( \rho \right) \frac{\ell \left( \ell +1\right) }{\rho ^{2}}, \end{aligned}$$
(26)

where \(\rho _{*}=\int \ g^{-1}\left( \rho \right) d\rho \) is the tortoise coordinate. By applying the coordinate transformation \(\rho \) to r such that \(S(\rho )=\left( 1+q\rho \right) ^{-2}\) and \(\rho =r\left( 1-qr\right) ^{-1}\), the effective potential (26) reduces to

$$\begin{aligned} V_{e}\left( r_{*}\right) =\ \left( 1-qr\right) ^{2}{\bar{f}}\left( r\right) \frac{\ell \left( \ell +1\right) }{r^{2}}, \end{aligned}$$
(27)

which \(\left( 1-qr\right) ^{2}{\bar{f}}\left( r\right) \) 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

$$\begin{aligned} \frac{d^{2}\Psi \left( r_{*}\right) }{dr_{*}^{2}}+\left[ \omega ^{2}-V_{\ell }\left( r_{*}\right) \right] \Psi \left( r_{*}\right) =0, \end{aligned}$$
(28)

with the potential

$$\begin{aligned}&V_{\ell }\left( r_{*}\right) =f\left( r\right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}+\left( 1-s^{2}\right) \right. \nonumber \\&\quad \times \left. \left( \frac{\left( 4-s^{2}\right) b}{2}-\frac{d}{r^{3}}-\frac{s\left( c-1\right) }{3r^{2}} \right) +\left( 1-s\right) ^{2}\frac{c^{2}-1}{3rd}\right] ,\nonumber \\ \end{aligned}$$
(29)

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=-\Lambda /3\) into ( 29), one can obtain

$$\begin{aligned} V_{\ell }\left( r_{*}\right) =f\left( r\right) \left[ \frac{\ell \left( \ell +1\right) }{r^{2}}+\left( 1-s^{2}\right) \left( \frac{2M}{r^{3}}-\frac{ \left( 4-s^{2}\right) \Lambda }{6}\right) \right] , \end{aligned}$$
(30)

for the Schwarzschild-(a)dS black holes [48].

6 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 \(\Psi \left( r_{*}\right) \) at \(r_{*}\rightarrow \pm \infty \) 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

$$\begin{aligned} \left\{ \begin{array}{cc} \Psi \left( r_{*}\right) \sim e^{-i\omega r_{*}}, &{} r_{*}\rightarrow -\infty \ \left( r\rightarrow r_{e}\right) \\ \Psi \left( r_{*}\right) \sim e^{i\omega r_{*}}, &{} r_{*}\rightarrow +\infty \ \left( r\rightarrow r_{c}\right) \end{array} \right. , \end{aligned}$$
(31)

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 \(\Psi \left( r_{*}\right) \) 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

$$\begin{aligned} \omega ^{2}=V_{0}+\sum _{j=1}\Omega _{2j}-i\left( n+\frac{1}{2}\right) \sqrt{ -2V_{0}^{\prime \prime }}\left( 1+\sum _{j=1}\Omega _{2j+1}\right) , \end{aligned}$$
(32)

where n is the overtone number, \(V_{0}\) is the value of the effective potential at its local maximum, and \(\Omega _{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

$$\begin{aligned} -4\frac{\partial ^{2}\Psi \left( u,v\right) }{\partial u\partial v}=V_{\ell }\left( u,v\right) \Psi \left( u,v\right) , \end{aligned}$$
(33)

in terms of the light-cone coordinates \(u=t-x\) and \(v=t+x\), and \(\Psi \) assumed to have time dependence \(e^{-i\omega 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 \(\Psi \left( u,v\right) \) and expanding this operator for sufficiently small grids

$$\begin{aligned}&\Psi \left( u+\Delta ,v+\Delta \right) =\Psi \left( u+\Delta ,v\right) \nonumber \\&\quad +\Psi \left( u,v+\Delta \right) -\Psi \left( u,v\right) \nonumber \\&\quad -\frac{\Delta ^{2}}{8}\left[ V_{\ell }\left( u+\Delta ,v\right) \Psi \left( u+\Delta ,v\right) \right. \nonumber \\&\quad \left. +V_{\ell }\left( u,v+\Delta \right) \Psi \left( u,v+\Delta \right) \right] , \end{aligned}$$
(34)

which \(\Delta \) 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 1, 2, 3 and 4. We calculated the lowest frequencies (Tables 1, 2) and the second overtone (Tables 3, 4) for some values of the free parameters b, c, d, and \(\ell \). The frequencies are written as \(\omega =\omega _{r}-i\omega _{i}\) where \(\omega _{r}\) (\(\omega _{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 \(\Psi \left( t,r\right) \) 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.

Table 1 The fundamental QNMs for \(\ell =2\), calculated by using the AIM method (first row) and WKB formula (second row)
Table 2 The fundamental QNMs for \(\ell =3\), calculated by using the AIM method (first row) and WKB formula (second row)
Table 3 The QNMs for \(n=1\) and \(\ell =2\), calculated by using the AIM method (first row) and WKB formula (second row)
Table 4 The QNMs for \(n=1\) and \(\ell =3\), calculated by using the AIM method (first row) and WKB formula (second row)
Fig. 1
figure 1

This figure evaluated at \(r=2\) (\(r_{e}<2<r_{c}\)) for \(\ell =2\) and \( c=0.3\). The left (right) panel indicates the (absolute) value of the wave function \(\Psi \left( t,r\right) \) of gravitational perturbations as a function of time

Fig. 2
figure 2

This figure evaluated at \(r=2\) for \(\ell =2\), \(b=-0.1\), and \(c=0.3\) . It shows the QN modes of scalar (\(s=0\)) and electromagnetic (\(s=1\)) perturbations alongside the gravitational ones (\(s=2\)). By comparing the left panel and right panel we find that the modes of all perturbations decay faster with more oscillations by increasing d

Fig. 3
figure 3

\(Re[\omega ]\) and \(Im[\omega ]\) as a function of c for \(\Lambda =0.02\) and unit mass calculated by using the sixth order WKB formula. The vertical line indicates the frequencies of the Schwarzschild-dS black hole. Both the real and imaginary parts of QN frequencies decrease whenever c deviates from one

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=-\Lambda /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

$$\begin{aligned} f\left( r\right) =c-\frac{2M}{r}-\frac{c^{2}-1}{6M}r-\frac{\Lambda }{3}r^{2}, \end{aligned}$$
(35)

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\le -0.5\) and \(c\ge 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 (\(\ell \rightarrow \infty \)) can be obtained by the first order WKB formula (32)

$$\begin{aligned} \omega ^{2}=V_{0}-i\left( n+\frac{1}{2}\right) \sqrt{-2V_{0}^{\prime \prime } }, \end{aligned}$$
(36)

and in this regime, the effective potential (29) is given by

$$\begin{aligned} V_{\ell }\left( r_{*}\right) \sim f\left( r\right) \frac{\ell ^{2}}{r^{2} }, \end{aligned}$$
(37)

which is still a function of c due to the presence of \(f\left( r\right) \). Interestingly, we find that these black holes, unlike the non-singular black holes in conformal gravity [57], can be distinguished from the Schwarzschild ones even in the eikonal limit.

6.1 Error estimation of QN frequencies

Although the WKB formula gives the best accuracy at \(\ell >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

$$\begin{aligned} \Delta _{k}=\frac{\left| \omega _{k+1}-\omega _{k-1}\right| }{2}, \end{aligned}$$
(38)

for the error estimation of \(\omega _{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]

$$\begin{aligned} \Delta _{k} > rsim \left| \omega -\omega _{k}\right| , \end{aligned}$$
(39)

where \(\omega \) 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 5 and 6, we summarize the error estimation of the WKB formula for the fundamental and first overtone QNMs of fixed \(\ell =2\), \( b=-0.08\), \(c=0.3\), and \(d=-1\) (related to the last line of Tables 1 and 3). From Table 5, 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 6 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 1 (3) are reliable. Our calculations based on (39), not shown here due to economic reason, show that the frequencies given in Tables 1, 2, 3 and 4 are reliable up to a minimum of 3 digits.

7 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.

Table 5 The fundamental modes of gravitational, electromagnetic, and scalar perturbations calculated with the WKB formula of different orders. The calculated frequency is related to the last line of Table 1
Table 6 The first overtone frequency of gravitational, electromagnetic, and scalar perturbations calculated with the WKB formula of different orders. The calculated frequency is related to the last line of Table 3

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.