A regular black hole as the final state of evolution of a singular black hole

We propose a novel black hole model in which singular and regular black holes are combined as a whole and more precisely singular and regular black holes are regarded as different states of parameter evolution. We refer to them as singular and regular states, respectively. Furthermore, the regular state is depicted by the final state of parameter evolution in the model. We also present the sources that can generate such a black hole spacetime in the framework of F(R) gravity. This theory of modified gravity is adopted because it offers a possible resolution to a tough issue in the thermodynamics of regular black holes, namely the discrepancy between the thermal entropy and Wald entropy. The dynamics and thermodynamics of the novel black hole model are also discussed when a singular state evolves into a regular state during the change of charge or horizon radius from its initial value to its extreme value.

The definition of an RBH is twofold.One is based on the coordinate-independent approach, known as finiteness of curvature invariants [1,10,11], i.e., the RBH is defined as a BH spacetime with finite curvature invariants.This definition is related to Markov's limited curvature conjecture [25].The other is based on the coordinate-dependent approach, called geodesic completeness [26,27], i.e., the RBH is defined as a BH spacetime with complete null and timelike geodesics.These two definitions are not equivalent generally, which means that some models satisfy one of the definitions but violate the other [28,29].
In addition, we have to clarify two concepts of geodesic completeness so that the reader can understand our subsequent discussions unambiguously.At first, the "completeness" implies that a test (null or timelike) particle arriving at any point of spacetime requires an infinite affine parameter.Secondly, the "coordinate dependence" means that it is possible for the affine parameter to be finite for a certain choice of coordinates even if this spacetime is regular everywhere.
In order to clarify specifically the "coordinate dependence" mentioned above, we take two examples.One is the Rindler spacetime which can be described [30,31] by the line element, ds 2 = −z 2 dt 2 + dx 2 + dy 2 + dz 2 .The affine parameter is finite as a test particle approaches to z = 0, i.e., it seems that the Rindler spacetime is "incomplete" at z = 0 in this choice of coordinates.Nevertheless, it is known that the Rindler spacetime is a part of the Minkowski spacetime, thus it cannot include any incomplete points.In fact, the Rindler spacetime will convert to the Minkowski spacetime via an appropriate choice of coordinates.The other example is the Hayward BH spacetime [32], where the affine parameter is finite at r = 0 in a certain choice of coordinates, seeming that the Hayward BH is "geodesic incomplete".However, the Hayward BH is regular if the Painlevé-Gullstrand coordinate is chosen.In Sec.2.2, we will demonstrate this point and note that such a coordinate should be adopted when we analyze the geodesic completeness at the center of RBHs with the spherical symmetry and one shape function.
The differences in dynamics and/or thermodynamics between RBHs and SBHs are a topic of great interest in the field of RBHs [33][34][35][36][37][38][39][40][41][42][43][44][45].By investigating these differences, one expects to study the effects of singularities on the dynamics and/or thermodynamics of BHs.In addition, RBHs are generally considered to be the product of quantum gravity.From this point of view, if a BH with a singularity at the center is regarded as a classical state, the BH evolves to its final state through self-heating radiation, where the final state is a quantum state without a singularity at the center.That is to say, RBHs can be regarded as the final state of parameter evolution of SBHs.To this end, we construct a novel model which combines RBHs and SBHs into one whole.This proposal can highlight the differences in dynamics and/or thermodynamics between RBHs and SBHs, where an RBH appears as the final state of parameter evolution of an SBH.It provides a basis for the study of singularities' effects when an SBH evolves into an RBH.The idea to combine RBHs and SBHs into one whole is essentially similar to the situation in Refs.[14][15][16][17].In these references, the BHs are in singular states when M ̸ = 0 and in regular states when M = 0, i.e., the entire contribution to the ADM mass of BHs comes from electromagnetic interactions.
The paper is organized as follows.At first, we prove in Sec. 2 that it is equivalent to determine RBHs by using complete geodesics or finite curvature invariants for the BHs that are spherically symmetric and have only one shape function.Then, in Sec. 3 we propose a spherically symmetric BH model with one shape function, which combines singular and regular BHs as a whole.In particular, an RBH can be regarded as the final state of an SBH when the charge or horizon radius goes to its extreme value.
The matter sources that generate such a BH are established in terms of modified gravity F (R) rather than Einstein's gravity in Sec. 4. The motivation comes from the need for the thermodynamics of RBHs, i.e., we want to make the Wald entropy in F (R) gravity consistent with the entropy calculated by the normal thermodynamic formula, dM/T , where M and T are BH mass and temperature, respectively.As physical applications, we investigate the dynamics and thermodynamics of our BH model in Sec. 5 and Sec.6, respectively.The purpose of these two sections is to analyze the essential changes in physical properties that take place during the evolution of an SBH into an RBH.In addition, as a by-product, we also discuss the differences in the interpretation of RBHs under F (R) gravity and Einstein's gravity.
The conclusions with future outlooks are summarized in Sec. 7, whereas App.A includes the calculations of monodromy that are applied in Sec.5.2 and App.B gives the actions of matter sources Eq.(31) explicitly.

Criterion of regular black holes
The condition of finite curvature invariants is widely applied to determine [1,10,11,13,14] whether a BH is regular or not, which relates to the Markov limiting curvature conjecture [7,25,46,47].Moreover, the other condition to ensure a non-singular spacetime is the completeness [30,48] of timelike and null geodesics.These two conditions are generally not equivalent [26][27][28][29]32].In this section, we shall show that the finiteness of curvature invariants and the completeness of geodesics are equivalent for the spherically symmetric BHs with one shape function.

Finite curvature invariants
Let us begin with the metric with a spherical symmetry, where the shape function reads Traditionally, the Kretschmann scalar K := R µναβ R µναβ , Ricci scalar R := g µν R µν , and contraction of two Ricci tensors R 2 := R µν R µν are applied to check [13,14] the regularity of BHs.These three invariants are connected by a syzegy [49] (also known as Ricci decomposition) as follows, where C := C µναβ C µναβ is contraction of two Weyl tensors.Alternatively, R 2 and K can be replaced [7] by S and C, where S is defined by S := S µν S µν and S µν by S µν := R µν − 1 4 g µν R. Fourteen independent curvature invariants can be established [50] based on the Riemann curvature in the 4-dimensional spacetime, while a complete set of curvature invariants includes seventeen elements which are called Zakhary-Mcintosh (ZM) invariants [49,[51][52][53].Therefore, it is natural to wonder whether the set of three invariants {K, R 2 , R} or {C, S, R} can represent the finiteness of all curvature invariants.If not, what elements should be included in the minimum set in order to determine an RBH?To answer this question for the metric Eq. ( 1), we calculate the seventeen invariants, where six of them vanish.The non-vanishing eleven invariants can be separated into three groups.
1. Ricci type constructed by Ricci tensors: Weyl type constructed by Weyl tensors: 3. Mixed type constructed by both Ricci and Weyl tensors: where we have defined The three quantities A 1 , A 2 and A 3 in Eq. (7) determine the finiteness of all non-vanishing invariants.Furthermore, I 6 , as a specific representation, is necessary to determine if a BH is regular when σ(r) is not a constant.The reason is that the squares of A 1 and A 2 in I 6 guarantee that no divergence terms appear.For more general cases, we can combine another invariant, for instance, I 1 with I 6 .Since A 1 and A 2 are finite, A 3 is also finite if 6M σ/r 3 is finite.In summary, I 1 and I 6 , as representatives, can fully determine the regularity of the BHs with the spherical symmetry and one shape function.We note that the number of elements in the complete set of spherically symmetric BHs is less than that of rotating RBHs [54].
A similar result can be obtained if the Kretschmann scalar If I 6 is finite, the regularity can be determined fully by K, i.e., {K, R 2 } is a complete set to determine whether a spherically symmetric BH with one shape function is regular.Furthermore, the requirement of a BH being regular at center r = 0 can be obtained In other words, σ(r) approaches zero not slower than r 3 .Alternatively, a BH is regular at the center if we have lim This behavior of σ around r = 0 will help us to construct the RBHs in the next section.

Complete geodesics
In order to determine the condition that guarantees the geodesics being complete at r = 0,1 we start with the effective Lagrangian of null/timelike geodesics [55], where τ is an affine parameter along a geodesic.For the spacetime given by metric Eq. ( 1), the Lagrangian reads which does not depend on t and φ explicitly, i.e., t and φ are cyclic coordinates.Thus, the canonical momenta conjugate to t and φ are conserved quantities, where E and L are the total energy and total angular momentum, respectively.To simplify the following discussion, we just concentrate on the motions on an equatorial plane, θ = π/2, and set the initial angular momentum to be zero, L = 0. Next, we introduce the Painlevé-Gullstrand time [56], and recast the Lagrangian Eq. ( 12), Moreover, the definition of the Painlevé-Gullstrand time provides where we have replaced ṫ by E via Eq.( 13).Thus, we can obtain the equation of radial motion by using Eqs.( 15) and ( 16) together with 2L * = k, where k = 0, −1 corresponds to null and timelike geodesics, respectively, and we consider only the ingoing motion.
For timelike geodesics with k = −1, if we suppose that the test particle starts with dr/dt * = 0, i.e., E = 1, Eq. ( 17) can be reduced [57] to then the Painlevé-Gullstrand time can be calculated which is divergent at r = 0 if σ has the same asymptotic behavior around r = 0 as Eq. ( 9).In other words, the timelike geodesics are complete at r = 0 if σ ∼ O(r n ) with n ≥ 3, which is consistent with the requirement of finite curvatures.
For null geodesics with k = 0, the initial velocity should be dr/dr * = −1 at infinity, which also corresponds to E = 1.Thus, Eq. ( 17) becomes The Painlevé-Gullstrand time of null geodesics is divergent if σ ∼ 1 as r → ∞, which also meets the requirement of asymptotic flatness, i.e.
The null geodesics are complete for a massless test particle starting at infinity.However, this cannot provide any constraints to σ(r) in the vicinity of a BH center.In summary, we have proved that the finiteness of curvature invariants and the completeness of geodesics are equivalent for the spherically symmetric BHs with one shape function.

Construction of our model
We propose a shape function according to the asymptotic behavior given by Eq. ( 9), where δ(Q) is a function of parameter Q and vanishes when Q reduces to an extreme value Q = Q ext which corresponds to the extreme horizon r ext , meanwhile we have performed a transformation, r → r/(2M ) and Q → Q/(2M ), such that all variables are dimensionless.The existence of an extreme horizon indicates that this BH model has two horizons at least.Furthermore, from f (r ext , Q ext ) = 0 and δ(Q ext ) = 0, we derive the extreme horizon radius and charge, respectively.Thus, we select δ to be Actually, δ(Q) as a function of Q is not unique, and it is a relatively simple one we have selected.Alternatively, δ is a function depending on the difference between outer and inner horizons, r + −r − , and vanishes when the BH reduces to its extreme configuration, r + = r − = r ext , e.g., δ = r + − r − .
The horizon curve f (r H , Q) = 0 is shown in Fig. 1, where it clearly shows that there are two horizons for a given value Q ≤ Q ext .The maximum value of outer horizon r + is which can be normalized to unity by a certain transformation.It is not necessary to make such a normalization because Eq. ( 22) does not reduce to the shape function of Schwarzschild BHs when Q = 0.  22)- (24).The corresponding curvature invariants around r = 0 are which are divergent if Q ̸ = Q ext .Therefore, the BH that we constructed has curvature singularity at r = 0, that is, this BH is located at its singular state when Q ̸ = Q ext .As Q approaches the extreme value Q ext , the behaviors of curvature invariants around r = 0 become which implies that the curvature invariants are finite and positive at r = 0. Thus, the BH is located at its regular state when Q = Q ext .The transition of BHs due to the jump in the parameter Q in our model is similar to the changes between BHs and wormholes discussed in Refs.[58,59].In these two references, however, the change in the parameter a or m does not have to correspond to a physical process.In other words, the BHs and wormholes are more likely to be regarded as independent physical objects that are described by the same metric.The following contexts are to investigate the changes of physical properties for our BH model as it evolves from its singular to regular state.

Sources and energy conditions
In this section, we discuss the matter sources that may generate the BH model Eqs.( 1) and ( 22) in the framework of F (R) gravity and give the energy conditions of the matters sources.The motivation for interpreting our BH model in the context of F (R) gravity is that the entropy of RBHs has a deviation from the entropy-area law of Einstein's gravity if it is calculated by S = dM/T , where T is BH temperature.In other words, either the entropy-area theorem or this thermodynamic formulas is no longer valid.We anticipate that the modified gravity may shed some light on resolving this problem.

Matter sources in our model
We suppose that the action is of the following form, where I M represents the action of the matter sources that will be determined below.The equation of motion with respect to metric reads [60] M where we have defined which leads to a minus sign in the right hand side of Eq. ( 29) Now we analyze the algebraic structure of M µν in order to establish the action of matter sources.Since the metric Eq. ( 1) is spherically symmetric and has one shape function Eq. ( 22), the Ricci tensor R µν has the algebraic structure [(1, 1)(11)], so does F ′ (R)R µν .Here we adopt the Segré notation [61,62], where symbol 1 corresponds to one component of Ricci tensors, all components in parentheses are equal, and a comma separates timelike and spacelike components.The remaining part of M µν is of [1, 1(11)] because the structure of g µν is [1, 1(11)] and R depends only on radial coordinate r.As a result, the algebraic structure of M µν is [1, 1(11)].Therefore, the energy-momentum tensor T µν must have the same structure, otherwise, the equation of motion Eq. ( 29) would not be satisfied.
Let us turn to the structure of the energy-momentum tensor of nonlinear electrodynamics.It is [(1, 1) (11)], such that it can generate spherically symmetric RBHs with one shape function.However, only nonlinear electrodynamics is not enough [8] for its energy-momentum tensor to have the same algebraic structure as M µν .To achieve our purpose, we introduce a scalar field following the method of Ref. [18] because the energy-momentum tensor of a scalar field has the algebraic structure [1(11, 1)].
The combination of nonlinear electromagnetic and scalar fields provides the algebraic structure [1, 1(11)].Thus, we write down the action for the matter source, where V (ϕ) is potential of scalar fields and F = F µν F µν is contraction of strength tensors.The introduction of W (ϕ) is to guarantee that the kinetic term ∂ µ ϕ∂ µ ϕ is positive definite, but this also leads to a redundant degree of freedom.It seems a disadvantage but actually provides more possibilities for us to construct reasonable scalar fields.The energy-momentum tensor of the scalar field reads After substituting the metric Eq. ( 1) into the above equation, we obtain which has the structure [1 (11,1)].The energy-momentum tensor of nonlinear electrodynamics takes the form, where L F denotes the derivative of L with respect to F. The explicit form of where the contraction of strength tensors reads and Q is magnetic charge [14,18].The algebraic structure of ] as we expected.To give T ν µ [ϕ]+T ν µ [A] explicitly, we suppose that F (R) has the same form with the Starobinsky inflation model [63], but with parameter α that will be determined via the formula of entropy below.Then, we solve L, V , and W ϕ ′2 as functions of radial coordinate r by substituting Eqs. ( 22), ( 32), ( 34) and (37) into Eq.( 29), see App.B for the derivations and results. 2 We verify that the following Klein-Gordon equation holds automatically, when the metric Eq. ( 1) with the shape function Eq. ( 22), together with the relevant L, V , W and ϕ, is substituted into it, which means that the metric Eq. ( 1) with the shape function Eq. ( 22) is indeed the solution of F (R) gravity with the matter source described by Eqs. ( 28) and (31).We then draw L, V , and W ϕ ′2 in the density plots of full two-dimensional parameter space, {r, Q}, see Fig. 2. We make two comments.The first is that L, V and W ϕ ′2 may change their signs which can be verified via the Klein-Gordon equation.In other words, L, V and W φ 2 satisfy the Klein-Gordon equation automatically.The results are shown in the density plots of full two-dimensional parameter space, {r, Q}, see Fig. 2.
-5000  We make two comments.The first is that L, V and W φ 2 may change their signs along the r-axis.This signifies that it is appropriate and necessary to introduce W , otherwise a negative φ 2 , φ 2 < 0, leads to a contradiction.The second comment is that these solutions may diverge 10 along the r-axis.This signifies that it is appropriate and necessary to introduce W , otherwise a negative ϕ ′2 , ϕ ′2 < 0, would lead to a contradiction.The second comment is that these solutions may diverge at the BH center.To exhibit this divergent nature clearly, we plot the solutions with specific values of Q in Figs. 3 and 4, where L, V and W ϕ ′2 are divergent at r = 0 for a singular state but approach finite values for a regular state.
at the BH center.To exhibit this divergent nature clearly, we plot the solutions with specific values of Q in Figs. 3 and 4, where L, V and W φ 2 are divergent at r = 0 for a singular state but approach finite values for a regular state.
Figure 3: Nonlinear electrodynamic source.The orange curves correspond to a regular state, while the gray curves a singular state.In addition, the blue dotted line in Fig. 3b depicts Maxwell's asymptotic, L(F) = F. α = 1 is set.
Figure 3: Nonlinear electrodynamic source.The orange curves correspond to a regular state, while the gray curves are a singular state.In addition, the blue dotted line depicts Maxwell's asymptotic, L(F) = F. α = 1 is set.
We further list the asymptotic behaviors of L, V , and W ϕ ′2 for both singular and regular states in order to illustrate the properties of the two states clearly.
We consider a singular state at first.As r → ∞, we have Figure 3: Nonlinear electrodynamic source.The orange curves correspond to a regular state, while the gray curves a singular state.In addition, the blue dotted line in Fig. 3b depicts Maxwell's asymptotic, L(F) = F. α = 1 is set.
Figure 4: Scalar field source.The orange curves correspond to a regular state, while the gray curves a singular state, where α = 1 is set.
We further list the asymptotic behaviors of L, V and W φ 2 for both a singular and regular states in order to illustrate the properties of the two states clearly.
We consider a singular state at first.As r → ∞, we have where L 0 is an integration constant which can be fixed by L → 0 as r approaches infinity [59], 11 where L 0 is an integration constant which can be fixed by L → 0 as r approaches infinity [64], Correspondingly, we obtain and As r → 0, we have which are divergent at r = 0 in the form of r −3 .Now we turn to the asymptotic behaviors of a regular state.As r → ∞, we derive which is different from Eq. ( 40) when compared with the Lagrangian of a singular state.But the asymptotic behaviors of V and W ϕ ′2 remain unchanged when compared with those of a singular state, see Eqs. ( 42) and ( 43), just with the replacement of Q by Q ext .On the other hand, as r → 0, the obvious differences emerge between singular and regular states, i.e., L, V and W ϕ ′2 do not diverge any longer for a regular state, In order to determine W and ϕ separately, we introduce the following ansatz for a scalar field, Thus, W and ϕ ′2 can be separated as shown in Fig. 5, where ϕ ′2 is positive definite although it diverges in a certain region of the parameter space, {r, Q}, see Fig. 5b.Furthermore, we provide the numeric solution of ϕ in Fig. 6.We note that ϕ is a monotone decreasing function of radial coordinate r and vanishes as r approaches infinity.Moreover, ϕ diverges at the BH center except for Q = Q ext , i.e., when the BH stays at its regular state, the scalar field is bounded, max{ϕ(r, The dependence of V and W on ϕ is exhibited in Fig. 7.It is worth emphasizing that V and Figure 7: Dependence of V and W on ϕ.The orange curves (regular states) end up at blue points, while the gray curves (singular states) do not.α = 1 is set.
W will end up in the extreme value of ϕ, see Eq. ( 52), when the BH is in its regular state (orange curves in Fig. 7).Further, we note that W (ϕ) vanishes at ϕ max ≈ 3.3386.This implies W ϕ ′2 = 0 at the BH center.In other words, the scalar field does not have kinetic energy at r = 0 when the BH is in its regular state.

Energy conditions of matter sources
It is known that RBHs can bypass [5,65] the Penrose singularity theorem because they violate the strong energy condition (SEC).In addition, the other energy conditions, like the null energy condition (NEC), the weak energy condition (WEC), and the dominant energy condition (DEC), are applied to examine the pathological behaviors of spacetime, see, e.g., Ref. [66] and the references within.Therefore, the energy conditions play an important role in the study of RBHs [58,[67][68][69].
In the current subsection, we investigate the energy conditions in our model Eq. ( 22) from two aspects: One is the effective energy conditions of spacetime and the other aspect is the energy conditions of each matter source.
For the matter sources Eq. ( 31), we define the energy densities and pressures of scalar and monopole fields, respectively, using our conventions in Sec.4.1, Thus, the effective energy density and pressure of spacetime in F (R) gravity can be introduced directly by the contributions of scalar field Eq. (54a) and monopole field Eq. (54b), where M within a square bracket denotes the matter consisting of ϕ and A. The effective energy conditions of spacetime are shown in Fig. 8, where we use the blue area to mark the validity of energy conditions and the orange curve to mark the horizon, f (r H , Q) = 0.  We note that the effective energy conditions of spacetime are badly damaged, especially the DEC.In addition, the NEC, the WEC, and the SEC are disrupted outside the horizon when Q becomes small.To figure out which ingredient of matter is responsible for such damage, we draw the energy conditions associated with the scalar and monopole fields in Figs. 9 and 10, respectively.
It can be seen that the DEC of scalar and monopole fields is destroyed outside the horizon, whereas the NEC, WEC, and SEC are damaged outside the horizon by the monopole field when the charge is small.In addition, Fig. 8c demonstrates that the SEC is broken near r = 0 even in the singular state, and such a violation is caused by the scalar field, which can be deduced when we compare Fig. 9c with Fig. 10c.As a result, it is not possible to give the reason for the existence of singularity just from the point of view of energy conditions.Perhaps we need to start with the complete theory of geodesic convergence and then analyze Raychaudhuri's equation [42].

Quasinormal frequencies in our model
In this section, we consider the dynamic properties of our model by investigating the frequencies of quasinormal modes (QNMs) and their damping limit for a massless scalar field perturbation, i.e., the so-called asymptotic quasinormal frequencies (QNFs).The calculations we shall make are based on the 13th-order WKB approach improved by the Padé approximation [70,71], the light ring/QNMs correspondence [72,73], and the monodromy method [74][75][76].

Scalar field perturbation
We write down the equation of a massless scalar field [77][78][79], To separate the variables, we take the ansatz, and then obtain the equation satisfied by the radial component, where r * := dr/f (r) denotes the tortoise coordinate and the effective potential reads Substituting the shape function Eq. ( 22) into the above equation, we plot the effective potential for three Q values in Fig. 11, where Fig. 11a is for V (r) and Fig. 11b for V (r * ).The shape of the effective potential implies that we are dealing with a scattering problem.The so-called quasinormal modes are defined as the solution of Eq. ( 58) satisfying the particular boundary conditions, Ψ ∼ e ±iωr * , r * → ±∞.
To compute the QNFs, we use the 13th-order WKB approach improved by the Padé approximation [70,71].For a specific value of l = 6, we compute the QNFs with different values of Q.The relation is shown in Fig. 12, where we find the curves that fit the data best via the polynomials of Q.The linear relation of the real part of the QNFs with charge is obvious in Fig. 12a, where the fit line is of the form, with c 1 ≈ 0.8989 and c 0 ≈ 2.1675.However, the imaginary part shows the nonlinear relation with respect to the charge, where the fit coefficients take the values as follows: Fig. 13 depicts the relationship between ℜ(ω) and −ℑ(ω) with more values of charge Q ∈ [0, Q ext ] and multipole number l ∈ [1,10], where the frequencies with the same multipole number l are joined into a curve.Meanwhile, all curves have a similar shape: One global minimum is located at Q ext and one global maximum at a critical value Q C .We calculate these critical values Q C for a varying multipole number l ∈ [1, 10] and their corresponding QNFs, see Tab. 1.The data circled by the red boxes in Tab. 1 show a possible phenomenon that Q C and −ℑ(ω) converge to constants in the eikonal limit (l ≫ 1).To provide an analytic result for the QNFs, we apply the light ring/QNMs correspondence [72], according to which the QNFs can be cast in the following form, where Ω c denotes the angular velocity of a test scalar particle moving along an unstable null geodesic and λ c is the Lyapunov exponent.The applicability of this correspondence is limited, as evidenced by the case of the Einstein-Lovelock theory [80].Furthermore, if the magnetic field's influence on photon spheres is disregarded in our model, the correspondence is valid.However, if the magnetic field's influence on photon spheres is considered, the correspondence is invalid, as in the cases discussed in Refs.[39,81].Note that the subscript c means that the angular velocity and Lyapunov exponent are calculated at the radius of a circular null geodesic, r c , determined by where f c := f (r c ) and a prime means the derivative with respect to r.Thus we have [73] For the shape function Eq. ( 22), we obtain the equation of a photon sphere, the angular velocity, and the Lyapunov exponent, where k 0 = 0.3339 and k 1 = 0.1360.When the above equation is compared with Eq. ( 61), the former differs from the latter by a factor of l = 6 due to ℜ(ω) = Ω c l.We can see from Fig. 14b that the Lyapunov exponent takes its maximum value, λ max c ≈ 0.3753, when Q C ≈ 0.2689, which corresponds to the maximum of the minus imaginary part due to −ℑ(ω) = n + 1 2 |λ c |.These data are consistent with those in Tab. 1 for the case of l = 6.
We end this subsection by summarizing the dynamic evolution in our model: As the charge Q goes to its extreme value Q ext , the minus imaginary part approaches to its minimum.In other words, the decay amplitude becomes minimum.As a result, the regular state corresponds to the most stable state of our BH model.

Asymptotic quasinormal frequencies
In this subsection, we calculate asymptotic quasinormal frequencies using the monodromy method [74,75].This method is based on the Stokes portrait [76] which consists of Stokes lines ℜ[r * ] = 0, Stokes fields V st = {ℑ [1/f ] , ℜ [1/f ]}, and complex horizons, i.e., the complex roots of the equation, f (r H ) = 0.In fact, the Stokes portrait shows a branch of Riemann surfaces described by r * (r).The operations, ℜ and ℑ, appear because we analytically continue the radial coordinate into the complex plane, r → x + iy or r → ρ e iφ , where x, y, ρ, φ ∈ R.
We start with the analysis of singularities of differential equations.The master equation, Eq. ( 58), has four regular singular points3 due to the behavior of effective potential Eq. ( 59), where only r 1 , r 2 and r 3 are "singular points" of Stokes lines if Q ̸ = 0, see Figs. 15, 17 and 19.Here, the "singular point" refers to the point where the curve (Stokes line) has self-intersection in the sense of algebraic geometry [82].Thus, when using the monodromy method, r = r 0 should be treated as an ordinary point except Q = 0.If Q equals zero, the four regular singular points, see Eq. ( 70), merge into one.In this case, r 0 reduces to one singular point of Stokes lines, see Fig. 15c.
The Stokes portraits of the model Eq. ( 22) can be separated into three groups according to the ways how the Stokes lines cross over a regular singular point.
1. Stokes lines cross over the point r = r 0 , see Figs. 15 and 16, where the Stokes lines in Fig. 15 are built in the Cartesian complex plane, while those in Fig. 16 are in the polar complex plane, i.e. r = ρe iφ .However, only in Fig. 15c, the monodromy of asymptotic solutions of Eq. ( 58) along the Stokes lines is non-trivial.The reason is that r = r 0 is the branch point in Fig. 15c, but not in Figs.15a and 15b.Therefore, we ignore the cases in Figs.15a and  15b.
2. Stokes lines cross over the point r = r 2 or r = r 3 , which is the second non-trivial case, see the Cartesian Stokes portraits in Fig. 17 and polar ones in Fig. 18, respectively.
3. Stokes lines cross over the point r = r 1 , see the Cartesian Stokes portraits in Fig. 19 and polar ones in Fig. 20.Although it is non-trivial, the monodromy is associated with the non-physical horizon at the negative axis of x.Thus, we will not take into account this case, either.At first, we calculate the asymptotic quasinormal frequencies in the case of Q = 0 based on the Stokes portraits, see Figs. 15c and 16c.The tortoise coordinate has the following asymptotic behavior in the limit of r → 0, where z is the analytical continuation of r * into a complex plane.The leading order of the effective potential reads Thus, the master equation becomes which does not contain the multipole number l. From this equation, we solve one asymptotic solution via the combination of Bessel functions, where A 1 and A 2 are two integration constants.
Then we try to find the formula of asymptotic quasinormal frequencies by matching the monodromy of the above asymptotic solution along two contours we choose depending on the Stokes lines in Fig. 15c.The first circles the outer horizon.The second contour starts at infinity in the lower left corner, rotates at the merged regular singular point, and then goes out to infinity in the upper left corner, and finally goes along a large arc to infinity in the lower left corner.The formula of asymptotic quasinormal frequencies takes the form, where T + H (T − H ) is BH temperature at the outer (inner) horizon, n > 0 and n ∈ Z.The details of the calculation can be found in App. A. Comparing with the case in the Schwarzschild BH [74], we note that the real part of ω/T + H is weakened.Now let us turn to the Stokes portrait Fig. 17b, where Q ∈ (0, Q ext ).Similarly, we can obtain the asymptotic quasinormal frequency by the monodromy method, which has the same form as that of Hayward BHs [76].Here ν = √ 1 − 4V 0 with V 0 defined by Eq. ( 106), see App. A. When Q goes to Q ext , the inner and outer horizons are close to each other.When Q is equal to Q ext , the inner and outer horizons merge into the extreme horizon, and the Stokes line seems to cross over the extreme horizon, see Fig. 17a, but this is not the case.
To clarify whether the Stokes line crosses over the extreme horizon or not, we write down the equations of Stokes lines describing this case, The extreme horizon corresponds to the point {x, y} = {0, 2/3}, which is not the solution of Eq. ( 77).In other words, the Stokes line does not cross over the extreme horizon.More precisely, no Stokes line can be defined at that point in the complex plane.The monodromy method is not available for the extreme case, and it is not feasible in the limit of Q → Q ext for Eq. ( 76), either.The reason is that the exponential function e ω/T + H blows up when T + H goes to zero.

Black hole thermodynamics
The thermodynamics of RBHs is associated currently with many inconclusive issues, where two of the most important ones are entropy and the first law of thermodynamics.On the one hand, the entropy of RBHs is problematic because a simple thermodynamic formula, S + = dM/T , yields an inconsistency with that obtained by Wald's Noether-charge method [83] or Hawking's path integral method [84].On the other hand, the problem is that there are ambiguous terms in the first law of thermodynamics for RBHs if the first law is constructed in the thermodynamic phase space without redundant degrees of freedom.
In this section, we propose an alternative method to resolve the inconsistency between the thermodynamic formula, S + = dM/T , and Wald's Noether-charge method.In addition, using a self-consistent entropy, we also examine the thermodynamics of our BH model described by Eq. ( 22) and the relationship between the divergence of heat capacity and the divergence of thermodynamic curvatures at the end of this section.

Thermal entropy and Wald entropy
At first, we denote S + as the entropy obtained by the thermal formula, S + := dM/T , and S W as the Wald entropy.Secondly, we recover the mass parameter in Eq. ( 22) by the replacements of r → r/(2M ), Q → Q/(2M ), and where r, Q, and M have the same dimension, . Next, we represent the outer horizon r + via the following form, and express the temperature in terms of r + , Q, and Q ext , which has an inverse dimension of mass, [T ] = [M ] −1 .Let us compute S + and S W , respectively.If M is regarded as internal energy, the entropy S + can be calculated by the usual thermodynamic formula, The dimension of S + is the square of mass, [ S + ] = [M ] 2 .The first term in Eq. ( 81) corresponds to the entropy-area law, which equals A/4 in Einstein's gravity.The second term is a deviation, which inspires us to interpret our model in terms of F (R) gravity.Moreover, the Wald entropy of F (R) theory is of the form [85], According to our proposal in Sec. 4, if the the Starobinsky form is chosen, F (R) = R + αR 2 , see Eq. ( 37), we have where the dimension of Ricci scalars is the square of the inverse of mass, [R] = [M ] −2 , and α has the dimension of square of mass, [α] = [M ] 2 .The magnitude of R at the outer horizon reads If α is chosen to be the special form, S + equals S W .We note that many possible forms of F (R) give the equality, S + = S W .For instance, if we take rather than the Starobinsky form, where ∆(R) is a non-trivial function of Ricci scalars, e.g., ∆(R) = 1/R, etc., we still obtain this equality for our model.Such an observation can be understood when the different Lagrangians of F (R) gravity lead to the same equations of motion.

Phase transition
In this subsection, we first analyze whether there are Davies points [86] in our model, i.e., the divergent points of heat capacity, corresponding to the second-order phase transition points of our BH model.The motivation comes from the fact that our model, Eq. ( 78), reduces to the Hayward model when Q equals Q ext .The heat capacity of the Hayward BH has a single Davies point [87], thus our model may have the same thermodynamic structure as the Hayward BH.Further, we consider the correspondence [88,89] between the divergences of heat capacity and the divergences of different thermodynamic curvatures in our model.
To this end, we calculate the heat capacity when Q is constant, i.e., which gives an analytical result, where the divergent point (Davies point) is unique and approximately located [90] at the zeros of ∂ 2 S + M = 0, i.e., Q c /(2M ) ≈ 0.428243 that is less than the extreme charge Q ext , see Fig. 21.In other words, the singular state with Q/2M < 0.428243 must undergo a second-order phase transition in order to evolve to the regular state.Now let us turn to the singularities of thermal curvatures.According to Ruppeiner's geometry, the thermodynamic metric is cast in the following form, where ) and Q 2 /(2M ) comes from another factor in the denominator of R R , see Tab. 2 below.However, neither of these two singularities corresponds to Q c , i.e., the singularities of heat capacity and Ruppeiner curvature do not have any evident relation.
A similar situation of singularity will happen in Weinhold's geometry, where the metric reads because Eq. ( 91) and Eq. ( 89) are related to each other via a conformal factor [93].Therefore, R W has the same singularities with those of R R , see Fig. 22b, where Q 1 /(2M ) now corresponds to the zero of (∂ In addition, the singularity, Q 1 /(2M ), coincides [89] with the critical point of first-order phase transition, which corresponds to the zeros of ( Another way to introduce a thermal metric was proposed in Refs.[94,95], which is called Geometrothermodynamics (GTD).This theory is invariant under the Legendre transformation, which leads to three metrics.One of them is [96] which belongs to the type-II of Ref. [89] and the singularities of corresponding thermal curvatures are expected to coincide with the critical points of the second-order phase transition.The thermal curvature, R GTD , is exhibited in Fig. 23.We note that there are three singularities, where Q 4 /(2M ) coincides with the singularity of heat capacity Q c /(2M ) and is the zero of ∂ 2 S + M , whereas Q 3 /(2M ) and Q 5 /(2M ) are zeros of ∂ 2 Q M , which have no counterparts in heat capacity.In a word, R GTD provides more singularities than heat capacity.This phenomenon is the same for some other thermodynamic metrics, such as in Ref. [90].As long as there are multiple types of singularities in the thermodynamic curvature that correspond to the zeros of formulas rather than ∂ 2 S + M (see Tab. 2), we can always reach the same conclusion as above.We end this subsection with a comment on the correspondence between the divergence of heat capacity and the divergence of GTD curvature invariants.Unlike previous works by others, we interpret the RBH in terms of F (R) gravity rather than Einstein's gravity.This interpretation has its advantage because it guarantees the consistency of the Wald entropy S W and the entropy calculated by the usual thermodynamic formula, S + = dM/T , but it explains only partially the connection between the divergence of heat capacity and the divergence of curvature invariants in the GTD.The reason comes from the complicated nature of RBH thermodynamics, which is

Conclusion and outlook
In this paper, we propose an RBH model that can be used to simulate the final state of an SBH.Our motivation is to investigate the physical differences between RBHs and SBHs as a result of singularities' presence or absence.We also interpret this RBH model in terms of F (R) gravity, where the matter sources are a nonlinear scalar field and a magnetic monopole in nonlinear electrodynamics.
The interpretation by F (R) gravity resolves the well-known problem of RBHs, namely the inconsistency between the Wald entropy and the entropy calculated by the usual thermodynamic formula.Moreover, we find that the four energy conditions of two matter constituents are violated in varying degrees, especially the DEC which is completely violated outside the outer horizon.
Further, we investigate the dynamics and thermodynamics of our model in order to examine its physical properties.On the side of dynamics, we study the perturbation of a massless scalar field on this BH background.We compute the corresponding QNFs by the 13th-order WKB approach with a Padé approximation, and their asymptotic limit by the monodromy method.We note from the spectrum of QNFs that the RBH as the final state is also the most stable dynamically, and has the maximum oscillation frequency.As to the asymptotic QNFs, the monodromy method cannot give any valid predictions.The reason for this is not the disappearance of singularity, but the merging of the inner and outer horizons on the Stokes portraits.
On the side of thermodynamics, we investigate the phase transition in addition to explaining the entropy shift by F (R) gravity.We are attempting to discover a connection between the phase transition predicted by heat capacity and three types of thermodynamic curvatures, particularly based on self-consistent entropy.We draw two conclusions here.The first is that the divergence of heat capacity does not coincide with the divergences of both Ruppeiner and Weinhold curvatures in our model, but that it appears as one of the divergent points in the GTD curvature.In other words, the GTD curvature predicts more critical points of second-order phase transitions than the heat capacity.The second conclusion is that the regular state appears as a divergent point in both Ruppeiner and Weinhold curvatures, indicating that it is also a critical point of first-order phase transitions according to the current study on the GTD.
Finally, we mention that the present paper focuses just on a static and spherically symmetric BH that evolves to its regular state when the charge or horizon radius goes to its extreme value -a regular but non-rotating BH.Next, we plan to investigate the evolution of rotating BHs [97], especially the issue of how a singular rotating BH evolves to its regular counterpart, together with the relevant topics, such as the superradiance of rotating BHs [98].To this end, the Newman-Janis algorithm (NJA) [99,100] will be applied to Eqs. ( 1) and ( 22) in order to construct the metrics of rotating BHs.This work is in progress.
we obtain the asymptotic solution after such a rotation, Then, closing the contour and comparing it with the rotation around the outer horizon, we arrive at At last, we derive the formula of asymptotic quasinormal frequencies, Eq. ( 75) associated with the Stokes lines in Fig. 15c, from Eqs. ( 98) and ( 101).Now we calculate the asymptotic quasinormal frequencies by using the monodromy method, see Fig. 17b.At first, we compute the tortoise coordinate at r = r 3 , where see also Fig. 24.Then, we expand the effective potential at r = r 3 and obtain  V eff = − (−1) Thus the master equation becomes where the prefactor takes the form, By imitating the Hayward BH [76], we obtain the leading order of the perturbative solution, where ν = √ 1 − 4V 0 .The distance between r = r 2 and r = r 3 is ℑ[r 2 − r 3 ] = √ 3Q, i.e., the shift of z is z → z + √ 3Qi.After repeating a similar procedure, we finally derive Eq. ( 76) that is associated with the Stokes portrait Fig. 17b.
Taking radical coordinate r as a parameter, we solve the dependences of V and W on ϕ numerically and draw them in Fig. 7.

Figure 1 :
Figure 1: The orange curve depicts the equation f (r H , Q) = 0 and the blue line Q = Q ext .

Figure 2 :
Figure 2: Density plots of matter sources.The orange curves correspond to the horizons, f (r H , Q) = 0, where α = 1 is set.

Figure 2 :
Figure 2: Density plots of matter sources.The orange curves correspond to the horizons, f (r H , Q) = 0, where α = 1 is set.

Figure 4 :
Figure 4: Scalar field source.The orange curves correspond to a regular state, while the gray curves a singular state, where α = 1 is set.

Figure 8 :
Figure 8: Effective energy conditions of spacetime in F (R) gravity.α = 1 is set.

Figure 9 :Figure 10 :
Figure 9: Energy conditions of the scalar field

Figure 11 :
Figure 11: Effective potential with respect to r and r * , respectively.α = l = 1 is set.

Figure 12 :
Figure 12: Dependence of QNFs on the charge Q with l = 6.

Figure 13 :
Figure 13: Quasinormal frequencies in the complex frequency plane

Figure 15 :
Figure 15: Stokes portraits with r 0 = 0 as the central point over which the Stokes lines (orange curves) must cross.The blue vector field corresponds to the Stokes field and the purple points to the complex horizons.

Figure 24 :
Figure 24: Complex function C(Q).The blue point denotes the place where Q = Q ext .

Table 1 :
[1,10]the real part, and the maximum of the minus imaginary part with respect to the multipole number, l ∈[1,10].

Table 2 :
428243 0.191286 2 2/3 /3 2 2/3 /3 The first row displays the formulas that appear in the denominators of curvature invariants in different thermodynamic geometries.The second and third rows give the zeros of different formulas, where the stars in the third row imply that no real roots exist.related to various factors, such as the interpretation and parameterization of spacetime metrics, as well as the construction of thermal metrics.Our model suggests that it is important to consider the poles of curvature invariants, i.e., the zeros in the denominators of curvature invariants, in relation to heat capacity when constructing a thermal metric.