Singularities of regular black holes and the monodromy method for asymptotic quasinormal modes

Regular black holes (RBHs) are a collection of black holes (BHs) that do not have singularities in spacetime, especially at their centers [1–3], that is, all the curvature invariants of BHs, such as the Kretschmann scalar, are finite everywhere. This definition is related to Markov's limiting curvature conjecture [4, 5], which states that the curvature invariants should be restricted by a universal value. There is also a well-known definition of the regularity of spacetime through the completeness of null and timelike geodesics [6, 7]. Unfortunately, these two definitions are not equivalent. Some BHs have finite curvature everywhere, but the spacetime contains incomplete geodesics, for example, the Taub–NUT BH [8, 9]. Others are the opposite, namely, the geodesics in spacetime are complete everywhere, but the curvatures can be divergent in a certain area. For example, the wormhole model in Ref. [10].

The study of RBHs dates back to Sakharov and Gliner [11, 12] and is interesting because RBHs are not restricted by the Penrose singularity theorem [13]. The spacetime singularity can be avoided by replacing the center with a de Sitter core [14]. The first RBH model was implemented by Bardeen [15], in which a RBH was obtained by a formal modification of the Schwarzschild BH. A significant breakthrough came in 1998, when Ayon-Beato and Garcia offered the first interpretation of matter, in terms of classical fields, that generates RBHs [1]. In their follow-up research [16], Ayon-Beato and Garcia also showed that the Bardeen BH can be interpreted by the nonlinear magnetic monopole. This idea was generalized further by Fan and Wang [17]. Now, all spherically symmetric RBHs with singular shape functions can be interpreted in the context of nonlinear electrodynamics, that is, given the metric of a RBH, the Lagrangian of matter generation can be derived in terms of electrodynamic field strength. These results stimulated further research on RBHs, including not only their construction [2, 18–23], but also the interpretation of RBH models [24–27] and various physical phenomena in RBHs spacetime, such as BH thermodynamics [28, 29], BH shadows [30, 31], and quasinormal modes (QNMs) [32–34].

The research program for RBHs may differ from that of traditional BHs with singularities. In the latter, the full action of the system is first given, and the BH solution can then be solved using the Euler–Lagrange equation obtained from the variational principle. In the former, the metric of an RBH can be constructed from the aspect of finite curvature invariants. Subsequently, information on the matter source can be deduced within a specific context that facilitates the solution, such as nonlinear electrodynamics [35–37], in the presence of a phantom scalar field [25, 38] or modified gravity [39, 40].

The area/entropy spectrum of BHs is known to be related to their QNMs [38, 41] and, specifically, to the imaginary component of QNMs in the large damping limit, that is, asymptotic QNMs (AQNMs) [42]. In other words, the area/entropy spectrum is obtained via the Bohr–Sommerfeld quantization of an adiabatic invariant [43] constructed using the AQNMs, whereas the AQNMs of BHs can be determined by the so-called monodromy approach [44], in which singularities and Stokes lines play crucial roles [45]. As a result, studying the quantum entropy spectrum of RBHs requires a thorough understanding of AQNMs. The goal of this study is to investigate the AQNMs of RBHs in greater detail using a clear Stokes diagram.

For a traditional BH with a singularity at the center, such as the Schwarzschild BH, the Stokes lines emit from the central singularity. Thus, when the asymptotic solutions of the Regge–Wheeler master equation rotate from one Stokes line to another around the zero point of the radial coordinate (that is, the center of the BH), the change is not trivial, and the monodromy along a closed Stokes line will give an analytical expression for AQNMs (see, for example, [46, 47]).

For an RBH [3], the zero point of the metric is no longer singular in most cases, and the corresponding Stokes lines do not converge to that point, which implies that the rotation of the asymptotic solution around the zero point is trivial. Several attempts have been made by ignoring this problem [32, 48], where the Stokes lines of RBHs are borrowed from Reissner–Nordström (RN) or RN-like BHs based on the similarity of the metrics at infinity, and the asymptotic frequencies of these models have a universal form (see Eq. (12) in Ref. [32]). In other words, to obtain the monodromy, the RBH is regarded as a singular BH. Therefore, one may ask what the true Stokes lines of RBHs are and what the monodromy based on these is. In the current study, we give a detailed discussion.

This paper is organized as follows: In Sec. II, we illustrate and classify the singularities of RBHs from the perspective of complex analysis, which can be regarded as a prelude for the calculation of AQNMs in the subsequent sections.

After introducing an alternative approach to displaying the Stokes diagram in Sec. III, we calculate the AQNMs of RBHs that still have complex singularities as the poles of the shape functions based on the explicit Stokes portrait in Sec. IV. Sec. V contains other examples in which the RBHs have essential singularities or no singularities at all.

We also provide a brief comparison of the monodromy method and complex WKB method from the perspective of the Stokes portrait in Sec. VI, and a conclusion is given in Sec. VII.

Let us start with spherically symmetric BHs with a single shape function $ f(r) $ in the metric

where we suppose that $ f(r) $ is of the following form:

Here, $ M $ is the BH mass, and $ \sigma(r) $ is a dimensionless function of the radial coordinate. A BH with the metric in Eq. (1) is regular^① if all its curvature invariants referring to the Riemann tensor are finite for $ r\in [0,\infty) $. For example, the Ricci curvature $ R:= g^{\mu\nu} R_{\mu\nu} $ and Weyl curvature $W:= W_{\mu\nu\alpha\beta}W^{\mu\nu\alpha\beta}$, where $ R_{\mu\nu} $ and $ W_{\mu\nu\alpha\beta} $ are the Ricci and Weyl tensors, respectively. From a physical perspective, "regularity" implies that there are no gravitational singularities within the BH horizons. Furthermore, all RBHs with the Eq. (1) metric can be interpreted as magnetic monopoles [17].

For the Eq. (1) metric, we can compute the Ricci and Weyl curvatures

from which we note that the singularities of the curvature invariants are unavoidable if the function $ f(r) $ has singularities in $ r\in [0,\infty) $. Additional restricted conditions for $ f(r) $ that guarantee the finiteness of the curvature invariants at $ r=0 $ have been discussed, for example, in Refs. [5, 17, 21, 29].

Furthermore, the Ricci and Weyl curvatures can be simplified if we replace $ \sigma\to s_1(r)/(2 M r) $ for $ R $ and $ \sigma \to \sqrt{3} r^2 s_2(r)/(2M) $ for $ W $, yielding

which provides convenience when constructing RBHs with various types of mathematical singularities.

In contrast, according to the monodromy method, as we compute the AQNMs for BHs with gravitational singularities at $ r=0 $, such as Schwarzschild and RN BHs, the radial coordinate $ r $ in the Regge–Wheeler master equation is analytically continued into the complex plane, and the physical singularities (gravitational singularities and BH horizons)^② play a key role [44]. Therefore, it is natural to wonder what to do when dealing with the AQNMs of RBHs if there are no gravitational singularities.

In fact, the mechanisms for removing the singularities of BHs are considerably different from each other from the perspective of complex analysis, and these differences are reflected in the calculations of AQNMs. Therefore, we first divide the RBHs into three classes according to the types of singularities of the curvature invariants after analytically continuing the radial coordinate into the complex plane.

1. The curvature invariants are entire functions of $ r $, that is, they have either no singularities on the entire complex plane $ \mathbb{C} $ or removable singularities at $ r=0 $.

2. The curvature invariants are meromorphic functions on $ \mathbb{C} $. In other words, the point $ r=0 $ is a regular point; however, the singularities (as poles) still exist and are located in the ranges beyond the non-negative axis of $ r $, that is, the singularities of the curvature invariants are dragged into the non-physical domain.

3. $ r=0 $ is an essential singularity of the curvature invariants. Via Picard's great theorem [50], RBHs with an essential singularity at $ r=0 $ must have infinitely many complex horizons on any punctured neighborhood of $ r=0 $.

Several BH examples can be used to help understand this classification. First, we propose an RBH model with^③

where $ Q $ is a parameter introduced to balance the dimensions.

To reduce the parameters, we rescale $ r\to l^{2/3} (2M)^{1/3} z $ and $ Q\to 2^{7/6} l^{1/3} M^{2/3} $, where $ z $ is a dimensionless radial coordinate, and $ l $ is a parameter of length dimension and is expected to be the Planck length or of the same order. Then, the horizons $ z_{\rm{H}} $ can be determined using the equation

Because sine is a transcendental function, there must be infinitely many complex horizons for any $ \lambda\in\mathbb{R} $ and no Picard exceptional values according to Picard's little theorem [50, 51]; however, the number of complex horizons is finite in any circle around $ z=0 $.

We also note that the number of physical horizons increases discretely with increasing $ \lambda $ in the interval $ (0.238, 1/4] $. The model has only one horizon when $ \lambda\lesssim 0.238 $, two horizons when $ \lambda\approx 0.238 $, and an infinite number of horizons when $ \lambda=1/4 $. The Ricci and Weyl curvatures can be computed as

where $ z=0 $ is a removable singularity belonging to the first type above. Moreover, it is easy to verify that this model meets the dominant energy condition (DEC) [52].

To understand another case in the first type, we propose a model of a RBH with

Under the rescalings $ r\to l^{2/3} (2 M)^{1/3} z, l \to l^{2/3} (2 M)^{1/3} $, we can fix the horizons via

where $ \lambda $ is formally the same as defined in Eq. (6). When $ \lambda\lesssim 0.173 $, there are two physical horizons, whereas there are infinitely many complex horizons. In addition, the DEC of this model also holds. The typical curvature invariants are then

which are surely finite in the physical domain. In contrast, as we analytically continue $ z $ into the entire complex plane, although there are no singularities in the finite range, divergence occurs at complex infinity owing to the Stokes phenomenon of the function $ \mathrm{e}^{-z} $. In other words, as $ z $ approaches infinity along a path restricted in the sector

the curvature invariants become divergent, namely, complex infinity is an essential singularity of the model.

Now, we turn to the Hayward BH [2] with $ \sigma=r^3/\left(r^3+ 2Ml^2\right) $, where $ l/(2M)\ll 1 $. The curvature invariants of this model are finite in the domain $ r\in [0,\infty) $. However, after analytically continuing $ r $ into the complex plane, the Weyl curvature

starts having three singularities in the non-physical domain, which are determined using the algebraic equation $ 2 l^2 M+r^3=0 $. More precisely, the Hayward BH has three poles after the analytical continuation. This is what we have in the second type, that is, the singularities are hidden in the non-physical domain.

Furthermore, we can determine that all three singularities are located in the outermost horizons. The proof is a typical application of Rouché's theorem [50]. First, we note from $ 2 l^2 M+r^3=0 $ that all the singularities are located on a circle with the radius $ |{r}|= (2M l^2)^{1/3} $. Second, we can rewrite $ f(r)=0 $ as $ r^3-2M r^2+2M l^2=0 $, that is, we obtain three complex horizons using the fundamental theorem of algebra. Third, on the circle $ |{r}|= (2M l^2)^{1/3}+\epsilon $, where $ \epsilon\ll(2M l^2)^{1/3} $, we have $ |{2M r^2}|>|{r^3}|+|{2M l^2}| $ if $ l/2M<1/(2\sqrt{2})\approx 0.35 $ and $ \epsilon/(2M) < \left(\sqrt{5}-1\right)/4\approx 0.31 $. Finally, because $ l/(2M)\ll 1 $, we can conclude that there are two horizons in the circle $ |{r}|= (2M l^2)^{1/3}+\epsilon $ by applying Rouché's theorem. In other words, the three complex horizons are sandwiched between the outermost and next to outermost horizons.

At last, the widely discussed model [21] with $ \sigma(r)= \exp[-P^2/(2Mr)] $ belongs to the third type. After rescalings $r\to l^{2/3} \sqrt[3]{2 M} z,~P\to \sqrt[3]{l} (2 M)^{2/3}$, its curvature invariants can be represented by

where $ z=0 $ becomes an essential singularity after analytical continuation into the complex plane. A detailed discussion on this model is given in Sec. V.A.

We end this section with a comment on constructing physical RBHs, Eqs. (5) and (8). Here, "physical" implies that the RBHs must obey the DEC [52]. Instead of the traditional method of selecting well-behaved $ \sigma $ [5, 36], we start with the curvature invariants Eqs. (3) and (4) directly, and choose strictly positive and finite functions as curvatures, solving $ \sigma $-functions inversely from the differential equations. The "strict positivity" of Ricci curvature is a minimum requirement for satisfying the DEC, whereas the "finiteness" is the regularity in $ r\in[0,\infty) $ and asymptotic flatness at infinity in their given sense. This process apparently provides an efficient approach to establishing physical RBHs, and it has indeed helped us create a series of models.

However, several problems in the process must be treated carefully. For example, the existence of horizons and the integrality of $ \sigma $. In other words, the $ \sigma $-functions obtained from this method do not always provide horizons for the models. Because these problems go beyond the scope of the current study, we leave them for future research.

Furthermore, in this study, we concentrate on not only physical RBHs, but also non-physical models within the framework of Einstein's gravity to focus on how the singularities of RBHs affect the spectrum of ANMQs.

The perturbation of a BH can be realized by laying a probing field onto the BH spacetime [47]. For a massless scalar field without backreaction, the perturbation reduces to the propagation equation of a scalar field,

In the spherically symmetric spacetime given by Eq. (1), we can substitute the ansatz $ \Psi = \mathrm{e}^{- \mathrm{i}\omega t} Y_l(\theta,\phi) r^{-1} \psi(r) $ and then select the radial component

which is a second order differential equation. Or in the canonical form, we have

where

Then, substituting the shape function, Eq. (2), the effective potential $ V $ can be represented in terms of $ \sigma $ and its derivatives

where

is a regular function at $ r=0 $ when $ r\in[0,\infty) $, because for RBHs, $ \sigma\sim O(r^3) $ as $ r\to 0 $ [5, 17, 21, 29]. Thus, for any BHs without degenerate horizons (that is, the horizons are the simple roots of the algebraic equation $ f(r)=0 $), one may find that $ r=0 $ and horizons are regular singular points of the master equation (Eq. (16)).

Furthermore, to diagonalize Eq. (16), we can define the tortoise coordinate by an indefinite integral,

or the solution of the differential equation

The master equation, Eq. (16), then becomes

The QNMs are the solutions of Eq. (22) under the following boundary conditions:

Thus, the high damping limit of QNMs will be defined as AQNMs, that is, $ |{\text{Im}(\omega)}|\gg |{\text{Re}(\omega)}| $ as the overtone number approaches infinity for BHs with asymptotic flatness (see, for example, Refs. [46, 47]).

Our goal in the current section is to find the AQNMs of Eq. (22) using an asymptotic analysis known as the monodromy method, in which the Stokes lines play a pivotal role. Therefore, before studying the AQNMs of specific models, let us take a close look at the geometry of Stokes lines. A detailed illustration of the monodromy method can be found in Refs. [44, 45, 53].

To understand the geometric aspects of Stokes lines, we first extend the radial coordinate $ r $ into the complex plane and explicitly rewrite it with its real and imaginary components $ x $ and $ y $. We can then recast Eq. (20) as

where the tortoise coordinate $ z $ is also complex in general. Based on this representation, the real and imaginary parts of $ z(x, y) $ can be regarded as two surfaces, and $ z(x, y)=C $ are two groups of contours on each surface, that is,

The Stokes lines are defined as the zero-contours, $ \text{Re}[z(x, y)]=0 $, on the first surface [44]. The Stokes lines must be symmetric with respect to the $ y $-axis, that is,

The one-form $ \mathrm{d} \text{Re}[z] $ corresponds to a vector field,

whereas its Hodge dual (with respect to the Euclidean metric) or the Stokes field [54]

is tangent to the contours at each point. In the following, we also use the component notation to represent the Stokes field, that is,

Moreover, it is easy to verify that both the divergence and curl of the Stokes fields are constantly zero, that is, $ \star \mathrm{d}\star (\star \mathrm{d} \text{Re}[z])=0=\star \mathrm{d}\,(\star \mathrm{d} \text{Re}[z]) $. In other words, the Stokes field $ \star \mathrm{d} \text{Re}[z] $ is a harmonic one-form [55] because $ z(r) $ is a holomorphic function, and its real and imaginary parts satisfy the Cauchy–Riemann equations.

The critical points of the Stokes field are defined as a set of points $ \{(x,y)|\star \mathrm{d} \text{Re}[z(x,y)]=0\} $, that is, the zeros of $ 1/f(r)=0 $ on the complex plane of $ r\in\mathbb{C} $. Meanwhile, the critical points of the field are also the singular points [56] of the Stokes lines. Then, one can immediately deduce a fact about the point $ r=0 $: the origin $ r=0 $ is a critical point for singular BHs but not for RBHs, which implies that the Stokes lines of singular BHs converge to (or emit from) zero, whereas for RBHs, the convergence points may not include $ r=0 $.

Let us first review two integrable examples^④, Schwarzschild and RN BHs. The two plots in Fig. 1 show the Stokes surfaces $ \text{Re}[z(x, y)] $, where the peaks on the surfaces correspond to the horizons. Fig. 2 is a Stokes portrait/diagram, which consists of the vector field $ \star \mathrm{d} \text{Re}[z] $ (blue), Stokes lines $ \text{Re}[z] $=0 (yellow-brown), and complex horizons (purple). The Stokes lines can then be understood as the integral curves passing through the critical points in the Stokes portrait. Meanwhile, the complex horizons must be located inside the closed parts of the Stokes lines.

Figure 1.  (color online) Stokes surfaces $ \text{Re}[z(x, y)] $.

Figure 2.  (color online) Stokes diagrams in Cartesian coordinates; $ \star \mathrm{d}\text{Re}[z(x, y)] $–blue vector fields, and the yellow-brown stream corresponds to the Stokes curves $ \text{Re}[z(x, y)]=0 $, where we set $ k=1/2 $ for RN BHs.

For the Schwarzschild BH with the shape function $ f=1-1/r $, the Stokes vector field reads

The origin $ r=(0,0) $ on the complex plane is a critical point, that is, this Stokes line has a self-intersection at the origin. Near this point $ (0,0) $, there is an approximation, ${\boldsymbol V}_{\rm{st}}\sim\left\{-y, -x\right\}$. The Stokes line of the Schwarzschild BH can be integrated out analytically,

which is a transcendental curve. Moreover, around the origin $ r=(0,0) $, Eq. (31) reduces to $ y=\pm x $, that is, there are $ 2\times 2=4 $ lines emitting from the origin, and the angle between any two adjacent lines is $ \pi /2 $.

For the RN BH with the shape function $f= (r-1)(r-k)/r^2$, where $ 1\ge k >0 $, the Stokes field is

The origin $ r=(0, 0) $ is also a self-intersecting point, near which there is an asymptotic behavior, ${\boldsymbol V}_{\rm{st}}\sim\left\{2 x y/k, (x^2-y^2)/k\right\}$. The Stokes line for this case is also a transcendental curve,

Around the origin $ r=(0, 0) $, $ x=0 $ and $ x=\pm \sqrt{3} \,y $, that is, there are $ 3\times 2=6 $ lines emitting from the origin, and the angle between any two adjacent lines is $ \pi /3 $.

Alternatively, to fix the angles between any two adjacent Stokes lines emitting from the critical point, we can also rewrite Eq. (28) in the polar coordinates $ \{\rho,\phi\} $

The angle relationship can then be directly visualized from the diagram in Fig. 3, that is, the intersection of the Stokes lines with $ \rho = 0 $ corresponds to the incident/exit angle of each line.

Figure 3.  (color online) Stokes diagrams in polar coordinates; $ \star \mathrm{d}\text{Re}[z(\rho, \phi)] $–blue vector fields, and the yellow-brown stream corresponds to the exact Stokes lines. The parameter $ k=1/2 $ is adopted in the RN BH.

From an analytical point of view, the Stokes line of the Schwarzschild BH in polar coordinates becomes

where $ \rho\ge 0 $, and $ 2 \pi \ge\phi\ge 0 $ (see the plot in Fig. 3(a)). To fix the angles of the emitting lines from origin, let us consider the corresponding vector field approaching $ \rho\sim0 $, which gives

The stability condition $ \mathrm{d}\rho/ \mathrm{d} \phi\equiv 0 $ implies that $ \phi= \pi n/2 $ with $ n\in \mathbb{Z} $. The angles of the emitting lines from origin can then be calculated as

Similarly, we can compute the Stokes line and emitting angle for the RN BH, arriving at

and $ \phi_{\rm{st}}= \pi n/3+ \pi /6 $, $ n\in \mathbb{Z} $ (see the plot in Fig. 3(b)).

Now, let us start with our first example of RBHs, the Bardeen BH [15], which can be generated by a nonlinear magnetic monopole. The Lagrangian of matter generation can be found in Refs. [17, 35],

where $ \mathcal{F}=F_{\mu\nu}F^{\mu\nu} $ is the field strength, and $ P $ and $ M $ are the magnetic charge and mass parameters, respectively. We first rescale the shape function by

such that it becomes

where all variables and parameters are dimensionless. The Stokes diagrams are displayed in Fig. 4.

Figure 4.  (color online) Stokes diagrams for the Bardeen BH; $ \star \mathrm{d}\text{Re}[z(x, y)] $–blue vector fields, and the yellow-brown curves correspond to the Stokes lines, where $ M=4 $.

Note that there are two critical points $ \pm \mathrm{i} $, that is, the singularities of the Bardeen BH are stretched from zero to the pure imaginary axis compared with the RN BH (see Fig. 2(b)). This can also be understood from the Weyl curvature $ W $, which is regular at zero but divergent at $ \pm \mathrm{i} $, that is, $ \pm \mathrm{i} $ are the poles of the seventh-order,

Therefore, $ \pm \mathrm{i} $ are the singular points of the master equation according to the effective potential, Eq. (18). Meanwhile, they are regular singular points because

are regular at the corresponding points, where $ p $ and $ q $ are the coefficients in Eq. (17) of the master equation for the Bardeen BH. Furthermore, around these two critical points, there is an approximation of the tortoise coordinate

and thus, there are five lines emitting from each critical point [57]. The polar diagram in Fig. 4(b) shows the angles of the Stokes lines passing through the zero point. Because $ r=0 $ is no longer a critical point, the angles at this point are trivial $ \pi $, the monodromy of the asymptotic solutions around $ r=0 $ is trivial.

To obtain information on the angles around the critical points $ \pm \mathrm{i} $, we should apply the polar coordinates with the starting point at each critical point separately. Polar diagrams with starting points at $ r= \mathrm{i} $ and $ r=- \mathrm{i} $ are shown in Fig. 5. The angle between adjacent lines can be computed as $ 2 \pi /5 $. Meanwhile, it must be emphasized that the line $ \phi= \pi /2 $ emits from the lower point, whereas $ \phi=3 \pi /2 $ emits from the upper point. This information will help us determine a closed contour, which is applied to calculate the monodromy.

Figure 5.  (color online) Stokes diagrams for the Bardeen BH in polar coordinates with starting points at each critical point.

For convenience of later discussions, let us denote the origin as $ \mathcal{O} $, the upper and lower critical points as $ \mathcal{A} $ and $ \mathcal{B} $, respectively, and a point close to the right inner horizon as $ \mathcal{C} $.

Now, let us start with the master equation at zero located on the Stokes line to calculate the AQNMs. Because the tortoise coordinate is approximately $ z\sim r $ at $ r=0 $, we have

which gives the solution around $ \mathcal{O} $

where $ \nu=2l+1 $, and $ Q_\pm $ are two arbitrary constants. Then, considering the asymptotics of the Bessel functions,

with $ \alpha_\pm= \pi \left(1\pm \nu\right)/4 $, we find that at $ - \mathrm{i} \infty $ ,

and the condition

where the boundary condition $ \psi\to \mathrm{e}^{- \mathrm{i} z\omega} $ as $ z\to \infty $ is used. To pass through $ \mathcal{O} $, we must rotate $ \pi $ and hence obtain

which will be matched with the asymptotics at point $ \mathcal{A} $

Here, we use $ \hat{\mathcal{O}} $ to denote the position after rotation around $ \mathcal{O} $. This notation also applies to the following discussion. The matching $ \psi_{\mathcal{A}}\leftrightarrow\psi_{\mathcal{\hat{O}}} $ provides two more conditions:

To pass through $ \mathcal{A} $, we also rotate $ \pi $ instead of $ 2 \pi $, which provides

and the asymptotics as $ z $ approaches the inner horizon are

where $ \delta = \mathrm{i} /(2 T_{\rm{H}}^-) $, and $ T_{\rm{H}}^- $ is the "temperature" of the inner horizon. We use

because $ \omega (z-\delta) $ is negative on this branch. Then, the matching $ \psi_{\mathcal{\hat A}}\leftrightarrow\psi_{\mathcal{C}} $ gives

Subsequently, we rotate to the branch containing $ \mathcal{B} $. We find

which should be matched with the asymptotics at $ \mathcal{B} $

The matching $ \psi_{\mathcal{\hat C}}\leftrightarrow\psi_{\mathcal{B}} $ leads to another two conditions:

We then pass to the origin $ \mathcal{O} $ again but in a different branch. To distinguish the solution of the master equation with the original one, Eq. (46), let us denote it as $ \widetilde \psi_{\mathcal{ O}} $, that is,

which should be matched with the asymptotics at $ \mathcal{B} $, $ \psi_{\mathcal{ B}}\leftrightarrow\psi_{\mathcal{\widetilde O}} $

Finally, after rotating $ \pi $ around the origin, we arrive at

Then, as we close the counter and compare it with the change around the outer horizon, we obtain the condition

where $ T_{\rm{H}}^+ $ is the temperature of the outer horizon. The condition of asymptotic QNMs is computed by combining Eqs. (49), (52), (56), (59), (61), and (63), which gives

or

where $ n\in \mathbb{Z} $. This result is completely different from that obtained in Ref. [32] for two reasons. First, the origin $ r=0 $ is no longer a singular point of the Stokes lines in our calculation, and second, the trajectory of the asymptotic solution along the Stokes lines is different from that of the RN BH.

Before we start a new model, let us comment on the above process. First, the origin point $ r=0 $ is a common point of the Stokes lines for the Bardeen BH (that is, there is only one incoming line and one outgoing line), although it is a regular singular point of the master equation; therefore, when we follow the contour around it, the rotation of the solution is trivial, and thus the multipole number $ l $ does not contribute to the spectrum of AQNMs.

Second, to calculate the monodromy of the solution, we start the asymptotics, Eq. (46), at the point $ r=0 $. In fact, we can start with any point on the Stokes lines, and the result will be the same. Nevertheless, in practice, if one begins the analysis at $ r=+ \mathrm{i} $, the master equation becomes

where $ V_0=3 \sqrt[5]{-1} M^{2/5}/5^{8/5} $, and $ z\sim -4\left(1- \mathrm{i}\right) (r- \mathrm{i})^{5/2}/ (5M) $. This equation is not solvable because the recursion relation obtained using the Frobenius method is not solvable as a difference equation. Therefore, the asymptotic analysis of the solution becomes complicated.

To provide an effective analysis, one can use the perturbation method by considering the following equation with a perturbative parameter $ \epsilon $

the zero-order $ O(\epsilon^0) $ of which can be applied in the calculation of the monodromy. Owing to this approximation, a deviation from the result in Eq. (64) is allowed.

Third, there is no linear dependence of $ \omega $ on $ n $, which is considerably different from the case of traditional BHs with singularities at the centre, such as Schwarzschild and RN BHs (see Ref. [44]).

The second example refers to the Hayward BH [2]. Though the original model is proposed following Markov's idea based on the limiting curvature conjecture [4], it can also be reparameterized and interpreted in the context of nonlinear electrodynamics[17]. We apply the latter in our discussion. The matter Lagrangian reads as [17]

With the help of the coordinate transformation

the shape function becomes

After analytically continuing $ r $ into the complex plane, we find that the Hayward BH has three horizons. The three critical points are $ -1 $, $ (-1)^{1/3} $, and $ -(-1)^{2/3} $, and the Stokes diagrams are shown in Fig. 6, where the Stokes lines are separated into two groups. One crosses the critical point $ -1 $, and the other crosses both $ (-1)^{1/3} $ and $ -(-1)^{2/3} $ simultaneously; however, neither crosses the zero point. Furthermore, because $ r=-1 $ corresponds to a temperature with a negative BH radius, we calculate the monodromy based on Fig. 6(b). In other words, we follow the Stokes line crossing the points $ (-1)^{1/3} $ and $ -(-1)^{2/3} $.

Figure 6.  (color online) Stokes diagrams for the Hayward BH; $ \star \mathrm{d}\text{Re}[z(x, y)] $–blue vector fields, and the yellow-brown stream corresponds to the exact Stokes curves, where $ M=4 $.

Information on the angles at each critical point is shown in Fig. 7, that is, the angle between two adjacent lines is $ \pi /2 $. The tortoise coordinate around $ (-1)^{1/3} $ and $ -(-1)^{2/3} $ are, respectively,

Figure 7.  (color online) Stokes diagrams of the Hayward BH in polar coordinates, with starting points at each critical point.

Because $ r=0 $ is not a point on the Stokes line, we start with the lower critical point $ r=(-1)^{1/3} $. The master equation then becomes

where $ V_0=-\sqrt[6]{-1} \sqrt{M}/(2 \sqrt{6}) $, and its coefficient multivalues owing to $ z^{3/2} $. To perform the asymptotic analysis, we apply the perturbation method based on a Bessel-type equation, as mentioned at the end of previous subsection, that is, we recast the master equation with

and expand the wave function as a power series of $ \epsilon $, $ \psi(\epsilon, z)=\psi_0(z)+\epsilon \psi_1(z)+\dots $, where $ \epsilon $ is the perturbative parameter. Then, we can obtain the perturbative equations order by order

Next, to calculate the monodromy, we only use the zero-order equation and drop the order subscript of the wave function to simplify the notation. In other words, our starting point is the solution of the zero-order perturbative equation

where $ \nu=\sqrt{1+4V_0} $ is a complex number. Moreover, as with the Bardeen BH, we denote the upper and lower critical points as $ \mathcal{A} $ and $ \mathcal{B} $, respectively, and a point close to the middle horizon as $ \mathcal{C} $.

From the lower critical point, we can obtain an asymptotic behavior of the wave function in Eq. (75),

and the condition

where $ \alpha_\pm= \pi (1\pm \nu)/4 $. Then, with the matching $ \psi_{\mathcal{\hat B}}\leftrightarrow\psi_{\mathcal{C}} $, the asymptotics after rotating around $ \mathcal{B} $ with that around $ \mathcal{C} $ give another two conditions:

where $ \delta = \mathrm{i} /(2 T_{\rm{H}}^-) $, and $ T_{\rm{H}}^- $ is the "temperature" of the inner horizon.

With the matching $ \psi_{\mathcal{\hat C}}\leftrightarrow\psi_{\mathcal{A}} $, the asymptotics of the wave function at $ \mathcal{C} $ after rotating around the upper critical point with that at $ \mathcal{A} $ provide two more conditions:

where $ \sqrt{3} $ originates from a shift in the variable in Eq. (75) from the lower critical point to the upper, that is, $ z\to z+\sqrt{3} $, because the difference between the upper critical point and the lower is $\text{Im}\left[(-1)^{1/3}\right] - \text{Im}\left[-(-1)^{2/3}\right]=\sqrt{3}$.

Finally, after closing the contour and comparing it with the change around the outer horizon, we arrive at

where $ T_{\rm{H}}^+ $ denotes the temperature of the outer horizon. Furthermore, combining all the above conditions, we obtain an analytical expression for the AQMNs,

which has a similar form to the result in Ref. [32] because the track along the Stokes lines is similar to that of the RN BH. However, there is an essential difference: in our result, the AQNMs should not depend on the multipole number $ l $ because $ r=0 $ is not a singular point of the Stokes lines, and thus rotation of the asymptotic solution around $ r=0 $ is trivial, and $ l $ does not appear in the formula of the AQNMs.

For the third class of RBHs generalized by the nonlinear electrodynamic source in Ref. [17], we simply choose $ \mu=3 $, yielding the matter Lagrangian [17]

and the shape function

and use the transformation $ r\to P r $ and $ M\to M P/2 $ to rescale the variables, such that all the parameters appearing in the shape function are dimensionless. Furthermore, there is one critical point $ r=-1 $ and three horizons in this model (see Fig. 8), and from the critical point, there are eight Stokes lines emitting with adjacent angles of $ \pi /4 $.

Figure 8.  (color online) Stokes diagrams for the third class of RBHs; $ \star \mathrm{d}\text{Re}[F(x, y)] $–blue vector fields, and the yellow-brown curves correspond to the Stokes lines, where $ M=7 $.

Approaching the critical point, the leading term of the tortoise coordinate is

and thus the master equation becomes

where $ V_0= 3 \sqrt[4]{-M}/(8 \sqrt{2}) $. The same as before, we apply perturbation theory and take the zero-order solution as the asymptotics of the master equation at the critical point. This procedure gives us a similar Bessel-type function.

On the other side, the closed contour starts from the bottom left of Fig. 8(a) and approaches the critical point before rotating $ \pi /2 $ in $ r $ (that is, $ 2 \pi $ in $ z $). Subsequently, it moves toward the middle horizon, returns to the critical point, and then finally reaches the top left corner by rotating $ \pi /2 $ in $ r $.

The expression for the AQNMs of this third class is formally the same as that of the Hayward BH (Eq. (81)), but with a different value of $ \nu $. This implies that the matching of $ \psi_{\mathcal{\hat C}}\leftrightarrow\psi_{\mathcal{A}} $ in the case of the Hayward BH can be ignored. Moreover, this reflects the fact that a transformation, such as Eq. (69), should not affect the physical results, even though the numerical values of the upper and lower critical points depend on the scale transformations.

From the above examples, we clearly see that the AQNMs have the same form as long as the traces of the asymptotic solutions along the Stokes lines are similar. Let us analyze one more integrable example.

Inspired by [58], we construct the following shape function:

where $ P $ is interpreted as magnetic charge, as before. The complete roots of $ f(r)=0 $ include $ \pm P $ and $ \pm\sqrt{P^2+4M^2} $, among which only the positive are physical. The AQNMs of a similar type of RBH have been considered in Refs. [59, 60].

The physical radius of this BH cannot be less than $ P $, otherwise the metric becomes complex. Even though the curvature invariants are mathematically divergent at the zero point, the physical condition $ r\ge P $ forbids any test particle from passing into the inner horizon $ r_-=P $. Therefore, all the curvature invariants are finite in the physical domain $ r\in [P,\infty) $. Moreover, the origin is the only critical point, and the tortoise coordinate can be obtained analytically,

where $ z_0=-1/2+\ln (4)+ \mathrm{i} \pi $, and we use the normalization $ r\to 2 M r $ and $ P\to 2 M P $ and set $ P=1 $. Thus,

which implies that there are twelve Stokes lines emitting from the origin, and the angle between any two adjacent lines is $ \pi /6 $ (see Fig. 9). The master equation is

Figure 9.  (color online) Stokes diagrams for the PK BH; $ \star \mathrm{d}\text{Re}[F(x, y)] $–blue vector fields, and the yellow-brown curves correspond to the Stokes lines, where $ M=1/2 $ and $ P=1 $.

Thus, the monodromy relation is similar to that for the RN BH,

where $ \nu=\sqrt{14}/3 $. Substituting the explicit forms of the temperatures of the inner and outer horizons,

we arrive at

As a step summary, we have shown that the analytical form of AQNMs may be universal if the trajectories of the asymptotic solutions along the Stokes lines are similar; however, they should not depend on the multipole number $ l $.

The models considered in Sec. III have a common aspect: the singularities of RBHs still exist (defined as the second type), but they are simply moved to the classically forbidden region. In the current section, we investigate several exotic examples that belong to the first and third types discussed in Sec. II.

Let us express the shape function with an essential singularity at $ r=0 $ as [21, 61]

where $ P $ is interpreted as magnetic charge. The Lagrangian of the matter source as a nonlinear magnetic monopole can be found as

To perform an appropriate analysis, we use a dimensionless representation with the help of the rescaling transformation [28]

The shape function then becomes

When $ r\in \mathbb{R} $, there are two real horizons $r_{\rm{H}}= -1/W_n\left(-P^{-2}\right)$ with $ n=0, -1 $ if $ P^2> \mathrm{e} $. Here, $ W_n(z) $ is the Lambert W function [62]. When $ r $ is analytically continued into the complex plane, $ n $ is enlarged to all natural numbers $ \mathbb{N} $, that is, the roots of $ f(r_{\rm{H}})=0 $ are infinitely many according to Picard's great theorem [50] because $ r=0 $ is an essential singularity of $ f(r) $. Meanwhile, all of these roots except the physical horizons are bounded by $ |{r_{\rm{H}}}|\le |{-1/W_0\left(-P^{-2}\right)}| $ and located symmetrically with respect to the real axis because $ W^*_n\left(-P^{-2}\right)=W_{-n-1}\left(-P^{-2}\right) $ if $ n\neq 0, -1 $.

To find the critical points, it is equivalent to calculate the zeros of $ r \mathrm{e}^{1/r} $ because

but $ r \mathrm{e}^{1/r} $ does not have any zeros on the complex plane, that is, the zero is a Picard exceptional value of $ 1/f $. Thus, there are no critical points for this model. However, the origin $ r=0 $ is a jump discontinuity,

which signifies that it may play a special role in the Stokes diagram, even though it is not a critical point. The superscript in $ 0^+ $ ($ 0^- $) indicates that the path of the limit to $ 0 $ is taken in the right (left) half-plane. From the perspective of the Weyl curvature

the jump discontinuity of $ r=0 $ implies that

that is, the Weyl curvature is divergent at the essential singularity $ r=0 $ as $ r $ approaches zero from the left half-plane. Thus, we dub $ r=0 $ as the semi-critical point because it is different from regular points and the critical points considered in Sec. III.

To perform a quantitative analysis of the Stokes lines, we apply the asymptotic relation $ 1/f(r)\sim - r \mathrm{e}^{1/r}/P^2 $ as $ r\to0^- $ and rewrite Eq. (21) as

which provides the solution

where $ {\rm{Ei}} $ denotes the exponential integral [62]. Thus, the Stokes lines $ \text{Re}[z]=0 $ around $ 0^- $ are approximately

Along the longitudinal direction, $ x=0 $, this reduces to

As a result, we find that there are infinite Stokes lines quantified by $ 1/y=n \pi $ with $ n\in\mathbb{Z}/\{0\} $ owing to the periodicity of the sine function (see Fig. 10). This phenomenon of infinitely many Stokes lines around the essential singularity $ r=0 $ is also a reflection of Picard's great theorem.

Figure 10.  (color online) Stokes diagrams in Cartesian coordinates for the Balart–Vagenas BH close to $ 0^- $. The yellow-brown curves are Stokes lines, whereas the blue points are their intersections with the longitudinal axis $ y $.

To fix the angle between two adjacent Stokes lines at $ 0^- $, we rewrite Eq. (103) in the polar system $ \{\rho, \phi\} $, which gives

When $ \rho $ becomes small, we find that

That is, $ \phi= \pi $ as $ r\to 0^- $. In other words, the angles of all Stokes lines emitting from $ 0^- $ are zero (see Fig. 11).

Figure 11.  (color online) Stokes diagrams in polar coordinates for the Balart–Vagenas BH close to $ 0^- $.

The full Stokes lines are computed using a numeric integral and are depicted in Fig. 12. Because of the extreme characteristic of the integrand near the singularities, the lines are not smooth at all, as shown in the plot. Alternatively, we develop a formula based on the Cauchy–Riemann equations [50] to estimate the global feature of Stokes lines in the range far from the origin,

Figure 12.  (color online) Stokes diagrams for the Balart–Vagenas BH; $ \star \mathrm{d}\text{Re}[z(x, y)] $–blue vector fields, and the yellow-brown curves correspond to the Stokes lines, where $ P=2 $.

which loses its validity as $ x\to 0^- $ because an infinite number of complex horizons are located there (see Fig. 13).

Figure 13.  (color online) Improved Stokes lines for the Balart–Vagenas BH based on the numeric integral in Eq. (107) with $ P=2 $.

Now, we turn to the master equation. First, we obtain the asymptotic form of the potential in Eq. (18) as $ r\to 0^- $

and then with the help of the asymptotic relation in Eq. (102), we can rewrite the master equation as

However, because the angles of the Stokes lines emitting from the origin are zero, a trivial monodromy relation is obtained, that is, $ 1=\exp(\omega/T^+_{\rm{H}}) $, which gives a pure imaginary spectrum of the AQNMs, $ \omega=2 \pi \mathrm{i} n T^+_{\rm{H}} $.

Next, we turn to an RBH model inspired by noncommutative geometry [63],

where $ \epsilon $ denotes the noncommutative parameter, and $ \gamma\left(\dfrac{3}{2},\dfrac{r^2}{4\epsilon}\right) $ is the lower incomplete gamma function [62]. The original model was proposed by considering the existence of a minimum distance scale. However, we use the monopole interpretation instead of the original model, that is, we replace the parameter $ \epsilon $ with $ P^2 $. Thus, the Lagrangian can be obtained as

Then, with the help of the following transformation [28]:

the shape function becomes

where the parameter $ P $ is bounded by $ 0.93\gtrsim P>0 $, otherwise there will be no physical horizons.

The complex extended horizons can be calculated numerically using $ f(r)=0 $. The distribution of complex horizons depends on the value of $ |{r}| $, which can be illustrated by the argument principle^⑤; see Fig. 14, where we show winding indices for the images of circles with different moduli $ |{r}| $. The number of complex horizons increases with increasing modulus $ |{r}| $. Moreover, we can find a property of the horizons, that is, if $ r $ is a horizon, its conjugate $ r^* $ is also a horizon because $\gamma^*(3/2, z^2)= \gamma(3/2, (z^*)^2)$, and thus, $ 0=f^*(r)=f(r^*) $ (see Fig. 15(a)).

Figure 14.  (color online) Winding indices. The images (black closed curves) of circles with radii $ |{r}|=1/2, 1, 3, 4 $ from left to right under the map of Eq. (113), where $ P=1/2 $. The yellow-brown rays from zero to infinity are auxiliary lines for calculating the winding indices. The signs "$ \pm $" correspond to the counter-clockwise and clockwise intersections of the curve with the ray, respectively.

Figure 15.  (color online) Stokes diagrams for the Eq. (110) BH; $ \star \mathrm{d}\text{Re}[F(x, y)] $–blue vector fields with $ P=1/2 $.

As the modulus of $ r=x+ \mathrm{i} y $ becomes large, the complex horizons can be approximately computed by solving the system

which gives

with $ n\in\mathbb{Z} $. We can also verify that there are neither critical nor semi-critical points in this model because $ \gamma\left(3/2, r^2\right)/r $ is an entire function. In particular, $ 1/f $ approaches $ 1 $ as $ r\to 0 $. Therefore, the monodromy of the asymptotic solutions along any closed curve is trivial, which leads to a similar result to that of the above example, $ 1=\exp(\omega/T^+_{\rm{H}}) $.

However, if we consider the small-charge approximation $ P\sim 0 $, that is, if we treat the RBH model as a singular one, we obtain

which has a critical point at zero and one real horizon. In other words, the Stokes diagram of this model in the small-charge approximation is similar that of the Schwarzschild BH. The tortoise coordinate and effective potential close to zero are, respectively,

The master equation becomes

and the monodromy relation is then

with

As we find in this section, the analytical spectrum of the AQNMs reduces to an extreme simple form because the monodromy of the asymptotic solutions along any closed Stokes lines is trivial.

To study the AQNMs using the complex WKB approach, Andersson and Howls used another method to diagonalize the differential operator in the master equation, where the radial coordinate, unlike the tortoise coordinate, is no longer mutivalued [64]. The diagonalized master equation in Andersson and Howls' approach is

with

where $ V(r) $ is the same as in Eq. (18), and $ \Phi $ connects with $ \psi $ via $ \Phi=f^{1/2}\psi $. The Stokes lines are defined by^⑥

Generally, the integral in Eq. (123) cannot be calculated analytically for RBHs. Therefore, we use the WKB Stokes field $ \{\text{Re}[Q],-\text{Im}[Q]\} $ to depict the Stokes portrait; see Fig. 16, where the Stokes field for Schwarzschild and RN BHs are shown.

Figure 16.  (color online) Stokes field $ \{\text{Re}[Q],-\text{Im}[Q]\} $. The purple points denote horizons, whereas the red points are zeros of $ Q $.

We use the term "WKB Stokes field" to distinguish it from the one we discuss in the monodromy method. However, for a given BH model, these two fields are related. To see it clearly, let us take the damping limit $ \omega\to- \mathrm{i} \infty $ in the WKB Stokes field, which gives us $ Q\sim - \mathrm{i} |{\text{Im}(\omega)}|/f(r) $ for $ r\ge r_{\rm{H}} $. Thus, the WKB Stokes field in the damping limit becomes

which is the scaled Stokes field in the monodromy method, Eq. (28). Furthermore, for singular BHs, the critical points (zeros) of $ Q $ converge to the origin (essential singularities) of BHs as $ \text{Im}(\omega)\to \infty $, and the WKB Stokes field becomes the Stokes field in the monodromy method.

The situation for RBHs is slightly different from the above case for singular BHs. $ Q_0^2 $ of RBHs has two zeros around the origin $ r=0 $ owing to $ V\sim r^{-2} $ (see Fig. 17), even though $ r=0 $ is not a zero of $ f^{-2} $. However, these two zeros do not converge to the origin but disappear as $ \omega \to - \mathrm{i} \infty $ because $ \omega $ dominates in $ Q_0^2 $ at that moment. As a result, $ r=0 $ becomes a regular point in the damping limit. In other words, the AQNMs via the complex WKB approach, like in the monodromy method, are not dependent on the behavior around the BH center.

Figure 17.  (color online) Stokes field $ \{\text{Re}[Q],-\text{Im}[Q]\} $. The purple points denote horizons, whereas the red points are zeros of $ Q $.

In this study, to calculate the AQNMs of RBHs, we apply the Stokes field to extensively investigate the local characteristics of the Stokes lines and classify RBHs based on the types of their complex singularities, which arise as a result of the analytical continuation of the radial coordinate into the complex plane.

From the calculation, we find several novel aspects of the AQNMs of RBHs, that is, the analytical forms of the asymptotic frequency spectrum are not universal for spherically symmetric RBHs with single shape functions and do not depend on the multipole number $ l $ because $ r=0 $ is not a point on the Stokes lines; even if $ r=0 $ is a point on the Stokes lines, it must not be an accumulation point, such that the rotation of the asymptotic solution around $ r=0 $ is trivial and $ l $ does not appear in the formula of the AQNM.

The forms of asymptotic frequency may depend on the structures of the Stokes portraits.

1. The existence of singularities. The absence of singularities leads to a trivial monodromy, such that the asymptotic frequency is purely imaginary.

2. The rotation angle of the asymptotic solutions around singularities. If the rotation angle is trivial, the asymptotic frequency also has only the imaginary part.

3. The trajectory of the asymptotic solutions along the Stokes lines. It is not the topology of the Stokes line that plays a decisive role, but the way of bypassing the trajectory.

In a broader sense, specific research on AQNMs relates to several mathematical topics, such as transcendental curves and the value distribution of holomorphic functions. As shown in this study, even the simplest Stokes lines obtained from the Schwarzschild BH cannot be depicted by a polynomial, whereas the Stokes lines for RBHs are usually not integrable and have no analytical expressions.

On the one hand, to provide aspects of the Stokes lines, we attempt several numerical methods in this study to remove the integral of the tortoise coordinate, for example, Newton–Cotes quadrature rules. Nevertheless, the results obtained from these numerical methods more or less lose some important information on the Stokes lines.

On the other hand, the situation of the complex singularities (curvature and coordinate singularities) of RBHs becomes intricate because the Stokes lines must have self-intersections at curvature singularities, and their closed parts must surround the coordinate singularities. Thus, the distribution of the zeros and poles of the shape functions directly affect the appearance of the Stokes lines. Clarification of the value distribution may help us construct information on the Stokes lines. Therefore, we will devote ourselves to developing more effective approaches to obtaining the features of Stokes lines in future studies.

Finally, because our main motivation for studying the AQNMs is to obtain the quantum entropy spectrum of RBHs, our subsequent work will focus on how to derive the correct quantum entropy spectrum based on the AQNMs obtained in this study.