Detectable universes inside regular black holes

While spacetime in the vicinity outside astrophysical black holes is believed to be well understood, the event horizon and the interior remain elusive. Here, we discover a degenerate infinite spectrum of novel general relativity solutions with the same mass-energy and entropy that describe a dark energy universe inside an astrophysical black hole. This regular cosmological black hole is stabilized by a finite tangential pressure applied on the dual cosmological-black hole event horizon, localized up to a quantum indeterminacy. We recover the Bekenstein-Hawking entropy formula from the classical fluid entropy, calculated at a Tolman temperature equal to the cosmological horizon temperature. We further calculate its gravitational quasi-normal modes. We find that cosmological black holes are detectable by gravitational-wave experiments operating within the $\mu{\rm Hz}-{\rm Hz}$ range, like LISA space-interferometer.


Introduction
As early as 1966, Sakharov [1] proposed that the proper equation of state of matter and energy at very high densities is that of a dark energy fluid P = −ρc 2 . About the same time Gliner [2] suggested that a spacetime filled with vacuum could provide a proper description of the final stage of gravitational collapse, replacing the future singularity [2]. Black hole solutions where the singularity is avoided are called regular black holes [3,4,5,6,7,8,9,10] and may or may not involve a de Sitter core. The, so called, dark energy stars or gravastars [11,12,13,14,15,16,17,18,19] generally do not predict the presence of an event horizon.
The idea that a new universe is generated inside a black hole has been put forward in [20,21,22,23,24]. Gonzalez-Diaz [5] was, to our knowledge, the first to explicitly propose that a de Sitter space may complete an exterior Schwartzschild metric with the presence of a kind of cosmological black hole horizon in-between. Later, it was realized by Poisson & Israel that in this case a singular tangential pressure will be exerted on the horizon [6]. The present work elaborates on the Poisson-Israel solution regularizing the horizon by considering that the quantum uncertainty principle applies.
We discover an infinite spectrum of solutions, which describe the fluid shell which matches the interior de Sitter core with the exterior Schwarzschild spacetime. The metric's derivatives are continuous up to any required order. All states of the spectrum have the same energy and entropy. This fluid entropy of the dual horizon recovers the Bekenstein-Hawking black hole entropy for a fluid temperature equal to the cosmological horizon temperature. This spacetime spectrum describes a novel kind of regular black hole we shall call the "cosmological black hole" for brevity.
Gravitational-wave astronomy has opened up the possibility to detect such objects. The cosmological black holes may exist independently than singular or other types of regular black holes, or may describe the state of all detected black holes. The detectability of cosmological black holes is founded on the fact that the fundamental quasi-normal mode, calculated here, is distinctively different than the one of Schwarzschild black holes for any mass. These modes are closely related to the ringdown phase of a post-merger object. This phase follows the inspiral phase of a binary black hole merger. The ringdown phase is dominated by the natural frequencies of black hole spacetime, like a ringing bell. We argue that LIGO-Virgos's detections could involve cosmological black holes, because LIGO-Virgo is not able to discriminate between cosmological and singular black holes, due to the well-known "mode camouflage" mechanism [25] and the inadequate frequency sensitivity. On the other hand, the frequency spectrum of quasi-normal modes of the cosmological black hole interior lie within the detectability range frequencies of the planned LISA space interferometer.

The cosmological black hole solution spectrum
A black hole is formed from material that crossed the Schwartzschild horizon. Thus, inside the horizon it is proper to use the full Einstein equations G µν = (8πG/c 4 )T µν instead of the field vacuum equation G µν = 0. Outside the horizon the Schwartzschild metric should apply assuming that all of material has crossed the horizon. One such solution requires the interior to be a de Sitter vacuum P = −ρ 0 c 2 = const., with a tangential pressure P T = −ρ 0 c 2 θ 1 − r r H + 1 2 ρ 0 c 2 δ r r H − 1 being applied on the horizon r H , where θ and δ denote the Heaviside and Dirac functions respectively. This expression was mentioned (without a derivation) for the first time, to our knowledge, by Poisson & Israel [6]. We derive in detail this solution in Appendix A. Poisson & Israel remarked that this tangential pressure diverges for an observer at some proper distance outside the horizon (see equation (45) of Appendix A). Nevertheless, we see little qualitative difference regarding the physical problems encountered by the Poisson-Israel solution and the Schwarzschild black hole solution, which does present a curvature singularity in the centre. It is only that the problem in the former case is transferred from the center to the horizon, having the curvature singularity of the Schwarzschild solution replaced by a pressure singularity in the Poisson-Israel solution.
However, assuming that within r H there is distributed a mass M • in some non-singular way up to r H , quantum physics suggests that the boundary r H of M • cannot be localized with accuracy greater than the Compton wavelength It is justified therefore to assume there exists a length-scale α that specifies the quantum fuzziness of the horizon ∆r H = α (2) in this case. For an astrophysical black hole with mass M • = O(M ⊙ ), α may equal the Compton wavelength or the Planck scale or a few times the latter, so that in each of these cases. Lacking a quantum theory of gravity we cannot know its precise value, still we shall be able here to reach definite quantitative results, irrespective from the value of α. As we shall now show, the Poisson-Israel solution gets regularized by an infinite spectrum of solutions with the same energy and entropy. Let us assume the static, spherically symmetric ansatz for the metric and the following components of an anisotropic diagonal energy-momentum tensor T 0 0 = −ρ(r)c 2 , T 1 1 = P r (r), T 2 2 = T 3 3 = P T (r), where ρ(r) is the mass density and P r (r), P T (r) the radial and tangential pressures. It is straightforward to show (see Appendix A) that the Einstein equations admit the following formulation and P r (r) = −ρ(r)c 2 (6) Note that another density distribution function, besides the Poisson-Israel solution (40)-(43), that solves this system was identified in Ref. [7]. We discover here a new infinite spectrum of solutions that regularizes the Poisson-Israel solution within the fuzziness α of the horizon  The cosmological black hole spectrum of density (8) and tangential pressure (7) inside the α-shell, that designates the horizon. The values of P T on the boundaries are mentioned separately in brackets, because the scaling does not allow them to be distinguished optically. The value of the radial pressure is P r (r) = −ρ(r)c 2 everywhere. Exterior to the horizon, ρ, P T are zero and in the interior P T = −ρc 2 = −ρ 0 c 2 = const. Each curve in all panels represents a solution of the Einstein equations with the same energy and entropy. The solid curves represent the limiting allowed curves for which ρ ′ ≤ 0 everywhere, for the respective maximum order of continuous metric derivatives.
Upper panels: The metric, its first and second derivatives are continuous everywhere, that is K = 1. We considered N = 3 and each curve corresponds to a different value of A 2 ∈ [− 4,8]. Lower panels: The metric, its first, second and third derivatives are continuous everywhere, that is K = 2. We considered N = 5 and each curve corresponds to a different value of A and ε is given in (3). Proper choices of A (±) n (ε) ensure that the density and consequently the metric through (5) are continuous and have continuous derivatives. The maximum order of the continuous derivatives can be arbitrarily high. This is ensured by demanding to hold the following conditions The condition K ≥ 1 ensures that the metric and its first and second derivatives are continuous, (over)satisfying Lichatowich junction conditions, as well as that the tangential pressure (7) is contin-   Figure 1, where here is shown the solution with minimum order N of each density polynomial (8) for various values of K which is the maximum order of the derivative of density that is continuous. We observe that density variation is constrained for large K within much smaller region than α.
uous. The condition (12) certifies that the total mass of the system up to the radius R = r H + α 2 is equal to and therefore that This means that at r = r H coincide a cosmological and a black hole event horizon if the quantum indeterminacy is ∆r H = α. This renders our solution a regular, free of singularities, type of black hole, which we call the cosmological black hole solution. In the interior the scalar curvature R = 6/r 2 H ∝ M −2 • is finite. We shall see in the next section that the entropy of the cosmological black hole equals the Bekenstein-Hawking entropy at the cosmological horizon temperature. Note that expressing ρ 0 with respect to the horizon radius we get precisely the definition of the critical energy density in a Friedmann cosmological model for r H = c/H 0 , where H 0 denotes the Hubble parameter. The order N of the density polynomial (8) can be arbitrarily high independently from the order K (11) that designates the maximum order of continuous density derivatives. There can always be found solutions as long as N ≥ (3K + 2)/2. The maximum order of continuous metric derivatives is K + 1. In the Appendix B we provide the exact K = 1 spectrum for N = 3 and the exact K = 2 spectrum for N = 5 requiring that density is a decreasing function of radius This condition, along with conditions (10)- (12), certify that the equation of state (6), (7) for the solution (8) constrained by these conditions satisfies the Weak Energy Condition T µν ξ µ ξ ν ≥ 0 for any time-like ξ µ ; namely that ρ ′ ≤ 0, P r + ρ ≥ 0, P T + ρ ≥ 0. We plot the spectrum K = 1 for N = 3 and the spectrum K = 2 for N = 5 in Figure 1. Note that there exist also solutions, that are not symmetrical about the r = r H vertical axis. Furthermore, we remark that the fuzziness α does not necessirily constrain maximally the density variations. For large K, the density variation is localized within a region smaller than α. This is depicted in Figure  2.

Fluid entropy
The work performed by a fluid with stress tensor T ij , i, j = 1, 2, 3 may be written with respect to the strain tensor σ i j as The strain tensor can be decomposed as the sum of a pure shear (shape deformations) and a hydrostatic compression (volume deformations) [26] 1 The trace of strain expresses relative volume change, so that considering unit volume deformations we get in general dV = dTr(σ i j ) [26]. For a spherical anisotropic fluid the space component of the energy-momentum tensor may be written in spherical coordinates as T i j = diag(P r , P T , P T ). The work of the gravitational force under a spherical deformation (no pure shear) is therefore It is P that contributes to the work and not only P r despite the deformation being isotropic. There is an additional, to P r , contribution coming from 2 3 (P T − P r ) due to the stretching forces on the fluid sphere during any spherical deformation (during a spherical expansion/contraction the area of the sphere increases/decreases, therefore there are applied tangential forces). The relativistic thermodynamic Euler relation should involve the pressure that contributes to the work. Thus, for zero chemical potential, inside the α-shell the thermodynamic Euler relation is properly written as where s = s(r) is the total entropy density, including the tangential contribution. We denote T = T (r) the local temperature and we have used equation (7). Local temperature obeys the Tolman law where the constant T 0 is called the Tolman temperature and corresponds to the temperature of the fluid measured by an oberver at infinity. The interior, excluding the horizon, does not contribute to the fluid entropy since T s interior = ρ 0 c 2 + P interior = ρ 0 c 2 − ρ 0 c 2 = 0. Thus, the total fluid entropy (integrating the local entropy over the proper volume in General Relativity) equals the fluid entropy of the event horizon, which using (20), (21), is equal to We get consecutively For all solutions (8) the second term is a constant given in equation (12). We finally get where we used equation (14) that identifies the coincidence of cosmological and black hole event horizons. Note that this result holds for any choice of α and for all solutions of the cosmological black hole spectrum (8). Therefore, irrespectively from the exact value of the temperature T 0 , equation (24) shows that all solutions (8) with the same total mass-energy, correspond also to the same entropy.
Provided α accounts for the quantum indeterminacy of the event horizon (2), this Tolman temperature may be identified with the cosmological temperature T dS . In such a case, by direct substitution of the cosmological temperature in entropy (24), we reach the intriguing conclusion that the fluid entropy equals the Bekenstein-Hawking entropy S BH if where T BH denotes the Bekenstein-Hawking temperature. Equation (25) suggests the interpretation of the Bekenstein-Hawking entropy as the entropy of the horizon realized as a fuzzy fluid shell.
Assuming that the Bekenstein-Hawking entropy is a universal maximum bound, the temperature should be equal to the de Sitter temperature (to maximize the entropy) and the width α of the shell should not exceed the appropriate quantum fuzziness, not specified in this work. In this sense, the equation of state of matter (6), (7) for the cosmological black hole (8) is dictated by maximum entropy and General Relativity. It describes the smooth connection between two different vaccua.
It seems remarkable that we recover the black hole entropy for a Tolman temperature that equals the cosmological temperature, fusing effectively de Sitter and Schwarzschild horizons, considering classical relativistic fluid considerations. In this respect, Figures 1, 2 describe a new type of event horizon, we call a dual horizon, that is a fusion of cosmological and black hole event horizons. Note that the sole quantum assumption in this calculation is that the horizon's width is fuzzy. The exact measure of quantum fuzziness α does not affect the result.

Quasi-normal modes
A linear perturbation analysis about the static equilibrium (4)- (7), performed in Appendix C, shows that a radial perturbation cannot develop unstable radial modes. Unless it is identical to another static equilibrium, it may however develop non-radial oscillation modes. For any compact object, the latter are categorized in two types; polar and axial [27]. In Ref. [28] was argued that, similarly to the Schwarzschild black hole case, polar and axial perturbations are isospectral for an ultra-compact object with a de Sitter core and an ultra-thin shell (limit 2GM → Rc 2 as in our case ε ≪ 1). That is because the master equation for polar perturbations is continuous across the shell. Here we shall calculate the quasi-normal modes of axial perturbations.
Axial perturbations ψ ℓ (r, t) = e −iωt φ ℓ (r) about the static spacetime (4), (5) for any metric function h(r) are described by the Regee-Wheeler type of equation [29] where the scattering potential is and r * is the, so-called, tortoise coordinate defined as We determine the constants of integration (see Appendix D) by setting r * = r + r H ln(r/r H − 1) for r ≥ r H + α/2 and requiring that r * is continuous at r = r H ± α/2. The scattering potential, plotted in Figure 3, is strictly positive certifying the stability of the solutions (8) against axial perturbations. The Sturm-Liouville problem (27) can by solved (calculation of the quasi-normal modes ω n ) by use of the generalized Bohr-Sommerfeld method [30]. It has been shown to be sufficiently accurate  in the case of ultra-compact gravastars [31]. This method dictates identifying E R,n for some n such that (see Appendix E) The quasi-normal mode frequencies are then specified directly as ω n = ω R,n + iω I,n = c E R,n + iE I,n , where The quantities r * 0 , r * 1 corresponding to some r 0 , r 1 , respectively, are the roots of the equation E R − V (r) = 0 that define the bounding region E R ≥ V (r). The r 2 along with r 1 define the reflecting      Figure 1, since any differences are negligible in the scales used. Any larger value of K we tried (≤ 15) gives also identical results. Differences in the precise value of ε are also negligible. There are also depicted the values of Schwarzschild black holes for comparison. region E R ≤ V (r). This is depicted in Figure 3. We were able to calculate the modes ω n,ℓ by use of mixed analytical and numerical calculations, which we describe in detail in Appendix E.
In Figure 4 we depict, for black hole masses M • ∈ [10, 10 9 ]M ⊙ , the frequencies of the mode n = 0, ℓ = 2 for the two limiting solutions of Figures 1a, 1b. These correspond to A ± 3 = −4 and A ± 3 = 8 for {K = 1, N = 3}. The differences between the two solutions in the values of ω R , ω I are very small and decrease with increasing black hole mass. We find numerically that the two solutions bound the values of the modes for all solutions with K = 1, K = 2 and we find strong numerical evidence that this is true at least up to K = 15. Our results are identical for α ranging from 10 9 times the Planck length down to α equal to the Compton wavelength. Note that in contrast to the case of a Shwarzschild black hole, the value of (2GM • /c 3 )ω depends on the cosmological black hole mass, although this dependence is very small.
In Table 1 we list the values of the quasi-normal mode frequencies for 0 ≤ n ≤ 7 (the importance of the first seven quasi-normal modes has been emphasized recently [32,33,34,35,36]) and 2 ≤ ℓ ≤ 4 for a cosmological black hole with M • = 10M ⊙ . For black hole masses M • ∈ [10, 10 9 ]M ⊙ , the fundamental frequency lies in the range 10 −6 Hz ω R,0 50Hz and the n = 5 overtone is about ten times bigger, as depicted in Figure 5. The damping time of the overtones are on the other hand drastically lower with respect to the fundamental mode.
The fundametal mode alone may be used to reconstruct the full inspiral merger ringdown waveform of a binary black hole merger signal (e.g. GW150914 [37,38]). The highest possible fundamental mode of an astrophysical cosmological black hole (corresponding to the minimum possible mass M • = 10M ⊙ , that is a merger of two 5M • black holes) is 63Hz ( Figure 5) which lies outside the detection range of LIGO-Virgo. This does not mean that LIGO-Virgo observations exclude the possibility that they involved cosmological black holes, but only that LIGO-Virgo cannot discrimininate between a singular (Schwarzschild or Kerr) and a cosmological black hole. The reason is the, so called, "mode camouflage" mechanism [25]. The ringdown modes of a black hole (regular or not) are determined by the external null geodesic and not by interior fluctuations [39]. Fluctuations generated inside our cosmological black hole will dominate after the exterior perturbations are damped. Thus, LIGO-Virgo cannot discriminate between a Schwarzschild (or Kerr) black hole and a cosmological black hole. Following the full damping of the external light-ring modes, the internal fluctuation frequencies lie outside the frequency range detectability of LIGO-Virgo.
However, it is evident from Figure 5 that for M • 10 4 M ⊙ the fundamental mode of cosmological black hole fluctuations lies within the frequency detectability range (∼ 10 −1 − 10 −5 Hz) of the LISA space interferometer. It sounds in particular intriguing that LISA may be able to detect an intermediate mass cosmological black hole through its postmerger ringdown phase, even if the binary inspiral phase cannot be detected. Still, in order to estimate the minimum possible amplitude sensitivity of an interferometer so as to detect a cosmological black hole ringdown, the excitation factors of its quasi-normal modes, following a binary merger, have to be calculated. This is an involved task, that this work urges the community to perform. In every case, our results as in Figure 5, clearly suggest that despite the well-known mode camouflage mechanism of ultra-compact objects [25] mentioned above, cosmological black holes particularly are in principle detectable and distinguishable from singular black holes. If already LISA is not amplitude-wise sensitive enough, it is a matter of developing the appropriate technology to detect cosmological black holes provided they exist.
Finally, let us remark that regular black holes, like the one we propopse, should not suffer from the instabilities, such as the light-ring instability [40,25], the ergosphere instability [41] and the accretion instability [42,43,44], that have been argued to occur in gravastars and dark energy stars. The absence of a horizon is a key assumption that drives the appearance of these instabilities [45].

Discussion
We discover that in General Relativity, regular black holes containing a de Sitter core correspond to a spectrum of spacetime solutions assuming quantum indeterminacy of the localization of the horizon, which behaves as an anisotropic fluid shell. All spacetime states of the cosmological black hole spectrum have the same energy and entropy, resembling a quantum degeneracy. This is a fluid entropy. It recovers the Bekenstein-Hawking black hole entropy if the Tolman temperature of the fluid is identified with the temperature of the cosmological horizon, fusing the cosmological and black hole horizons in a single dual horizon.
The quasi-normal modes of cosmological black holes are distinctively different than the ones of Schwarzschild black holes. Still, LIGO-Virgo cannot disciminate between cosmological and Schwarzschild black holes, because of the well-known mode camouflage mechanism -ringdown waveform is dominated initially by spacetime fluctuations in the region of the external null geodesic, that is common in regular and singular black holes-and the fact that the mode frequencies of astrophysical cosmological black holes, namely 10 −6 Hz ω R 10Hz, lie outside the frequency's range detectability of LIGO-Virgo. Therefore, it remains open the possibility that LIGO-Virgo's detections are cosmological black holes. Most importantly, the quasi-normal frequency range of astrophysical cosmological black holes lies inside the detectability frequency range of the planned space interferometer LISA. Thus, this work urges the community to investigate further the properties of cosmological black holes, proposed here, and especially their inspiral and ringdown waveforms. There arises the fascinating possibility that black hole detections are also detections of dark energy universes. If these may evolve to inflationary universes similar to our own and if the latter is itself such an object remain open possibilities that beg for further investigation. [43] Baoyi Chen, Yanbei Chen, Yiqiu Ma, Ka-Lok R. Lo, and Ling Sun. Instability of exotic compact objects and its implications for gravitational-wave echoes. arXiv e-prints, page arXiv:1902.08180, February 2019.

A Derivation of the Poisson-Israel solution
Assuming the static, spherically symmetric ansatz (4), the Einstein equations give 8πG 8πG where prime denotes differentiation with respect to r. Let us denote T 0 0 = −ρ(r)c 2 , T 1 1 = P r (r), T 2 2 = T 3 3 = P T (r), where ρ(r) has dimensions of mass density and P r (r), P T (r) dimensions of pressure.
Assuming further a function m(r) such that the Einstein equations (34), (35) give One solution of these equations is where ρ 0 = M • /(4πr 3 H /3), r H = 2GM • /c 2 and M • = m (PI) (r H ). We denote δ(x) the Dirac δ-function and θ(x) the Heaviside step function The superscript (PI) is an acronym for "Poisson-Israel", because the expression (39) has appeared for the first time, to our knowledge, in Ref. [6]. Poisson & Israel remarked that an observer at a proper distance ∆s outside the horizon will perceive an infinite tangential pressure on the horizon Let us now prove that equations (40)-(43) satisfy the Einstein equations (37)- (39). We shall use the dimensionless variable and the dimensionless quantities The dimensionless functionm (PI) (u) is written as Using the property θ ′ (u) = δ(u) we have that which, considering the identity uδ(u) = 0, gives Comparing with (37) we getρ that is equation (41). Now, we have that where we have used the transformation x = u − 1. Thus, the tangential pressure (39) is that is the Poisson-Israel solution (43).

B Analytical expressions of the cosmological black hole spectrum
We shall use in the followings the dimensionless quantities (9), (47). The spectrum (8) is written in these dimensionless variables as For K = 1 and N = 3 the conditions (10)- (12) give The additional requirement ρ ′ ≤ 0 (62) imposes the constraint, to zero-th order of ε, For B 1 = B 2 = B, the requirement ρ ′ ≤ 0 imposes the constraint, to zero-th order of ε, Note that if we identify α with Compton wavelength h/(M • c) then and if we identify it with the Planck scale, then The parameter α expresses the quantum fuzziness of the horizon. It is a free parameter within our framework and may be bigger than the Planck length.
In our convention, the Einstein equations read Differentiating G 1 1 , substituting it into G 2 2 and using Einstein equations to substitute G µ ν for T µ ν we get the equation of hydrodynamic equilibrium We may reach to the same equation by using the continuity equation Let us now consider the perturbations λ = λ eq + δλ, ν = ν eq + δν, ρ = ρ eq + δρ, P r = P r,eq + δP r , P T = P T,eq + δP T , u = u eq + δu. (93) Subscript "eq" denotes equilibrium quantities. It is e νeq = e −λeq = h, P r,eq = −ρ eq c 2 , P T,eq = −ρ eq c 2 − 1 2 rρ ′ eq , u eq = 0, so that u = δu. To the first order we get δ (1) T 1 0 = −(ρ eq c 2 + P r,eq )δu = 0, 8πG where Γ denotes the Gamma function. Here x 0 , x 1 denote the roots of E n − V (x) = 0 that define the bounding region E n ≥ V (x) and x 2 is the upper limit of the reflecting region E n ≤ V (x), as in Figure 3. This expression can be simplified further for sufficiently low modes [47] as x 1 x 0 The imaginary part is a measure of the barrier penetrability, that is absent in the normal Bohr-Sommerfeld rule (113) since the barrier is infinite in the latter case. We shall denote E n = E R,n + iE I,n .
In case the imaginary part is negligible with respect to the real part E I,n ≪ E R,n the generalized Bohr-Sommerfeld rule (116) may be decomposed as follows which correspond to (30), (32). We used this generalized Bohr-Sommerfeld formulas to solve the Sturm-Liouville problem (27). The accuracy of the method in the calculation of the quasi-normal modes of gravastars has been verified in Ref. [31]. We calculated the integrals (30), (32) by a combination of analytical and numerical techniques. For ε ≪ 1 the point r 0 lies within the de Sitter core r 0 /r H < 1 − ε/2 and the points r 1 < r 2 lie outside the black hole within the Schwarzschild spacetime r 1 /r H > 1 + ε/2. In particular, we have and we use the dimensionless variable u = r/r H . It is This diverges logarithmically for ε ≪ 1 and therefore can be calculated analytically for a Compton or Planck or smaller ε (see (74), (75)) atanh In the interval r H − α/2 ≤ r ≤ r H + α/2 it applies the density profile (8) and we use the variable x = (r − r H )/α = (u − 1)/ε = O(1). The metric function h ± is of order ε, h ± = O(ε). Thus we have The function h ± is Taylor expanded about ε = 0 and the integral is calculated numerically. We find numerically that for K = 1 and K = 2 all quasi-normal modes are constrained (as in Figure 4) by the two limiting solutions of the K = 1, N = 3 solutions that correspond to A and ρ II (±) = −4x 3 ± 6x 2 − 3x + h II (±) = ε 3x 4 ∓ 6x 3 + 9 2 and the integral (123) is calculated numerically for both solutions I, II, decomposing it to regions −1/2 ≤ x ≤ 0 where h − applies and 0 ≤ x ≤ 1/2 where h + applies. There remains the region 1/2 ≤ x ≤ x 1 corresponding to r H + α/2 ≤ r ≤ r 1 , equivalently 1 + ε/2 ≤ u ≤ u 1 . The radius r 1 is the smallest root within r ≥ r H + α/2 of We write the integral as where P (u 1 ) = 0. Considering that (ln(u − 1)) ′ = 1/(u − 1) we get This expression is calculated numerically in a straightforward manner. The integral in the numerator of equation (32) is calculated directly numerically. The integral in the denominator of equation (32) is calculated by using analogous treatments as above, except in the Schwarzschild region where there appears the additional pole at u 1 . We approximated this non-regular integral by use of the transformation u = y + 1 and expanding P (y + 1) with respect to y. We get