Stability and Quasinormal Modes of Black Holes in Conformal Weyl Gravity

In this paper, we first investigate the thermal stability of black holes in conformal Weyl gravity with a comparison with the Schwarzschild black holes. Then, we consider a minimally coupled massive scalar perturbation and calculate the quasinormal modes in asymptotically dS spacetime by employing the sixth order WKB approximation and asymptotic iteration method. The deviations from those of the Schwarzschild-dS solutions are obtained and the possibility of the presence of quasi-resonance modes for Weyl black hole solutions is investigated. Finally, we consider a massless scalar perturbation in the background of asymptotically AdS solutions and calculate the quasinormal modes by using the pseudospectral method. The effects of the free parameter of the theory on the quasinormal modes are studied and deviations from those of the Schwarzschild-AdS black holes are investigated. The imaginary part of quasinormal frequencies in AdS spacetime is the time scale of a thermal state (in the conformal field theory) to approach thermal equilibrium.

black holes attracted attention during the past three years (for instance, see an incomplete list [22][23][24][25][26] and references therein). In the case of conformal gravity, the gravitational waves of small perturbations in Minkowski background have been investigated and the effective energy-momentum tensor of the gravitational radiation is calculated [27]. The astrophysical gravitational waves of inspiralling compact binaries have been also investigated [28]. More recently, the electromagnetic and gravitational perturbations of black holes in conformal Weyl gravity have been studied [29], and also, the QNMs of these black holes in the nearly extreme regime have been obtained [30]. In this paper, we consider black hole solutions in conformal Weyl gravity and investigate their thermal stability and QNMs of scalar perturbations in asymptotically (A)dS spacetimes. We also investigate deviations from the Schwarzschild-(A)dS black holes due to an additional linear r-term in conformal Weyl solutions.
The outline of this paper is as follows. In the next section, we give a brief review of 4-dimensional black holes in conformal Weyl gravity. Then, thermodynamic properties of the Weyl solutions are investigated and the thermal stability analysis of the solutions is obtained. We also consider a minimally coupled massive scalar perturbation in the background spacetime of the black holes and calculate the related QNMs by using the sixth order WKB approximation and the asymptotic iteration method (AIM). These calculations are done in asymptotically dS spacetime, and the conditions of possible quasi-resonance modes and anomalous decay rate of QNMs are studied. Finally, the QNMs of the black holes in asymptotically AdS spacetime are obtained and some relations for calculating the QNMs of large black holes are found. For both dS and AdS cases, deviation of results from those of the Schwarzschild-(A)dS black holes is discussed. We finish our paper with some concluding remarks.

II. REVIEW OF SOLUTIONS
The four-dimensional conformal Weyl action reads where the Weyl conformal tensor is Variation of action (1) with respect to the metric tensor leads to the following equation of motion W µν = 2C µλνκ ;λκ − C µλνκ R λκ = 1 2 g µν (R α α ) ;β ;β + R µν;β ;β − R µβ;ν ;β − R νβ;µ ;β − 2R µβ R ν It was shown that the static spherically symmetric solution of Weyl gravity in four dimensions can be written as where dΩ 2 is the line element of a 2−sphere and the metric function is [14] f (r) = c + d r where b, c, and d are integration constants. It is clear that for non-vanishing b that plays the role of the cosmological constant, Eq. (5) is not a vacuum solution of Einstein gravity with or without the cosmological constant, while as long as b = 0, c = 1, and d = −2M , the metric becomes identical to the Schwarzschild solution of Einstein gravity. It is worth mentioning that although b plays the role of the cosmological constant, it is arisen purely as an integration constant and is not put in the action by hand. Such a constant cannot be added to the action of WG because it would introduce a length scale and hence breaks the conformal invariance. For future comparison with the Schwarzschild-(A)dS black holes, we set b = −Λ/3 and d = −2M , and consider the following metric function from now where for c = 1 reduces to the Schwarzschild-(A)dS solutions. Fig. 1 shows the behavior of the metric function (6) for different c. The singularity of Weyl solution can be covered with an event horizon, and therefore, we can interpret the solution as a black hole.

A. Thermodynamic quantities
Here, we calculate temperature and entropy of the Weyl black holes by using the surface gravity at the event horizon and the Wald entropy formula, respectively [31]. After that, we investigate thermal stability of the black holes.
By calculating the surface gravity, κ = − (∇ µ χ ν ) (∇ µ χ ν ) /2 (χ = ∂ t is the null Killing vector of the horizon), we can obtain the Hawking temperature of the black hole at the outermost (event) horizon, r + . The temperature of the black hole is simplified as in which M + can be obtained by f (r + ) = 0 as and we consider the lower sign that is compatible with the Schwarzschild-(A)dS black hole at c → 1 limit. Thus, the temperature converts to Besides, the entropy of the black hole in higher derivative theories can be obtained by Wald's formula [31,32] which makes the dependence of entropy on gravitational action where α is a constant with dimension (length) 2 , L is the Lagrangian density of the theory, ξ µν is the binormal to the (arbitrary) cross-section M of the horizon, and h is the determinant of induced metric on M. Therefore, the entropy of our case study black hole takes the following form Now, we use the first law of thermodynamics (δM tot = T δS) to calculate the total mass of the black holes as Notably, for general c (c = 1), one cannot obtain the expression for the finite mass (12) by using either the Deser-Tekin [33] or the AMD [34,35] methods (see [36] for more details). Nevertheless, it is worth mentioning that the validity of Eq. (12) is explicitly checked in the appendix D of Ref. [36] through the Noether approach.

B. Thermal stability
Here, we explore thermal stability and the possibility of a phase transition of the obtained black hole solutions. For an ordinary thermodynamical system, thermal stability criteria can be governed by the sign of heat capacity and compressibility. In other words, the following conditions [37] guarantee the stability, and thus, heat capacities (at constant pressure c p , or constant volume c v ) and compressibilities (isothermal κ T , or isentropic κ S ) must be positive in a stable system. Now, we investigate the stability conditions of constructed black hole solutions to find the thermally stable criteria. To do so, we should only investigate a heat capacity since this system just has one intensive parameter, T . The heat capacity of the solutions has the following explicit form where we set the constant α = 1 without loss of generality. The sign of heat capacity help us to find the stability information. The positive sign shows stable solutions whereas the negative sign indicates unstable ones. The heat capacity changes sign whenever it meets root or divergence points. The root of heat capacity indicates a bound point which separates the physical black holes (positive temperature) from unstable ones (negative temperature). In addition, divergence points may separate stable and unstable regions and may signal the existence of phase transition. Figures 2 and 3 show the behavior of the heat capacity with different c for small black holes (SBHs) and large black holes (LBHs). We fixed the cosmological constant since we are going to investigate the effect of c which is a new parameter compared with the Schwarzschild-(A)dS solutions. From these figures we find that the new parameter c affects the stability conditions. Besides, the behavior of the heat capacity is different in the cases of SBHs and LBHs for c 1 and c 1. One may note that the SBHs of Einstein-Λ gravity (Schwarzschild-(A)dS solutions) are unstable for all values of Λ. But in the case of Weyl solutions, there are stable SBHs for c 1 (c 1) and positive Λ (negative Λ). In contrast, the LBHs of Einstein-Λ gravity are stable everywhere. However, we have large stable black holes in Weyl gravity just for negative Λ and −1 c 1. Thus, the impact of new parameter c is significant regarding the stability conditions, especially for SBHs.
From this section, we found that the free parameter c changes the stability conditions of black holes in the Weyl gravity, significantly, compared with the Schwarzschild-(A)dS black holes, i.e., stable black holes become unstable and vice versa. Thus, we can expect to see interesting and notable effects of c on the QNMs of Weyl solutions, as we will show in the following sections.

IV. QUASINORMAL MODES IN DS SPACETIME
Scalar fields have been extensively investigated in the context of cosmology as inflatons [38], dark energy [39], and dark matter [40]. They can be considered as candidates for constructing quantum gravity [41,42], and also, can be considered as fields in the strong field regime to modify the black hole background [43,44]. Scalar fields can be regarded to form clouds through accretion or instabilities around black holes [45,46]. In this regard, one can consider non-minimal interaction between the scalar field and gravity which leads to black hole hair [43,47]. In different models with non-minimal interaction of scalar fields with the spacetime metric, the gravitational waves h µν of the metric will be a linear combination of gravitational waves in general relativity, h where Φ is the scalar field, g µν is the background metric, and β is the non-minimal coupling parameter. The QNMs h µν , that could potentially be observed, would have components oscillating with a combination of the general relativity and scalar QNMs. Therefore, the signature of scalar fields on gravitational waves could be observed by future gravitational wave detectors. However, a minimally coupled scalar field describes the QNMs in a range of scalar-tensor theories. In this paper, we focus on minimal coupled scalar fields to investigate the effects of the new parameter c on the scalar QNM spectrum in asymptotically dS spacetime and find deviations from those of the Schwarzschild-dS black holes.

A. Setup
Here, we consider a massive scalar perturbation in the background of the black hole spacetime and obtain the QNFs by using two independent methods of calculations; the sixth order WKB approximation [48][49][50] and the asymptotic iteration method (AIM) [51]. Furthermore, we concentrate our attention on the asymptotically dS black holes (Λ > 0) with the obtained metric function (6). The asymptotic flat solutions (c = 1 and Λ = 0) reduce to the Schwarzschild black holes which have been investigated extensively.
The equation of motion for a minimally coupled massive scalar field is given by in which µ is the mass of the scalar field Φ and = ∇ ν ∇ ν . Now, we consider the following expansion of modes where Y l,m (θ, ϕ) denotes the spherical harmonics on 2-sphere. Here and in the rest of this paper, we omit the integral over frequency ω in the Fourier transform for notational simplicity. Substituting the decomposition into Eq. (16), the equation of motion (16) reduces to a wavelike equation for the radial part Ψ l (r) in the following way In this equation, r * is the known tortoise coordinate with the definition and the effective potential V l (r * ) is given by where l is the angular quantum number and it is notable that r in the right-hand side is a function of r * by (19). Figure 4 shows the behavior of this effective potential (20) versus the tortoise coordinate for different values of c and µ.
The spectrum of QNMs for a perturbed black hole spacetime is the solution of the wave equation (18). However, we have to impose some proper boundary conditions to obtain its solutions. The quasinormal boundary conditions imply that the wave at the event (cosmological) horizon is purely incoming (outgoing) where r e and r c are, respectively, the radius of event horizon and the cosmological horizon. One should consider the mentioned boundary conditions in order to obtain the QNFs.

B. WKB approximation
The method is based on the matching of WKB expansion of the wave function Ψ l (r * ) at the event horizon and cosmological horizon with the Taylor expansion near the peak of the potential barrier through two turning points. Therefore, this method can be used for an effective potential that forms a potential barrier and takes constant values (or zero) at the event horizon (r * → −∞) and cosmological horizon (r * → ∞) (like Fig. 4). The WKB approximation was first applied to the problem of scattering around black holes [48], and then extended to the third order [49], 6th order [50] and recently to the 13th order [52]. The 6th order of WKB formula is given by where V 0 is the value of the effective potential at its local maximum, the correction terms Γ j 's correspond to the jth order and depend on the value of the effective potential and its derivatives at the local maximum, and n is the overtone number. The explicit form of the WKB corrections is given in [49] (for Γ 2 and Γ 3 ) and [50] (for Γ 4 , Γ 5 , and Γ 6 ). It is worthwhile to mention that the WKB approximation does not give reliable frequencies for n ≥ l. We use this formula up to the sixth order as a semi-analytical approach to obtain the QNFs of perturbations.

C. AIM
The AIM has been employed to solve the eigenvalue problems and solving second-order differential equations [53,54], and then it was shown that it is an accurate technique for calculating QNMs [30,51,55]. If one wants to employ the AIM, it is convenient to use the independent variable ξ = 1/r, and rewrite the wave equation (18) as where P , P ′ , and V l (ξ) are given by In order to choose the appropriate scaling behavior for quasinormal boundary conditions, one may define [55,56] in which κ j is the surface gravity at ξ j with f (ξ = ξ j ) = 0. This equation scale out the divergent behavior at some boundary ξ j and applies the boundary conditions (21) to the solution. Now, based on Eqs. (21) and (28), an appropriate choice for QNMs to scale out the divergent behaviour at the cosmological horizon is as follows so that the equation (23) converts to Besides, by considering the equations (28) and (29), the correct quasinormal condition at the event horizon, ξ e , is where U l (ξ) should be a finite and convergent function, and κ e is the surface gravity at the horizon By inserting (31) into (30), one can find the standard AIM form as follows so that λ 0 (ξ) and s 0 (ξ) are Once the standard AIM form of the master wave equation is obtained, we can differentiate it and apply the quantization condition to calculate the QN frequencies (see [51,55] for details of calculations). From table I, we find that although the WKB approximation does not give reliable frequencies for n ≥ l, for higher c, say c ≥ 1.15, the WKB formula gives better results for n = 0 = l. Also, most of the results of WKB approximation are in a good agreement with numeric results (tables II − IV ), and results get better for the higher multipole number, l > n, for all values of c as we expected. On the other hand, both the real and imaginary parts of the QNFs decrease with increasing c which shows that there are fewer oscillations for higher c at the ringdown stage and the modes decay faster for lower c. The tables I-IV also show deviations of the QNMs of Weyl solutions from those of the Schwarzschild-dS black holes. For lower c (c < 1), there are more oscillations and long life modes compared to higher values of c (c > 1). Thus, the QNMs spectrum of Weyl solutions deviates from those of the Schwarzschild-dS black holes and these deviations could potentially be observed by using future gravitational wave detectors for non-minimal interactions of scalar fields with the spacetime metric.
On the other hand, as one can see from Fig. 4 that the effective potential forms a potential barrier which is positive everywhere and vanishes at the event horizon and spatial infinity. This shows that we can find dynamically stable black holes undergoing massive scalar perturbations. It is worthwhile to mention that although the conformal Weyl solutions are dynamically stable for c > 1, they are thermally unstable (table I). In order to find black holes that are both thermally and dynamically stable for c > 1, one can consider larger black holes with r + > 15 (see the right panel of Fig. 2).     Besides, one of the motivations for considering the test massive fields comes from the fact that there are some QNMs with arbitrarily long life (purely real) modes called quasi-resonance modes [57]. For the quasi-resonance modes, the oscillations do not decay and the situation is similar to the standing waves on a string that is fixed at its both ends. The quasi-resonance modes occur for special values of the field mass and the QNMs disappear when the field mass takes higher values. However, this happens just for lower overtones whenever the effective potential is non-zero at least at one of the boundaries (the event horizon r * → −∞ or cosmological horizon r * → ∞). Now, let us investigate the possibility of the quasi-resonance modes for Weyl black hole solutions. The effective potential (20) vanishes at both infinities for all possible values of different parameters, and therefore, there are no quasi-resonance modes for (18). In addition, if one sets the integration constant Λ equals to zero, the effective potential still vanishes at both infinities. There is only one case so that the effective potential can be non-zero at spatial infinity and that is zero cosmological constant (Λ = 0) with c = 1. In this case, the effective potential reduces to the Schwarzschild black holes in flat spacetime which its quasi-resonance modes have been investigated before  [58]. Therefore, our black hole case study has no quasi-resonant oscillations in general and the imaginary part of the frequencies never vanishes (regardless of the mentioned trivial case). Figure 5 shows the behavior of QNFs with increasing in µ and confirms the above discussion. As µ increases, the real part of frequencies increases too, whereas the imaginary part first decreases rapidly and then takes a constant value.

E. Anomalous decay rate of QNMs in dS spacetime
It was recently shown that the decay timescales of the QNMs of a massive scalar field perturbations in the Schwarzschild background either grow or decay with increasing multipole number l, depending on whether the mass of the scalar field is small or large [59]. This anomalous behavior is due to an additional µ 2 -term in the sub-leading term in the eikonal expression for ω i . In this scenario, there is a critical massμ so that the imaginary part of the QNFs increases (decreases) with increasing in l for µ >μ (µ <μ). For low-l values, the value of the critical massμ decreases when l increases, but there is one fixed critical mass for large-l values.
Here, we investigate the possibility of this anomalous behavior for the Weyl solutions and Schwarzschild black holes in dS spacetime numerically. Fig. 6 shows the numerical results for ω i as a function of µ for different values of l and c. From the middle panel of this figure, we find that the curves cross over at a special mass, and thus the QNM spectrum of Schwarzschild black holes in dS spacetime contains this anomaly as for the flat case [59]. Also, one can see the same behavior for the Weyl solutions (left and right panels of Fig. 6). It is worthwhile to note that the free parameter c affects the value of the critical mass andμ increases with increasing in c.

V. QUASINORMAL MODES IN ADS SPACETIME
Here, we want to point out an application of the conformal Weyl black hole solutions in the context of AdS/CFT correspondence. The AdS/CFT correspondence describes a relation between the string theory on asymptotically AdS spacetime and conformal field theory on the boundary [60]. Some phenomena like the Nernst effect [61], Hall effect [62], superconductivity [63], and the decaying time scale of perturbations of a thermal state [17] in the field theory have dual descriptions in gravitational theory. Besides, this correspondence between a gravitational theory and a CFT can be extended to explain some aspects of nuclear physics [64].
In terms of the AdS/CFT correspondence, a large AdS black hole corresponds to an approximately thermal state in conformal field theory and scalar perturbations of the black hole correspond to perturbations of this thermal state. Therefore, the decay rate of the scalar field perturbations describes the decay of perturbations of this thermal state. In this scenario, we can calculate the QNMs of a large static black hole in asymptotically AdS spacetime to obtain the time scale of the thermal state to approach thermal equilibrium. Here, we shall obtain the QNMs of Weyl solutions with AdS asymptote to find the stability time scale of the corresponding thermal state. Investigating the dynamical stability of Weyl solutions undergoing scalar perturbations is another advantage of calculating the QNMs. We can follow either Horowitz-Hubeny approach [17] or pseudospectral method [65] to calculate the QNMs of asymptotically AdS black holes (Λ < 0). In this paper, we follow the pseudospectral method and use a public code given in [66] to obtain the QNMs. The pseudospectral method replaces the continuous variable by a discrete set of points and solves the resulting eigenvalue equation. Now, we consider the fluctuations of a massless scalar field in the background of Weyl solutions. It is convenient to obtain the master wave equation in Eddington-Finkelstein coordinates whenever we are going to use the pseudospectral method. In Eddington-Finkelstein coordinates, the background line element (6) is given by where u = 1/r and L is the AdS radius related to the cosmological constant by Λ = −3/L 2 . Thus, u = 0 represents the boundary and u = 1/r + corresponds to the event horizon. The equation of motion for a minimally coupled massless scalar field is governed by Eq. (16) with µ = 0. By considering (36) as the spacetime background and expanding the scalar field eigenfunction Φ versus the spherical harmonics as (17) and substituting the scalar field decomposition (17) into (16), we can find the following second-order differential equation for the radial part in which l is the multipole number and ω = ω r − iω i is the quasinormal frequency with an imaginary part ω i giving the damping of the perturbations and a real part ω r giving the actual oscillations. Thus, in terms of the AdS/CFT correspondence, τ = 1/ω i is the time scale that the thermal state should pass to meet the thermal equilibrium. On the other hand, the negativity of the imaginary part of the frequencies guarantees the dynamical stability of the black hole [17]. Otherwise, the perturbations increase as time increases and the spacetime becomes unstable. For the boundary conditions in AdS spacetime, causality requires ingoing modes at the event horizon and finite modes at the boundary that results in a discrete spectrum of QNFs ω. To analyze the behavior of modes Ψ (u) near the horizon and the boundary to see how to deal with the boundary conditions, without loss of generality, we consider r + = 1 and replace the value of M from the lower sign of Eq. (8). Therefore, u = 0 represents the spatial infinity and u = 1 corresponds to the event horizon radius. Starting with the horizon, by substituting an ansatz Ψ (u) = (1 − u) p in Eq. (38), one can obtain two solutions Ψ in (u) ∝ Const.
where λ is a constant with the following explicit form By considering the time dependence e −iωt of the modes, the Ψ out (t, u) behaves as In order to keep a constant phase, (1 − u) has to increase as t increases, and thus u should decrease which means that this solution is outgoing. Therefore, we must consider just the ingoing solution Ψ in (u) ∝ Const and discard the outgoing one. In addition, the ingoing mode is perfectly smooth near the horizon while as we approach the horizon, the outgoing mode oscillates more and more rapidly. Fortunately, the boundary condition of ingoing modes at the horizon is enforced automatically since we are working in the ingoing Eddington-Finkelstein coordinates.
On the other hand, there are two solutions near the boundary; a normalizable mode Ψ (u) ∝ u 3 and a nonnormalizable one Ψ (u) ∝ Const which we must discard. If we redefine Ψ (u) = u 2Ψ (u), then the normalizable mode tends to zero linearly whereas the non-normalizable mode diverges as ∼ u −2 . Doing this rescaling, the wave equation (38) converts to where in which M + is given in (8). Now, the normalizable solution behaves smoothly at the boundary and it should be considered, while we discard the other non-normalizable solution. The wave equation (41) is an input for the code and one can fix the free parameters L and c, and also, the event horizon radius r + or u + (presented in M + ) to calculate the QNMs. We recall that the large black holes correspond to the thermal states in CFT. Therefore, we shall focus on the QNMs of LBHs (r + >> L) and discard the small ones (r + << L). In addition, we set L = 1 as the AdS radius to compare the Weyl solutions in asymptotically AdS spacetime with the Schwarzschild-AdS black holes investigated in [17,67].
In table V , we list the QNFs for the fundamental mode (n = 0) and the first overtone (n = 1) of intermediate black holes (r + = 1, 10) and large ones (r + = 50, 100) for l = 0. From this table, one can see that as the overtone number and the event horizon radius increase, both the real and imaginary parts of frequencies increase as well. However, the free parameter c affects the frequencies differently based on the size of the event horizon radius. For intermediate black holes, the QNFs first increase, and then decrease when c increases. In the case of large black holes, the QNFs increase linearly with an increase in c. Thus, for black holes with r + = 100, the QNMs can be obtained by using the following equations and note that some similar linear equations can be found for different values of r + and overtone number. In addition, by considering table V , we find that the deviation between the QNMs of two consecutive values of c decreases when r + increases, and therefore, changing in c does not affect the QNMs significantly in the case of LBHs. Before studying this behavior and the reason for the linear relation between ω and c, let us, first, investigate some limitations on Eqs. (44). It is worthwhile to mention that the equation (8) puts two bounds on the lower and upper values of c. Therefore, c cannot take an arbitrary negative value to obtain zero ω r or ω i . As a result, there are bounds on the lower and higher values of the QNMs. For example, for black holes with r + = 100, c must obey −4641 c 64641 (the heat capacity for this range is illustrated in Fig. 7) We should mention that as the imaginary part of QNFs increases, the corresponds thermal state meets the stability faster. Therefore, the thermal state with c ≈ 64641 enjoys the fastest decay rate in its perturbations. In addition, the Weyl AdS solutions are dynamically stable under massless scalar perturbations since all the frequencies have a negative imaginary part, but they are thermally unstable in some areas (see Fig. 7).   Table V : The fundamental mode (first line) and the first overtone (second line) of the QNFs for l = 0 and different values of c and r + . The fundamental frequency and the first overtone of the Schwarzschild black holes (c = 1) agree with previous results [17,67] calculated by using Horowitz-Hubeny method.
In order to explain the linear relation between ω and c, and also, the weak effect of c on the QNMs of LBHs, we consider the temperature (9) for large black holes and then look at the relation between the QNMs and this temperature illustrated in Fig. 8. As one can see, both the real and imaginary parts of frequencies increase linearly with an increase in the temperature (46). On the other hand, for fixed r + , the temperature (46) is a linear function of c. Therefore, we can expect a linear relation between ω and c, as was obtained in (44) for r + = 100. More interestingly, since r + is present in the denominator of the second term in (46), changing in c does not affect the QNMs significantly when r + increases. Thus, the free parameter c has a weak effect on the QNMs of LBHs. This behavior was also found numerically from table V . The points in Fig. 8, representing the QNMs, lie on straight lines that their real part are given by ω r = 0.0261128 + 23.2396T, n = 0 ω r = 0.0427779 + 39.7233T, n = 1 for c = 0.5, ω r = 0.0829902 + 23.2372T, n = 0 ω r = 0.139923 + 39.7200T, n = 1 for c = 1.5, Re Ω , n 0 Im Ω , n 0 Re Ω , n 1 Im Ω , n 1 Re Ω , n 0 Im Ω , n 0 Re Ω , n 1 Im Ω , n 1 According to the AdS/CFT correspondence, τ = 1/ω i is the time scale of a thermal state to approach thermal equilibrium. Therefore, Eqs. (44), (49), and (50) give the value of τ and these equations are the main results of this section. In addition, the QNFs are linear functions of r + since the temperature of large black holes is a linear function of r + (46). Interestingly, some similar results were found for the Schwarzschild-AdS black hole [17] and AdS black holes in Born-Infeld massive gravity [68]. However, there is some differences between the Schwarzschild-AdS black hole and Weyl solutions; for the Schwarzschild case, the QNMs lie on straight lines through the origin [17] whereas this is not the case for Weyl black holes. Besides, the QNMs increase when the new parameter c increases and the effect of c on the QNMs decrease as r + increases.
As the final remark, we should mention that by considering the results of this section and previous results [17,30,68], the temperature has a significant role in determining the QNMs of large black holes in AdS spacetime and nearly extreme black holes in dS/flat spacetimes. In all cases, the QNM spectrum is a linear function of temperature (surface gravity for near-extremal black holes).

VI. CONCLUSIONS
We have studied the thermal stability of conformal Weyl black holes in the canonical ensemble. We observed that the presence of an additional linear r-term in the Weyl solutions modifies the stability conditions in comparison with the Schwarzschild-AdS black hole and this linear r-term changed the stability conditions significantly. Although the SBHs of Schwarzschild-(A)dS solutions are unconditionally unstable for all values of the cosmological constant, it was possible to find stable SBHs in the conformal Weyl gravity. In addition, we know that the LBHs of Schwarzschild-(A)dS solutions are always stable whereas we could find unstable LBHs in Weyl gravity. Furthermore, we have considered a minimally coupled massive scalar perturbation in the background spacetime of the conformal Weyl black holes and calculated the QNFs by using the sixth order WKB approximation and the AIM after 15 iterations. We have shown that although the WKB approximation does not give reliable frequencies for n ≥ l, this approximation gives better results for higher c when n = 0 = l. It was shown that most of the results of WKB approximation are in good agreement with AIM and results get better for the higher multipole number. Besides, we observed that there were more oscillations and the modes lived longer for lower c. We have shown that although it is possible to find dynamically stable black holes in the conformal Weyl gravity, they can be thermally unstable.
In addition, we argued that the black holes have no quasi-resonant oscillations even when one sets the integration constants Λ = 0 and c = 1. This happens due to the presence of the linear r-term. Therefore, the imaginary part of the frequencies never vanishes and there are always damping modes (QNFs are always complex). Moreover, there was a critical mass for the scalar field such that the decay timescales of the QNMs of the Schwarzschild-dS black holes and Weyl solutions either grow or decay with increasing multipole number l, depending on whether the mass of the scalar field was small or large. We have numerically shown that the free parameter c affects the value of the critical mass and the critical mass increases with increasing in c.
Also, the QNMs of the Weyl solutions in asymptotically AdS spacetime were obtained and deviations from those of the Schwarzschild-adS black holes investigated. The QNFs were calculated by using the pseudospectral method to investigate the dynamical stability of the black holes, the effects of the free parameter c on the QNMs, and obtain the time scale of the thermal state to approach thermal equilibrium in conformal field theory. It was seen that the Weyl solutions are dynamically stable. Besides, it was shown that the free parameter c affects the QNMs of the Weyl solutions differently based on the size of the black hole. As c increased, the QNFs of intermediate black holes first increased and then decreased. However, the QNMs of LBHs increased linearly with an increase of c. In addition, we have found that the QNMs are a linear function of the free parameter c and the event horizon radius (or temperature). We also have shown that as the size of the black hole increases, the effect of c on the QNMs decreases, and thus the effect of c on the QNMs is negligible in the case of very large black holes. Since a static large black hole in AdS spacetime corresponds to an approximately thermal state in CFT, the free parameter c has no effect on the time scale of the thermal state in the case of very large black holes. As a result, two thermal states correspond to extremely large black holes in Einstein and Weyl gravities are identical.