Black holes in Starobinsky-Bel-Robinson Gravity and the breakdown of quasinormal modes/null geodesics correspondence

We show that perturbations of a scalar field in the background of the black hole obtained with the Starobinsky-Bel-Robinson Gravity is unstable unless the dimensionless coupling $\beta$ describing the compactification of M-theory is small enough. In the sector of stability quasinormal spectrum show peculiar behavior both in the frequency and time domains: the ringing consists of two stages where two different modes dominate. The WKB method does not reproduce part of the spectrum including the fundamental mode, which is responsible for the first stage of the ringing. As a result, the correspondence between the high frequency quasinormal modes and characteristics of the null geodesics reproduces only one branch of the eikonal spectrum. The frequencies are obtained with the help of three methods (Frobenius, WKB and time-domain integration) with excellent agreement among them.


Introduction
Quasinormal modes (QNMs) of black holes are their fundamental characteristics, proper frequencies of oscillations which depend on the parameters of the black hole and not on the way in which a black was perturbed [27,51,36,5].The data from gravitational waves together with observations in the electromagnetic spectra [3,1,16] in principle should allow one in the future to determine the black hole's mass, angular momentum, and charge and thereby to test gravitational theory in the strong field limit.Nevertheless, by now the uncertainty in the mass and angular momentum of a black hole is large, what allows for big room for alternative theories of gravity [28,60].
Among various alternative/modified theories of gravity the Starobinsky inflation model plays an important role.It uses the idea that the early universe went through an inflationary de Sitter era [55] and was conceived by Alexei A. Starobinsky in [56].The model incorporates modifications to General Relativity by introducing a quadratic curvature term.This framework allows for an inflationary phase driven by a scalar field, addressing cosmological issues like the horizon and flatness problems.It generates a scale-invariant density perturbation spectrum, consistent with cosmic microwave background and large-scale structure observations, making it a significant inflationary model in modern cosmology.
Recently a new gravitational theory in four dimensions has been proposed in the form of a sum of the modified R + αR 2  Starobinsky inflation allowing for the leading Bel-Robinsontensor-squared term.These corrections are inspired by the gravitational effective action of superstrings/M-theory compactified to four dimensions [20,24,7,25].The Starobinsky-Bel-Robinson Lagrangian possesses only two free parameters, which makes it attractive for studying various physical phenomena which could potentially be observed and, thereby, the theory could be tested.Consequently, the corrected Schwarzschild-like spacetime was perturbatively obtained in [7], while cosmological applications were considered in [25,15,53].Various effects around such black holes have been recently considered in [4,2] Here we propose the first study of the quasinormal frequencies in the background of black holes in the Starobinsky-Bel-Robinson gravity.We will consider a test scalar field, which is the simplest case frequently sharing the properties of other spin fields.Apparently the most important aspect of the quasinormal spectrum is its (in)stability [59,14,19,26,9,57,13,63], because usually it is difficult to prove the stability analytically, and the thorough analysis of quasinormal modes allows one to judge about (in)stability.The onset of instability in the scalar sector may indicate the possible instability in the gravitational one.The neutral test scalar field in the Schwarzschild, Kerr and Kerr-Newman spacetimes are known to be stable under the quasinormal mode boundary conditions, because the corresponding quasinormal frequencies are damped [52,35,37].
Here we will consider the evolution of perturbations of the scalar field in the black hole background in the Starobinsky-Bel-Robinson gravity.We will show that the scalar field is unstable unless the higher curvature correction parameter β is small enough.We find the critical value of the parameter β at the threshold of instability and analyze the quasinormal spectrum in the stable sector.The latter has distinctive features: the ringdown consists of the two stages at which two different modes dominate.One of them cannot be reproduced by the WKB method what undermines the correspondence between eikonal quasinormal modes and null goedesics.
The paper is organized as follows.In sec.II we briefly describe the black hole spacetime and wave-like equation.Sec.III summarizes the methods used for the analysis of the quasinormal spectrum, while in sec.IV we present the obtained results for the instability threshold and quasinormal modes in the range of parameters corresponding to the stability.In the Conclusions we review the obtained results and mention some points for future investigation.

The metric
Using the compactification coupled with the presence of stringy fluxes, the appropriate 4D gravity models may be described by the following action: Here g is the metric determinant, R is the Ricci scalar; m is a free mass parameter, β > 0 is a dimensionless coupling describing the compactification of M-theory, P 2 4 and E 2 4 are the Pontryagin and Euler topological densities which are related to the Bel-Robinson tensor T µνλρ in four dimensions as follows: The Bel-Robinson tensor has the following form: It can be seen from Eq.(1) that the gravity action depends only on the two parameters, m and β.This allows for various applications such as Hawking radiation, entropy, inflation, optical phenomena etc. [20,24,7,25,2] The black hole solution in SBR gravity depends on the value of β corrected to first order perturbations.
The line element of the spherically symmetric metric is where the metric function has the form [7]: Here r s = 2M is the Schwarzschild radius and M is the mass parameter.The tidal force effects and shadows for this metric were considered in [2].Coefficient (8 √ 2π) 3 ≈ 44901.9 which means that β ∼ 3•10 −5 (at M = 1) must already be considered as large deformation of the Schwarzschild spacetime.As we are interested here also in the testing of the null geodesics/eikonal quasinormal modes correspondence correspondence for generic black-hole spacetimes, we will not be limited by the regime of tiny β only and include consideration of relatively large deformations of the Schwarzschild spacetime which already cannot be described as a small correction.

Wave-like equation and boundary conditions
We will consider the test neutral massless scalar field Φ which obeys the general covariant Klein-Gordon equation, The separation of variables could be performed using the stationary ansatz, which transforms the Klein-Gordon equation to the wave-like form: The effective potential is where the tortoise coordinate x can be expressed as follows: The effective potential is shown on fig 1.There one can see that the when β is increasing, the additional gap appears, which at sufficiently large values of β becomes negative.This is an indication that a bound state with negative energy may be formed, which means the onset of instability.Therefore, thorough analysis of quasinormal modes is necessary.The quasinormal modes are eigenvalues of the above second order differential equation under the requirements of purely ingoing waves at the event horizon and purely outgoing waves at infinity.

Higher order WKB method
For an additional check of accurate results obtained by the Frobenius method in the range of stability we will also use the WKB method.It was applied for the first time to the problem of finding of quasinormal modes and grey-body factors of black holes by B. Schutz and C. Will [54].The first order WKB formula was afterwords generalized to higher orders up to the 13-th order [21,29,49].Recently it has been further improved by using the Padé approximants [49], so that the relative error of the WKB formula with Padé approximants is usually from quite few times to a few orders smaller than without them.The WKB formula can be expressed in the following general form [41]: where for the quasinormal frequencies we have  Here n is the overtone number, U 0 is the value of the potential's maximum, U 2 is the value of the second derivative of the potential in this point, and A 2 , A 3 , A 4 , . . .depend on higher derivatives of the potential in its maximum.In the present work we apply the sixth order WKB method of [29], but with Padé approximants m = 4, as prescribed in [49,41].This kind of WKB approach is frequently and effectively used for calculations of quasinormal frequencies and greybody factors and usually shows good agreement with the precise methods (see, for example [48,18,47,30]).

The method of integration in time-domain
The integration in the time domain that we employ as the main method here is founded on the wavelike equation expressed in terms of the light-cone variables u = t − r * and v = t + r * .We implement the discretization approach introduced by Gundlach, Price and Pullin [17]: Here, the notations for the points are the following: and S ≡ (u, v).We assume that an initial Gaussian wave package on the surfaces u = u 0 and v = v 0 model the perturbation.The obtained quasinormal frequencies are known to be independent on the parameters of this package, such as its height and location of the the center.The quasinormal frequencies can then be found from the time-domain profiles with not very high but guaranteed accuracy using the representation of the signal as a sum of exponents which is the Prony method, as detailed, for instance, in [36].As a rule one can find quasinormal modes with the relative error of less than one percent at the first or higher multipoles [6,44], but not always at the ℓ = 0 case, because the time of the ringing is very short for that case.The challenging point of the method is the construction of the effective potential as a function of the tortoise coordinate r * which can be done in the following way.One can numerically integrate the equation dr * = dr/ f (r) and find r * (r) then using its connection with the coordinates u, v to built U(r * ).An important feature of the time-domain integration which we will use here is that it includes contributions of all overtones and, is, thereby, efficient for detecting the instability.

Frobenius method
The Frobenius or Leaver method is based on the convergent procedure and is therefore allows us to find quasinormal frequencies with great accuracy [45,46].As it is also broadly discussed in the literature, here we will only briefly mention the key points of the method.The differential master wave-like equation has a regular singular point at the event horizon r = r g and the irregular singular point at r = ∞.Using the new function of the radial coordinate P(r), we then choose the factor P such that it provides regularity if y(r) in the range r 0 ≤ r < ∞, for the outgoing wave at infinity and ingoing wave at the horizon.Therefore, we can write down the Frobenius series in the following way: Then, the Guassian eliminations reduce the problem to an algebraic equation.In addition, for further increasing of the convergence we will use the Nollert improvement [50], which was generalized to an arbitrary number of terms recurrence relations by A. Zhidenko in [62].The Frobenius method, however, requires an initial guess for the value of quasinormal frequency, when finding the roots of the corresponding algebraic equations and if the corresponding zeros are highly localized in the parametric space, it might be difficult to find such an initial guess for every mode, so that there is a risk of missing some modes.Therefore the WKB and time-domain integration methods are useful for getting such an initial guess.

The threshold of instability
From the form of the effective potentials it is evident that once β is increased, the negative gap outside the black hole becomes deeper (fig.1), so that one should check whether the scalar field is stable.Time-domain integration at various values of β and ℓ = 0 shows that, indeed, the instability occurs at some critical value: Although usually perturbations at higher multipole numbers are more stable, because the effective potential get higher, there are examples of instability developing at high ℓ as well [57,38], so that higher multipoles must be investigated for possible (in)stability as well.The modes with ℓ = 1, 2 and higher multipoles (see, for example, table II) are stable at least for those values of β for which the stability of the ℓ = 0 perturbations is guaranteed.The unstable mode has non-oscillatory exponential form as can be seen from fig. 2. This non-oscillatory character of the unstable mode was rigorously proved in [32] for spherically symmetric black holes.
Because of the large numerical pre-factor, the threshold of instability corresponds to a very large deformation of the Schwarzschild spacetime, so that, obviously, the regime of small β-corrections is safe as to the possible development of instability.

Quasinormal frequencies in the stable sector
When the multipole number ℓ is large, one can find quasinormal frequencies in the analytic form.It is convenient to use κ = ℓ + 1/2 instead of ℓ.We find the location of the peak of the effective potential expanded in powers of κ and, afterwards, to use the higher order WKB technique as prescribed in [40] in order to obtain an analytic expression for the quasinormal modes in the eikonal approximation and at the first order beyond it in β Time-domain WKB 10 −6 0.110022 -0.105870 i 0.109460 -0.103417 i 5 • 10 −6 0.110119-0.105857i 0.109426 -0.103701 i 10 −5 0.110244 -0.105843 i 0.109436 -0.103135 i 3 • 10 −5 0.110766 -0.105798 i 0.083291 -0.094837 i 5 • 10 −5 0.111322 -0.105778 i 0.091824 -0.129149 i 7 • 10 −5 0.111892 -0.105784 i 0.086486 -0.138424 i 10 −4 0.112733 -0.105896 i 0.032663 -0.143389 i 10 −3 0.121337 -0.115261 i -Table 1: Quasinormal modes obtained by the time-domain integration and subsequent extraction of frequencies with the Prony method for various values of β, ℓ = 0, M = 1.The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when β grows.
β Time-domain WKB 10 −6 0.293003 -0.097665 i 0.292990 -0.097690 i 5 • 10 −6 0.293090 -0.097568 i 0.293246 -0.097705 i 10 −5 0.293202-0.097447i 0.293569 -0.097658 i 2 • 10 −5 0.293430-0.097203i 0.294184 -0.097428 i 3 • 10 −5 0.293667-0.096956i 0.294743 -0.097104 i 4 • 10 −5 0.293912 -0.096706 i 0.295266 -0.096741 i 5 • 10 −5 0.294166 -0.096453 i 0.295776 -0.096360 i 10 −4 0.295543 -0.095136 i 0.299002 -0.094337 i 3 • 10 −4 0.304555 -0.091160 i 0.279144 -0.059728 i 5 • 10 −4 0.315643 -0.093030 i 0.278248 -0.138408 i 7 • 10 −4 0.323746 -0.099017 i 0.006171 -0.042030 i 10 −3 0.330159 -0.109229 i 0.082300 -0.106786 i Table 2: Quasinormal modes obtained by the time-domain integration and subsequent extraction of frequencies with the Prony method for various values of β, ℓ = 1, M = 1.The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when β grows.the regime of small β and 1/ℓ: Finally, applying the above expansion for r max and using higher order WKB formula we obtain Neglecting orders ∼ 1/κ in the above equation, we obtain the eikonal formula which is exact in the regime ℓ → ∞.The eikonal formula is interesting, because, as was stated in [10], parameters connected to unstable circular null geodesics around a spherically symmetric, and asymptotically flat or de Sitter black hole are dual to the quasinormal modes in the ℓ ≫ n limit.Thus, real and imaginary parts of the ℓ ≫ n quasinormal mode are proportional to the frequency and instability timescale of the circular null geodesics, respectively: Here Ω is the angular velocity at the unstable null geodesics, and λ is the Lyapunov exponent.An intricate aspect of this correspondence lies in its applicability, which extends to the Schwarzschild black hole and several other scenarios.However, in [33], it was demonstrated that this correspondence falters when the wave-like equation's dominant eikonal centrifugal term deviates from f (r)ℓ(ℓ+1)/r 2 .This divergence occurs, for instance, in the context of gravitational perturbations within the Einstein-Gauss-Bonnet [33,42] or Einstein-dilaton-Gauss-Bonnet theories [43].Broadly speaking, the aforementioned correspondence holds true when the first-order WKB formula [54] is applicable, signifying the presence of an effective potential with a singular maximum and two turning points.
Furthermore, even within these scenarios, the WKB formula might prove inadequate in reproducing the complete spectrum of the eikonal regime.This limitation becomes evident in the Schwarzschild-de Sitter spacetime, where two categories of modes emerge: Schwarzschild modes that are influenced by the Λ-term [61,34], and de Sitter modes that are modified by the presence of a black hole [39,8].Consequently, as elucidated in [31], although the correspondence suggested by equation (19) holds formally, the parameter n no longer accurately represents the overtone number.
Here we observe, in a sense, similar situation: as can be seen from fig. 3, there are two stages of the ringing on which two different modes dominate.The mode at the first stage is slower decaying and, formally, is the fundamental one, while the first overtone dominates at the second longer stage.However, the fundamental mode which is reproduced with the timedomain integration and Frobenius method with very good concordance between them, cannot be found by the WKB method which gives only the first overtone.Thus, while the correspondence works for the part of the eikonal spectrum, it does not allow one to reproduce the whole spectrum in the eikonal regime.It is essential that the effective potential at such small values of β appears as a WKB-good single peak barrier.It is worth mentioning that consequently the break-down of the null geodesics/eikonal quasinormal modes correspondence for our case must also violate the correspondence between the quasinormal modes and radius of the shadow considered in [23,22].One should also remember that null geodesics not always determine the trajectories for photon movement.Thus in the nonlinear electrodynamics, photons do not move along the null geodesics [11,12,58].
Lower modes are shown in tables I and II, where we can see that the WKB approach cannot be trusted unless the parameter β is tiny.At the same time, the time-domain integration always provides relatively good accuracy for the fundamental mode with relative error less than one percent at least for ℓ = 1 when there are sufficient number of oscillations in the profile.Nevertheless, it is difficult to achieve great accuracy with the time-domain integration, because of arbitrariness of the fitting of the profile with a sum of exponents within the Prony method.

Conclusions
This work is the first analysis of quasinormal modes of the black hole in the Starobinsky-Bel Robinson gravity suggested recently in [7].It proved out that the black hole possesses a rich quasinormal spectrum with the following features: 1.At unphysically large values of the parameter of compactification β it produces instability, while the regime of small β is free from instability 2. In the eikonal limit ℓ → ∞ the correspondence breaks down, because it reproduces only part of the eikonal spectrum, but, unlike [31], this breakdown is for an asymptotically flat space.
3. The ringing consists of the two stages at each of which different modes dominate.This could be well seen in the regime of large ℓ.
Our work could be extended to fields of other spin and considerations of higher overtones, which could be done with the help of the Frobenius method.