Quasibound States, Stability and Wave Functions of the Test Fields in the Consistent 4D Einstein–Gauss–Bonnet Gravity

Abstract

We examine the interaction between quantum test particles and the gravitational field generated by a black hole solution that was recently obtained in the consistent 4-dimensional Einstein–Gauss–Bonnet gravity. While quasinormal modes of scalar, electromagnetic, and Dirac fields have been recently studied in this theory, there is no such study for the quasibound states. Here, we calculate the spectrum of quasibound states for the test fields in a spherically symmetric and asymptotically flat black hole solution in the consistent 4-dimensional Einstein–Gauss–Bonnet gravity. The quasispectrum of resonant frequencies is obtained by using the polynomial condition associated to the general Heun functions. We also discuss the stability of the systems for some values of the Gauss-Bonnet coupling constant.

1. Introduction

The quasibound states, quasinormal modes, and shadows of black holes are among the most interesting characteristic of such an astrophysical objects in the observational (measurable) spectra. In addition, some physical phenomena have been observed in modern experiments involving condensate matter and optical systems. In order to interpret these data, it is used the general theory of relativity. However, since some fundamental questions, as for example, the issues related to the quantum gravity phenomenon [1,2,3,4,5], cannot be solved with the Einstein’s theory, some alternative theories of gravity have been proposed. Among these alternative approaches, we can mention the f(R), the Lovelock, and the Einstein–Gauss–Bonnet theories of gravity, where the last two deals with higher curvature corrections [6,7].

The Einstein–Gauss–Bonnet (EGB) theory is one of the most promising approaches developed to deal with the higher curvature corrections that appears in the standard Einstein’s theory (see Ref. [8] and references therein). The EGB theory is quadratic in the curvature and leads to non-trivial corrections of the equation of motion when the Gauss-Bonnet (GB) term is coupled to a matter field, which can be, for instance, a dilaton. On the other hand, Aoki–Gorji–Mukohyama (AGM) [9] developed the so-called consistent theory of 4-dimensional (4D) EGB gravity, where they used the Arnowitt–Deser–Misner (ADM) decomposition [10] to construct the Hamiltonian and then this theory has not infinite coupling.

In the present paper, we will use a black hole metric that is an exact solution of the (well-defined truly) 4DAGM theory. However, it is worth emphasizing that this black hole solution is also a (particular) result of the dimensional regularization suggested in the (controversial) work by Glavan and Lin [11]. Henceforth, the black hole metric under consideration can be safely used, since it is considered as a consistent solution [12]. Here, we will investigate the test fields, namely, the scalar, electromagnetic, and Dirac fields, which means that we will solve the Klein-Gordon, Maxwell, and Dirac equations in the background under consideration.

In this sense, we can also mention some works on the regularized/consistent 4DEGB gravity, as follows. Motivated by the recent results on searching for the metric solutions of the 4DEGB gravity in the spherically symmetric spacetime, as well as in the Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime, it was proved that, in a cylindrically symmetric spacetime, the novel 4DEGB gravity always has the same solutions as the ones in its regularized counterpart. Thus, these solutions satisfy the equations of motion in the novel 4DEGB gravity. In these works, the nontrivial effect of the GB coupling constant on the energy scale of the second order phase transitions (which give rise to cosmic strings) is also shown [13,14,15,16,17,18]. In addition, the innermost stable circular orbit and shadow of the 4DEGB black hole were studied by Guo and Li [19], where was pointed out that a negative GB coupling constant is allowed to retain a black hole in 4DEGB gravity and gave the allowed range.

In this paper, we calculate the spectrum of quasibound states, and their corresponding radial and angular wave eigenfunctions, for scalar, electromagnetic, and Dirac particles in an asymptotically flat 4-dimensional Einstein–Gauss–Bonnet black hole (4DEGBBH) by using the polynomial condition of the general Heun function. We show that the quasibound states depend on the GB coupling constant, a, and that all the test fields constitute stable systems, when $0<a<1/2$. In addition, we also compute both radial and angular wave eigenfunctions of the test fields in the 4DEGBBH spacetime.

The paper is organized as follows. In Section 2, we introduce the metric corresponding to the 4DEGBBH spacetime. In Section 3, we solve the master wave equations in the background under consideration. Section 4 is devoted to the quasibound states of the test fields. In Section 5, we provide both radial and angular wave eigenfunctions, by using some properties of the general Heun functions. Finally, in Section 6, we summarize the obtained results. Here we adopt the natural units where $G\equiv c\equiv \hslash \equiv 1$.

2. The Consistent 4-Dimensional Einstein–Gauss–Bonnet Black Hole Spacetime

A crucial aspect for our studies of black hole radiation (be it an emission, transmission and/or reflection of quasinormal modes and/or quasibound states) is such that the black hole solution should be a solution of the truly 4DAGM theory [9], as well as of the theories with extra scalar degrees of freedom [20,21,22,23]. Thus, we will use the description given by Churilova (see Ref. [12] and references therein), in which all of the aforementioned approaches are taken into account to construct the novel consistent 4DEGB theory, and its black hole solution as well. In what follows, we briefly review the basic ideas behind this approach.

Let us consider a 4-dimensional spacetime in which is valid the following gauge condition [24]

$${}^3\mathcal{G}=\sqrt{\gamma }{D}_k{D}^k({\pi }^{ij}{\gamma }_{ij}/\sqrt{\gamma })\approx 0,$$

(1)

where $D_k$ is the covariant derivative and ${\pi }^{ij}$ is the canonical momentum conjugate to ${\gamma }_{ij}$. Then, it can be shown that the unique gravitational action is given by

$$S_{\mathrm{AGM}}=\int dt{d}^{3}xN\sqrt{\gamma }{\mathcal{L}}_{\mathrm{EGB}}^{4\mathrm{D}},$$

(2)

where N is a lapse function. The 4DEGB Lagrangian ${\mathcal{L}}_{\mathrm{EGB}}^{4\mathrm{D}}$ is defined as

$${\mathcal{L}}_{\mathrm{EGB}}^{4\mathrm{D}}=\frac{M_{\mathrm{Pl}}^2}{2}\left\{2R-\mathcal{M}+\frac{\tilde{a}}{2}\left[8R^2-4R\mathcal{M}-{\mathcal{M}}^2-\frac{8}{3}(8R_{ij}{R}^{ij}-4R_{ij}{\mathcal{M}}^{ij}-{\mathcal{M}}_{ij}{\mathcal{M}}^{ij})\right]\right\},$$

(3)

where $M_{\mathrm{Pl}}$ is the reduced Planck mass characterizing the gravitational coupling strength, $\tilde{a}$ is the GB coupling constant, $R_{ij}$ is the Ricci tensor and R is the Ricci scalar. The ${\mathcal{M}}_{ij}$ tensor is given by

$${\mathcal{M}}_{ij}=R_{ij}+{\mathcal{K}}_k^k{\mathcal{K}}_{ij}-{\mathcal{K}}_{ik}{\mathcal{K}}_j^k,$$

(4)

with

$${\mathcal{K}}_{ij}=({\dot{\gamma }}_{ij}-2D_{(i}{N}_{j)}-{\gamma }_{ij}{D}^2{\lambda }_{\mathrm{GF}})/2N,$$

(5)

where the dot denotes the derivative with respect to the time t, and ${\lambda }_{\mathrm{GF}}$ is a gauge-fixing parameter. It worth noticing that this consistent 4DEGB theory has a spatial diffeomorphism invariance and a time reparametrization symmetry given by $t\to t=t\left(t^{\prime }\right)$. Finally, by doing an appropriate rescaling of the GB coupling constant, an exact solution describing the 4DEGBBH spacetime has the following form

$$d{s}^2=-f\left(r\right)d{t}^2+\frac{1}{f\left(r\right)}d{r}^2+r^{2}d{\theta }^2+r^2{sin}^2\theta d{\varphi }^2,$$

(6)

where the metric function, $f\left(r\right)$, is given by

$$f_{\pm }\left(r\right)=1+\frac{r^2}{a}\left(1\pm \sqrt{1+\frac{4aM}{r^3}}\right).$$

(7)

Note that we have chosen the gauge-fixing parameter as ${\lambda }_{\mathrm{GF}}=0$. The parameter M is the total mass centered at the origin of the system of coordinates.

The metric function $f_{+}\left(r\right)$ corresponds to an asymptotically de Sitter spacetime, while $f_{-}\left(r\right)$ corresponds to an asymptotically flat spacetime. Here, we will focus on the “minus” case and hence we will set $f\left(r\right)\equiv f_{-}\left(r\right)$. In this case, there are two solutions when $a>0$; otherwise, for $a<0$ (and $M>0$) the metric function is not real for sufficiently small values of r, that is, for $r^3<-4aM$. Here, we will focus on the positive values of the GB coupling constant and hence the surface equation can be parametrized as

$$f\left(r\right)=0=(r-r_{+})(r-r_{-}).$$

(8)

Thus, the two solutions of this surface Equation (8) are the exterior, $r_{+}$, and interior, $r_{-}$, event horizons, and given by

$$r_{\pm }=M\pm \sqrt{M^2-\frac{a}{2}}.$$

(9)

For simplicity, and without loss of generality, from now on we set $M=1/2$, which means that the GB coupling constant is such that $0<a<1/2$; that is the stability region. The behavior of the metric function $f\left(r\right)$, as well as the event horizons, is shown in Figure 1. Thus, the 4DEGBBH metric is well-behaved outside the exterior event horizon, and then we will obtain some results on quasibound states that are valid for these values of a.

In what follows, we will solve the equations of motion for the test fields and then discuss the quasibound states, as well as the radial and angular wave eigenfunctions, in the background under consideration.

3. Master Wave Equations

The general covariant equations of motion for the scalar, $\Phi$, electromagnetic, $F_{\mu \nu }$, and Dirac, ${\rm Y}$, massless test fields can be written as

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }(g^{\mu \nu }\sqrt{-g}{\partial }_{\nu }\Phi )=0,$$

(10)

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

(11)

$${\gamma }^{\alpha }(\frac{\partial }{\partial x^{\alpha }}-{\Gamma }_{\alpha }){\rm Y}=0,$$

(12)

where $F^{\mu \nu }(={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu })$ is the electromagnetic field tensor, $A_{\mu }$ is the 4-vector electromagnetic potential, ${\gamma }^{\alpha }$ are the noncommutative gamma matrices, and ${\Gamma }_{\alpha }$ are the spin connections in the tetrad formalism.

In order to solve these master equations, first we need to separate their angular and radial parts. To do this, we will follow the approaches described by Konoplya et al. [25], Churilova [26], and Aragón et al. [27] to deal with the Maxwell, Klein-Gordon, and Dirac equations, respectively. These approaches have the advantage of avoiding square roots of the metric function $f\left(r\right)$, so that it will be possible to find analytical solutions for the radial parts. However, these approaches involve too much algebra, so let us summarize it as follows.

For the scalar field, we will use the subscript “_B” (Bosonic, for simplicity but without loss of generality) and set the following ansatz for the wave function: $\Phi =u_{\mathrm{B}}\left(r\right)P_{\nu }^m(cos\theta )e^{im\varphi }{e}^{-i{\omega }_{\mathrm{B}}t}$. Here, $u_{\mathrm{B}}\left(r\right)=R_{\mathrm{B}}\left(r\right)/r$ is the scalar radial function, $P_{\nu }^m(cos\theta )$ are the associated Legendre functions, ${\omega }_{\mathrm{B}}$ is the scalar frequency, $m\le 0$ ($\in \mathbb{Z}$) is the magnetic quantum number, and $\nu$ ($\in \mathbb{C}$) is the degree (or azimuthal quantum number). Thus, the bosonic radial function $R_{\mathrm{B}}\left(r\right)$ satisfies the following equation,

$$\frac{d^2{R}_{\mathrm{B}}\left(r\right)}{d{r}^2}+\frac{1}{f\left(r\right)}\frac{df\left(r\right)}{dr}\frac{d{R}_{\mathrm{B}}\left(r\right)}{dr}+\frac{1}{f^2\left(r\right)}\left\{{\omega }_{\mathrm{B}}^2-f\left(r\right)\left[\frac{\lambda }{r^2}+\frac{1}{r}\frac{df\left(r\right)}{dr}\right]\right\}{R}_{\mathrm{B}}\left(r\right)=0,$$

(13)

where $\lambda$ is the separation constant.

For the Maxwell field, we will use the subscript “_E” (Electromagnetic, for simplicity but without loss of generality) and expand the four-vector electromagnetic potential $A_{\mu }$ in terms of the vector associated Legendre functions, as well as assume the time dependence as $e^{-i{\omega }_{\mathrm{E}}t}$, where ${\omega }_{\mathrm{E}}$ is the electromagnetic frequency. Here, $u_{\mathrm{E}}\left(r\right)=R_{\mathrm{E}}\left(r\right)/r$ is the electromagnetic radial function. Thus, the electromagnetic radial function $R_{\mathrm{E}}\left(r\right)$ satisfies the following equation,

$$\frac{d^2{R}_{\mathrm{E}}\left(r\right)}{d{r}^2}+\frac{1}{f\left(r\right)}\frac{df\left(r\right)}{dr}\frac{d{R}_{\mathrm{E}}\left(r\right)}{dr}+\frac{1}{f^2\left(r\right)}\left\{{\omega }_{\mathrm{E}}^2-f\left(r\right)\left[\frac{\lambda }{r^2}\right]\right\}{R}_{\mathrm{E}}\left(r\right)=0.$$

(14)

For the Dirac field, we will use the subscript “_F” (Fermionic, for simplicity but without loss of generality) and set the following ansatz for the wave function: ${\rm Y}=u_{\pm \mathrm{F}}\left(r\right)\otimes \varsigma (\theta ,\varphi )e^{-i{\omega }_{\mathrm{F}}t}$. Here, $u_{\pm \mathrm{F}}\left(r\right)=R_{\pm \mathrm{F}}\left(r\right)/r{f}^{1/4}\left(r\right)$ is the fermionic radial function, the signs ± label the spins, ${\omega }_{\mathrm{F}}$ is the fermionic frequency, and $\varsigma (\theta ,\varphi )$ is a two-components fermion. Thus, the fermionic radial function $R_{\pm \mathrm{F}}\left(r\right)$ satisfies the following equation,

$$\frac{d^2{R}_{\pm \mathrm{F}}\left(r\right)}{d{r}^2}+\left[\frac{1}{2f\left(r\right)}\frac{df\left(r\right)}{dr}+\frac{1}{r}\right]\frac{d{R}_{\pm \mathrm{F}}\left(r\right)}{dr}+\frac{1}{f^2\left(r\right)}\left\{{\omega }_{\mathrm{F}}^2-f\left(r\right)\left[\frac{\lambda }{r^2}\mp \frac{i{\omega }_{\mathrm{F}}}{r}\pm \frac{i{\omega }_{\mathrm{F}}}{2f\left(r\right)}\frac{df\left(r\right)}{dr}\right]\right\}{R}_{\pm \mathrm{F}}\left(r\right)=0.$$

(15)

Indeed, these three radial equations can be condensed in a single one, which is labeled according to the spin, s, of the test fields, such that $s=(0,1,1/2)=({}_{\mathrm{B}},{}_{\mathrm{E}},{}_{\mathrm{F}})$, where we are choosing the “plus” spin. Thus, Equations (13)–(15) can be rewritten as

$$\frac{d^2{R}_s\left(r\right)}{d{r}^2}+\left[{\left(\frac{1}{2}\right)}^{4s(1-s)}\frac{1}{f\left(r\right)}\frac{df\left(r\right)}{dr}+\frac{4s(1-s)}{r}\right]\frac{d{R}_s\left(r\right)}{dr}+\frac{1}{f^2\left(r\right)}[{\omega }_s^2-V_s(r,{\omega }_s)]R_s\left(r\right)=0,$$

(16)

where the effective potential, $V_s(r,{\omega }_s)$, is given by

$$V_s(r,{\omega }_s)=f\left(r\right)\left\{\frac{\lambda }{r^2}-\frac{4is(1-s){\omega }_s}{r}+\left[\frac{2(s-1)(s-1/2)}{r}+\frac{4is(1-s){\omega }_s}{2f\left(r\right)}\right]\frac{df\left(r\right)}{dr}\right\}.$$

(17)

In fact, the effective potentials $V_{+\frac{1}{2}}$ and $V_{-\frac{1}{2}}$ can be transformed one into another by using the Darboux transformation, which means that both potentials will give the same spectrum of resonant frequencies.

Now, by following the Vieira-Bezerra-Kokkotas (VBK) approach [28,29], we can provide an analytical solution for the radial Equation (16) in terms of the general Heun functions [30]. It is given by

$$R_s^j\left(x\right)=x^{\frac{1}{2}[\gamma -2^{4s(s-1)}]}{(x-1)}^{\frac{1}{2}[\delta -2^{4s(s-1)}]}{(x-b)}^{\frac{1}{2}[ϵ-4s(1-s)]}[C_1^j{y}_1^j\left(x\right)+C_2^j{y}_2^j\left(x\right)],$$

(18)

where $C_1^j$ and $C_2^j$ are constants to be determined, and $j=\{0,1,b,\infty \}$ labels the solution at each singular point. The new radial coordinate, x, and the singularity parameter, b, are defined as

$$\begin{array}{ccc}\hfill x& =& \frac{r-r_{-}}{r_{+}-r_{-}},\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill b& =& -\frac{r_{-}}{r_{+}-r_{-}}.\hfill \end{array}$$

(20)

Thus, the pair of linearly independent solutions at $x=0$ ($r=r_{-}$) is given by

$$\begin{array}{ccc}\hfill y_1^0& =& \mathrm{HeunG}(b,q;\alpha ,\beta ,\gamma ,\delta ;x),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill y_2^0& =& x^{1-\gamma }\mathrm{HeunG}(b,(b\delta +ϵ)(1-\gamma )+q;\alpha +1-\gamma ,\beta +1-\gamma ,2-\gamma ,\delta ;x).\hfill \end{array}$$

(22)

Similarly, the pair of linearly independent solutions corresponding to the exponents 0 and $1-\delta$ at $x=1$ ($r=r_{+}$) is given by

$$\begin{array}{ccc}\hfill y_1^1& =& \mathrm{HeunG}(1-b,\alpha \beta -q;\alpha ,\beta ,\delta ,\gamma ;1-x),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill y_2^1& =& {(1-x)}^{1-\delta }\mathrm{HeunG}(1-b,((1-b)\gamma +ϵ)(1-\delta )+\alpha \beta -q;\alpha +1-\delta ,\beta +1-\delta ,2-\delta ,\gamma ;1-x).\hfill \end{array}$$

(24)

The pair of linearly independent solutions corresponding to the exponents 0 and $1-ϵ$ at $x=b$ ($\left|b\right|<1$) is given by

$$\begin{array}{ccc}\hfill y_1^b& =& \mathrm{HeunG}(\frac{b}{b-1},\frac{\alpha \beta b-q}{b-1};\alpha ,\beta ,ϵ,\delta ;\frac{b-x}{b-1}),\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill y_2^b& =& \left(\frac{b-x}{b-1}{)}^{1-ϵ}\mathrm{HeunG}\right(\frac{b}{b-1},\frac{\left(b\right(\delta +\gamma )-\gamma )(1-ϵ)}{b-1}+\frac{\alpha \beta b-q}{b-1};\alpha +1-ϵ,\beta +1-ϵ,2-ϵ,\delta ;\frac{b-x}{b-1}).\hfill \end{array}$$

(26)

Finally, the pair of linearly independent solutions corresponding to the exponents $\alpha$ and $\beta$ at $x=\infty$ is given by

$$\begin{array}{ccc}\hfill y_1^{\infty }& =& x^{-\alpha }\mathrm{HeunG}(\frac{1}{b},\alpha (\beta -ϵ)+\frac{\alpha }{b}(\beta -\delta )-\frac{q}{b};\alpha ,\alpha -\gamma +1,\alpha -\beta +1,\delta ;\frac{1}{x}),\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill y_2^{\infty }& =& x^{-\beta }\mathrm{HeunG}(\frac{1}{b},\beta (\alpha -ϵ)+\frac{\beta }{b}(\alpha -\delta )-\frac{q}{b};\beta ,\beta -\gamma +1,\beta -\alpha +1,\delta ;\frac{1}{x}).\hfill \end{array}$$

(28)

Here, $\mathrm{HeunG}(b,q;\alpha ,\beta ,\gamma ,\delta ;x)$ denotes a general Heun function, which is analytic in the disk $\left|x\right|<1$, and has a Maclaurin expansion given by

$$\mathrm{HeunG}(b,q;\alpha ,\beta ,\gamma ,\delta ;x)=\sum _{n=0}^{\infty }{c}_n{x}^n,$$

(29)

with

$$\begin{array}{ccc}\hfill -q{c}_0+b\gamma c_1& =& 0(c_0=1),\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill P_n{c}_{n-1}-(Q_n+q)c_n+X_n{c}_{n+1}& =& 0(n\ge 1),\hfill \end{array}$$

(31)

and

$$\begin{array}{ccc}\hfill P_n& =& (n-1+\alpha )(n-1+\beta ),\hfill \\ \hfill Q_n& =& n\left[\right(n-1+\gamma \left)\right(1+b)+b\delta +ϵ],\hfill \\ \hfill X_n& =& (n+1)(n+\gamma )b.\hfill \end{array}$$

(32)

In these solutions, the parameters $\alpha$, $\beta$, $\gamma$, $\delta$, $ϵ$, and q are given according to each value of the spin s, as follows. For scalar fields, we have

$$\begin{array}{ccc}\hfill {\alpha }_{\mathrm{B}}& =& \frac{5}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{B}}}{r_{+}-r_{-}},\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill {\beta }_{\mathrm{B}}& =& -\frac{1}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{B}}}{r_{+}-r_{-}},\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{B}}& =& 1-\frac{2i{\omega }_{\mathrm{B}}}{r_{+}-r_{-}},\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill {\delta }_{\mathrm{B}}& =& 1-\frac{2i{\omega }_{\mathrm{B}}}{r_{+}-r_{-}},\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill {ϵ}_{\mathrm{B}}& =& 1-\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}},\hfill \\ \hfill q_{\mathrm{B}}& =& \frac{1}{2r_{-}{(r_{+}-r_{-})}^4}\{2i{\omega }_{\mathrm{B}}[\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}]\hfill \\ & & +r_{-}[\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}]\hfill \\ & & -r_{+}\{\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}-r_{-}^4\left(\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}+4\right)\hfill \\ & & +r_{-}^2[3\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}-8{\omega }_{\mathrm{B}}^2-6\lambda ]+6i{r}_{-}{\omega }_{\mathrm{B}}\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}+14i{r}_{-}^3{\omega }_{\mathrm{B}}\}\hfill \\ & & -r_{+}^4{r}_{-}+2r_{+}^3(\lambda +2r_{-}^2-i{r}_{-}{\omega }_{\mathrm{B}})+r_{+}^2{r}_{-}[3\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}-6r_{-}^2+10i{r}_{-}{\omega }_{\mathrm{B}}-6\lambda ]\hfill \end{array}$$

(37)

$$\begin{array}{ccc}& & -r_{-}^5\left(\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}+1\right)-2i{r}_{-}^4{\omega }_{\mathrm{B}}\left(\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-3\right)-2r_{-}^3(\lambda +4{\omega }_{\mathrm{B}}^2)\}.\hfill \end{array}$$

(38)

For electromagnetic fields, we have

$$\begin{array}{ccc}\hfill {\alpha }_{\mathrm{E}}& =& \frac{3}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{E}}}{r_{+}-r_{-}},\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill {\beta }_{\mathrm{E}}& =& \frac{1}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{E}}}{r_{+}-r_{-}},\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{E}}& =& 1-\frac{2i{\omega }_{\mathrm{E}}}{r_{+}-r_{-}},\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill {\delta }_{\mathrm{E}}& =& 1-\frac{2i{\omega }_{\mathrm{E}}}{r_{+}-r_{-}},\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill {ϵ}_{\mathrm{E}}& =& 1-\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}},\hfill \\ \hfill q_{\mathrm{E}}& =& \frac{1}{2r_{-}{(r_{+}-r_{-})}^4}\{2i{\omega }_{\mathrm{E}}[\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}]\hfill \\ & & +r_{-}[\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}]\hfill \\ & & -r_{+}\{\sqrt{r_{+}^5{r}_{-}(4\lambda +r_{+}{r}_{-})}+3\sqrt{r_{+}{r}_{-}^5(4\lambda +r_{+}{r}_{-})}+r_{-}^4\left(4-\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}\right)\hfill \\ & & +r_{-}^2[3\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}-8{\omega }_{\mathrm{E}}^2-6\lambda ]+6i{r}_{-}{\omega }_{\mathrm{E}}\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}+14i{r}_{-}^3{\omega }_{\mathrm{E}}\}\hfill \\ & & +r_{+}^4{r}_{-}+2r_{+}^3(\lambda -2r_{-}^2-i{r}_{-}{\omega }_{\mathrm{E}})+r_{+}^2{r}_{-}[3\sqrt{r_{+}{r}_{-}(4\lambda +r_{+}{r}_{-})}+6r_{-}^2+10i{r}_{-}{\omega }_{\mathrm{E}}-6\lambda ]\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& & r_{-}^5\left(\sqrt{1-\frac{4\lambda }{r_{+}{r}_{-}}}+1\right)-2i{r}_{-}^4{\omega }_{\mathrm{E}}\left(\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-3\right)-2r_{-}^3(\lambda +4{\omega }_{\mathrm{E}}^2)\}.\hfill \end{array}$$

(44)

Finally, for Dirac fields, we have

$$\begin{array}{ccc}\hfill {\alpha }_{\mathrm{F}}& =& \frac{3}{2}-\sqrt{\frac{\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{F}}}{r_{+}-r_{-}},\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill {\beta }_{\mathrm{F}}& =& \frac{1}{2}-\sqrt{\frac{\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{F}}}{r_{+}-r_{-}},\hfill \end{array}$$

(46)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{F}}& =& \frac{3}{2}-\frac{2i{\omega }_{\mathrm{F}}}{r_{+}-r_{-}},\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill {\delta }_{\mathrm{F}}& =& \frac{1}{2}-\frac{2i{\omega }_{\mathrm{F}}}{r_{+}-r_{-}},\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill {ϵ}_{\mathrm{F}}& =& 1-2\sqrt{\frac{\lambda }{r_{+}{r}_{-}}},\hfill \\ \hfill q_{\mathrm{F}}& =& \frac{1}{4r_{-}{(r_{+}-r_{-})}^4}\{24\sqrt{\lambda r_{+}^5{r}_{-}^3}+8i{\omega }_{\mathrm{F}}\sqrt{\lambda r_{+}^5{r}_{-}}+2r_{+}^4{r}_{-}-36\sqrt{\lambda r_{+}^3{r}_{-}^5}-24i{\omega }_{\mathrm{F}}\sqrt{\lambda r_{+}^3{r}_{-}^3}-6\sqrt{\lambda r_{+}^7{r}_{-}}\hfill \\ & & +r_{+}^3(4\lambda -9r_{-}^2-8i{r}_{-}{\omega }_{\mathrm{F}})+r_{+}^2{r}_{-}\left(15r_{-}^2+32i{r}_{-}{\omega }_{\mathrm{F}}-12\lambda \right)-6\sqrt{\frac{\lambda r_{-}^9}{r_{+}}}-8i{\omega }_{\mathrm{F}}\sqrt{\frac{\lambda r_{-}^7}{r_{+}}}+24\sqrt{\lambda r_{+}{r}_{-}^7}\hfill \end{array}$$

(49)

$$\begin{array}{ccc}& & +24i{\omega }_{\mathrm{F}}\sqrt{\lambda r_{+}{r}_{-}^5}+r_{+}{r}_{-}^2(12\lambda -11r_{-}^2-40i{r}_{-}{\omega }_{\mathrm{F}}+16{\omega }_{\mathrm{F}}^2)+3r_{-}^5+16i{r}_{-}^4{\omega }_{\mathrm{F}}-4r_{-}^3(\lambda +4{\omega }_{\mathrm{F}}^2)\}.\hfill \end{array}$$

(50)

Next, we will use these analytical solutions of the radial equations in the 4DEGBBH spacetime, and the properties of the general Heun functions, to compute the spectrum of resonant frequencies related to the quasibound states.

4. Quasibound States

The quasibound states, also known as quasistationary levels or resonance spectra, are a kind of wave phenomena occurring near the black hole exterior event horizon, that is, they are localized in the black hole potential well. The definition of quasibound states means that there exist a flux of particles crossing into the black hole, so that the spectrum of quasibound states is constituted by complex frequencies, which can be expressed as $\omega ={\omega }_R+i{\omega }_I$, where ${\omega }_R$ and ${\omega }_I$ are the real and imaginary parts, respectively. By using these results, it is possible, in principle, to get some information about the physics of black holes as well as to validate some alternative/modified theories of gravity.

The real part of the resonant frequency is the oscillation frequency, while the imaginary part determines the stability of the system. Thus, the wave solution is said to be stable when the imaginary part of the resonant frequency is negative (${\omega }_I<0$), which means a decay rate with the time. Otherwise, the wave solution is unstable when the imaginary part of the resonant frequency is positive (${\omega }_I>0$), which means a growth rate with the time.

In order to compute the spectrum of quasibound states, we need to impose two boundary conditions, which are related to the asymptotic behavior of the radial solution. First, it should describe an ingoing wave at the exterior event horizon. Second, it should tend to zero far from the black hole at asymptotic infinity, that is, the probability to find such a particle in the spatial infinity must be zero.

To derive the characteristic resonance equation, many authors have been using a method which consists of solving the radial equation in two different asymptotic regions and then matching these two radial solutions in their common overlap region. However, in this work, we will use the VBK approach (for details, see Refs. [28,29]) to derive the characteristic resonance equation and then find the spectrum of resonant frequencies related to quasibound states. This approach consists in (i) obtaining the ingoing wave solution at the exterior event horizon, (ii) imposing the polynomial condition of the general Heun functions, and then (iii) computing the resonant frequencies.

Thus, in the limit when $r\to r_{+}$, which implies that $x\to 1$, the radial solution given by Equation (18) behaves as

$$\underset{r\to r_{+}}{lim}{R}_s^1\left(r\right)\sim C_1^1{(r-r_{+})}^{\frac{1}{2}[\delta -2^{4s(s-1)}]}+C_2^1{(r-r_{+})}^{-\frac{1}{2}[\delta -2^{4s(s-1)}]},$$

(51)

where $C_1^1$ and $C_2^1$ include all remaining constants. On the other hand, from Equations (36), (42) and (48), we get

$$\frac{1}{2}[\delta -2^{4s(s-1)}]=-\frac{i{\omega }_s}{2{\kappa }_{+}},$$

(52)

where ${\kappa }_{+}$ is the gravitational acceleration on the exterior event horizon and given by

$${\kappa }_{+}\equiv \frac{1}{2}{\left.\frac{df\left(r\right)}{dr}\right|}_{r=r_{+}}=\frac{r_{+}-r_{-}}{2}.$$

(53)

Then, we can rewrite Equation (51) as

$$\underset{r\to r_{+}}{lim}{R}_s^1\left(r\right)\sim C_1^1{R}_{\mathrm{in}}^1\left(r\right)+C_2^1{R}_{\mathrm{out}}^1\left(r\right),$$

(54)

where the ingoing, $R_{\mathrm{in}}^1\left(r\right)$, and outgoing, $R_{\mathrm{out}}^1\left(r\right)$, wave solutions are given by

$$\begin{array}{ccc}\hfill R_{\mathrm{in}}^1(r>r_{+})& =& {(r-r_{+})}^{-\frac{i{\omega }_s}{2{\kappa }_{+}}},\hfill \end{array}$$

(55)

$$\begin{array}{ccc}\hfill R_{\mathrm{out}}^1(r>r_{+})& =& {(r-r_{+})}^{\frac{i{\omega }_s}{2{\kappa }_{+}}}.\hfill \end{array}$$

(56)

Therefore, in order to fully satisfy the first boundary condition, we must impose that $C_2^1=0$ in Equation (54), as well as in Equation (18). Finally, we obtain the radial wave solution describing ingoing test particles at the 4DEGBBH exterior event horizon, which is given by

$$\underset{r\to r_{+}}{lim}{R}_s^1\left(r\right)\sim C_1^1{R}_{\mathrm{in}}^1(r>r_{+}).$$

(57)

On the other hand, in the limit when $r\to \infty$, which implies that $x\to \infty$, the radial solution given by Equation (18) behaves as

$$\underset{r\to \infty }{lim}{R}_s^{\infty }\left(r\right)\sim C_1^{\infty }{r}^{{\sigma }_s},$$

(58)

with ${\sigma }_s=D_s-{\alpha }_s$, where the coefficients $D_s$ are given by

$$\begin{array}{ccc}\hfill D_{\mathrm{B}}& =& -\frac{1}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{B}}}{r_{+}-r_{-}},\hfill \end{array}$$

(59)

$$\begin{array}{ccc}\hfill D_{\mathrm{E}}& =& -\frac{1}{2}-\frac{1}{2}\sqrt{1+\frac{4\lambda }{r_{+}{r}_{-}}}-\frac{2i{\omega }_{\mathrm{E}}}{r_{+}-r_{-}},\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill D_{\mathrm{F}}& =& -1+\frac{1}{r_{+}-r_{-}}\left(\sqrt{\frac{r_{-}\lambda }{r_{+}}}-\sqrt{\frac{r_{+}\lambda }{r_{-}}}-2i{\omega }_{\mathrm{F}}\right).\hfill \end{array}$$

(61)

Thus, the sign of the real part of ${\sigma }_s$ determines the asymptotic behavior of the radial solution far from the black hole at asymptotic infinity. Therefore, the radial solution tends to zero if $\mathrm{Re}\left[{\sigma }_s\right]<0$ and then it will fully satisfy the second boundary condition, which describes the quasibound states. Otherwise, if $\mathrm{Re}\left[{\sigma }_s\right]>0$, the radial solution diverges. The final asymptotic behavior of the radial solution will be determined when we know the values of the frequencies ${\omega }_s$.

Now, to find the characteristic resonance equation, we match these two asymptotic radial solutions in their common overlap region by imposing the polynomial conditions of the general Heun functions, which are given by

$$\begin{array}{ccc}\hfill \alpha & =& -n,\hfill \end{array}$$

(62)

$$\begin{array}{ccc}\hfill c_{n+1}\left(q\right)& =& 0,\hfill \end{array}$$

(63)

where $n=0,1,2,\dots$ is the overtone number, which can be, without loss of generality, called the principal quantum number. From the first polynomial condition, given by Equation (62), we can find the frequency eigenvalues. On the other hand, from the second polynomial condition, given by Equation (63), we can determine the values of the separation constant $\lambda$ for each value (mode) n, and then write the radial and angular wave eigenfunctions. Thus, by imposing Equation (62) on each parameter ${\alpha }_s$, we obtain the exact analytical spectrum of quasibound states for massless test fields in the 4DEGBBH spacetime. They are given by

$$\begin{array}{ccc}\hfill {\omega }_{\mathrm{B}n}& =& -\frac{1}{4}i\sqrt{1-2a}\left(5+2n-\sqrt{\frac{a+8\lambda }{a}}\right),\hfill \end{array}$$

(64)

$$\begin{array}{ccc}\hfill {\omega }_{\mathrm{E}n}& =& -\frac{1}{4}i\sqrt{1-2a}\left(3+2n-\sqrt{\frac{a+8\lambda }{a}}\right),\hfill \end{array}$$

(65)

$$\begin{array}{ccc}\hfill {\omega }_{\mathrm{F}n}& =& -\frac{1}{4}i\sqrt{1-2a}\left(3+2n-\sqrt{\frac{8\lambda }{a}}\right).\hfill \end{array}$$

(66)

On the other hand, by imposing Equation (63) on each parameter $q_s$, we determine the exact analytical value of the separation constant ${\lambda }_{sn}$ (for details, see Refs. [31,32,33]). However, the final expressions are quite long, and for this reason no insight is gained by writing them out. Thus, instead of doing this, we will present some of their features in what follows.

In

Table 1. The quasibound states ${\omega }_{sn}$, the separation constants ${\lambda }_{sn}$, and the real part of coefficients ${\sigma }_{sn;\eta }$. We focus on the fundamental mode $n=0$.

a	${\mathit{\omega }}_{\mathbf{B}0}$	${\mathit{\lambda }}_{\mathbf{B}0;+}$	$\mathbf{Re}\left[{\mathit{\sigma }}_{\mathbf{B}0}\right]$	${\mathit{\omega }}_{\mathbf{E}0}$	${\mathit{\lambda }}_{\mathbf{E}0;+}$	$\mathbf{Re}\left[{\mathit{\sigma }}_{\mathbf{E}0}\right]$	${\mathit{\omega }}_{\mathbf{F}0}$	${\mathit{\lambda }}_{\mathbf{F}0;+}$	$\mathbf{Re}\left[{\mathit{\sigma }}_{\mathbf{F}0}\right]$
$0.05$	$-0.711512i$	$0.018750$	$-3$	$-0.474342i$	0	$-2$	$-0.249342i$	$0.023734$	$-2.500000$
$0.10$	$-0.670820i$	$0.037500$	$-3$	$-0.447214i$	0	$-2$	$-0.247214i$	$0.044861$	$-2.500000$
$0.15$	$-0.627495i$	$0.056250$	$-3$	$-0.418330i$	0	$-2$	$-0.243330i$	$0.063250$	$-2.500000$
$0.20$	$-0.580948i$	$0.075000$	$-3$	$-0.387298i$	0	$-2$	$-0.237298i$	$0.078730$	$-2.500000$
$0.25$	$-0.530330i$	$0.093750$	$-3$	$-0.353553i$	0	$-2$	$-0.228553i$	$0.091069$	$-2.500000$
$0.30$	$-0.474342i$	$0.112500$	$-3$	$-0.316228i$	0	$-2$	$-0.216228i$	$0.099934$	$-2.500000$
$0.35$	$-0.410792i$	$0.131250$	$-3$	$-0.273861i$	0	$-2$	$-0.198861i$	$0.104801$	$-2.500000$
$0.40$	$-0.335410i$	$0.150000$	$-3$	$-0.223607i$	0	$-2$	$-0.173607i$	$0.104721$	$-2.500000$
$0.45$	$-0.237171i$	$0.168750$	$-3$	$-0.158114i$	0	$-2$	$-0.133114i$	$0.097451$	$-2.500000$

, we present some values of the quasibound state ${\omega }_{sn}$, as well as the corresponding separation constants ${\lambda }_{sn}$ and coefficients ${\sigma }_{sn}=D_{sn}+n$, as functions of the GB coupling constant a. We also show the behavior of the quasibound states ${\omega }_{sn}$ in Figure 2, as functions of the GB coupling constant a. From these results, we can conclude that all solution for the massless resonant frequencies ${\omega }_{sn}$ are physically admissible in the fundamental mode, which represents the spectrum of quasibound states in the 4DEGBBH spacetime. In this scenario, the radial solution given by Equation (18) tends to zero far from the 4DEGBBHs at asymptotic infinity, since $\mathrm{Re}\left[{\sigma }_{sn}\right]<0$. In addition, we can see that the motion of all the test fields are overdamped (with purely imaginary frequencies). Furthermore, all the test fields may describe stable systems, since there is no change in the sign of the imaginary part of their massless resonant frequencies in the stability region, that is, when $0<a<1/2$. We can also see that the imaginary part of the massless resonant frequencies ${\omega }_{sn}$ (for all the test fields in these modes) decreases (in modulus) as the GB coupling constant a approaches to 1/2.

It is worth noticing that we cannot determine the quasibound states ${\omega }_{sn}$, nor the separation constants ${\lambda }_{sn}$ and the coefficients ${\sigma }_{sn}$, for the particular case when $a=0$, which corresponds to the General Relativity (GR) limit. This is due to the fact that when we set $a=0$, the singularity parameter b goes to 0 and hence the radial solution must be recalculated in this scenario, since it will occur a confluent process involving the already existing singularity at the point $x=0$ with the ones in the point $x=b\to 0$, which will lead to a new radial solution given in terms of the confluent Heun functions, with totally new parameters. Henceforth, for the case of scalar fields, we can compare this result with the one obtained by Muniz et al. [34], where the spectrum of quasibound states for the Schwarzschild black hole ($a=0$, standard GR case) is such that ${\omega }_n=-i(n+1)/2$, for $M=1/2$; we see that there exist a peculiar difference between the found solutions.

Finally, we can also discuss the time dependence of the wave eigenfunctions, ${\Psi }_{sn}\left(t\right)={\mathrm{e}}^{-i{\omega }_{sn}t}$, which is shown in Figure 3. From this plot, we can realize the “final flight” of the test particles when crossing into the 4DEGBBH exterior event horizon.

Transition Frequency

In addition, we can calculate the transition frequency, $\Delta \omega$, between two highly damped ($n\to \infty$) neighboring states [35]. It is given by

$$\Delta \omega \approx \mathrm{Im}\left[{\omega }_{s(n-1)}\right]-\mathrm{Im}\left[{\omega }_{sn}\right]={\kappa }_{+}=\frac{1}{2}\sqrt{1-2a}.$$

(67)

On the other hand, the natural adiabatic invariant quantity, $I_{\mathrm{ad}}$, for a test field-black hole system with total energy E, is given by

$$I_{\mathrm{ad}}=\int \frac{dE}{\Delta \omega }=\int \frac{T_{\mathrm{H}}d{S}_{\mathrm{BH}}}{\Delta \omega }=\frac{\hslash S_{\mathrm{BH}}}{2\pi },$$

(68)

where $T_{\mathrm{H}}$ is the Hawking temperature and $S_{\mathrm{BH}}(={\mathcal{A}}_{+}/4\hslash )$ is the Bekenstein-Hawking entropy, with ${\mathcal{A}}_{+}$ being the surface area of the exterior event horizon. In this limit, the Bohr-Sommerfeld quantization condition applies and hence $I_{\mathrm{ad}}$ is a quantized quantity, namely, $I_{\mathrm{ad}}=n\hslash$. Thus, from Equation (68) we get

$$S_{\mathrm{BH}n}=2\pi n.$$

(69)

Therefore, the area spectrum is given by

$${\mathcal{A}}_{+n}=8\pi n\hslash ,$$

(70)

so that its minimum change becomes

$$\Delta {\mathcal{A}}_{+}^{\min }=8\pi \hslash .$$

(71)

From Equations (69) and (70), we can conclude that both entropy and area spectra are equally spaced. Furthermore, they do not depend on the black hole parameters. Finally, Equation (71) shows that the 4DEGBBH exterior event horizon is made by patches with equal area.

Next, we will use the polynomial conditions of the general Heun functions, as well as the frequency eigenvalues, to discuss both radial and angular wave eigenfunctions.

5. Wave Eigenfunctions

In this section, we continue to use the VBK approach (for details, see Refs. [31,32,33]) in order to present both radial and angular wave eigenfunctions describing the quasibound states of massless test particles propagating in the 4DEGBBH spacetime.

5.1. Radial Eigenfunctions

The quasibound state radial wave eigenfunctions are given in terms of the general Heun polynomials, which can be denoted as ${\mathrm{HeunGp}}_{n;\eta }(b,q;-n,\beta ,\gamma ,\delta ;x)$, where the parameter $\eta$ is related to the appropriate determination of the accessory parameter $q=q_{n;\eta }$ via the second polynomial condition given by Equation (63), which are the solutions that cut (in a certain order) the power series describing the general Heun functions; it is such that $\eta =0,1,2,\dots ,n$. In our case, the general Heun polynomials depend also on the spin s and hence they should be denoted as ${\mathrm{HeunGp}}_{sn;\eta }(b,q_{sn;\eta };-n,\beta ,\gamma ,\delta ;x)$.

Therefore, the quasibound state radial wave eigenfunctions for massless test particles propagating in the 4DEGBBH spacetime are given by

$$R_{sn;\eta }\left(x\right)=C_{sn;\eta }{x}^{\frac{1}{2}[\gamma -2^{4s(s-1)}]}{(x-1)}^{\frac{1}{2}[\delta -2^{4s(s-1)}]}{(x-b)}^{\frac{1}{2}[ϵ-4s(1-s)]}{\mathrm{HeunGp}}_{sn;\eta }(b,q_{sn;\eta };-n,\beta ,\gamma ,\delta ;x),$$

(72)

where $C_{sn;\eta }$ is a constant to be determined. Thus, the full quasibound state radial wave eigenfunctions, $u_{sn;\eta }\left(r\right)$, are given by

$$\begin{array}{ccc}\hfill u_{\mathrm{B}n;\eta }\left(r\right)& =& \frac{R_{\mathrm{B}n;\eta }\left(r\right)}{r},\hfill \end{array}$$

(73)

$$\begin{array}{ccc}\hfill u_{\mathrm{E}n;\eta }\left(r\right)& =& \frac{R_{\mathrm{E}n;\eta }\left(r\right)}{r},\hfill \end{array}$$

(74)

$$\begin{array}{ccc}\hfill u_{\mathrm{F}n;\eta }\left(r\right)& =& \frac{R_{\mathrm{F}n;\eta }\left(r\right)}{r{f}^{1/4}\left(r\right)}.\hfill \end{array}$$

(75)

The squared full quasibound state radial wave eigenfunctions are presented in Figure 4, as functions of the radial coordinate r for some values of the GB coupling constant a. From these results, we can conclude that all the resonant frequencies ${\omega }_{sn}$, given by Equations (64)–(66), describe massless quasibound states in the 4DEGBBH spacetime, since the radial solutions tend to zero at asymptotic infinity and diverge at the exterior event horizon. Note that these eigenfunctions reach a maximum value (at the exterior event horizon $r_{+}=0.994975$ for $a=0.01$, $r_{+}=0.853553$ for $a=0.25$, and $r_{+}=0.570711$ for $a=0.49$) and then cross into the 4DEGBBH, as shown in the log-scale plots.

5.2. Angular Eigenfunctions

It is known that the angular solution must be regular at its two boundaries, namely, when $\theta =0$ and $\theta =\pi$. This requirements single out a discrete set of angular eigenvalues $\lambda$, which couples the angular and radial equations. In our case, we can obtain the values for the angular eigenvalues ${\lambda }_{sn;\pm }$ by substituting the first polynomial condition, given by Equation (62), into the second polynomial condition, given by Equation (63). The ±signs denote the positive and negative solutions; Here, we choose the positive values, in order to express all the angular solutions in terms of the associated Legendre functions.

For the fundamental mode $n=0$, the eigenvalues $q_{s0;\eta }$ must obey the relation $c_1=0$, where $c_1=q/b\gamma$. Thus, we have $q_{s0;\eta }=0$, which implies that there exist only on solution for q, namely, $q_{s0;0}=0$, and hence we can obtain the values for the angular eigenvalues ${\lambda }_{s0;+}$. They are presented in Table 1. As it was expected, the angular eigenvalues ${\lambda }_{s0;+}$ are in the set of real numbers $\mathbb{R}$, and hence the eigenvalues for the associated Legendre functions can be written as $\lambda =\nu (\nu +1)$. Therefore, we present the behavior of the quasibound state angular wave eigenfunctions in Figure 5, Figure 6 and Figure 7, for the scalar, electromagnetic, and Dirac fields, respectively, as functions of the new angular coordinate $z=cos\theta$ for some values of the GB coupling constant a. It is worth calling attention to fact that (i) the numerically satisfactory solutions of the associated Legendre equation of general (includng complex) degree $\nu$ are given in terms of the Ferrers functions of the first kind $P_{\nu }^{-m}(-z)$ and $P_{\nu }^{-m}\left(z\right)$, in the interval $-1<z<1$, and (ii) these solutions are regular at the two boundaries $\theta =0$ ($z=1$) and $\theta =\pi$ ($z=-1$).

6. Final Remarks

In this work, we obtained exact analytical solutions of the master wave equations for the test fields in a 4-dimensional Einstein–Gauss–Bonnet black hole spacetime, where the radial solution is given in terms of the general Heun functions, and the angular solution is given in terms of the associated Legendre (Ferrers) functions (of the first kind).

We imposed two boundary conditions on the radial solution in order to study its asymptotic behaviors, which led to the quasibound state phenomena. Near the exterior event horizon, the radial solution describes ingoing waves, which reachs a maximum value and then crosses into the black hole. On the other hand, far from the black hole at asymptotic infinity, the radial solution tends to zero, that is, the probability of finding any particles in the spatial infinity is null.

The spectrum of quasibound states for the test fields was obtained by using the polynomial condition of the general Heun functions. In fact, that is a new (analytical) approach developed by Vieira, Bezerra, and Kokkotas [28,29]. It is worth pointing out that these massless resonant frequencies were obtained directly from the general Heun functions, and, to our knowledge, there is no similar result in the literature for the background under consideration. In addition, it is worth emphasizing that all the numerical/graphical computations performed in this work were carried out by using the standard package for the general Heun functions installed in Wolfram Mathematica 12.3.

Finally, we have discussed the stability of the system. All the systems are stables in the fundamental mode, and present an overdamped motion, in the range $0<a<1/2$. We hope that our results, which describe an unquestionably phenomenon associated with purely quantum effects in gravity, may be used to fit some astrophysical data in the near future, as for example, the ones related to the observations of some spectrum of thin accretion disks with present and future X-ray facilities [36], as well as from the secondary object in GW190814 which is compatible with being a slowly-rotating neutron star in EGB theory of gravity [37], and hence shed some light on the physics of black holes and compact objects.

As a future perspective, it is possible to extend our results for the case of an asymptotically de Sitter spacetime, as well as for negative values of the GB coupling constant. In addition, we can obtain a new acoustic curved black hole embedded in the 4DEGB spacetime, which could be a very interesting framework within the analog models of gravity.