YITP-21-02

# Inner horizon instability and the unstable cores of regular black holes

Raúl Carballo-Rubio  Florida Space Institute, University of Central Florida,
12354 Research Parkway, Partnership 1, Orlando, FL, USA  Francesco Di Filippo
Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University, Kyoto 606-8502, Japan  SISSA - International School for
Advanced Studies, Via Bonomea 265, 34136 Trieste, Italy  Stefano Liberati
SISSA - International School for Advanced Studies, Via Bonomea 265, 34136
Trieste, Italy  IFPU - Institute for Fundamental Physics of the Universe, Via
Beirut 2, 34014 Trieste, Italy  INFN Sezione di Trieste, Via Valerio 2, 34127
Trieste, Italy  Costantino Pacilio  Dipartimento di Fisica, “Sapienza”
Università di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 5, 00185, Roma,
Italy  Matt Visser  School of Mathematics and Statistics, Victoria University
of Wellington, PO Box 600, Wellington 6140, New Zealand

###### Abstract

We discuss the generalization of the Ori model to black holes with generic
inner horizons (with arbitrary location and surface gravity), analyzing the
behavior of the metric around these inner horizons under the influence of both
a decaying ingoing flux of energy satisfying Price’s law and an outgoing flux
described by a null shell. We show that the exact solution for the late-time
evolution obtained by Ori is valid for all spacetimes in which the
Misner–Sharp mass is assumed to be linearly proportional to the time-dependent
metric perturbations introduced. We then consider generic situations in which
this linear proportionality is not assumed, obtaining asymptotic analytical
expressions and performing numerical integrations that show how the
exponential divergence associated with mass inflation is generically replaced
at late times by a power-law divergence. We emphasize that all these
geometries initially experience a first phase in which mass inflation proceeds
exponentially, and describe the physical implications that follow for generic
regular black holes. The formalism used also allows us to make some remarks
regarding the early-time transients associated with a positive cosmological
constant, known to modify the late-time behavior of ingoing perturbations from
Price’s law to an exponential decay. Finally we compare our analysis with that
in arXiv:2010.04226v1, and illustrate specific shortcomings in the latter work
that explain the differing results and, in particular, make it impossible for
that analysis to recover the well-known Ori solution for Reissner–Nordström
backgrounds.

##  I Introduction

Singularity theorems Penrose1964 ; Hawking1967 ; Hawking:1970  demonstrate
that, within the framework of standard general relativity, singular black
holes are unavoidably formed as the end state of gravitational collapse.
Observational tests of black hole spacetimes coming from the detection of
gravitational waves emitted by the merger of two black holes Abbott2016 ;
Abbott2016b ; Abbott2017 ; Abbott2017c ; TheLIGOScientific2017  are so far in
perfect agreement with the predictions of general relativity. However, there
is still room for alternatives to classical black holes Carballo-Rubio2018a ;
Cardoso:2019rvt , a consideration that becomes more pertinent given that it is
reasonable to assume that singularities will be tamed once quantum gravity
effects are taken into consideration, which naturally leads to some of these
alternatives Carballo-Rubio2019b ; Carballo-Rubio2019a .

Among the different classes of black hole mimickers, a very simple approach
consists of replacing the singularity with a regular core without necessarily
introducing substantial long-range modifications to the geometry. The price to
pay for such a conservative approach is the introduction of an inner horizon
Carballo-Rubio2019b ; Carballo-Rubio2019a  whose stability is not guaranteed.
In fact, it is well known that inner horizons in general relativity are
generically unstable under small perturbations (see e.g. Poisson1989 ;
Ori:1991zz ; Hamilton:2008zz ). Such instability has indeed been found in the
context of regular black holes, first for a particular model Brown2011  and
then for a generic regular black hole metric Carballo-Rubio:2018pmi . While
these analyses capture the main physical ingredients, they consider a somewhat
idealized scenario consisting of two non-interacting null shells on top of a
static background.

Recent work Bonanno:2020fgp , which aimed to extend the result of Ori’s work
for Reissner–Nordström black holes Ori:1991zz  to regular black hole
spacetimes, claims that a less idealized configuration of perturbations would
lead to a stable inner horizon. However, we point out a technical issue in the
analysis of Bonanno:2020fgp , and that the inner horizon of a regular black
hole is indeed generically unstable. To this end, in sec. II we extend Ori’s
calculation to a generic spherically symmetric metric. In sec. III we repeat
the analysis using the formalism of Bonanno:2020fgp , as this formalism may
provide a more straightforward physical interpretation and it is better suited
to numerically solve the exact time evolution of the geometry. In sec. IV we
consider the role of the cosmological constant, showing that also in this case
the instability is present, and that it could become sub-exponential at most
at late times (when the backreaction would be already non-negligible anyway).
Finally, in sec. V we solve numerically the differential equations governing
the growth of instability for many specific choices of spacetimes, and we show
that the solutions are in agreement with the analytical treatment of the
previous sections. Our conclusions are summarized in sec. VI.

##  II Generalized Ori model

In reference Ori:1991zz  Ori developed a model to describe the dynamical
evolution of the neighborhood of the inner horizon of a Reissner–Nordström
black hole under perturbations, in order to explicitly analyze the singularity
generated by the mass-inflation instability. The Ori model can be easily
generalized to any (singular or non-singular) spacetime with a metric of the
form

$ds^{2}=-f(v,r)dv^{2}+{2}dvdr+r^{2}d\Omega^{2}\,.$ (1)

Most of the derivation can be worked out only using the parameterization
$f(v,r)=1-2M(v,r)/r$ in terms of the Misner–Sharp mass $M(v,r)$. At this point
in the discussion we do not need to assume any specific form of the
Misner–Sharp mass, but only that its radial dependence is such that the
geometry contains an inner horizon or, equivalently, that $f(v,r)$ has en even
number of zeroes, and that all its dependence on $v$ enters through a single-
variable function $m(v)$ which coincides with the value of the Misner–Sharp
mass at large values of $r$, namely

$M(v,r)=M(m(v),r),$ (2)

with

$m(v)=\lim_{r\rightarrow+\infty}M(v,r).$ (3)

Well-known examples include the Vaidya–Reissner–Nordström spacetime in which
$M(m,r)=m-e^{2}/2r$, or the Hayward spacetime in which
$M(m,r)=mr^{3}/(r^{3}+2m\ell^{2})$ Hayward2006 .

Following reference Ori:1991zz , we now consider a dynamical model of a null
shell $\Sigma$ which divides the black hole into two subregions
$\mathcal{R}_{-}$ and $\mathcal{R}_{+}$, such that $\mathcal{R}_{-}$ is on the
same side of $\mathscr{I}^{-}$, as shown by the Penrose diagram in fig. 1.

Figure 1: Relevant sections of the Penrose diagram of a regular black hole.
The shell $\Sigma$ divides the spacetime in two regions $\mathscr{R}_{-}$ and
$\mathscr{R}_{+}$. $\mathscr{H}$ denotes the event horizon while $\mathscr{C}$
denotes the Cauchy horizon.

Therefore we can parametrize $\mathcal{R}_{-}$ by the same advanced null
coordinate $v_{-}$ of $\mathscr{I}^{-}$, while we shall choose a distinct null
coordinate $v_{+}$ to parametrize $\mathcal{R}_{+}$. By continuity, the radial
coordinate $r$ is the same on both regions, while the mass function $m(v)$
depends on the region, so that we shall distinguish $m_{-}(v_{-})$ and
$m_{+}(v_{+})$. On the other hand, the matching conditions imply that both
$v_{-}$ and $v_{+}$ can be expressed _on $\Sigma$_ in terms of a single affine
coordinate $\lambda$ Barrabes1991 . We will define $\lambda$ in such a way
that it is negative and $\lambda\to 0$ when $v_{-}\to\infty$. Additionally,
the radius of the shell can be expressed as a function $R(\lambda)$ of the
affine parameter.

The functional form of $m_{-}(v_{-})$ is dictated by Price’s law Price1971 ;
Price1972  to be

$m_{-}(v_{-})=m_{0}+\delta m=m_{0}-\frac{\beta}{v_{-}^{p}}\,.$ (4)

When $v_{-}\to\infty$, the shell crosses the corresponding inner horizon
$r_{0}$, which is implicitely defined by

$\displaystyle\lim_{\lambda\to 0}R=r_{0};$ (5a) $\displaystyle\lim_{\lambda\to
0}f_{-}(\lambda,R)=0;$ (5b)

where $R\equiv R(\lambda)$ is the radius of the shell and
$f_{-}(\lambda,R)=f(M(m_{-}(\lambda),R),R)$. It is convenient to define the
surface gravity in $\mathcal{R}_{-}$ at $\lambda=0$ as

$\kappa_{0}=\frac{1}{2}\lim_{\lambda\to 0}\left.\frac{\partial
f_{-}(\lambda,r)}{\partial r}\right|_{r=R(\lambda)}\,.$ (6)

We now set up the basic ingredients for the discussion of the generalized Ori
model. It can be shown that the matching conditions on $\Sigma$ in eqs.
(7)–(9) of Ori:1991zz  generalize, respectively, to

$\displaystyle z{{}_{i}}R^{\prime}=\frac{R}{2}f{{}_{i}}(v{{}_{i}},R);$ (7a)
$\displaystyle v{{}_{i}}(\lambda)=\int^{\lambda}d\lambda\frac{R}{z{{}_{i}}};$
(7b) $\displaystyle
z{{}_{i}}(\lambda)=Z{{}_{i}}+\frac{1}{2}\int^{\lambda}d\lambda\left(f{{}_{i}}(\lambda,R)+R\frac{\partial
f{{}_{i}}}{\partial r}(\lambda,R)\right);$ (7c)

where $z{{}_{i}}=R/v{{}_{i}}^{\prime}$, the index $i$ takes two values $+$ and
$-$ (for the two regions $\mathcal{R}_{+}$ and $\mathcal{R}_{-}$,
respectively), and a prime denotes differentiation with respect to $\lambda$.
The quantity $Z{{}_{i}}$ is an integration constant, while we omitted a
similar integration constant in (7b) because it is irrelevant. These equations
are valid on both sides of the ingoing null shell. Now, in region
$\mathcal{R}_{-}$, we can expand eq. (7c) around $\lambda=0$ to obtain

$z_{-}(\lambda)=Z_{-}-r_{0}|\kappa_{0}|\lambda+\mathcal{O}(\lambda^{2})\,.$
(8)

However, given that $\lim_{\lambda\to 0}v_{-}=\infty$, it follows from eq.
(7b) that $Z_{-}=0$, otherwise $v_{-}$ would tend to a constant proportional
to $1/Z_{-}$, and therefore (hereafter, for simplicity, we will use $v$ in
place of $v_{-}$)

$v(\lambda)=-\frac{1}{|\kappa_{0}|}\ln{|\lambda|}+\mathcal{O}(\lambda^{0})\,.$
(9)

Notice that the derivation of (8) and (9) proceeded exactly as in reference
Ori:1991zz . Furthermore, eq. (7a) in region $\mathcal{R}_{-}$ can be written
as a differential equation in $v$ as

$\frac{dR(v)}{dv}=\frac{1}{2}f_{-}({m_{-}(v)},R(v)).$ (10)

We can perform a Taylor expansion of $f_{-}$ around the values $m_{-}=m_{0}$
and $R=r_{0}$ at $v=\infty$ to obtain

$\frac{d\delta R(v)}{dv}=-\frac{A}{r_{0}}\delta m(v)-|\kappa_{0}|\delta
R(v){+\mathcal{O}(\delta m^{2},\delta m\delta R,\delta R^{2})},$ (11)

where we have written ${R(v)=}r_{0}+\delta R(v)$ and $A=\partial
M_{-}/\partial m_{-}|_{{m_{-}=m_{0},r=r_{0}}}$. The general solution to the
differential equation above (neglecting subleading orders) is given by

$\delta
R(v)=c_{1}e^{-|\kappa_{0}|v}-\frac{A}{r_{0}}e^{-|\kappa_{0}|v}\int\text{d}\tilde{v}\,\delta
m(\tilde{v})e^{|\kappa_{0}|\tilde{v}},$ (12)

where $c_{1}$ is an integration constant. This equation is valid for a generic
$\delta m(v)$. If we further consider the functional form of $\delta m(v)$
from eq. (4) we obtain

$\begin{split}\delta
R(v)&={c_{1}}e^{-|\kappa_{0}|v}+\frac{A\beta}{|\kappa_{0}|r_{0}v^{p}}\sum_{k=0}^{\infty}\left[\frac{(p+k-1)!}{(p-1)!}\frac{1}{|\kappa_{0}|^{k}v^{k}}\right]\\\
&=\frac{A\beta}{|\kappa_{0}|r_{0}v^{p}}\left\\{1+\frac{p}{|\kappa_{0}|v}+\frac{p(p+1)}{|\kappa_{0}|^{2}v^{2}}+\mathcal{O}(v^{-3})\right\\}\,.\end{split}$
(13)

Let us now consider the evolution of the metric functions in region
$\mathcal{R}_{+}$ on its boundary $\Sigma$, hence as functions of $\lambda$
only (alternatively, $v$). One can manipulate eqs. (7a) and (7c) for $i=+$ to
obtain

$\frac{df_{+}{(\lambda,R)}}{d\lambda}=R^{\prime}\left.\frac{\partial
f_{+}{(\lambda,r)}}{\partial
r}\right|_{r=R}+\frac{R^{\prime\prime}}{R^{\prime}}f_{+}{(\lambda,R)},$ (14)

or equivalently, as a differential equation in $v$,

$\begin{split}\frac{dM_{+}{(v,R)}}{dv}&=R_{,v}\left.\frac{\partial
M_{+}{(v,r)}}{\partial
r}\right|_{r=R}-\left(\frac{\lambda_{,vv}}{\lambda_{,v}}-\frac{R_{,vv}}{R_{,v}}\right)\left(M_{+}{(v,R)}-\frac{R}{2}\right)\\\
&=R_{,v}\left.\frac{\partial M_{+}{(v,r)}}{\partial
r}\right|_{r=R}+\left(|\kappa_{0}|-\frac{p+1}{v}{+\mathcal{O}(v^{-2})}\right)\left(M_{+}{(v,R)}-\frac{R}{2}\right)\,.\end{split}$
(15)

In the equation above the $v$ subindices denote differentiation with respect
to this variable. In the case of Reissner–Nordström,

$\left.\frac{\partial M_{+}{(v,r)}}{\partial
r}\right|_{r=R}=\frac{e^{2}}{2{R}^{2}},$ (16)

so the first term of the right-hand side of eq. (15) does not depend on
$M_{+}$ and the equation displays (exponentially) inflating solutions with
leading behaviour

$M_{+}{(v,R(v))}\propto\frac{e^{|\kappa_{0}|v}}{v^{p+1}}\,.$ (17)

This is precisely the exact solution obtained by Ori Ori:1991zz . More
generally, we see that this result by Ori is recovered unchanged for a wider
family of geometries satisfying the linear ansatz
$M(v,r)=g_{1}(r)m(v)+g_{2}(r)$, which includes for instance the Bardeen metric
for a regular black hole.

Not all metrics potentially of interest are covered by this linear ansatz. For
instance, in the case of the regular black hole described by the Hayward
metric Hayward2006 , the functional form of $M(m,r)=mr^{3}/(r^{3}+2m\ell^{2})$
implies that

$\frac{\partial M}{\partial r}=\frac{6\ell^{2}}{r^{4}}M^{2}.$ (18)

While these more general situations can be analyzed using the equations above,
it is not difficult to see, either using eq. (15) or eq. (7) directly, that
the linear ansatz for $M(v,r)$ is necessary in order to derive Ori’s solution
in eq. (17). We will discuss these more general situations in more detail,
using both analytical and numerical methods, after an equivalent formalism is
described in the next section.

##  III Alternative derivation

It is instructive to derive the same results with a slightly different,
although equivalent, formalism. Moreover, this will allow us to critically
analyze the results presented in reference Bonanno:2020fgp . The analysis of
this section will closely follow the discussion in that work up to the point
where we correct a technical error therein. We will then discuss how this
technical error invalidates the conclusions presented in reference
Bonanno:2020fgp .

The setting is the same as in the Ori model discussed in the previous section,
and we will therefore use the same notation. The main difference in this
section is the way in which the gluing of the two spacetime regions along
$\Sigma$ is dealt with. Instead of using the existence of a common (for both
regions) affine parameter $\lambda$, we will focus on the junction conditions
at the shell Barrabes1991

$\left[T_{\mu}\,^{\nu}s^{\mu}s_{\nu}\right]=0,$ (19)

where $T_{\mu}\,^{\nu}$ is the effective stress energy tensor obtained by
imposing Einstein’s equations on both regions of the spacetime, and

$s^{\mu}=\left(2/f_{\pm},1,0,0\right)$ (20)

is the outgoing null normal to the shell. The relevant components of the
stress energy tensor are

$T_{v}\,^{v}=-4\pi r^{2}\frac{\partial M}{\partial r},\qquad T_{v}\,^{r}=4\pi
r^{2}\frac{\partial M}{\partial v},\qquad T_{r}\,^{r}=-4\pi
r^{2}\frac{\partial M}{\partial v}.$ (21)

Eq. (19) then becomes

$\left.\frac{1}{f_{+}^{2}}\frac{\partial M_{+}(v_{+},r)}{\partial
v_{+}}\right|_{r=R(v_{+})}=\left.\frac{1}{f_{-}^{2}}\frac{\partial
M_{-}(v_{-},r)}{\partial v_{-}}\right|_{r=R(v_{-})}.$ (22)

We can eliminate the $v_{+}$ dependence from this equation noting that along a
null trajectory

$dr=\frac{1}{2}f_{\pm}dv_{\pm},$ (23)

or, equivalently,

$\frac{dv_{+}}{dv_{-}}=\frac{f_{-}}{f_{+}}.$ (24)

Using this equation in order to relate the total derivatives of $M_{+}$ with
respect to $v_{+}$ and $v_{-}$ allows us to write eq. (22) as

$\left.\frac{1}{f_{+}}\frac{\partial M_{+}}{\partial
v}\right|_{r=R(v)}=\left.\frac{1}{f_{-}}\frac{\partial M_{-}}{\partial
v}\right|_{r=R(v)}$ (25)

where, as in the previous section, we use $v$ in place of $v_{-}$. We
emphasize that, when evaluating eq. (25) on the shell $\Sigma$, the partial
derivatives with respect to $v$ must be taken before imposing $r=R(v)$ on
$\Sigma$ (one would be calculating total derivatives otherwise, thus changing
the content of the equation). Crucially, this observation lies at the core of
the technical flaw in Bonanno:2020fgp , as it will be discussed in more detail
at the end of this section: indeed, in deriving their eq. (13), the authors
imposed the condition $r=R(v)$ inside the partial derivative signs in eq.
(25).

Given that we are following the evolution of $M(v,r)$ on the trajectory
$r=R(v)$ or, in other words, we are solving a coupled system of differential
equations for $R(v)$ and $M(v,R(v))$, we cannot directly integrate eq. (25)
because it involves partial derivatives with respect to $v$. Instead we need
to consider an equivalent equation involving total derivatives. To this end,
let us note that

$\frac{dM_{+}\left(v,R(v)\right)}{dv}=\left.\frac{\partial
M_{+}{(v,r)}}{\partial v}\right|_{r=R(v)}+\left.\frac{\partial
M_{+}{(v,r)}}{\partial r}\right|_{r=R(v)}\frac{dR{(v)}}{dv}.$ (26)

Using eq. (25), the above equation can be rewritten as

$\begin{split}\frac{dM_{+}{\left(v,R(v)\right)}}{dv}=\frac{f_{+}}{f_{-}}\left.\frac{\partial
M_{-}{(v,r)}}{\partial v}\right|_{r=R(v)}+R_{,v}\left.\frac{\partial
M_{+}{(v,r)}}{\partial r}\right|_{r=R(v)}\,.\end{split}$ (27)

Consider now the identity

$\frac{1}{f_{-}}{\left.\frac{\partial M_{-}{(v,r)}}{\partial
v}\right|_{r=R(v)}}=\frac{1}{2R_{,v}}{\left.\frac{\partial
M_{-}(v,r)}{\partial m_{-}}\right|_{r=R(v)}}\frac{dm_{-}(v)}{dv}$ (28)

which holds true on $\Sigma$ due to the functional form of the Misner–Sharp
mass, $M(v,r)=M(m(v),r)$, and the equation of motion for the null shell
$\Sigma$, namely eq. (10). Moreover, in eq. (13) we obtained the solution of
the equation of motion for the null shell, which using the derivative of eq.
(4) can be written as

$R_{,v}=-\frac{A}{|\kappa_{0}|^{2}r_{0}}{\frac{dm_{-}(v)}{dv}}\left(|\kappa_{0}|+\frac{p+1}{v}\right)+\mathcal{O}(v^{-p-3}).$
(29)

It is then straightforward to see that the derivatives of the function
$m_{-}(v)$ cancel on the right-hand side of eq. (28), while the definition
$A=\partial M_{-}/\partial m_{-}|_{m_{-}=m_{0},r=r_{0}}$ allow us to further
simplify the resulting expression to

$\begin{split}\frac{1}{f_{-}}{\left.\frac{\partial M_{-}{(v,r)}}{\partial
v}\right|_{r=R(v)}}&=-\frac{r_{0}{|\kappa_{0}|^{2}}}{2}\left(|\kappa_{0}|+\frac{p+1}{v}\right)^{-1}+\mathcal{O}(v^{-2})\\\
&=-\frac{r_{0}}{2}\left(|\kappa_{0}|-\frac{p+1}{v}\right)+\mathcal{O}(v^{-2})\,.\end{split}$
(30)

As a consequence we obtain

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}{=}\left(|\kappa_{0}|-\frac{p+1}{v}{+\mathcal{O}\left(v^{-2}\right)}\right)\left\\{M_{+}{\left(v,R(v)\right)}-\frac{r_{0}}{2}\right\\}+R_{,v}\left.\frac{\partial
M_{+}}{\partial r}\right|_{r={R(v)}}.$ (31)

If $\partial M_{+}/\partial r$ is proportional to $M_{+}$, namely
$M_{+}(v,r)=g_{1}(r)m_{+}(v)+g_{2}(r)$, the last term becomes subdominant at
late times due to being multiplied by $R_{,v}=\mathcal{O}(v^{-p-1})$ terms and
we have

$M_{+}{\left(v,R(v)\right)}={c_{2}}\frac{e^{|\kappa{{}_{0}}|v}}{{v^{p+1}}}{\left[1+\mathcal{O}\left(v^{-1}\right)\right]+\mathcal{O}(v^{0})},$
(32)

where $c_{2}$ is an integration constant. Hence, we reproduce Ori’s result.
However, if $\partial M_{+}/\partial r$ is not linear in $M_{+}$, the
differential equation (31) becomes nonlinear and the last term on its right-
hand side may drive the evolution at late times, thus modifying the
exponential behavior associated with mass inflation. Hence, the analysis must
proceed in a case-by-case basis for these more general situations.

For instance let us consider the case where $\partial M_{+}/\partial r$ is
polynomial in $M_{+}$. We would get

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}{=}\left(|\kappa_{0}|-\frac{p+1}{v}{+\mathcal{O}\left(v^{-2}\right)}\right)\left\\{M_{+}{\left(v,R(v)\right)}-\frac{r_{0}}{2}\right\\}{-\frac{\gamma}{v^{p+1}}M_{+}^{n}\left(v,R(v)\right)}.$
(33)

for some given constants $\gamma$ and $n$ depending on the specific metric
under consideration. The behavior of the solution of this differential
equation depends on the sign of $\gamma$. It is possible to show, for
$\gamma<0$, that $M_{+}{\left(v,R(v)\right)}$ diverges for a finite value of
$v$. On the other hand, for $\gamma>0$ we can solve this differential equation
in different regimes. At early times the first term on the right hand side
dominates the dynamics, and this equation predicts an exponential growth of
$M_{+}{\left(v,R(v)\right)}$, up to a critical time $v_{0}$ for which

$M_{+}^{n-1}(v_{0},R(v_{0})){\approx}\frac{|\kappa_{0}|}{\gamma}v_{0}^{p+1}.$
(34)

From that point on the first and the second term of the right hand side assume
comparable values and the growth is polynomial rather than exponential,

$M_{+}^{n-1}{\left(v,R(v)\right)\propto}v^{p+1},\qquad{v\gg v_{0}}.$ (35)

As a concrete example, let us consider Hayward’s metric for which the relation
(18) allows us to write

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}\simeq|\kappa_{0}|M_{+}{\left(v,R(v)\right)}-\gamma\frac{M_{+}^{2}{\left(v,R(v)\right)}}{v^{p+1}}-\frac{r_{0}|\kappa_{0}|}{2},$
(36)

where the symbol $\simeq$ means that we are neglecting subdominant terms in
the $v\rightarrow\infty$ limit, and

$\gamma=\frac{6p\beta\ell^{2}}{|\kappa_{0}|}\frac{r_{0}}{\left(2m_{0}\ell^{2}+r_{0}^{3}\right)^{2}}.$
(37)

By neglecting the last term on the right hand side, the solution of (36) is
given by

$M_{+}(v,R(v)){\approx}\frac{e^{|\kappa_{0}|v}v^{p}}{c_{1}v^{p}-(-v)^{p}|\kappa|^{p-1}\gamma\Gamma(-p,-|\kappa_{0}|v)}.$
(38)

Expanding the incomplete Gamma function at late times we obtain

$M_{+}(v,R(v)){\approx}\frac{|\kappa_{0}|}{\gamma}v^{p+1},$ (39)

in perfect agreement with eq. (35). Therefore, the growth rate of the
instability is slowed down, becoming polynomial rather than exponential at
late times. However, it is reasonable to expect that this polynomial behaviour
will not have any physical consequence, as it occurs beyond the validity of
our linear analysis and after a certain timescale during which the
Misner–Sharp mass has already become substantially large in an exponential
fashion. Indeed, assuming that the regularization scale $\ell$ satisfies the
condition $\ell\ll m_{0}$, (as would be the case if, for instance, $\ell$ is
given by the Planck length), and substituting (37) into (34) we see that
$M_{+}$ grows exponentially up to the critical value

$M_{+}(v,R(v))\approx\frac{|\kappa_{0}|^{2}\left(2m_{0}\ell^{2}+r_{0}^{3}\right)^{2}}{r_{0}\ell^{2}}\frac{v_{0}^{p+1}}{6p\beta}\sim\frac{m_{0}^{2}}{\ell}\frac{v_{0}^{p+1}}{6p\beta},$
(40)

where we have used

$r_{0}\sim|\kappa{{}_{0}}|^{-1}\sim\ell.$ (41)

Eq. (40) implies

$\frac{M_{+}(v,R(v))}{m_{0}}\sim\frac{m_{0}}{\ell}\frac{v_{0}^{p+1}}{6p\beta}\gg
1.$ (42)

So the transition between exponential and polynomial regimes happens when the
backreaction of the perturbation on the background geometry is already very
large due to an initial phase of exponential growth of the Misner–Sharp mass,
and at this stage we expect nonlinear effects to become relevant and supersede
our linear treatment. It is also important to keep in mind that, in any case,
$M_{+}(v,R(v))$ is always divergent.

Let us now pause for a brief moment to understand why reference
Bonanno:2020fgp  fails to capture the divergent nature of $M_{+}(v,R(v))$ that
we have described above. As previously mentioned, in deriving eq. (13) of
reference Bonanno:2020fgp , the condition $r=R(v)$ was imposed inside the
derivatives of eq. (25). That is, the partial derivative is effectively
replaced with a total derivative. Therefore, eq. (27) is (incorrectly)
modified into

$\displaystyle\frac{dM_{+}\left(v,R(v)\right)}{dv}$
$\displaystyle=\frac{f_{+}}{f_{-}}\left(\left.\frac{\partial
M_{-}(v,r)}{\partial v}\right|_{r=R(v)}+\left.\frac{\partial
M_{-}(v,r)}{\partial r}\right|_{r=R(v)}\frac{dR(v)}{dv}\right)$
$\displaystyle=\frac{f_{+}}{f_{-}}\left(\left.\frac{\partial
M_{-}(v,r)}{\partial v}\right|_{r=R(v)}+\frac{f_{-}}{2}\left.\frac{\partial
M_{-}(v,r)}{\partial r}\right|_{r=R(v)}\right).$ (43)

In order to extract the leading behavior encapsulated in this equation let us
note that the definition of the surface gravity of the inner horizon,

${\left.\frac{\partial\left[M_{-}(v,r)/r\right]}{\partial
r}\right|_{r=R(v)}}=|\kappa_{0}|+\mathcal{O}\left(v^{-p}\right),$ (44)

implies that

${\left.\frac{\partial M_{-}(v,r)}{\partial
r}\right|_{r=R(v)}}=\frac{M_{-}{\left(v,R(v)\right)}}{r_{0}}+r_{0}|\kappa_{0}|{+\mathcal{O}\left(v^{-p}\right)}.$
(45)

Using the relation above, together with eqs. (28) and (29), we can therefore
write

$\displaystyle\left.\frac{\partial M_{-}(v,r)}{\partial
v}\right|_{r=R(v)}+\frac{f_{-}}{2}\left.\frac{\partial M_{-}(v,r)}{\partial
r}\right|_{r=R(v)}$ $\displaystyle=$ $\displaystyle\left.\frac{\partial
M_{-}(v,r)}{\partial
v}\right|_{r=R(v)}-\frac{1}{r_{0}|\kappa_{0}|}\left.\frac{\partial
M_{-}(v,r)}{\partial
v}\right|_{r=R(v)}\left\\{\frac{M_{-}\left(v,R(v)\right)}{r_{0}}+r_{0}|\kappa_{0}|+\mathcal{O}\left(v^{-p}\right)\right\\}$
$\displaystyle=$
$\displaystyle-\frac{1}{r_{0}^{2}|\kappa_{0}|}\left.\frac{\partial
M_{-}(v,r)}{\partial
v}\right|_{r=R(v)}\left\\{M_{-}\left(v,R(v)\right)+\mathcal{O}\left(v^{-p}\right)\right\\}.$
(46)

This implies that the leading behavior of the incorrect expression (III) is
encapsulated in the equation

$\frac{dM_{+}\left(v,R(v)\right)}{dv}\simeq\frac{M_{-}\left(v,R(v)\right)}{2r_{0}}f_{+}=\frac{M_{-}\left(v,R(v)\right)}{2r_{0}}-\frac{M_{-}\left(v,R(v)\right)}{r_{0}^{2}}M_{+}\left(v,R(v)\right).$
(47)

This differential equation does not have growing exponential solutions, but,
due to the sign of the second term, decaying ones. In fact, its general
solution at leading order is

$M_{+}\left(v,R(v)\right)\simeq\frac{r_{0}}{2}+c_{1}e^{-M_{0}v/r_{0}^{2}},$
(48)

where we have defined the constant $M_{0}=M_{-}(m_{0},r_{0})$. In particular,
$M_{+}\left(v,R(v)\right)$ tends to the finite value $r_{0}/2$ for
$v\gg{r_{0}^{2}/M_{0}}$. Finally, let us note, as a further evidence of
something wrong with eq. (III), that, exactly for the aforementioned reasons
(and also due to the fact that for the discussion above to work it is not
necessary to select an ansatz for the Misner–Sharp mass $M(m,r)$), it also
fails to reproduce the well known instability of the inner horizon of a
standard Reissner–Nordström black hole.

##  IV Role of the cosmological constant

In the discussion so far we have ignored the role of the cosmological
constant. There is a very simple and intuitive reason behind this decision.
The instability we have discussed originates close to the inner horizon, where
the energy density is Planckian and the role of the cosmological constant is
completely negligible. However, in the presence of a non-zero cosmological
constant, the late time behavior of the perturbation is no longer described by
Price’s law (4), but has an exponential decay barreto1997 ; dyatlov2012
characterized by

$m_{-}{(v)=m_{0}+\delta m(v)}=m_{0}-\alpha e^{-\omega_{I}v},$ (49)

where $\omega_{I}>0$ is the imaginary part of the least damped quasinormal
mode. This modification only affects the extremely late time behavior of the
perturbation, and the very small value of the cosmological constant of our
universe implies that the fall-off of the perturbations follows the polynomial
Price’s law of eq. (4) for a very long time. Therefore in our universe the
instability described in the previous sections will arise well before the
cosmological constant could play any significant role. However, it is
interesting from the theoretical side to investigate whether or not the values
of the cosmological constant and of the black hole parameters can be fine
tuned, in such a way that the inner horizon does not became unstable.

Using eq. (49), eq. (12) allows us to obtain

$\delta
R(v)=c_{1}e^{-|\kappa_{0}|v}+\frac{A\alpha}{r_{0}(|\kappa_{0}|-\omega_{I})}e^{-\omega_{I}v}$
(50)

where $c_{1}$ is an integration constant. This modifies the rate at which
$f_{-}$ approaches zero which, let us recall, due to eq. (10) is given by the
derivative of the expression right above. Eq. (27) can now be written as

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}\simeq\frac{f_{+}A\alpha\omega_{I}e^{-\omega_{I}v}}{-2|\kappa_{0}|c_{1}e^{-|\kappa_{0}|v}-2\omega_{I}A\alpha
e^{-\omega_{I}v}{/r_{0}(|\kappa_{0}|-\omega_{I})}}+\frac{f_{-}}{2}{\left.\frac{\partial
M_{+}(v,r)}{\partial r}\right|_{r=R(v)}}.$ (51)

We have to analyze the cases $|\kappa_{0}|>\omega_{I}$ and
$|\kappa_{0}|<\omega_{I}$ independently. For $|\kappa_{0}|>\omega_{I}$ we have

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}{\simeq}\frac{r_{0}(|\kappa_{0}|-\omega_{I})\left(2M_{+}{\left(v,R(v)\right)}/r_{0}-1\right)}{2}+\frac{f_{-}}{2}{\left.\frac{\partial
M_{+}(v,r)}{\partial r}\right|_{r=R(v)}}.$ (52)

If $\partial_{r}M_{+}$ is linear in $M_{+}$ the last term can be neglected
leading to

$M_{+}{\left(v,R(v)\right)}\propto e^{(|\kappa_{0}|-\omega_{I})v}.$ (53)

This result has a simple physical interpretation. We have the same blueshift
factor $e^{|\kappa_{0}|v}$ that we had in the absence of the cosmological
constant, while the polynomial decay is replaced by the exponential decay
described in equation (49).

Let us now study the case in which $\partial_{r}M_{+}$ is polynomial in
$M_{+}$. We have

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}\simeq(|\kappa_{0}|-\omega_{I})M_{+}{\left(v,R(v)\right)}-\gamma
e^{-\omega_{I}v}M_{+}^{n}{\left(v,R(v)\right)}-r_{0}(|\kappa_{0}|-\omega_{I}).$
(54)

The late-time behavior of the solution of this differential equation describes
an exponential instability

$M_{+}{\left(v,R(v)\right)}\propto e^{\omega_{I}v/(n-1)}.$ (55)

Finally, for $|\kappa_{0}|<\omega_{I}$ we have

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}{\simeq}\frac{\left(2M_{+}{\left(v,R(v)\right)}/r{{}_{0}}-1\right)}{2|\kappa_{0}|}A\alpha
e^{-(\omega_{I}-|\kappa|_{0})v}-2|\kappa_{0}|e^{-|\kappa_{0}|v}{\left.\frac{\partial
M_{+}(v,r)}{\partial r}\right|_{r=R(v)}}.$ (56)

The exact form of the solution depends on the specific choice of the metric.
However, the asymptotic behavior is generically regular (that is, the
Misner–Sharp mass does not diverge asymptotically). The only scenario left to
consider is the special case $|\kappa_{0}|=\omega_{I}$. While this scenario
requires some fine tuning to be realized, let us study it for the sake of
completeness. The general solution in eq. (12) now reads

$\delta R=c_{1}e^{-|\kappa_{0}|v}+\frac{{\alpha}}{r_{0}}e^{-|\kappa_{0}|v}v,$
(57)

leading to

$\frac{dM_{+}{\left(v,R(v)\right)}}{dv}\approx\frac{A\alpha}{v}M_{+}{\left(v,R(v)\right)}-\frac{A\alpha
r_{0}}{2v}+\frac{f_{{-}}}{2}{\left.\frac{\partial M_{+}(v,r)}{\partial
r}\right|_{r=R(v)}}$ (58)

If $\partial_{r}M_{+}$ is linear in $M_{+}$ we get

$M_{+}{\left(v,R(v)\right)\propto}v^{\alpha A}.$ (59)

Therefore, in this case $M_{+}{\left(v,R(v)\right)}$ exhibits a polynomial
divergence.

##  V Numerical analysis

In this section we provide some explicit numerical solutions of the equations
discussed above for some particular choices of the function $f(v,r)$. We can
use two equivalent methods for the numerical integrations below. The first one
is solving eqs. (10) and (25) for the variables $R(v)$ and $m_{+}(v)$. The
other method considers instead the variables $R(v)$ and
$M_{+}\left(r,R(v)\right)$, solving the system of differential equations given
by eqs. (10) and (27). The main advantage of the second method is that it uses
the Misner–Sharp mass as one of the variables, while the function $m_{+}(v)$
used in the first method has no clear physical meaning. We will nevertheless
use the first method in order to highlight some of the issues that can arise
in specific examples due to this choice of variables.

### Reissner–Nordström black hole

In order to connect our discussion with Ori’s original work, let us start by
analyzing the case of a Reissner-Nordström black hole. The metric function
reads

$f_{\pm}(v,r)=1-\frac{2m_{\pm}(v)}{r}+\frac{e^{2}}{r^{2}}.$ (60)

Then, eq. (25) becomes

$\frac{m_{+}^{\prime}(v)}{R^{2}-2R\,m_{+}(v)+q}=\frac{m_{-}^{\prime}(v)}{R^{2}-2R\,m_{-}(v)+q}.$
(61)

The asymptotic time evolution of the mass parameter in the past of the shell
is determined by Price’s law (4). Eq. (61) can now be solved numerically. The
results of the numerical integration are plotted in fig. 2. As expected, the
system develops a mass inflation instability and the Misner–Sharp quasilocal
mass diverges exponentially.

Figure 2: Left: Numerical evolution of the mass parameter $m_{+}$ and of the
metric functions $f_{\pm}$ of a Reissner–Nordström black hole. Right:
Comparison between the numerical and the analytical evolution of the
Misner–Sharp mass $M_{+}$. In both plots we have used the parameters
$\beta=1$, $p=12$, $m_{0}=10$, $q=5$ with the initial conditions $R(v=1)=5$
and $m_{+}(v=1)=m_{0}+1$.

### Regular black hole: Hayward

Having reproduced Ori’s results for the case of a Reissner–Norström black
hole, we are now confident in extending the computation for different choices
of regular black hole geometries. The first metric that we are going to
analyze is the Hayward regular black hole metric introduced in Hayward2006 ,

$f_{\pm}(v,r)=1-\frac{2m_{\pm}(v)r^{2}}{r^{3}+2m_{\pm}(v)\ell^{2}}.$ (62)

The analog of eq. (61) now reads (notice the difference with respect to the
incorrect eq. (13) in Bonanno:2020fgp )

$\frac{m_{+}^{\prime}(v)}{\left(2\ell^{2}m_{+}(v)+R^{3}\right)\left(2\left(\ell^{2}-R^{2}\right)m_{+}(v)+R^{3}\right)}=\frac{m_{-}^{\prime}(v)}{\left(2\ell^{2}m_{-}(v)+R^{3}\right)\left(2\left(\ell^{2}-R^{2}\right)m_{-}(v)+R^{3}\right)}$
(63)

We now impose Price’s law (4) for $m_{-}(v)$, and solve numerically for
$m_{+}(v)$. The result of the integration is shown in fig. 3. The mass
parameter $m_{+}(v)$ has a vertical asymptote and, therefore, diverges for
finite value of $v$. This mathematical divergence does not have any physical
implications as the Misner–Sharp mass stays finite. After the asymptote,
$m_{+}(v)$ changes sign and its late time limit is given by

$\lim_{v\rightarrow+\infty}m_{+}(v)=-\frac{r^{3}}{2\ell^{2}}.$ (64)

Figure 3: Left: Numerical evolution of the mass parameter $m_{+}$ and of the
metric functions $f_{\pm}$ of a Hayward regular black hole. Right: Numerical
evolution of the Misner–Sharp mass. In both plots we picked the parameters
$\beta=1$, $p=12$, $m_{0}=10$, $\ell=0.5$ with the initial conditions
$R(v=1)=5$ and $m_{+}(v=1)=m_{0}+1$.

As a consequence, we can see from eq. (62) that the Misner–Sharp mass diverges
for large values of $v$. However, as shown in the right panel of fig. 3 this
divergence is polynomial rather than exponential. As discussed above, an
alternative method that can be used to obtain the same results without dealing
with the (unphysical) divergence of $m_{+}(v)$ at finite $v$ is changing the
variables used for the numerical integration to the physical variables $R(v)$
and the Misner–Sharp mass $M_{+}\left(r,R(v)\right)$, and solving the system
of differential equations given by eqs. (10) and (27). This leads directly to
the right plot in fig. 3.

### Regular black hole: Bardeen

As our next example, let us analyze what happens if we consider the Bardeen
metric

$f_{\pm}=1-\frac{2m_{\pm}(v)r^{2}}{\left(r^{2}+\ell^{2}\right)^{3/2}}.$ (65)

With this choice we obtain

$\frac{m^{\prime}_{+}(v)}{\left(\ell^{2}+R^{2}\right)^{3/2}-2R^{2}m_{+}(v)}=\frac{m^{\prime}_{-}(v)}{\left(\ell^{2}+R^{2}\right)^{3/2}-2R^{2}m_{-}(v)}$
(66)

As before, we numerically integrate the system of equations for $m_{+}$ and
$R$, and we report the results in fig. 4. Once again, the numerical
integration is in perfect agreement with the analytic result: in this case,
the mass parameter $m_{+}(v)$ diverges exponentially leading to an exponential
divergence of the Misner–Sharp mass $M_{+}\left(v,R(v)\right)$. It is
important to remark that the different behaviors of $m_{+}(v)$, a variable
that has no direct physical meaning, in the Hayward and Bardeen cases is
simply caused by the different functional forms of the Misner–Sharp $M(m,r)$
as a function of $m(v)$. Indeed, in the Hayward case $m(v)$ appears both in
the numerator and denominator of $M(m,r)$, while in the Bardeen case $m(v)$
only appears in the numerator. The Misner–Sharp mass
$M_{+}\left(v,R(v)\right)$, which has a clear physical meaning, displays a
divergent behavior in both cases regardless of this difference; in terms of
this quantity, the difference between both cases manifests in its different
rate of divergence.

Figure 4: Left: Numerical evolution of the mass parameter $m_{+}$ and of the
metric functions $f_{\pm}$ of a Bardeen regular black hole. Right: Numerical
evolution of the Misner–Sharp mass. In both plots we picked the parameters
$\beta=1$, $p=12$, $m_{0}=10$, $\ell=1$ with the initial conditions $R(v=1)=5$
and $m_{+}(v=1)=m_{0}+1$.

### Regular black hole: Dymnikova

Next, let us consider a Dymnikova black hole Dymnikova1992

$f_{\pm}=1-\frac{2m_{\pm}\left(1-e^{-r^{3}/\ell^{2}m_{\pm}}\right)}{r}.$ (67)

This is a particularly interesting regular black hole metric to consider as
the derivative of the Misner–Sharp mass with respect to the radial coordinate
is not a polynomial. Therefore, we cannot apply the reasoning of the previous
section to analytically determine the late time behavior of
$M_{+}\left(r,R(v)\right)$. As depicted in fig. (5), the numerical integration
shows an exponential divergence also for this geometry.

Figure 5: Left: Numerical evolution of the mass parameter $m_{+}$ and of the
metric functions $f_{\pm}$ of a Dymnikova regular black hole. Right: Numerical
evolution of the Misner–Sharp mass. In both plots we picked the parameters
$\beta=1$, $p=12$, $m_{0}=10$, $\ell=0.5$ with the initial conditions
$R(v=1)=5$ and $m_{+}(v=1)=m_{0}+1$.

### Regular black hole: Cosmological constant

Finally, let us consider the presence of a non-zero cosmological constant. For
this example we consider again the Hayward metric (62) and we assume Price’s
law (4) up to a given time $v_{0}$, and an exponential decay after $v_{0}$.
The parameter $\alpha$ in eq. (49) is picked so that $m_{-}$ is continuous,
while we have considered two choices of $\omega_{I}$, one for which
$|\kappa_{0}|>\omega_{I}$, and one for which $|\kappa_{0}|<\omega_{I}$. As
shown in fig. 6 the numerical integration agrees with what we expect from the
analysis of sec. IV. For $|\kappa_{0}|>\omega_{I}$ the quasilocal mass grows
exponentially while, in the case $|\kappa_{0}|<\omega_{I}$ after an
exponential phase, the Misner–Sharp mass saturates to a constant value.

Figure 6: Comparison of the evolution of the Misner–Sharp mass for a Hayward
regular black hole, with and without a non-zero cosmological constant. In both
plots we picked the parameters $\beta=1$, $p=12$, $m_{0}=10$, $\ell=0.5$,
$R(v=1)=5$, and $m_{+}(v=1)=m_{0}+1$, leading to a surface gravity
$|\kappa|\approx 1.8$. Left: $\omega_{I}=1<|\kappa_{0}|$. Right:
$\omega_{I}=3>|\kappa_{0}|$

##  VI Conclusions

In this paper, we have analyzed the stability of inner horizons in both
singular and regular black holes by extending Ori’s model to a generic
spherically symmetric spacetime. We have reproduced Ori’s results and
recovered an exponential mass inflation in the case of Reissner–Norström black
holes, but we have also extended our analysis to several families of regular
black hole geometries, showing that mass inflation at the inner horizon is a
robust and very general prediction.

Our findings contradict the conclusions of Bonanno:2020fgp , the discrepancy
being due to a technical flaw in the derivation therein. In particular, while
in reference Bonanno:2020fgp  the Misner–Sharp mass was found to reach finite
values asymptotically in time, we have shown that, once this technical flaw is
accounted for, the Misner–Sharp mass becomes divergent for the all the cases
analyzed there.

A novel result of our analysis is that the late-time instability can become
polynomial instead of exponential. In fact, this happens for well-known
regular black hole metrics such as the Hayward metric. However, this regime is
always preceded by a period of exponential growth, and we quantified that the
onset of the polynomial instability likely occurs when the linear
approximation has already broken down. In all the cases that we have
considered, there is always an initial phase in which mass inflation proceeds
exponentially, and which persists for large enough times so that the
Misner–Sharp mass has already grown considerably. Therefore, the polynomial
regime is unlikely to be physically relevant. However, at the moment we cannot
completely discard the possibility that it might be possible to cook up a
metric for which this transition happens at very early stages, so that the
transition to the polynomial regime takes place before the initial exponential
growth of the Misner–Sharp becomes significant. The slower polynomial rate of
growth would therefore make the inner horizon stable for a longer time.

Similar conclusions apply when a nonvanishing cosmological constant is
allowed. While in principle a cosmological constant could stabilize the inner
horizon, it requires the imaginary part of the least damped quasinormal mode
to be larger than the surface gravity of the inner horizon for every type of
perturbation and for every choice of the mass. This is definitely not the case
for the Reissner–Norström–de Sitter case Cardoso2018b ; Dias2018 ; Hod2018 ;
Cardoso:2018nvb . Furthermore, even if this condition is realized, the
regularization of the instability will only arise after a very long transient
in which the curvature invariants grow exponentially. Hence, even if in these
hypothetical scenarios the Misner–Sharp mass would be finite, it would reach
substantially large values that make imprescindible to consider the
backreaction due to the initial exponential phase of mass inflation.

Finally, let us stress that the instability of the inner horizon does not
necessarily imply the formation of a singularity. We have proved that a small
perturbation has a huge effect on the geometry and leads to the unbounded
growth of the Misner–Sharp mass. While in general relativity this usually
implies the formation of a singularity, in a full theory of quantum gravity,
we might expect the backreaction (perhaps dominated by semiclassical effects
Barcelo:2020mjw ) to drive the geometry to a different class of non-singular
spacetime Carballo-Rubio2019a . However, in order to address this problem, a
geometrical analysis cannot be enough, and we would need to specify the
dynamical field equations of a specific theory. Possible quantum gravity
frameworks to understand this issue are given by asymptotic safety and loop
quantum gravity. In fact, within these theories, there are works predicting
regular black holes with an inner horizon Bonanno2000 ; Platania:2019 ;
Held:2019xde ; Modesto:2008im ; Alesci:2011wn ; Rovelli2014 . These works are
mathematically self-consistent as they only deal with vacuum spacetimes,
showing that quantum gravity effects can regularize the singularity. However,
we have shown in this paper that these analyses are incomplete as the presence
of perturbations cannot be ignored. Taking into account the perturbations in a
consistent way would be a very demanding but essential computation.

###### Acknowledgements.

FDF acknowledges financial support by Japan Society for the Promotion of
Science Grants-in-Aid for Scientific Research No. 17H06359
SL acknowledges funding from the Italian Ministry of Education and Scientific
Research (MIUR) under the grant PRIN MIUR 2017-MB8AEZ.
CP acknowledges: financial support provided under the European Union’s H2020
ERC, Starting Grant agreement no. DarkGRA–757480; support under the MIUR PRIN
and FARE programmes (GW-NEXT, CUP: B84I20000100001); support from the Amaldi
Research Center funded by the MIUR program “Dipartimento di Eccellenza” (CUP:
B81I18001170001).
MV was supported by the Marsden Fund, via a grant administered by the Royal
Society of New Zealand.

## References

  * (1) R. Penrose, “Gravitational collapse and space-time singularities,” Phys. Rev. Lett. 14 (1965) 57–59.
  * (2) S. Hawking, “The occurrence of singularities in cosmology. III. Causality and singularities,” Proc. Roy. Soc. Lond. A A300 (1967) 187–201.
  * (3) S. W. Hawking and R. Penrose, “The Singularities of gravitational collapse and cosmology,” Proc. Roy. Soc. Lond. A314 (1970) 529–548.
  * (4) Virgo, LIGO Scientific Collaboration, B. P. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116 no. 6, (2016) 061102, arXiv:1602.03837 [gr-qc].
  * (5) Virgo, LIGO Scientific Collaboration, B. P. Abbott et al., “GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence,” Phys. Rev. Lett. 116 no. 24, (2016) 241103, arXiv:1606.04855 [gr-qc].
  * (6) VIRGO, LIGO Scientific Collaboration, B. P. Abbott et al., “GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2” Phys. Rev. Lett. 118 no. 22, (2017) 221101, arXiv:1706.01812 [gr-qc].
  * (7) Virgo, LIGO Scientific Collaboration, B. P. Abbott et al., “GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence,” Phys. Rev. Lett. 119 no. 14, (2017) 141101, arXiv:1709.09660 [gr-qc].
  * (8) Virgo, LIGO Scientific Collaboration, B. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 no. 16, (2017) 161101, arXiv:1710.05832 [gr-qc].
  * (9) R. Carballo-Rubio, F. Di Filippo, S. Liberati, and M. Visser, “Phenomenological aspects of black holes beyond general relativity,” Phys. Rev. D 98 no. 12, (2018) 124009, arXiv:1809.08238 [gr-qc].
  * (10) V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Rev. Rel. 22 no. 1, (2019) 4, arXiv:1904.05363 [gr-qc].
  * (11) R. Carballo-Rubio, F. Di Filippo, S. Liberati, and M. Visser, “Opening the Pandora’s box at the core of black holes,” Class. Quant. Grav. 37 no. 14, (2020) 145005, arXiv:1908.03261 [gr-qc].
  * (12) R. Carballo-Rubio, F. Di Filippo, S. Liberati, and M. Visser, “Geodesically complete black holes,” Phys. Rev. D 101 (2020) 084047, arXiv:1911.11200 [gr-qc].
  * (13) E. Poisson and W. Israel, “Inner-horizon instability and mass inflation in black holes,” Phys. Rev. Lett. 63 (1989) 1663–1666.
  * (14) A. Ori, “Inner structure of a charged black hole: An exact mass-inflation solution,” Phys. Rev. Lett. 67 (1991) 789–792.
  * (15) A. J. Hamilton and P. P. Avelino, “The Physics of the relativistic counter-streaming instability that drives mass inflation inside black holes,” Phys. Rept. 495 (2010) 1–32, arXiv:0811.1926 [gr-qc].
  * (16) E. G. Brown, R. B. Mann, and L. Modesto, “Mass Inflation in the Loop Black Hole,” Phys. Rev. D 84 (2011) 104041, arXiv:1104.3126 [gr-qc].
  * (17) R. Carballo-Rubio, F. Di Filippo, S. Liberati, C. Pacilio, and M. Visser, “On the viability of regular black holes,” JHEP 07 (2018) 023, arXiv:1805.02675 [gr-qc].
  * (18) A. Bonanno, A.-P. Khosravi, and F. Saueressig, “Regular black holes have stable cores,” arXiv:2010.04226v1 [gr-qc].
  * (19) S. A. Hayward, “Formation and evaporation of regular black holes,” Phys. Rev. Lett. 96 (2006) 031103, arXiv:gr-qc/0506126 [gr-qc].
  * (20) C. Barrabès and W. Israel, “Thin shells in general relativity and cosmology: The lightlike limit,” Phys. Rev. D 43 (Feb, 1991) 1129–1142. https://link.aps.org/doi/10.1103/PhysRevD.43.1129.
  * (21) R. H. Price, “Nonspherical perturbations of relativistic gravitational collapse. 1. Scalar and gravitational perturbations,” Phys. Rev. D5 (1972) 2419–2438.
  * (22) R. H. Price, “Nonspherical Perturbations of Relativistic Gravitational Collapse. II. Integer-Spin, Zero-Rest-Mass Fields,” Phys. Rev. D5 (1972) 2439–2454.
  * (23) A. S. Barreto and M. Zworski, “Distribution of resonances for spherical black holes,” Mathematical Research Letters 4 no. 1, (1997) 103–121.
  * (24) S. Dyatlov, “Asymptotic distribution of quasi-normal modes for kerr–de sitter black holes,” in Annales Henri Poincaré, vol. 13, pp. 1101–1166, Springer.  2012\.
  * (25) I. Dymnikova, “Vacuum nonsingular black hole,” Gen. Rel. Grav. 24 (1992) 235–242.
  * (26) V. Cardoso, J. a. L. Costa, K. Destounis, P. Hintz, and A. Jansen, “Quasinormal modes and strong cosmic censorship,” Phys. Rev. Lett. 120 (Jan, 2018) 031103. https://link.aps.org/doi/10.1103/PhysRevLett.120.031103.
  * (27) O. J. Dias, F. C. Eperon, H. S. Reall, and J. E. Santos, “Strong cosmic censorship in de Sitter space,” Phys. Rev. D 97 no. 10, (2018) 104060, arXiv:1801.09694 [gr-qc].
  * (28) S. Hod, “Strong cosmic censorship in charged black-hole spacetimes: As strong as ever,” Nucl. Phys. B941 (2019) 636–645, arXiv:1801.07261 [gr-qc].
  * (29) V. Cardoso, J. L. Costa, K. Destounis, P. Hintz, and A. Jansen, “Strong cosmic censorship in charged black-hole spacetimes: still subtle,” Phys. Rev. D98 no. 10, (2018) 104007, arXiv:1808.03631 [gr-qc].
  * (30) C. Barceló, V. Boyanov, R. Carballo-Rubio, and L. J. Garay, “Black hole inner horizon evaporation in semiclassical gravity,” arXiv:2011.07331 [gr-qc].
  * (31) A. Bonanno and M. Reuter, “Renormalization group improved black hole spacetimes,” Phys. Rev. D 62 (Jul, 2000) 043008. https://link.aps.org/doi/10.1103/PhysRevD.62.043008.
  * (32) A. Platania, “Dynamical renormalization of black-hole spacetimes,” Eur. Phys. J. C79 no. 6, (2019) 470, arXiv:1903.10411 [gr-qc].
  * (33) A. Held, R. Gold, and A. Eichhorn, “Asymptotic safety casts its shadow,” JCAP 06 (2019) 029, arXiv:1904.07133 [gr-qc].
  * (34) L. Modesto, “Semiclassical loop quantum black hole,” Int. J. Theor. Phys. 49 (2010) 1649–1683, arXiv:0811.2196 [gr-qc].
  * (35) E. Alesci and L. Modesto, “Particle Creation by Loop Black Holes,” Gen. Rel. Grav. 46 (2014) 1656, arXiv:1101.5792 [gr-qc].
  * (36) C. Rovelli and F. Vidotto, “Planck stars,” Int. J. Mod. Phys. D23 no. 12, (2014) 1442026, arXiv:1401.6562 [gr-qc].