Features and stability analysis of non-Schwarzschild black hole in quadratic gravity

Black holes are found to exist in gravitational theories with the presence of quadratic curvature terms and behave differently from the Schwarzschild solution. We present an exhaustive analysis for determining the quasinormal modes of a test scalar field propagating in a new class of black hole backgrounds in the case of pure Einstein-Weyl gravity. Our result shows that the field decay of quasinormal modes in such a non-Schwarzschild black hole behaves similarly to the Schwarzschild one, but the decay slope becomes much smoother due to the appearance of the Weyl tensor square in the background theory. We also analyze the frequencies of the quasinormal modes in order to characterize the properties of new back holes, and thus, if these modes can be the source of gravitational waves, the underlying theories may be testable in future gravitational wave experiments. We briefly comment on the issue of quantum (in)stability in this theory at linear order.


Introduction
Since Einstein's proposal in 1915, General Relativity (GR) has been established as the standard theory of gravitation for one hundred years. As a pillar of modern science, the predictions of GR have been confirmed in all observations to date. However, one of the most challenging task that theoretical physicists are facing today is how GR can be reconciled with the laws of quantum physics to produce a consistent ultraviolet (UV) complete theory of quantum gravity.
To address this issue as well as to be of phenomenological interest, extensions of Einstein gravity with the presence of higher order derivative terms arise in fundamental theories, such as string theory, loop quantum gravity, asymptotically safe gravity and others. In particular, it was found in [1] that adding quadratic curvature terms could improve the renormalizability of the underlying gravitational theory although this theory would suffer from an instability of ghost degrees of freedom [2]. As shown in [3], however, if the path integral quantization is evaluated in Euclidean space and then Wick rotated to Lorentzian space, this path integral can yield a theory of quantum gravity without a negative norm state. This approach has been applied into inflationary cosmology and provided an interesting interpretation for primordial perturbations [4].
In addition, black hole solutions were obtained in a gravitational theory involving higher order terms within the scenario of asymptotical safety [32] and its stability issue was addressed in [33] by analyzing the so-called quasinormal modes. This study provides a representative example to show how the (in)stability issue of a black hole solution could be investigated via the method of analyzing quasinormal modes, of which the generation is due to the quasinormal ringing of the background spacetime under perturbations and hence is associated with the characteristics of black hole. The identification of the quasinormal modes is considered to possibly falsify various black hole solutions derived in a large class of gravity theories through the imminent gravitational wave surveys. Note that the investigation of black hole perturbations has drawn a lot of interest for decades since it is associated with black hole stability, gravitational wave detection as well as some fundamental symmetries such as the gauge/gravity duality. Analyses of these modes in theories of higher derivative gravity were performed in [34][35][36][37]. We refer to Refs. [38,39] and references therein for recent comprehensive reviews.
In the present paper we aim at examining the quasinormal modes seeded by a test scalar field propagating in black hole solutions of a pure Einstein-Weyl gravity as derived in ref. [16]. In Section 2 we briefly review the background gravitational theory and the associated black hole solutions. Then in Section 3 we perform a detailed analysis of quasinormal modes seeded by a massless test scalar field that is propagating freely within various black hole solutions beyond Schwarzschild and investigate their behaviors. Section 4 is devoted to a brief discussion of the field equations for gravitational waves, in which one obtains two poles on the dispersion relation and hence this implies an instability of quantum fluctuations. We eventually summarize our results with a discussion in Section 5. Throughout the whole paper we use geometrized units with G = c = 1 and the (−, +, +, +) convention for the metric.

Quadratic gravity and black holes beyond Schwarzschild
We start with a brief introduction to a general theory of quadratic gravity. Consider a four-dimensional Einstein gravity involving quadratic curvature terms of which the action is written as (following ref. [16]) is defined as the Weyl tensor in four dimension with R being the Ricci scalar. In the above action we have introduced two model parameters α and β, which are of area dimension: , and, describe the deviations from the regular Einstein gravity.

Spherically symmetric and static solutions
In order to find a black hole solution beyond Schwarzschild, we consider a static and spherically symmetric ansatz as follows, in which two dimensionless metric factors N and F have been introduced as functions of the radial coordinate. As was shown in [14] as well as argued in [16], the Ricci scalar vanishes for any static black hole solutions of the action (2.1). As a result, the general theory of quadratic gravity can reduce to the pure Einstein-Weyl gravity at the classical level by setting β = 0. We refer to the appendix of the present paper for the details of studying a black hole solution beyond Einstein analytically. In the main context, we simply summarize the steps of constructing such a solution as follows. Firstly, one can vary the action with respect to the metric factors and then derive the field equations for N and F . Secondly, in order to exhibit the difference between this solution and the Schwarzschild one, one can parameterize the leading term of the F factor (denoted by F 1 introduced in (A.7) in the appendix) as follows, where r H represents for the position of the black hole horizon and δ is the amount of deviation from the Schwarzschild solution since in GR we have δ = 0. Such a parametrization can provide a boundary condition for numerically solving the vacuum structure of the underlying gravity theory, which is the last step to implement in the whole construction. In the following subsection we numerically repeat the result of ref. [16] to demonstrate an existence of a black hole beyond Schwarzschild.

Numerical estimates
Note that, as has been observed in [16], for each given α, the viable value of r H is bounded. Specifically, if r H is too small there is no opportunity to form a non-Schwarzschild black hole; however, if r H is too large the black hole mass would become negative and hence leads to a quantum instability at high energy scales. For example, by setting α = 0.5, one numerically derives a bound: 0.876 < r H < 1.143. For any selected value of r H in the above bounded interval, there exists only one value of δ that allows for a healthy non-Schwarzschild black hole. This phenomenon is related to the fine-tuning of Γ + = 0 in the . Note that, as in the treatment of ref. [16], we impose the normalization factor to be 0.6 for N and unity for F at infinity, in order to avoid an asymptotic overlap.  weak field limit (see the second part of the appendix). In figure 1 we show an exact example of the numerical construction introduced in the present subsection. Afterwards, using the numerical fitting one can also relate the exact numerical solution to the weak field limit approximate solution in (A.10) as developed in the appendix. This shows that physically the weak field linearized theory can roughly describe the Einstein-Weyl black hole solution in the large radius regime. A specific fitting is provided in figure 2. One can read from the figure that the evolution of F and N in the weak field limit is in agreement with the exact numerical results at large radii but deviate badly when evolving to the regime near the horizon.

Quasinormal modes
In this section we study quasinormal perturbations of Einstein-Weyl black holes. Quasinormal modes have been shown to be very useful to uncover intrinsic properties of the geometry (for instance, see [40]). Consider a massless test scalar field ψ propagating in a black hole background governed by the Einstein-Weyl theory. Its evolution obeys the massless Klein-Gordon equation Plugging in the black hole metric into the above Klein-Gordon equation, one gets, In order to solve the above Klein-Gordon equation analytically, it is convenient to use the following standard separation of variables by making use of the spherical harmonic functions. Accordingly, the Klein-Gordon equation can be greatly simplified as: for each fixed value of l. By introducing the generalized tortoise coordinate we can obtain the Schrödinger-type equation for each value of l as follows, where the effective potential is given by This form can be treated systematically in the analysis of quasinormal modes, which will be given in the following subsections. We will apply the methods of the characteristic integration and the WKB analysis for quantitative estimation of numerical quasinormal modes, respectively, in the following up subsections.

Characteristic integration
A simple but efficient way of solving 1+1 dimensional d'Alembert equations is established in the pioneering work [41]. In this formalism, the standard (r * , t) coordinates are replaced by the light-cone variables, in terms of which all wave equations can have the same form.
Considering the fact that we do not have full analytic solutions for the time-evolving wave equation, one efficient approach is to discretize the mode function as and h is the discrete step size. Note that, in order to solve the mode function for a fixed l, one needs to impose the initial condition at the null boundary u = u 0 and v = v 0 . As will be confirmed by the following numerical simulations, however, the characteristics of the associated field decay are basically insensitive to the initial conditions imposed. From now on we would like to drop the subscript l from the mode function Φ but specify its value in detailed calculations. Through the difference equation mentioned above, one can get a series of time domain data Φ(t 0 ), Φ(t 0 + h), Φ(t 0 + 2h), etc, for all possible fixed r * , and the quasinormal vibrations can be read-off through the transformation from the time domain to the frequency domain.
The numerical computation of the evolution for a mode function |Φ| is shown in figure 3 on a logarithmic scale. In this plot one can see a representative field decay evolution in the time domain. We choose α = 0.7, which ensures that the conformal term contribution is smaller than the minimal coupling in the action, and r H = 1.1, where exists a numerical non-Schwarzschild solution far from the negative mass region. One can see that the scalar field evolution is firstly dominated by the quasinormal vibration, and then decays with a power-law tail. This is a standard scenario in the time domain profile of black holes in analogue with the case of the Schwarzschild solution (e.g. see the review [39]).
One can also investigate the field decays with different background parameters as shown in figure 4. The left panel of figure 4 shows the field decays by varying the value of the multipole number from 0 to 2. In this panel, one can read that for larger values of the multipole number l, the longer the power-law tail and the larger the slope can be in a given time region. This behavior is also similar to the Schwarzschild case. However, from the right panel of the figure, one can explicitly read that the slope of the field decay strongly depends on the newly introduced model parameter α. In the limit α → 0, it is expected that one can recover the field decay of a Schwarzschild black hole. Moreover, a larger value of α leads to a smoother slope for the field decay for a fixed value of l as shown in the lower panel of figure 4.
We recall that there exists a bound for the model parameter α that allows for a non-Schwarzschild black hole solution as has been addressed in section 2.2. In the theory of Einstein-Weyl gravity, this bound is stronger for a smaller value of α. Moreover, this bound interval would move along the negative direction of the real axis with α decreasing, so we cannot adjust r H as a fixed parameter. However, we find that the influence of changing r H (in the regular bound interval) on the properties of the field decay is not so dramatic as the variation of α. Also, with a smaller value of α the field decay evolution has a larger slope with a longer tail, which is similar to the tendency when decreasing l. The reason could be understood as follows. The value of α estimates the deviation from the standard Einstein gravity and thus from a standard Schwarzschild black hole, which is the highest symmetric solution (satisfying F = N and δ = 0). On the other hand, the multipole number l is also associated with the symmetry of a dynamical system, which is similar to the standard Hydrogen atom problem in quantum mechanics textbooks. Thus, the effects on the field decay by decreasing α and by decreasing l are similar from this perspective.

WKB analysis
In order to understand the quasinormal perturbations (semi-)analytically, it is useful to perform the WKB analysis of the mode function. Imposing s = iω on the Schrödinger-like equation (3.6), one gets which is very common in usual cases, one can get a discrete set of possible values for s(ω).
The WKB semi-analytic approach [42,43] is a very successful and efficient method to calculate quasinormal frequencies. In the present study we expand the computation up to the third order, and correspondingly, the square of those frequencies can be written as with ξ = n + 1/2. In addition, the superscript (i) denotes the i-th order differentiation with respect to r * of the potential V (r(r * )). The subscript 0 means that the potential and its derivatives are calculated at the point r * 0 , where V (r(r * )) is an extremum. The solution of ω can be determined when Thus we can use the formula (3.12)-(3.14) to find the quasinormal modes of the non-Schwarzschild black holes. The accuracy of the WKB method is sensitive to specific black hole solutions [44]. However, as has been observed in [36,45], for low overtones (l > n) the accuracy becomes better. In order to compare explicitly the method of the characteristic integration and the WKB analysis, we perform detailed analyses of two methods by varying the value of α from   Table 1. Values for the quasinormal frequencies for the mode function propagating in the non-Schwarzschild geometry based on the WKB method and the algorithm of characteristic integration, respectively. The value of model parameter α varies from 0.1 to 1.0 smoothly. The integer value of the multipole number varies from l = 0 to 2, respectively.
0.1 to 1.0 and the multipole number from l = 0 to 2, respectively. Our results are presented in table 1. From this table, we find that the results of the WKB analysis and the characteristic integration results are fairly consistent but not with high accuracy. Regarding this issue, we argue that there exists a limit on the accuracy of both two methods. For the method of the characteristic integration, the error by solving the differential equations could lead the wave function to deviate from the original form of the differential equations and the corresponding boundary conditions could be affected as well. To overcome this numerical deviation, one needs to improve the method of computer algorithm, which is a very detailed technical issue. Moreover, the accuracy of the results obtained from the WKB approach mainly rely on to which order one truncates the computation. Consider that the WKB calculations at higher order would become extremely lengthy and do not change the results qualitatively, we would like to simply take the third order truncation in our detailed analysis.
Some generic features can be concluded here. First, in the regular bounded interval of r H , all quasinormal mode frequencies have a negative imaginary part, which shows that scalar perturbations in the non-Schwarzschild black hole backgrounds is stable in the time evolution 1 . Second, it is observed that with increasing α, l and decreasing r H , the absolute value of the imaginary part of the frequency decreases. This observation is consistent with our previous argument on the tendency of the slope because the imaginary part of the frequency stands for the slope of the logarithmic time domain decay. Third, the real part of the frequency stands for the trigonometric vibration of the scalar field. This real part dramatically increases along with a larger value of the multipole number l.

Field equations for gravitational waves
In this section we briefly discuss tensor perturbations around a non-Schwarzschild black hole solution. Here we use Greek letters µ, ν, ... to denote the coordinates on the four-dimensional spacetime, Latin letters i, j, k... to denote the coordinates on the two-dimensional space submanifold S 2 (namely, (θ, φ) coordinates), and r, t to denote the radial and time coordinate respectively. And also, we use a comma to denote an ordinary derivative, while a semicolon denotes a covariant derivative.
Let us consider linear perturbations about the background metric, where the transverse and traceless tensor satisfies The first order perturbation theory will provide simple results for the perturbation of the Riemann tensor, Ricci tensor and Einstein tensor for spherically symmetric and static backgrounds [46]. We simply summarize them as follows. By writing the components of the Ricci and Riemann tensors as where G = − N (−2 + 2F + rF ,r ) + rF N ,r 2r 2 N , we can derive the perturbations of these geometric tensors Here the LHS has indexes (µ, ν), while the RHS may only have (i, j). This convention means that the related tensors are nontrivial only when (µ, ν) are on the submanifold S 2 .
As a result, one gets δR = 0. After a very lengthy computation, one can derive the field equation for tensor fluctuations at leading order, which is given by where we have introduced the Schouten tensor It is interesting to notice that, there exist the 2 operator which appears in the first term of the above field equation and cannot be canceled by any constraint equations. This implies that at high energy the tensor fluctuations would obtain an extra degree of freedom. Such a new degree of freedom often corresponds to a ghost mode that would spoil the stability of the vacuum state quantum mechanically, such as in the Lee-Wick theory of particle physics.
In this regard, the Einstein-Weyl theory of gravitation may still suffer from the quantum instability issue even though this theory is classically stable since the classical scalar perturbations can well behave as analyzed in the previous section. This instability was recently also addressed in ref. [47], where the authors applied the Stückelberg approach to show that the interplay between the ghost graviton and the healthy graviton allows the theory to evade the usual strong coupling issue widely existing in massive gravity theories and become renormalizable, at the expense of stability.

Conclusion
Recently, the theory of quadratic gravity has drawn the interest of theoretical physicists in the literature from various aspects [48][49][50][51]. In particular, it was found in [16] that the Einstein-Weyl gravity can allow for a black hole solution beyond Schwarzschild. In the present work we have revisited this type of new black hole solutions at both the background and perturbative levels. At the background level, we analyzed the solutions of the metric factors in the weak field limit and compared them with the exact numerical results. At the perturbative level, we have studied in detail the propagation of quasinormal perturbations seeded by a test scalar field within such a black hole background.
Specifically, we have analyzed the frequencies and time domain evolutions of quasinormal modes seeded by this massless scalar field in the exterior of such a non-Schwarzschild black hole as derived in the Einstein-Weyl theory of gravity. Our results show that the time domain evolution of the quasinormal modes is similar to that obtained in the regular Schwarzschild case where the mode functions decay in a power-law form. However, due to the existence of the higher derivative term, the slope of the field decay is generally smoother than that in the Schwarzschild one. In addition, we present a brief analysis of tensor perturbations, which are regarded as gravitational waves. We show explicitly that the linearized field equation for gravitational waves involves higher derivative operators that would bring the theory to be unstable quantum mechanically at high energy scales. However, it would be interesting to study in more depth this issue under non-perturbative approaches in order to reveal the relation between a ghost graviton mode and quantum renormalizability of gravity theories.

A Appendix
In the first part of this Appendix, we provide a detailed instruction to the background theory of quadratic gravity. In the second part, we show how a black hole solution beyond Schwarzschild can be obtained in this theory.

A.1 The theory of quadratic gravity
Varying the action (2.1) with respect to the metric yields the background field equation as follows, where we have introduced which is from the conformal gravity part, and also, the Bach tensor due to the Weyl-Eddington term.
In order to study the vacuum structure of a static spacetime satisfying spherical symmetry, we assume the absence of matter fields in the above system. As was shown in [14] as well as argued in [16], the Ricci scalar vanishes for any static black hole solutions of the action (2.1). That is, R = 0 in the above classical theory. It turns out that, at the classical level, one can greatly simplify the quadratic gravity action (2.1) by taking β = 0, and therefore, the generic action reduces to a theory of pure Einstein-Weyl gravity.
As was pointed out in [16], however, the requirement of R = 0 does not simply lead to the trivial solution of the Schwarzschild black hole. This is because, by setting R = 0 and integrating the trace of the field equation (A.1) over the spatial region could yield a nontrivial and non-vanishing Ricci tensor of the four-dimensional spacetime, although the surface term remains zero. This is also the key reason that there might exist static and spherically symmetric black holes over and above the Schwarzschild one as shown in [16].

A.2 Black holes beyond Schwarzschild
Plugging the spherically symmetric and static ansatz (2.3) into the background field equation (A.1) leads to the following two second order differential equations: where the subscript ,r represents a derivative with respect to r. To fully determine the solutions of these two differential equations, one also needs to impose the horizon condition that with r H being defined as the position of the black hole horizon.

A.2.1 Near horizon limit
Since the above two metric factors vanish at the horizon, in order to grasp the physics of the black hole solution near the horizon, it is convenient to make Taylor expansions as follows, Note that, among all the coefficients of the Taylor expansions, there exists at least one parameter that is a normalization factor and accordingly can be absorbed by a time rescaling. In the following we take N 1 to be the normalization factor without loss of generality. One can plug these expansions into the field equations (A.4) and solve for F i and N i order by order. In addition, in the parametrization (2.4) we have introduced δ to characterize the difference between this solution and the Schwarzschild one.

A.2.2 The linearized treatment in weak field limit
From the other side, it is well known that the metric factors should approach unity when far away from the black hole in order to be consistent with the boundary condition of Minkowski spacetime. Therefore, one can analyze these metric factors in the weak field limit at large scales. This method has been widely applied in the literature and turned out to be very useful in analyzing black hole systems in modified gravity theories, for instance, in massive gravity models [52,53]. In this limit, we can expand the metric factors around a Minkowski background as F (r) = 1 + f (r) , N (r) = 1 + n(r) , (A. 8) in the limit where r is much larger than the horizon. Keeping leading order terms in n and f , the field equations can then be greatly simplified, of which the forms are given by, r 2 f (r) − 4αf (r) + r 3 n ,r + 2r 2 αf ,rr = 0 , 2f (r) + 2rf ,r + 2rn ,r + r 2 n ,rr = 0 . (A.9) Consequently, the metric factors in the weak field limit can be approximately solved as n(r) = Γ 0 r + Γ − e −mr r + Γ + e +mr r , (A.10) f (r) = Γ 0 r + Γ − (1 + mr)e −mr 2r + Γ + (1 − mr)e +mr 2r , with m 2 ≡ 1/2α being introduced. Moreover, the coefficients Γ 0 , Γ ± are integral constants that can be determined by numerically evolving the metric factors from the horizon to large length scales. From the approximate solution in the weak field limit, one can immediately notice that the appearance of the e mr /r term would severely spoil the classical stability of the black hole system governed by the Einstein-Weyl theory. One possible way of circumventing this issue is to finely tune the value of δ introduced in the parametrization (2.4) to ensure Γ + = 0. Under this condition one gets n(r) = Γ 0 r + Γ − e −mr r , f (r) = Γ 0 r + Γ − (1 + mr)e −mr 2r . (A.11) This limit also shows that we cannot choose a negative α for a regular black hole solution, in which case m will be imaginary and cause a vibration at large radii for f (r). Since for a specific numerical solution as will be constructed in the next section one can numerically fits the values of Γ 0 and Γ − by matching the linearized solution with the exact one.