Regular black holes with improved energy conditions and their analogues in fluids

As is well-known, Einstein's general relativity lacks [1] the ultraviolet (UV) completeness that is reflected [2] in the singular solutions of Einstein's equations at the classical level and in the non-renormalizability at the quantum level. Regular black holes (RBHs) [3, 4], which have no curvature singularities at the centers, challenge the UV incompleteness at the classical level. This challenge originated from the change of vacuum [5, 6] and was implemented through various approaches, such as the introduction of nonlinear matter [7], the deformation of the commutative spaces [8], the regularization of singularities by quantum effects [9, 10], and the assistance of alternative theories of gravity [11–13]. Meanwhile, the apparent differences that can be tested between RBHs and singular BHs (SBHs) have motivated numerous studies [14–16]. Whether we can distinguish RBHs from SBHs by theoretical and experimental evidences is a critical point for research programs in the field of RBHs. Because RBHs are widely considered to be related to quantum physics, the discovery of RBHs in the universe will certainly provide a new hope to search for quantum gravity.

There are two ways to construct RBHs. The first starts with establishing BH metrics via certain mathematical rules [17–19] that guarantee the finite curvatures at BH centers, followed by enduing these metrics with physical meanings; e.g., the action of matter was provided, and then the theory of RBHs was established [20]. The second way draws support from physical theories or phenomena, e.g., the existence of a finite length scale [8] or the asymptotic safety [10], involves the derivation of the corrected metrics, which give rise to finite curvatures.

Nevertheless, among all the RBHs constructed in the above two ways, more than a few models break physical conditions or conjectures, in particular, the limited curvature conjecture [21], which states that the curvature invariants should be bounded by some universal value; the weak energy condition (WEC) that is associated with the second law of BH mechanics [22] or the dominant energy condition (DEC) that is related to the causal structure of spacetime [23]. The violation of these energy conditions motivates us to interpret RBHs from the perspective of quantum corrections.

If the violation is located inside regular black holes, the interpretation would be reasonable. The first reason is that an event horizon prevents any observer from observing a black-hole interior, and the second one is that some unknown effects that occur within event horizons or reachable microscopic scales may cause quantum corrections. However, if the violation occurs outside horizons, the problem immediately arises. For instance, if the DEC is violated outside horizons, the causal structure of spacetime will be broken, indicating that the observer outside horizons will encounter chaotic causal phenomena. This is not acceptable from the perspective of physics.

For more information about the energy conditions of RBHs, we refer the reader to Ref. [24], where the energy conditions of four well-known regular black holes are reviewed; in particular, the DEC is taken into consideration as a criterion to determine whether a regular black hole is realistic. In brief, the weak, null, and dominant energy conditions are the primary prerequisites for us to construct realistic RBHs, where the strong energy condition is an exception. Note that the violation of the DEC is pointed out for some RBHs in Ref. [24]; our aim is to further study how to achieve the recovery of the DEC for these RBHs.

Analogue gravity as a tool for gaining insights into general relativity has shown [25] its significance, representing a great leap from passively waiting for signals from external galaxies to actively studying BHs in ground laboratories. Among various manifestations of analogue gravity, acoustic BHs (ABHs) have not only a long history [26] but also an active status in current research [27–29]. Following our previous work [30], in which we proposed a new method to construct acoustic regular BHs (ARBHs), we explore ARBHs in terms of the energy conditions of their astronomical counterparts in the present work. Meanwhile, the simulation strategy used in the present work is different from our previous one; i.e., we adopt the approach proposed in Ref. [31], where the singular Schwarzschild and Reissner-Nordström spacetimes can be simulated in fluids. Our aim is to investigate realistic RBHs with the help of analogue gravity and try to find apparent evidences or phenomena for distinguishing RBHs from SBHs.

The remainder of this paper is organized as follows. In Sec. II, we clarify what a realistic RBH means by discussing the energy conditions. In Sec. III, we propose a remedy to those RBHs that break the dominant energy condition, including the Bardeen BH, the Hayward BH, and their extensions. Sec. IV is dedicated to conformally related RBHs, where we prove a no-go theorem under two general situations. In Sec. V, we simulate a realistic RBH in a fluid by using the properties of flows, where two specific models are discussed in terms of the asymptotic analysis associated with the numerical method. Inspired by the locally polytropic behaviors in the equations of state (EoSs) in the above section, we address the question of whether it is possible to obtain an RBH that possesses a globally polytropic EoS in Sec. VI. Secs. VII and VIII cover cylindrical RBHs and lower dimensional RBHs with polar symmetry, respectively. The conclusions, along with future outlooks, are summarized in Sec. IX. The appendices are dedicated to detailed discussions of the differential inequalities (Apps. A and B), the derivation of Eq. (8) (App. C), the regularity conditions of n-dimensional RBHs (App. D), and the asymptotic analysis for solving nonlinear differential equations (App. E).

Let us first consider the simplest case of RBHs whose metrics are spherically symmetric and of the following form:

where f is the shape function, ξ is the radial coordinate, and σ is dimensionless and may contain several parameters $ x_i $, $ i=1,\ldots,N $, such as mass and charge. Moreover, these parameters must appear in σ via the combinations $ \xi^{n_0} x_1^{n_1}\cdots x_N^{n_N} $, which are also dimensionless. If every combination includes a non-zero $ n_0 $, we can reduce one parameter and obtain $ N-1 $ independent dimensionless parameters by following the Buckingham π theorem [32].

From the mathematical perspective, the regularity of curvature invariants at BH centers demands the limit $ \sigma\sim O(\xi^n) $ with $ n\ge 3 $, and the asymptotic flatness requires the limit $ \sigma\sim O(\xi^m) $ with $ 0\le m< 1 $ [33]. When $ m=1 $, f may converge to a non-zero and non-unit constant, such that the spacetime is Ricci flat, with $ R=0 $, at infinity. The general properties of the shape function are illustrated in Ref. [18].

Given a Riemann tensor, one can construct 17 curvature invariants in total, which are called the Zakhary-Mcintosh (ZM) invariants [34]. For the metric Eq. (1), all 17 ZM invariants consist of the combinations of σ and its first and second derivatives, i.e., $ \sigma' $ and $ \sigma'' $. If certain conditions are further considered, e.g., σ is a positive definite or monotonic function of the radial coordinate, one can show that $ \sigma\sim O(\xi^n) $ with $ n\ge 3 $ and simultaneously guarantee that all the ZM invariants are finite at the center of the RBH described by Eq. (1).

From the physical perspective, the constructed RBHs should not violate the weak, null, and dominant energy conditions, which play important roles [35, 36]. The three energy conditions together with the strong energy condition can be formulated in three classes via different approaches: geometric, physical, and effective ones [37]. If the mechanism of constructing RBHs does not change the gravitational part of Einstein's equations, e.g., the gravitational field coupled with nonlinear electric fields or magnetic monopoles [7, 38], the three classes of definitions are equivalent. In this situation, the energy-momentum tensor can be represented via Eq. (1), and Einstein's equations read as follows:

where the prime denotes the derivative with respect to ξ. Because $ T^{t}_{\;t}=T^{\xi}_{\;\xi} $, there is no need [39] to distinguish the definitions of energy densities inside and outside the horizon; i.e., the energy density inside the horizon is the same as that outside the horizon for the metric we are considering. We define the energy density $ \epsilon $ and pressures $ p_\xi $ and $ p_\perp $ by the diagonal components of $ T^\mu_{\;\nu} $:

Thus, the four energy conditions can be cast [40] in terms of σ and its derivatives:

where WEC denotes the weak energy condition, NEC denotes the null energy condition, SEC denotes the strong energy condition, and DEC denotes the dominant energy condition.

It is not difficult from the four energy conditions in Eq. (4) to find that the NEC, i.e., $ \xi \sigma ''\leq 2 \sigma ' $, must be maintained; otherwise, the other three conditions will be broken. In other words, $ \xi \sigma ''\leq 2 \sigma ' $ is an inequality that ensures the four energy conditions; furthermore, this differential inequality can be solved using the Grönwall-Bellman lemma [41], and its solution reads $ \sigma\le \sigma_0 \xi^3 $, where $ \sigma_0:=\lim_{\xi\to 0}\sigma/\xi^3 $ is a positive constant (see App. A for details). As a counterexample, we consider the widely discussed model [19] with $ \sigma = \exp[-q^2/(2 M \xi)] $, where q represents the charge, in which this inequality is invalid because of $ \sigma_0=0 $. Therefore, the model of Ref. [19] is suggested to be ruled out from realistic RBHs, as are its extensions [42], because matter generating such RBHs breaks the four energy conditions.

Next, the first inequality in the WEC and DEC, i.e., $ \sigma'\ge 0 $, provides a solution, i.e., $ \sigma\ge 0 $, under the boundary condition $ \sigma(\xi)|_{\xi=0}=0 $; that is, σ is a non-negative and monotonically increasing function of ξ. It is known that the non-minimal Wu-Yang monopole [43–47] is a counterexample because its σ function, i.e., $\sigma=-\xi^3 \left(Q^2- 2 M \xi\right)/ [2 M \left(2 q Q^2 +\xi^4\right)]$, where Q is a charge parameter, is not monotonic and not strictly positive either. Thus, the WEC of the Wu-Yang monopole is broken, as is the DEC. In addition, the breaking of $ \sigma'\ge 0 $ may lead to other problems in the construction of RBHs. For instance, when σ is bell shaped, i.e., $\sigma =4\exp\left(-q^2/\xi^2\right)- \exp\left(-2 q^2/\xi^2\right)$, the corresponding BH has two horizons and all the curvature invariants are finite, but the extreme horizon radius is the maximum of the horizons, and the temperature is divergent as the radial coordinate approaches the extreme horizon radius.

The SEC implies an attractive interaction due to the Landau-Raychaudhuri equation [48]; i.e., when the affine parameter increases, the expansion scalar of a family of neighboring time-like geodesics decreases because of the condition $ -(\epsilon+p_\xi+2p_\perp)\propto \sigma''\le 0 $. Therefore, the violation of the SEC leads to a repulsive interaction. However, this violation is nothing to be concerned about, because the SEC of an RBH must be broken [39] near an RBH center. Moreover, the zeros of the equation, i.e., $ \sigma''(\xi_*)=0 $, separate the spacetime into different types of interactions, which will be discussed later with concrete examples.

Compared with the other energy conditions, the DEC has its particularity reflected in the inequality $ \xi\sigma''\ge-2\sigma' $, which can be visualized from the Ricci curvature $ R\propto \xi \sigma ''+2 \sigma ' $; i.e., the negative Ricci curvature violates [49] the DEC. However, the differential inequality $ \xi\sigma''\ge-2\sigma' $gives the solution $ \xi\sigma\ge 0 $under the boundary conditions $ \sigma(\xi)|_{\xi=0}=0=\sigma'(\xi)|_{\xi=0} $. This solution is trivial and provides no more constraints to the σ function. In practice, one does not need to verify the four energy conditions. If the DEC is valid, the WEC ad NEC are also valid. Therefore, checking the DEC is enough to guarantee the WEC and NEC. As for the SEC, we do not need to check it individually for RBHs, because it is valid beyond an RBH central region but invalid within an RBH central region because of a repulsive interaction. As a result, if the DEC is maintained for an RBH, the WEC and NEC are ensured automatically.

In summary, we list the requirements for a realistic RBH from the perspective of energy conditions. If an RBH with the metric Eq. (1) is realistic, its σ-function has the following behaviors:

● σ is a non-negative and monotonically increasing function of ξ, where $ \xi\in[0,\infty) $;

● σ must be bounded by $ \sigma_0 \xi^3 $ from above, i.e., $ \sigma\le \sigma_0 \xi^3 $, where $ \sigma_0:=\lim_{\xi\to0}\sigma/\xi^3 $ and $ \sigma_0 $ must be positive.

These two conditions are necessary but not sufficient for an RBH to be realistic; see App. A for a detailed explanation. In the next section, we show that some well-known examples, such as the Bardeen and Hayward BHs, comply with these two conditions, but their dominant energy conditions are broken. To solve this problem, we provide a phenomenological approach to restore their dominant energy conditions.

The problem of Hayward BHs depicted [17] by $ \sigma=\xi^3/(\xi^3+q^3) $, where q represents the charge, is the violation of the DEC in the region $ \xi> 2^{1/3} q $, even if this σ satisfies the two items above. The reason is explained in App. A.

As we mentioned above regarding the special status of the DEC, Hayward BHs also violate the WEC and NEC when $ \xi> 2^{1/3} q $. For constructing a Hayward-like BH that ensures the DEC in $ \xi\in[0,\infty) $, we propose the following σ function:

where $ M^{\mu-3} $ is introduced for balancing the dimension. The DEC requires $ 2<\mu \le (\sqrt{145}-7)/2\approx 2.52 $, under which the Hayward-like BH given by Eq. (5) would be realistic in the whole region of ξ, i.e., $ \xi\in[0,\infty) $. Alternatively, the dimensionless σ can be established via parameterization, i.e., $ \sigma=(\xi/l)^3/[1+(\xi/l)^{\mu }] $, where l is a parameter with the length dimension.

A similar procedure can be applied to the Bardeen BH, which gives rise to a Bardeen-like σ function:

For this model, the DEC gives rise to $4/3 < \mu \le\left(\sqrt{113}- 7\right)/2\approx 1.82$.

In fact, we can construct a general σ function

which satisfies the DEC if the parameters μ and ν take the values in the following regions (see App. C for the derivation):

It is not difficult to verify that the RBHs described by Eq. (7) are realistic because the numerator plays a decisive role, $ \sigma\sim O(\xi^3) $ when $ \xi\to 0 $, and the asymptotic flatness is maintained [33] simultaneously, i.e., $ f\to 1 $, because the power of the denominator, $ \mu\nu $, is greater than 2 when $ \xi\to \infty $.

Nevertheless, we cannot remedy all the RBHs that break the DEC by simply changing the power of radial coordinates. If σ is not a rational function, for instance, $ \sigma = \exp[-q^2/(2 M \xi)] $, this model cannot be repaired. On the other hand, although the RBH obtained via quantum corrections, e.g., Refs. [9, 10], can be remedied via the above phenomenological method, the remedied model will lose the original motivation of quantum corrections. Let us take the RG-improved Schwarzschild BH [50] as an example, which is motivated by the theory of gravitational asymptotic safety [51, 52]. The shape function reads

where $ G(r) $ is the running Newton constant, which plays a similar role to the σ-function, $ G_0 $ is identified with the experimentally observed value of Newton's constant, and ω and γ are two positive parameters. This RG-improved BH is regular from the perspective of finite curvatures, but it breaks the DEC because $ G(r) $ violates $ -2 G'\le r G'' $ beyond a certain value $ r_0 $, where $ r_0 $ is determined by a positive root of the algebraic equation

Thus, according to our remedy used above, we change $ r^3 $ to $ r^\mu $ and multiply by $ M^{3-\mu} $ for balancing the dimension in the denominator of $ G(r) $:

which reveals that the DEC requires $ 0\le \mu \le (\sqrt{145}-7)/ 2\approx 2.52 $. However, such a modification loses the original motivation of the RG-improvement, which can be understood from the distance scale λ that provides the relevant cutoff for the Newton constant. Using Eq. (10) and the formula [50] $ \widetilde G(r) =G_0 \lambda ^2/(G_0 \omega +\lambda ^2) $, we obtain

and give the asymptotic behaviors at zero and infinity, respectively:

where the second one violates the original asymptotic requirement, i.e., $\lambda^{ \undertilde {r \rightarrow \infty}}r $.

In the next section, we demonstrate that the conformally related RBHs cannot be repaired either, by proving a no-go theorem.

We discuss two classes of conformally related regular black holes: the conformally related Schwarzschild-type black holes and the astronomical counterparts of the ARBHs with the unit speed of sound.

We claim that one cannot establish a scale factor Ω that regularizes the Schwarzschild BH and makes the metric satisfy the DEC at the same time. To specify our statement, let us first express the metric of conformally related Schwarzschild BHs [53]:

where the scale factor is set to be $ \Omega=\exp [S(\xi)]>0 $ and $ 2M=1 $ is chosen for the discussions in this subsection. The metric being regularized implies that the corresponding curvature invariants are finite in the whole spacetime, particularly at the BH center. Next, instead of observing the Kretschmann scalar K, we concentrate on the contraction of two Weyl tensors $ W_{\mu\nu\alpha\beta} $ and $ W^{\mu\nu\alpha\beta} $, where $ W:= W_{\mu\nu\alpha\beta} W^{\mu\nu\alpha\beta} $, which is referred to as the Weyl curvature hereinafter. Because of the Ricci decomposition [54, 55], we obtain $ W=K - 2 R_2+R^2/3 $ and see that the Kretschmann scalar and Weyl curvature are equivalent for diagnosing the singularity in the four-dimensional spacetime, where $ R_2:= R_{\mu\nu} R^{\mu\nu} $ is the contraction of two Ricci tensors and $ R:= g^{\mu\nu}R_{\mu\nu} $ is the Ricci scalar. The Weyl curvature corresponding to the metric Eq. (13) reads

which is finite at the BH center if $ \mathrm{e}^{- S(\xi)} $ converges to zero no slower than $ \xi^3 $. When $ \mathrm{e}^{- S(\xi)} $ converges to zero on the order of $ \xi^3 $, i.e., $ \mathrm{e}^{- S(\xi)}\sim O(\xi^3) $, $ S(\xi) $ diverges positively, and its first-order derivative must be negative. In contrast, the asymptotic flatness requires $ \Omega\to 1 $ as $ \xi\to \infty $; i.e., $ S(\xi) $ must converge to zero at infinity. Summarizing the above properties of $ S(\xi) $, we find that $ \Omega^{-1}= \mathrm{e}^{-S} $ is a bounded function on the whole non-negative axis of ξ.

The energy conditions of the conformally related Schwarzschild BH given by Eq. (13) should be investigated inside and outside the horizon, because $ T^t_{\; t} $ no longer equals $ T^\xi_{\;\xi} $. In other words, the constraint $ T^t_{\; t}=T^\xi_{\;\xi} $ breaks the finiteness of the Weyl curvature. This can be understood easily by solving $ G^t_{\; t}= G^\xi_{\;\xi} $ as a differential equation of $ S(\xi) $, which provides a solution $ \mathrm{e}^{-S}=c_2 (\xi+2 c_1)^2 $ that converges to zero slower than $ \xi^3 $, where $ c_1 $ and $ c_2 $ are two integration constants. Consequently, the energy conditions inside and outside the horizon are different and should be treated separately. The energy density and pressures are defined inside the horizon ($ \xi<1 $) as

and are defined outside the horizon ($ \xi>1 $) as

where $ G^t_{\; t} $, $ G^\xi_{\; \xi} $, and $ G^\theta_{\; \theta} $ are components of the Einstein tensor calculated using the metric of Eq. (13). Thus, the DEC is reduced to four differential inequalities in terms of $ S(\xi) $ and its derivatives $ S'(\xi) $ and $ S''(\xi) $ in the range of $ \xi<1 $ or $ \xi>1 $. Among all the differential inequalities, $ \epsilon^{\rm in}+p_\xi^{\rm in}\ge 0 $ and $ \epsilon^{\rm out}+p_\xi^{\rm out}\ge 0 $ provide the same differential inequality:

Multiplying both sides of this inequality by a non-negative factor $ \mathrm{e}^{-S/2} $, we arrive at

from which we can conclude that $ \mathrm{e}^{-S/2} $ is a convex function in the range of $ \xi\in[0,1)\cup (1,\infty) $. However, the finiteness of curvature invariants and asymptotic flatness of the metric demand that $ \mathrm{e}^{-S} $ is bounded; therefore, $ \mathrm{e}^{-S/2} $ must be a constant,^① which is obviously contradictory to the asymptotic behavior of the Weyl curvature at $ \xi\to 0 $, $ \mathrm{e}^{-S}\sim O(\xi^n) $ with $ n\ge 3 $.

In other words, there exists no such a conformal factor Ω that can regularize the Schwarzschild BH and guarantee the DEC simultaneously. This conclusion can be extended to the conformally related Schwarzschild-type BHs with singularity at $ \xi=0 $. For such a BH with the metric

where $ \sigma(\xi)/\xi $ is of a unique pole at $ \xi=0 $ and goes to zero as $ \xi\to\infty $, there exists no conformal factor Ω that satisfies the following two conditions simultaneously:

● The Weyl curvature of metric $ \Omega\, g_{\mu\nu} $ is finite in $ \xi\in [0,\infty) $, where $ \Omega\to 1 $ at $ \xi\to\infty $;

● The DEC based on $ \Omega\, g_{\mu\nu} $ is valid.

This is the so-called no-go theorem for the conformally related Schwarzschild-type BHs that belong to conformally related RBHs.

The proof exactly follows the case of conformally related Schwarzschild BHs. At first, the Weyl curvature of metric Eq. (19) reads

Assuming that $ \xi=0 $ is the d-th order pole of σ, we can obtain an asymptotic relation:

which ensures that the Weyl curvature is finite. When $ \xi\to\infty $, the asymptotic flatness demands $ \mathrm{e}^{- S} \to 1 $. Moreover, the conditions $ \epsilon^{\rm in}+p_\xi^{\rm in}\ge 0 $ and $ \epsilon^{\rm out}+p_\xi^{\rm out}\ge 0 $ provide exactly the same differential inequality as Eq. (17), from which we find that $ \mathrm{e}^{- S} $ is a convex function. Nevertheless, the combination of the convexity and boundness of $ \mathrm{e}^{- S} $ leads to a contradiction with the asymptotic relation of Eq. (21). Therefore, our statement is proved.

On the premise that the speed of sound is set to be unity, we proposed [30] a general method to construct ARBHs in a fluid, where the metrics are similar to those of conformally related BHs [53]. However, as we noted in Ref. [30], the astronomical counterparts of ARBHs under a certain parameterization violate the DEC; i.e., the ARBHs we constructed hardly have any physical counterpart in the universe. Thus, it is natural to ask if we can find a way that the DEC for the astronomical counterparts of the ARBHs can be repaired and consequently the astronomical counterparts of the ARBHs can be detected in the universe.

The present case differs from that in the above subsection, as the unregularized metric $ g_{\mu\nu} $ has no singularity. Moreover, the ARBHs are different from the conformally related BHs because their conformal factors are proportional to the energy density of fluids and not constrained by any dynamical equations.

We discuss this in detail by following the strategy used in Ref. [30], which is opposite to that of Sec. IV.A. We express [30] the metric of acoustic RBHs with spherical symmetry:

where r is the radial coordinate in fluids,^② ρ represents the mass density, and v represents the radial velocity of fluids. The density and velocity are related [30] by $ v=A/(\rho r^2) $, where A is a positive constant. The radial velocity v is supposed to be positive; otherwise, the density will be negative. Moreover, because we set the speed of sound to be unity, v is dimensionless and A has the same dimension as $ \rho r^2 $.

By considering the similarity between Eqs. (13) and (22) and dealing with ρ as a scale factor, i.e., $ \rho:= \mathrm{e}^{S(r)} $, we can express the Weyl curvature as

where the prime denotes the derivative with respect to r. The regularity at $ r=0 $ demands

Namely, ρ diverges because $ \rho\propto r^{-n} $, and v converges to a constant [30] in order to make sure that the Weyl curvature is finite when $ r\to 0 $. As mentioned in the second paragraph of this subsection, if ρ was removed from Eq. (22), the remaining metric would still give rise to finite curvature invariants everywhere, which is different from the situation in Sec. IV.A. In addition, the asymptotic flatness requires $ \rho\to 1 $ or $ v\sim O(r^{-2}) $ as $ r\to \infty $; i.e., v is a monotonically decreasing function at infinity.

Similar to the discussion in Sec. IV.A, we suppose that the acoustic metric Eq. (22) directly corresponds to a spacetime metric. Thus, the energy density and pressures corresponding to the astronomical matter generating astronomical BHs must be defined inside and outside the horizon separately; otherwise, the equation $ G^{t}_{\;t}=G^{r}_{\;r} $ gives a false solution, i.e., $ S(r)=c_4-2 \ln(r+2 c_3) $, where $ c_3 $ and $ c_4 $ are integration constants. If this solution is consistent with the regularity Eq. (24), we have $ c_3=0 $, which leads to the result that f degenerates to a constant, i.e., $ 1-A^2c_4^2 $, such that the corresponding metric is no longer a BH solution. We then deduce that the forms of density and pressure inside the horizon must be different from those outside. Next, following the proof process in Sec. IV.A, we derive the inequality from two similar inequalities, i.e., $ \epsilon^{\rm in}+p_r^{\rm in}\ge 0 $ and $ \epsilon^{\rm out}+p_r^{\rm out}\ge 0 $:

which is similar to Eq. (17). Therefore, we conclude that $ \mathrm{e}^{-S/2} $ is a convex function in the region of $r\in \{r| r > 0, v\ne 1\}$, which contradicts the regular condition of Eq. (24) and asymptotic flatness. In other words, even if $ \widetilde g_{\mu\nu}/\rho $ is regular in the sense of finite curvatures, the regularized metric $ \widetilde g_{\mu\nu} $ cannot satisfy the DEC. That is, the DEC is violated in the astronomical counterparts of the ARBHs with the unit speed of sound.

Now, we discuss the simulation of realistic RBHs in a fluid. Our aim is to construct the spacetime of realistic RBHs with spherical symmetry using acoustic waves and verify the conditions under which the realistic RBHs can be simulated in the fluid, i.e., find the equations of state. Our result may have guiding significance for experiments.

We start with the general stationary acoustic metric [25]

where ρ and $ v^i $ represent the mass density and velocity of fluids, respectively, $ c:= \sqrt{|\partial p/\partial \rho|} $ represents the local speed of sound, and p represents the pressure. We shall use Eq. (26) to simulate realistic RBHs by providing the equations of state.

First, we suppose that the fluid is spherically symmetric and its velocity contains only a radial component, i.e., $ v^i=\{v_r(r), 0, 0\} $; thus, Eq. (26) is reduced to the following form in spherical coordinates:

where $ \mathrm{d} \Omega^2:= \mathrm{d} \theta^2+\sin^2\theta \mathrm{d} \phi^2 $. By using the solution of the continuity equation

to replace $ v_r $ and defining a new variable

we rewrite the acoustic metric as

which is supposed to equal the astronomical metric formally, where F and H are defined as follows:

Here, the prime denotes the derivative with respect to r. Note that our F and H are slightly different from those in Ref. [31], where they were represented in terms of ξ and v. The purpose of our expression is to derive analytical expressions for the density ρ and pressure p of the fluid.

Second, we impose the condition $ F H=1 $, i.e., the simulated metric has only one shape function, which leads to the differential equation

Now, we solve c using the first equation of Eq. (31):

where the negative root has been ignored owing to the positive c and ρ. By substituting $ c=\beta/\rho $ into Eq. (32), we obtain

where there are no rules for selecting any one of the two solutions at this moment. Then, we derive the density analytically:

where $ \rho_0 $ is an integration constant. In addition, using Eqs. (33) and (35) together with the definition of c, we compute the pressure:

where $ p_0 $ is an integration constant. Eqs. (35) and (36) give the equation of state for the fluid.

In practice, F as a function of ξ corresponds to the shape function of the realistic RBH that we will to simulate in the fluid. The relationship between ξ and r, i.e., Eq. (29), can be represented [31] by the nonlinear differential equation

or equivalently by

which does not have analytical solutions generally.^③ Furthermore, the variables associated with the fluid can be written in terms of the following functions of r:

We note from Eq. (39) that there exists a special position $ r_c $ that is a stationary point of both ρ and p, and this point is determined by $ \xi''(r_c)=0 $ because $ \rho' $ and $ p' $ are proportional to $ \xi'' $. Detailed numerical analyses with concrete examples are presented below.

Before studying specific models, we perform an asymptotic analysis for Eq. (38) and provide general properties of solutions at $ r\to 0 $ and $ r\to\infty $.

For the simulated RBH at $ r\to 0 $, we obtain an asymptotic F with the help of Ref. [33], i.e.,

when ξ approaches 0. Here $ R(0) $ is the limit of Ricci scalars at $ \xi=0 $. Using the dominant balance [57, 58] and the boundary condition $ \xi(r)|_{r=0}=0 $, we find the asymptotic solution of Eq. (38) when ξ approaches 0 (see App. E for details):

where $ c_6 $ is an integration constant. From Eqs. (33), (40), and (41) we derive the asymptotic forms of β, the density, and the pressure, respectively:

which do not depend on specific RBH models at the leading order. However, if we substitute Eq. (41) directly into Eq. (39), we are able to fix $ c_6 $ and rule out the redundant root through comparison with Eq. (42). As a result, we obtain $ c_6 = \rho_0 $ and know that the solution with subscript "$ + $'' is physical; see App. E for details. Thus, we omit the subscript "$ + $'' in ρ and p for simplifying the notation in the following.

We observe from Eq. (42) that the pressure of fluids at $ r\to 0 $ must be divergent, while the density converges to zero. Furthermore, according to the first equation of Eq. (39), we obtain the speed of sound:

which is divergent at $ r\to0 $. In other words, if the maximum speed of sound exists [59], there must be a cutoff $ r_0 $, such that the speed of sound is regularized. This will be particularly important in the numerical calculation later.

Combining the density and pressure in Eq. (42), we give the equation of state around the BH center:

where $ W_0(\cdot) $ is the Lambert W function. To make Eq. (44) more intuitive, we perform the asymptotic expansion and give the leading order for $ \rho\to 0 $:

Here, the relation between ρ and p around $ r=0 $ is linear.

In contrast, for the simulated RBH at $ r\to\infty $, we have $ F\sim 1-2M/ \xi^{n} $, $ 0<n\le 1 $ from Eq. (31), and the asymptotic solution $ \xi= c_7r $, where $ c_7 $ is constant (see App. E). Thus, we obtain $ \beta\sim 1-2M/(c_7\, r)^n $ using Eq. (33). The density and pressure can be solved via Eqs. (35) and (36) with the plus subscript:

where we have kept the leading term of p valid by adding $ n\neq 1/4 $. The limits of the density and pressure are $ \rho\to \rho_0 $ and $ p\to p_0-(1-4n)/\rho_0 $, respectively. In particular, the approximate equation of state at infinity reads

which is polytropic; more precisely, it describes a thermal process similar to that of the Chaplygin gas [60]. If $ n=1/4 $, the equation of state becomes approximately $ p= p_0-1/(2 \rho) $.

Thus far, the discussions in this section have focused on the simulation of the RBHs given by Eq. (1) in terms of acoustic analogy. In Secs. V.A and V.B, we simulate two realistic RBHs whose DEC is valid. Details can be found in Sec. III for the first model and Ref. [20] for the second one.

We simulate the remedied Bardeen BH by taking $ \mu=3/2 $, $ M=1 $, and $ q=1 $ in Eq. (6). The solution of Eq. (37) obtained numerically with the boundary condition $ \xi(0.2)=0.2 $ is shown in Fig. 1, where the two horizons are $ r_-\approx 0.35 $ and $ r_+ \approx 1.39 $ or equivalently $ \xi_-\approx1.69 $ and $ \xi_+ \approx 14.34 $. The initial point starting at $ r_0=0.2 $ instead of $ 0 $ is based on the possibility of the existence of the maximum speed of sound, and such a setting can avoid dealing with the divergent speed of sound, velocity of fluid, and pressure in numerical calculations. Furthermore, the critical point that separates the spacetime into different types of interactions is determined by $ \sigma''(\xi_*)=0 $; see Eq. (4) and the following discussions about the SEC, i.e., $ \xi_*=2^{2/3}\approx 1.59 < \xi_- $, which is located inside the inner horizon.

Figure 1.  (color online) $ \xi(r) $.

Figure 2 shows the speed of sound and velocity of fluid and highlights their difference by using the Mach number, i.e., $ \mathcal{M}:= v/c $. We note that the Mach number is located in the range of $ \mathcal{M}\in[0.8, 1.2] $ between the inner and outer horizons, which indicates that the transonic phenomenon occurs. A similar phenomenon was observed [31] for SBHs. As a matter of fact, the existence of horizons for the acoustic model described by Eq. (27) separates the spacetime into different regions according to the signs of $ c^2-v^2 $. For the simulated RBHs with one horizon, the fluid inside the horizon flows with the transonic phenomenon. For the simulated RBHs with two horizons, the transonic flow is sandwiched between the two horizons.

Figure 2.  (color online) $ v(r) $ and $ c(r) $. The two curves are almost overlapped. For details, see the inset graph of $ \mathcal{M}(r) $.

Generally, the Mach number can be computed with the help of Eq. (39):

Because $ F(\xi)<0 $ as long as $ \xi\in(\xi_-,\xi_+) $, we have $ z<0 $ between the two horizons; meanwhile, we find that the Mach number is constrained by the following inequality:^④

For the remedied Bardeen model, the minimum of z can be calculated numerically, i.e., $ z_{\rm min}\approx -0.033 $, under our settings, i.e., $ A=M=q=1 $, $ \rho_0=1 $, and $ p_0=0 $, and this minimum corresponds to $ r\approx 0.773 $, which is located between the two horizons. Thus, the Mach number is in the range $ 1<\mathcal{M}<1.034 $.

We also provide the numerical calculations of Eqs. (33), (35), and (36) with subscript "+" for the density and pressure in Figs. 3(a) and 3(b), respectively, and the EoS in Fig. 3(c). In Figs. 3(a) and 3(b), there are a global maximum of ρ and a global maximum of p located between the two horizons, i.e., $ r_c=0.489 $. This point plays a special role in Fig. 3(c) because it shows a sharp discontinuity of the EoS. In addition, the curve of the EoS ends at the green dot, where its values of ρ and p, i.e., $ \rho\approx0.509 $ and $ p\approx 49240.914 $, are estimated numerically when r approaches$ 500 $ as the infinity of our numerical calculations.

Figure 3.  (color online) Numerical solutions for the remedied Bardeen model, where $ A=M=q=1 $, $ \rho_0=1 $, and $ p_0=0 $.

Now, we turn to the model associated with nonlinear electromagnetic fields [20], whose σ is a rational and sigmoid function:

where ξ denotes the radial coordinate in the RBH. Note that σ is non-negative and monotonically increasing in the whole region of ξ, i.e.,$ \sigma\ge 0 $ and $ \sigma'=q \xi^2/(q+\xi)^2\ge 0 $; meanwhile, σ is bounded by $ \sigma\le \xi^3/q^3 $ because of $ \sigma_0= 1/q^3 $. The critical point of this model can be obtained by solving $ \sigma''(\xi_*)=0 $, which gives $ \xi_*=q $. Moreover, the existence of horizons demands $ q\le 4 \xi_{\rm Sch} /27 $; i.e., the critical point $ \xi_* $ is not greater than $ 4 \xi_{\rm Sch} /27 $, where $ \xi_{\rm Sch}=2M $ is the Schwarzschild horizon radius.

The numerical results are shown in Fig. 4, where $ M=1/2 $, $ q=0.1 $, and $ A=1 $ are set. Because p and c are divergent at $ r=0 $, see Eqs. (42) and (43). We have performed a cutoff for the lower boundary by setting $ r_0=0.2 $ as we did for the remedied Bardeen model.

Figure 4.  (color online) Numerical solutions for the model associated with nonlinear electromagnetic fields, where $ q=0.1 $, $ A=1 $, $ \rho_0=1 $, and $ p_0=0 $ are set.

Here are four points that need to be demonstrated.

● The oscillation of Mach numbers at the left tail in Fig. 4(b) arises from our computational accuracy;

● The critical point $ r_c $at which ρ and p are maximized is no longer located between the two horizons (see Figs. 4(c) and 4(d)), which is different from the case of the remedied Bardeen model;

● The green dot in Fig. 4(e) denotes the end of the EoS and corresponds to $ \rho\approx 13.322 $ and $ p\approx 48416.215 $;

● The upper boundary of Mach numbers can be estimated via a formula similar to Eq. (49), i.e., $ \mathcal{M}<1.010 $, where z reaches its global minimum $ z_{\rm min}=-0.010 $ at $ r=0.209 $.

Motivated by the polytropic behaviors of the equations of state in Sec. V, we try to determine whether the fluid with polytropic equations of state can be used to simulate realistic RBHs. For this purpose, we suppose that the fluid is barotropic and has a polytropic equation of state $ p=\tilde B\rho^\gamma $, where γ is a real and nonzero number. Its valid region is determined below, and $ \tilde B $ is constant. The local speed of sound can be calculated as

Meanwhile, the fluid should satisfy the continuity equation, which provides the same relationship between ρ and v as Eq. (28). Then, by substituting Eqs. (51) and (28) into Eq. (27), we give the metric in terms of ρ as follows:

where f is the shape function.

Now, we investigate whether the above metric can mimic the RHBs with the metric of Eq. (1). To this end, we use the new variable ξ defined by Eq. (29), which takes the following form when Eq. (51) is considered:

When we replace r with ξ by using Eq. (53), the metric of Eq. (52) becomes

where $ \xi':= \mathrm{d} \xi/ \mathrm{d} r $. Therefore, the condition $ g_{tt} g_{\xi\xi}=-1 $, leads to the following differential equation of ρ:

whose general solution is

where $ c_8 $ is an integration constant. Furthermore, we note that the asymptotic behavior of ρ is $ \rho\to B^{-2/(\gamma +1)} $ as $ r\to0 $ if $ -1<\gamma <3 $, while for $ \gamma <-1\cup\gamma >3 $, it is $ \rho\to c_8^{-4/(\gamma +1)} r^{-4/(3-\gamma)} $. Thus, the asymptotic behaviors of Weyl curvatures for these two cases are

The asymptotic relation associated with $ -1<\gamma<3 $ shows that the Weyl curvature inevitably has a singular point at $ r=0 $; as for the case of $ \gamma <-1\cup \gamma >3 $, the regularity requires $ 1\leq \gamma <3 $, which contradicts $ \gamma <-1\cup \gamma >3 $. As a result, the fluid with polytropic equations of state cannot simulate the RBHs with the metric of Eq. (1).

Let us investigate under what conditions the metric of Eq. (52) describes an RBH solution in the whole spacetime. We calculate the Weyl curvature of the metric of Eq. (52) directly,

To give the conditions just mentioned, we make an asymptotic ansatz [1] $ \rho\sim r^{-n} $ as $ r\to 0 $, and substitute it into Eq. (58), where n is a real and positive number. We find that the square root of Weyl curvatures consists of the following two terms:

where we have omitted the relevant coefficients. The regularity at $ r=0 $ provides two inequalities:

On the other hand, the asymptotic flatness requires

when $ r\to\infty $. In summary, we conclude that the metric of Eq. (52) describes an RBH if Eqs. (60) and (61) are satisfied.

We take the Chaplygin gas as an example, where $ \gamma=-1 $ in Eq. (52) and the density has the form

where l is introduced for balancing the length dimension. The corresponding Weyl curvature is

which has an asymptotic relation

when $ r\to 0 $. Thus, the Weyl curvature is regular. Moreover, the bracket in the denominator of Eq. (63) is an algebraic quadratic function of r, but it has no real roots because $ l\in \mathbb{R} $ and $ B\ne 0 $, which consequently indicates that the Weyl curvature is finite on the non-negative axis.

It is time for us to investigate the energy conditions for the astronomical counterpart of the metric Eq. (52). We study the vacuum equation $ T^{t}_{\; t}=T^{r}_{\; r} $ to clarify whether we have to define the energy density and pressure inside and outside the horizon, respectively. The vacuum equation leads to a second-order nonlinear differential equation of ρ:

whose general solution is

where $ c_9 $ and $ c_{10} $ are integration constants. Furthermore, the asymptotic analysis at $ r\to 0 $ gives us two situations:

None of them can realize a regular Weyl curvature at $ r=0 $; i.e., we have to discuss the energy conditions inside and outside a horizon separately.

The definitions of the energy density $ \epsilon $ and radial pressure $ p_r $ depend on the number of horizons $ \mathfrak{n} $; meanwhile, the horizons separate the spacetime into $ \mathfrak{n}+1 $ regions. If we start from the region outside the outermost horizon and denote that area as $ 1 $, then for the region with odd number $ \mathfrak{n}\in2\mathbb{N}+1 $, the energy density and radial pressure are defined by

while for even $ \mathfrak{n}\in2\mathbb{N} $, they are defined by

It can be verified that the model of Eq. (62) violates the DEC because there is no intersection between $ \epsilon\ge |{p_r}| $ and $ \epsilon\ge |{p_t}| $, regardless of whether the number of horizons is odd or even. Thus, the Chaplygin gas cannot be used to mimic an astronomical counterpart. In fact, the inverse problem, i.e., constructing ρ from the energy conditions, is rather complicated because the DEC leads to four second-order nonlinear differential inequalities, which are difficult to deal with. Therefore, we stop searching for the models with a polytropic EoS and leave this for future studies.

In this section, we study the RBHs with cylindrical symmetry. The metric can be cast as follows:

where ζ is a real function of ξ. Instead of analyzing the Weyl curvature, we analyze the Ricci scalar because of its simplicity:

where the prime denotes the derivative with respect to ξ. The Ricci scalar is regular at $ \xi=0 $ if σ and ζ have the asymptotic forms $ \sigma\sim O(\xi^m) $ and $ \zeta\sim O(\xi^n) $, where $ m\ge 3 $ and $ n\ge 2 $ or $ n=0 $. In contrast, the asymptotic flatness demands $ \sigma\sim O(\xi^{\tilde m}) $ with $ \tilde m<1 $ and $ \zeta\sim O(1) $.

Now, let us consider the energy conditions of this type of black hole. When $ \zeta=0 $, $ G^t_{\; t}=G^\xi_{\; \xi} $ is valid in the regions inside and outside a horizon. The energy density and pressures can be calculated only in terms of σ and its derivatives:

Thus, the DEC is

Note that the equality comes from $ \epsilon\ge p_\phi \cap p_z\ge -\epsilon $ and can be used to fix σ, i.e., $ \sigma = c_{11} \xi+c_{12} \xi \ln (\xi) $, where $ c_{11} $ and $ c_{12} $ are integration constants. However, the corresponding Ricci scalar, i.e., $ R=2 c_{12}M/\xi^2 $, is singular unless the metric is trivial, $ c_{12}=0 $. In other words, the metric Eq. (70) with $ \zeta=0 $ can never represent a realistic RBH.

By omitting the flow along the z direction, i.e., considering only the equatorial section (ES), we reduce the DEC to

where the contribution from $ p_z $ is ignored. Therefore, we have the following three differential inequalities:

We take the modified Hayward BH as an example (see Eq. (5)) but discuss its cylindrical counterpart:

By substituting Eq. (76) into Eq. (75), we obtain

which hold for $ \xi\in [0,\infty) $ if $ -4 (\sqrt{2}+1)\le \alpha \leq 4 (\sqrt{2}-1) $.

To simulate the equatorial sections of cylindrical regular black holes, we take $ \alpha = 3/2 $ in Eq. (76). Then, we can obtain the nonlinear differential equation (see Eq. (12) of Ref. [31]) that describes the relationship between the radial coordinate of black holes (ξ) and that of the simulation in fluids (r):

where the prime denotes the derivative with respect to r. From Eq. (78), we can obtain the asymptotic relations of ξ and r:

and

where $ c_{13} $ and $ c_{14} $ are integration constants. Meanwhile, because the physical variables of fluids can be represented by ξ and its derivatives, see Eq. (39),

their asymptotic behaviors for the case of the positive sign are

when $ r\to 0 $ and

when $ r\to \infty $. The phenomenon of transonic flows occurs outside the horizon, and the Mach number converges to 1, i.e., $ \mathcal{M}\to1 $, as r approaches infinity. The corresponding numerical analysis is shown in Fig. 5, where we have adopted the settings $ M=1/2 $, $ q=1/2 $, and $ A=1 $ and chosen $ r_{\rm H}=0 $, which corresponds to $ \xi(0)\approx -3.837 $. The upper boundary of the Mach numbers is determined numerically, i.e., $ \mathcal{M}\le 1.124 $, and the maximum is reached at $ r\approx0.796 $.

Figure 5.  (color online) Numerical solutions for the ES of repaired cylindrical Hayward-like models with $ M=q=0.5 $ and $ A=1 $, where $ r_c\approx -1.352 $ and $ r_0\approx-3.675 $, which is determined by $ \xi(r_0)=0 $, and $ r_\infty\approx 1.684 $, which is determined by a very large value of ξ in the numerical calculation.

For the case of $ \zeta\ne 0 $, the vacuum equation $ G^t_{\;t}=G^\xi_{\; \xi} $ gives rise to

which is valid for an arbitrary $ \xi\in [0,\infty) $. The solution $ \zeta =\ln (\xi-c_{15})+c_{16} $, where $ c_{15} $ and $ c_{16} $ are integration constants, does not satisfy the regular condition. Thus, we have to separate the discussion for the inside of the horizon from that for the outside of the horizon. However, $ p_\xi^{\rm out}\ge -\epsilon^{\rm out} $ and $ p_\xi^{\rm in}\ge -\epsilon^{\rm in} $ provide the same inequality

which can be solved by multiplying both sides by $ \exp(\zeta) $:

This inequality indicates that $ \exp(\zeta) $ is a concave function in $ \xi\in [0, \xi_{\rm H})\cup (\xi_{\rm H}, \infty) $; thus, ζ is a concave function, because the logarithm of a non-negative and concave function is concave.^⑤ This can also be seen from Eq. (85); $ \zeta'' $ is nonpositive because $ (\zeta')^2 $ is nonnegative. Nevertheless, ζ is bounded to satisfy the condition of finite curvatures, which contradicts the concavity of ζ. In other words, the metric Eq. (70) can in no way represent a realistic RBH. In contrast, the $ 2 $D BHs with polar symmetry can be regarded as cylindrical BHs with the z direction being suppressed. In the following section, we investigate the $ 2 $D polar-symmetric RBHs and their simulations in fluids.

The $ (2+1) $D simulation is rather different from the $ (3+1) $D case discussed in the previous section. First, the Weyl curvature tensor vanishes identically for any $ (2+1) $D spacetime, thus the Weyl scalar is no longer an appropriate candidate for analyzing the curvature divergence. Second, the regularity condition is closely related to the dimension of spacetime. When one studies the RBHs in a $ (2+1) $D spacetime, the criteria for shape functions will change, see App. D. Thirdly, the form of acoustic line elements also depends on dimension, in particular the prefactor [25]. Moreover, the $ (2+1) $D analogue BHs in a fluid are rather interesting and relatively easy to be realized in a laboratory.

Let us start considering a circularly symmetric BH,

where μ is mass-like parameter. According to App. D, such a metric is curvature regular at the BH center if $ \sigma\sim O(\xi^n) $, $ n\ge 2 $, as $ \xi\to 0 $, and asymptotic flat if $ \sigma \to 0 $ as $ \xi\to \infty $. In order to construct an example that satisfies these conditions, we make an ansatz,

and substitute it into Kretschmann scalars. We find that the Kretschmann scalar is regular at the BH center if $ \alpha\ge 2 $ and that the asymptotic flatness is satisfied if $ \beta>\alpha $. Next, given a particular α we are going to fix β by the DEC. For instance, taking $ \alpha=2 $, we reduce the DEC to the two inequalities,

Because the left hand sides of the two inequalities are positive, the two inequalities hold for all non-negative ξ and q if their right hand sides are non-positive, i.e., the DEC is satisfied in the whole spacetime. Following this idea, we find $ 0<\beta \le 2 \sqrt{3}-2\approx 1.46 $. However, this result contradicts to the asymptotic flatness. In fact, for the case of $ \alpha\ge 2 $, the positive energy density $ \epsilon \ge 0 $ leads to $ \alpha\ge \beta $ in $ \xi\in[0,\infty) $, while the asymptotic flatness requires $ \beta>\alpha $. No intersections exist.

If relaxing the asymptotic flatness, we replace it with the Ricci flatness, $ R=0 $, at infinity. Then substituting the ansatz Eq. (88) into the Ricci scalar, we obtain

Since β is non-negative, the power of ξ in the denominator is larger than that in the numerator, thus R vanishes as ξ approaches infinity. In other words, if $ \alpha=2 $, the metric with Eq. (88) automatically satisfies the condition of Ricci flatness. As a result, the model Eq. (88) together with the Ricci flatness is regular and satisfies the DEC, that is, it is a realistic RBH.

Furthermore, before we focus on the analogue in a fluid, we make a note on the toy model we just constructed, see Eqs. (87) and (88). The causal structure of the $ (2+1) $D RBH with Eq. (88) is exotic. Since the power of ξ in the numerator of σ is larger than that of the denominator, i.e., $ \alpha=2 $ and $ 0< \beta \le 2 \sqrt{3}-2\approx 1.46 $, σ is an increasing function with respect to ξ. Thus, the shape function f is greater than zero inside the horizon but less than zero outside the horizon. This indicates that this $ (2+1) $D RBH has an opposite structure of lightcones when compared with that of usual BHs, like the Schwarzschild BH.

Now let us turn to the simulation. From Eq. (37), we obtain [31] the relation between ξ and r,

For a specific case, $ \alpha=2 $ and $ \beta=1 $, i.e., $ \sigma= \xi^2/(\xi+q) $, we find the asymptotic solutions,

and

where $ k:=\sqrt{(\sqrt{4 A^2+\mu ^2}-\mu)/2} $, and $ c_{17} $ and $ c_{18} $ are integration constants. Since the situation with the positive sign corresponds to the positive correlation between ξ and r, we would like to select it as the candidate for the simulation. Meanwhile, we note from Eq. (93) that there is a movable singularity in the asymptotic solution as $ \xi\ \to \infty $. In other words, $ \xi^+_\infty(r) $ diverges at a finite $ r_\infty $, where $ r_\infty $ depends on the choice of $ \xi^+_\infty(0) $. To estimate the value of $ r_\infty $, we apply the condition $ \xi^+_\infty(0)=\xi_{\rm H} $ which determines the integration constant $ c_{18}=2 \sqrt{\sqrt{2}-1} $ and $ r_\infty = 2 \sqrt{(\sqrt{2}-1)/{(\sqrt{4 A^2+\mu ^2}-\mu)}} $. Here $ \xi_{\rm H} $ denotes the horizon radius of the astronomical counterpart depicted by Eq. (87).

Similarly, by applying Eq. (81) we find the asymptotic behaviors for the situation with the positive sign around $ r=0 $,

and in the limit of $ r\to\infty $,

The phenomenon of transonic flows occurs outside the horizon and the Mach number converges to a constant $ \mathcal{M}\to A/k^2 $ as r approaches infinity. The numerical analysis is shown in Fig. 6, where we have adopted the setting $ A=1 $, $ q=1/2 $, and $ \mu=1/2 $, together with the condition $ \xi^+(0)=1+\sqrt{2} $. The upper boundary of Mach numbers is determined numerically, $ \mathcal{M}\le (1+\sqrt{17})/4 $.

Figure 6.  (color online) Numerical solutions for the $ (2+1) $D repaired Hayward model with $ \mu=q=0.5 $ and $ A=1 $, where $ r_0\approx-2.203 $ which is determined by $ \xi(r_0)=0 $, and $ r_\infty\approx 1.030 $ which is determined by a very large value of ξ.

As an extension of our previous work [30], we analyze two specific questions on RBHs in the present paper: The first is how to remedy astronomical RBHs whose DEC is invalid; and the second is how to simulate realistic RBHs through acoustic gravity. We emphasize that the research strategies of the previous and present works are completely opposite. In the previous work [30], we construct an acoustic metric which is regular at first, then we investigate the energy conditions of the astronomical counterpart of the acoustic RBH. In the present work, we remedy an astronomical RBH to ensure the DEC and finite curvatures at first, and then we simulate it in a fluid.

The DEC of astronomical RBHs occupies a fundamental and decisive position in the research of RBHs. It determines whether an RBH is observable in the universe or is only theoretical. On the premise of the DEC, we have remedied several widely known RBHs whose DEC was broken [24]. In other words, the dominant energy condition, which is violated in the regular black holes listed in Ref. [24], can be recovered using our strategy given in Sec. III. Thus, these regular black holes revert back to the realistic. The procedure we proposed in the present paper works for a broad class of RBHs with the broken DEC, in particular, for those models whose shape functions are rational fraction functions of the radial coordinate. In addition, we have demonstrated that two types of conformally related RBHs can never meet the four energy conditions by proving a no-go theorem.

Although the analogue gravity is widely regarded as a tool of gaining insight into general relativity [25], the first simulation of Schwarzschild and Reissner-Nordström black holes was not realized until 2021 [31]. Prior to this simulation, the acoustic counterparts could not distinguish [61] those astronomical BHs that differ by a conformal factor, where the differences would be shown in the investigations of desired phenomena from the present point of view. For instance, the quasinormal modes (QNMs) of BHs would be affected [62] by conformal factors. We hope such an analogue made in a fluid would mimic those astronomic BHs just mentioned. To this end, starting with realistic RBHs, we have constructed their counterparts in acoustic gravity. Our ultimate goals focus on the guidance on simulation of realistic RBHs in a fluid and the possibility to distinguish RBHs from SBHs. Our guidance on simulation has been illustrated by the equations of state, see Figs. 3(c), 4(e), and 6(e), and by the speed of sound and velocity of flow, see Figs. 2, 4(b), and 6(b). Moreover, the similarities and differences are listed below when we compare our RBHs with the RN BH [31] in the aspects of velocities and Mach numbers, and in the aspects of densities and pressures for $ (3+1) $D models, their equatorial sections, and $ (2+1) $D models as follows.

● Velocity and Mach number in (3+1)D models

The velocity of flow is divergent at $ r=0 $ and tends to zero as $ r\to \infty $ for the RBHs considered in Sec. V, thus the Mach number converges to unity as $ r\to 0 $ and vanishes as $ r\to \infty $. For the RN BH [31], the velocity of flow is finite at $ r=0 $ and converges to a constant as $ r\to \infty $, thus the Mach number vanishes at $ r=0 $ and becomes a finite number as $ r\to \infty $.

● Density and pressure in (3+1)D models

$ \rho(r) $ and $ p(r) $ for our RBHs in Sec. V have maximums at the same value of r, which is related to the roots of $ \xi''(r)=0 $. This property leads to a sharp discontinuity in the plot of EoS, which does not occur for the RN BH. Meanwhile, the pressures of our RBHs are negatively infinite at $ r=0 $, while they are positively infinite at this point for the RN BH.

● Velocity, Mach number, density, and pressure in the equatorial sections of $ (3+1) $D models

The velocity, Mach number, density, and pressure for the equatorial section of the $ (3+1) $D RBH constructed in Sec. VII exhibit similar configurations to those of a RN BH [31], which implies that the regularity of our model does not appear in the simulation of equatorial sections. Note that we have applied the method [31] to turn a non-compact dimension into a compact one in the construction of the relation between ξ and r.

● Velocity and Mach number in $ (2+1) $D models

The velocity and Mach number of our $ (2+1) $D RBH constructed in Sec. VIII differ greatly from those of the $ (2+1) $D RN BH. At first, our RBH may have only one horizon; then, the transonic flow of our RBH occurs outside the horizon. For the $ (2+1) $D RN BH, it has two horizons and its transonic flow occurs between the inner and outer horizons.

● Density and pressure in $ (2+1) $D models

The differences between our $ (2+1) $D RBH constructed in Sec. VIII and the $ (2+1) $D RN BH [31] are obvious in the density and pressure. For our $ (2+1) $D RBH, $ \rho(r) $ and $ p(r) $ are no longer divergent at $ r=0 $ and the EoS becomes smooth in the whole region of r, see Figs. 6(c), 6(d), and 6(e). Meanwhile, the behaviors of the two variables in the $ (2+1) $D RBH differ from those in the equatorial sections of $ (3+1) $D models, despite the fact that we use the identical compactification in both situations. This difference indicates that the dimensions of RBHs affect the properties of flow simulations.

Finally, we summarize that the EoS of fluid can be used to simulate realistic RBHs. Meanwhile, we show that the acoustic analogues of RBHs have apparently different features from those of SBHs, such as the RN BH, and that the differences are indeed caused by singularities. In other words, the acoustic gravity can be applied as a tool to study astronomic RBHs, which also offers a theoretical basis to investigate more phenomena of astronomic RBHs in a fluid.

Here we solve the differential inequalities appeared in this paper, in particular, in Sec. II. We start with solving the differential inequality, $ \xi \sigma''\le 2 \sigma' $, with the boundary conditions, $ \sigma(0)=0=\sigma'(0) $. It can be rewritten as

After considering the boundary conditions, we obtain $ 3 \sigma -\xi \sigma '\ge 0 $. Next, multiplying $ \xi^{-4} $ on its both sides, we derive

If we define $ \sigma_0:=\lim_{\xi\to 0} \sigma/\xi^3 $, we arrive at the solution,

which can also be obtained when one directly uses the differential form of the Grönwall-Bellman lemma [41].

Similarly, we can obtain $ \sigma\ge 0 $ from $ \sigma'\ge 0 $ when the boundary condition, $ \sigma(0)=0 $, is considered.

We emphasize that the differential inequalities' solutions provide the necessary condition for an RBH to meet the energy conditions, but not the sufficient one. For example, taking a bell-shaped function,

we can see that the inequality's solution, $ \sigma\ge 0 $, is satisfied but the WEC and DEC are violated because of $ \sigma'\ngeq 0 $. However, if $ \sigma\ge 0 $ is broken, then the WEC and DEC must be violated. In addition, the Hayward BH gives us another example, that is, it satisfies $ 0\le\sigma\le \sigma_0 \xi^3 $, but breaks the DEC. The reason comes from the characteristics of differential inequalities, i.e., a differential inequality signifies that all functions satisfying this differential inequality must be bounded by its solution, while the functions bounded by the solution may not necessarily meet the original differential inequality.

Now we give an explanation from the energy conditions by analyzing the local properties of a realistic RBH at its center, that is, why a realistic RBH cannot have a flat or an AdS core around its center. The similar discussion can be found in Ref. [24].

Summarizing the energy conditions Eq. (4), we have four differential inequalities in total,

Supposing $ \xi\ll 1 $, we expand σ by an asymptotic series,

meanwhile, we have the property, $ |{a_n \xi^n}|\gg|{a_{n+1} \xi^{n+1}}| $, as ξ approaches 0. Then substituting the above series into Eq. (B1), we obtain

If $ a_0\neq 0 $ and $ a_1\neq 0 $, the leading terms of Eqs. (B3) and (B5) lead to $ a_0> 0 $, which inevitably violates Eq. (B6), and the leading term of Eq. (B4) gives $ a_1< 0 $. If $ a_0= 0 $, the DEC must be broken because Eq. (B3) and Eq. (B4) lead to the results that are contradictory to each other. If $ a_0\neq 0 $ and $ a_1= 0 $, the leading term of Eq. (B4) gives $ a_m< 0 $, $ m>1 $, where m is the ordinal number of the first non-zero term.

In contrast, based on Ref. [33] we know

as ξ approaches 0. Thus the case of $ a_0> 0 $ and $ a_1< 0 $ implies

For $ a_1=0 $ but $ a_2\neq 0 $, we obtain $ R(0)> 0 $ and $ R''(0)< 0 $, i.e. $ \xi=0 $ is a local maximum. Moreover, if an RBH is Ricci flat at its center, $ R(0)=0 $, or it has an AdS core, $ R(0)< 0 $, its energy conditions must be violated around the center.

In order to derive the conditions given by Eq. (8), we substitute Eq. (7) into the DEC of Eq. (4). Since the DEC holds for $ \xi\in[0,\infty) $ and $ q\in[0,\infty) $, we obtain

Meanwhile, the asymptotic flatness demands

Then, combining Eqs. (C1) and (C2) and considering the positiveness of all parameters, we finally arrive at Eq. (8).

Now we derive the regularity conditions which have been applied in Sec. VIII. We write down the d-dimensional metric with the spherical symmetry [63],

where f is shape function, $ f=1-\mu\sigma(\xi)/\xi^{d-3} $, and μ is mass-like parameter. As we did in Ref. [33], we are going to use the following three curvatures,

to represent σ and its derivatives. With the help of the relations,

we arrive at

where $ \mathfrak{s}_{1,2}=\pm $ are two signs which are not much important for the discussion of the finiteness of curvatures. Meanwhile, it is not difficult to see that the metric Eq. (D1) has finite curvatures if σ has asymptotic relation $ \sigma\lesssim O(\xi^{d-1}) $.

We analyze the local properties of the differential equation Eq. (38),

at the two boundaries, $ \xi \to 0 $ and $ \xi \to \infty $, by means of the dominant balance [57, 58]. Here $ \xi\in [0,\infty) $ is the radial coordinate of RBHs, while r is radial coordinate of fluids with the spherical symmetry.

First of all, we consider the asymptoticity of the shape function as $ \xi\to 0 $,

According to the dominant balance, we can separate the discussions into three situations.

In the first case, we have

where Eq. (E2) gives the solution,

and $ {\tilde c}_1 $ is an integration constant. This solution becomes complex as $ \xi \to 0 $, which contradicts the physical requirement that $ r(\xi) $ must be real.

In the second case, the asymptotic relations become

The solution of Eq. (E5) is

However, $ r_+ $ is not consistent with the asymptotic assumption depicted by Eq. (E6). For $ R(0)\neq 0 $, we have a divergent limit as $ \xi\to 0 $,

As to $ r_- $, we note that it becomes complex if $ R(0)< 0 $, thus this second case should also be ruled out.

In the third case, we suppose

Eq. (E9) provides the solution,

where $ {\tilde c}_2 $ and $ {\tilde c}_3 $ are integration constants. This solution is consistent with the the asymptotic assumption Eq. (E10) because we have

However, we exclude $ r_- $ from physical solutions due to $ r_-<0 $ when $ \xi\sim 0^+ $.

Next, we turn to the asymptotic solutions as $ \xi\to \infty $. The shape function then becomes

where $ 0<n \le 1 $. Similarly, the discussion can also be separated into three situations.

In the first case, we have

The asymptotic assumption Eq. (E13) gives two solutions,

where $ {\tilde c}_4 $ and $ {\tilde c}_5 $ are integration constants. It can be verified that $ r_- $ contradicts to the asymptotic assumption Eq. (E14), while $ r_+ $ does not.

In the second case, we suppose

The first relation leads to

where $ {\tilde c}_6 $ is an integration constant. We can see that $ r_\pm $ converge to zero as $ \xi\to\infty $. This implies that the transformation between r and ξ is not injective, which is obvious because both $ \xi=0 $ and $ \xi=\infty $ map to the single point $ r=0 $. As a result, we eliminate this case.

In the last case, the asymptotic assumption involves

whose solutions are inevitably complex. Thus, this case is not in our consideration.