Black Holes in Starobinsky-Bel-Robinson Gravity and the Breakdown of Quasinormal Modes/Null Geodesics Correspondence

1. Introduction

Quasinormal modes (QNMs) of black holes are proper oscillation frequencies which depend on the parameters of the black hole and not on the way of perturbations [1,2,3,4]. The data from gravitational waves together with observations in the electromagnetic spectra [5,6,7] in principle should allow 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 is large, what allows for big room for alternative theories of gravity [8,9].

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 [10] and was conceived by Alexei A. Starobinsky in [11]. 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+\alpha R^2$ Starobinsky inflation allowing for the leading Bel–Robinson-tensor-squared term. These corrections are inspired by the gravitational effective action of superstrings/M-theory compactified to four dimensions [12,13,14,15]. 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 [14], while cosmological applications were considered in [15,16,17]. Various effects around such black holes have been recently considered in [18,19]

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 [20,21,22,23,24,25,26,27], because usually to prove the stability analytically is difficult, and the thorough analysis of quasinormal modes allows 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 [28,29,30].

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 $\beta$ is small enough. We find the critical value of the parameter $\beta$ at the threshold of instability and analyze the quasinormal spectrum in the stable sector.

The paper is organized as follows. In sec. II we briefly describe the black hole spacetime and wave-like equation. Sec. III summarize 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.

2. The metric

Using the compactification coupled with the presence of stringy fluxes, the appropriate $4D$ gravity models may be described by the following action:

$$S=\frac{M_{pl}}{2}\int d^{4}x\sqrt{-g}\left(R+\frac{R^2}{6m^2}-\frac{\beta }{32m^6}\left(P_4^2-E_4^2\right)\right).$$

(1)

Here g is the metric determinant, R is the Ricci scalar. m is a free mass parameter, $\beta >0$ is a dimensionless coupling describing the compactification of M-theory, $P_4^2$ and $E_4^2$ are the Pontryagin and Euler topological densities which are related to the Bel-Robinson tensor $T_{\mu \nu \lambda \rho }$ in four dimensions as follows:

$$T^{\mu \nu \delta \rho }{T}_{\mu \nu \delta \rho }=\frac{1}{4}\left(P_4^2-E_4^2\right).$$

(2)

The Bel-Robinson tensor has the following form:

$$T^{\mu \nu \delta \sigma }=R^{\mu \rho \gamma \delta }{R}_{\rho \gamma }^{\nu \sigma }+R^{\mu \rho \gamma \sigma }{R}_{\rho \gamma }^{\nu \delta }-\frac{1}{2}{g}^{\mu \nu }{R}^{\rho \gamma \alpha \delta }{R}_{\rho \gamma \alpha }^{\sigma }.$$

(3)

It can be seen from Eq.(1) that the gravity action depends only on the two parameters, m and $\beta$. This allows for various applications such as Hawking radiation, entropy, inflation, optical phenomena etc. [12,13,14,15,19] The black hole solution in SBR gravity depends on the value of $\beta$ corrected to first order perturbations.

The line element of the spherically symmetric metric is

$$d{s}^2=-f\left(r\right)d{t}^2+\frac{1}{f\left(r\right)}d{r}^2+r^{2}d{\mathsf{\Omega }}^2.$$

(4)

where the metric function has the form [14]:

$$f\left(r\right)=1-\frac{2M}{r}+\beta M^6{\left(\frac{8\sqrt{2}\pi M}{r^3}\right)}^3\left(\frac{108r-194M}{5r}\right).$$

(5)

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 [19].

Coefficient ${\left(8\sqrt{2}\pi \right)}^3\approx 44901.9$ which means that $\beta \sim 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 $\beta$ only and include consideration of relatively large deformations of the Schwarzschild spacetime which already cannot be described as a small correction.

3. The methods

3.1. Wave-like equation and boundary conditions

We will consider the test neutral massless scalar field $\mathsf{\Phi }$ which obeys the general covariant Klein-Gordon equation,

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }\left(g^{\mu \nu }\sqrt{-g}{\partial }_{\nu }\mathsf{\Phi }\right)=0.$$

(6)

After the separation of variable,

$$\mathsf{\Phi }=\frac{{\psi }_{\ell }\left(r\right)}{r}{Y}_{\ell m}(\theta ,\varphi )e^{-i\omega t},$$

(7)

the equation for the scalar field can be reduced to the wave-like form:

$$\frac{d^2{\psi }_l}{d{r}_{*}^2}+\left[{\omega }^2-U\left(x\right)\right]{\psi }_{\ell }=0.$$

(8)

The effective potential is

$$U\left(r\right)=f\left(r\right)\left(\frac{\ell (\ell +1)}{r^2}+\frac{f^{\prime }\left(r\right)}{r}\right),$$

(9)

where the tortoise coordinate x can be expressed as follows:

$${\displaystyle \frac{d{r}_{*}}{dr}}=\frac{1}{f\left(r\right)}.$$

(10)

The effective potential is shown on fig Figure 1. There one can see that the when $\beta$ is increasing, the additional gap appears, which at sufficiently large values of $\beta$ 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.

3.2. 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 [31]. The first order WKB formula was afterwords generalized to higher orders up to the 13-th order [32,33,34]. Recently it shas been further improved by using the Padé approximants [34], 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 [35]:

$$\begin{array}{ccc}\hfill {\omega }^2& =& U_0+A_2\left({\mathcal{K}}^2\right)+A_4\left({\mathcal{K}}^2\right)+A_6\left({\mathcal{K}}^2\right)+\dots \hfill \\ \hfill & -& i\mathcal{K}\sqrt{-2V_2}\left(1+A_3\left({\mathcal{K}}^2\right)+A_5\left({\mathcal{K}}^2\right)+A_7\left({\mathcal{K}}^2\right)\dots \right),\hfill \end{array}$$

(11)

where for the quasinormal frequencies we have

$$\mathcal{K}=n+\frac{1}{2},n=0,1,2,\dots$$

(12)

Here n is the overtone number, $U_0$ is the value of the effective potential in its maximum, $U_2$ is the value of the second derivative of the potential in this point, and $A_2,A_3,A_4,\dots$ depend on higher derivatives of the potential in its maximum. We will use the sixth order WKB method obtained in [33] with Padé approximants $\tilde{m}=4$ and 5, where the definition of $\tilde{m}$ and further details on this method can be found in [34,35]. The higher order WKB method has been used for finding quasinormal modes and grey-body factors in numerous works, showing very good concordance with the result obtained by other methods (see, for example [36,37,38,39] and references therein).

3.3. Time-domain integration

The integration in the time domain that we employ as the main method is founded on the wavelike equation expressed using the light-cone variables $u=t-r_{*}$ and $v=t+r_{*}$. We implement the discretization approach introduced by Gundlach-Price-Pullin [40] as follows:

$$\mathsf{\Psi }\left(N\right)=\mathsf{\Psi }\left(W\right)+\mathsf{\Psi }\left(E\right)-\mathsf{\Psi }\left(S\right)-$$

$${\mathsf{\Delta }}^{2}V\left(S\right)\frac{\mathsf{\Psi }\left(W\right)+\mathsf{\Psi }\left(E\right)}{4}+\mathcal{O}\left({\mathsf{\Delta }}^4\right),$$

(13)

Here, the notations for the points are the following: $N\equiv \left(u+\mathsf{\Delta },v+\mathsf{\Delta }\right)$, $W\equiv \left(u+\mathsf{\Delta },v\right)$, $E\equiv \left(u,v+\mathsf{\Delta }\right)$, and $S\equiv \left(u,v\right)$. We impose the initial Gaussian wave package on the null surfaces $u=u_0$ and $v=v_0$. The quasinormal frequencies can then be extracted from the time-domain profiles reasonably accurately using the Prony method, as detailed in, for instance, [3]. Although the time-domain integration method is indeed an accurate and convergent technique for generating time-domain profiles, the extraction of frequencies with high precision can be challenging. Generally, achieving a convergence with the exact results from the Leaver method, with a discrepancy of less than one percent, is attainable for the first and higher multipoles [41,42]. 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 $d{r}^{*}=dr/f\left(r\right)$ and find $r^{*}\left(r\right)$ 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.

3.4. 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 [43,44]. 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=\infty$. Using the new function of the radial coordinate $P\left(r\right)$,

$$\mathsf{\Psi }\left(r\right)=P(r,\omega ){\left(1-\frac{r_g}{r}\right)}^{-i\omega /F^{\prime }\left(r_g\right)}y\left(r\right),$$

(14)

we then choose the factor P such that it provides regularity if $y\left(r\right)$ in the range $r_0\le r<\infty$, for the outgoing wave at infinity and ingoing wave at the horizon. Therefore, we can write down the Frobenius series in the following way:

$$y\left(r\right)=\sum _{k=0}^{\infty }{a}_k{\left(1-\frac{r_g}{r}\right)}^k.$$

(15)

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 [45], which was generalized to an arbitrary number of terms recurrence relations by A. Zhidenko in [46]. 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.

4. The threshold of instability

From the form of the effective potentials it is evident that once $\beta$ is increased, the negative gap outside the black hole becomes deeper (Figure 1), so that one should check whether the scalar field is stable. Time-domain integration at various values of $\beta$ and $\ell =0$ shows that, indeed, the instability occurs at some critical value:

$${\beta }_{cr}\approx 0.1364.$$

(16)

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 [25,47], so that higher multipoles must be investigated for possible (in)stability as well. The modes with $\ell =1,2$ and higher multipoles (see, for example, table II) are stable at least for those values of $\beta$ for which the stability of the $\ell =0$ perturbations is guaranteed. The unstable mode has non-oscillatory exponential form as can be seen from Figure 2. This non-oscillatory character of the unstable mode was rigorously proved in [48] 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 $\beta$-corrections is safe as to the possible development of instability.

5. Quasinormal frequencies in the stable sector

Table 1. Quasinormal modes obtained by the time-domain integration and subsequent extraction of frequencies with the Prony method for various values of $\beta$, $\ell =0$, $M=1$. The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when $\beta$ grows.

$\mathit{\beta }$	Time-domain	WKB
${10}^{-6}$	0.110022 - 0.105870 i	0.109460 - 0.103417 i
$5·{10}^{-6}$	0.110119- 0.105857 i	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 $\beta$, $\ell =0$, $M=1$. The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when $\beta$ grows.

$\mathit{\beta }$	Time-domain	WKB
${10}^{-6}$	0.110022 - 0.105870 i	0.109460 - 0.103417 i
$5·{10}^{-6}$	0.110119- 0.105857 i	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 2. Quasinormal modes obtained by the time-domain integration and subsequent extraction of frequencies with the Prony method for various values of $\beta$, $\ell =1$, $M=1$. The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when $\beta$ grows.

$\mathit{\beta }$	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.097447 i	0.293569 - 0.097658 i
$2·{10}^{-5}$	0.293430- 0.097203 i	0.294184 - 0.097428 i
$3·{10}^{-5}$	0.293667- 0.096956 i	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 $\beta$, $\ell =1$, $M=1$. The higher order WKB data with Padé approximants are evidently diverge from the correct time-domain integration results when $\beta$ grows.

$\mathit{\beta }$	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.097447 i	0.293569 - 0.097658 i
$2·{10}^{-5}$	0.293430- 0.097203 i	0.294184 - 0.097428 i
$3·{10}^{-5}$	0.293667- 0.096956 i	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

When the multipole number ℓ is large, one can find quasinormal frequencies in the analytic form. It is convenient to use $\kappa =\ell +1/2$ instead of ℓ. We find the location of the peak of the effective potential expanded in powers of $\kappa$ and, afterwards, to use the higher order WKB technique as prescribed in [49] in order to obtain an analytic expression for the quasinormal modes in the eikonal approximation and at the first order beyond it in the regime of small $\beta$ and $1/\ell$:

$$r_{max}=3M-\frac{M}{3{\kappa }^2}+\beta \left(\frac{2502656\sqrt{2}{\pi }^3}{295245M^5{\kappa }^2}-\frac{210944\sqrt{2}{\pi }^3}{32805M^5}\right)+O\left(\frac{1}{{\kappa }^4},{\beta }^2\right)$$

(17)

Finally, applying the above expansion for $r_{max}$ and using higher order WKB formula we obtain

$$\omega =\frac{\kappa }{3\sqrt{3}M}-\frac{i\mathcal{K}}{3\sqrt{3}M}+\frac{29-60{\mathcal{K}}^2}{1296\sqrt{3}M\kappa }+\left(\frac{13312\sqrt{\frac{2}{3}}{\pi }^3\kappa }{59049M^7}+\frac{115712i\sqrt{\frac{2}{3}}{\pi }^3\mathcal{K}}{59049M^7}+\frac{64\sqrt{\frac{2}{3}}{\pi }^3\left(91860{\mathcal{K}}^2-11783\right)}{7971615M^7\kappa }\right)\beta +O\left(\frac{1}{{\kappa }^2},{\beta }^2\right).$$

(18)

Neglecting orders $\sim 1/\kappa$ in the above equation, we obtain the eikonal formula which is exact in the regime $\ell \to \infty$. The eikonal formula is interesting, because, as was stated in [50], 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 $\ell \gg n$ limit. Thus, real and imaginary parts of the $\ell \gg n$ quasinormal mode are proportional to the frequency and instability timescale of the circular null geodesics, respectively:

$${\omega }_n=\mathsf{\Omega }\ell -i(n+1/2)\left|\lambda \right|,\ell \gg n.$$

(19)

Here $\mathsf{\Omega }$ is the angular velocity at the unstable null geodesics, and $\lambda$ 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 [51], it was demonstrated that this correspondence falters when the wave-like equation’s dominant eikonal centrifugal term deviates from $f\left(r\right)\ell (\ell +1)/r^2$. This divergence occurs, for instance, in the context of gravitational perturbations within the Einstein-Gauss-Bonnet [51,52] or Einstein-dilaton-Gauss-Bonnet theories [53]. Broadly speaking, the aforementioned correspondence holds true when the first-order WKB formula [31] 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 $\mathsf{\Lambda }$-term [54,55], and de Sitter modes that are modified by the presence of a black hole [56,57]. Consequently, as elucidated in [58], 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 Figure 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 time-domain 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 $\beta$ 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 [59,60].

Lower modes are shown in tables I and II, where we can see that the WKB approach cannot be trusted unless the parameter $\beta$ 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 $\ell =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.

6. Conclusions

This work is the first analysis of quasinormal modes of the black hole in the Starobinsky-Bel Robinson gravity suggested recently in [14]. It proved out that the black hole possesses a rich quasinormal spectrum with the following features:

• At unphysically large values of the parameter of compactification $\beta$ it produces instability, while the regime of small $\beta$ is free from instability
• In the eikonal limit $\ell \to \infty$ the correspondence breaks down, because it reproduces only part of the eikonal spectrum, but, unlike [58], this breakdown is for an asymptotically flat space.
• 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.