From regular black holes to horizonless objects: quasi-normal modes, instabilities and spectroscopy

We study gravitational and test-field perturbations for the two possible families of spherically symmetric black-hole mimickers that smoothly interpolate between regular black holes and horizonless compact objects accordingly to the value of a regularization parameter. One family can be described by the Bardeen-like metrics, and the other by the Simpson-Visser metric. We compute the spectrum of quasi-normal modes (QNMs) of these spacetimes enlightening a common misunderstanding regarding this computation present in the recent literature. In both families, we observe long-living modes for values of the regularization parameter corresponding to ultracompact, horizonless configurations. Such modes appear to be associated with the presence of a stable photon sphere and are indicative of potential non-linear instabilities. In general, the QNM spectra of both families display deviations from the standard spectrum of GR singular BHs. In order to address the future detectability of such deviations in the gravitational-wave ringdown signal, we perform a preliminary study, finding that third generation ground-based detectors might be sensible to macroscopic values of the regularization parameter.


I. INTRODUCTION
It is widely believed that quantum-gravity effects change the internal structure of black holes (BHs) at some scale ℓ and cure the central singularity.Without specifying the actual theory responsible for these effects, the possible regular spherically symmetric spacetimes can be classified into two families [1] in which the singularity can be either replaced by a global or local minimum radius hypersurface, i.e. a spacelike wormhole throat hidden inside a trapping horizon, as in the Simpson-Visser (SV) spacetime [2], or by an inner horizon shielding a non-singular core [3][4][5][6], as the Bardeen-like regular black hole (RBH) models.
For both families, depending on the value of the regularization parameter ℓ, regular, horizonless configurations can appear [7].The SV metric can actually describe a traversable wormhole connecting two symmetric regions of our universe or two universes, while Bardeen-like RBHs can be continuously deformed into horizonless objects.Remarkably, in addition to the usual unstable photon sphere, these objects may also possess a stable photon sphere, whose existence and position depend on the values of the model parameters.
It has been observed that generally linearized perturbations of ultracompact objects with stable photon spheres, as gravastars or constant density stars, decay extremely slowly [8][9][10].This strongly suggests that the presence of stable photon spheres can lead to a non-linear instability.This link between long living modes and non-linear instabilities has been recently confirmed by a pseudo-spectrum analysis [11].
To investigate whether such instabilities can also be present in the above mentioned horizonless configurations, and to explore the possibility to discriminate the regular from the standard general relativistic BHs through observations, here we study testfield and linear gravitational perturbations in such spacetimes, varying the regularization parameters so to pass smoothly from RBHs to the ultracompact horizonless objects. 1 We stress that regular models are not vacuum solutions of general relativity (GR) and they are often proposed as effective metrics useful for phenomenological analyses.Nonetheless, in several cases, reverse engineering techniques allow for interpreting these regular models as exact solutions of GR coupled to some suitable matter source [13][14][15][16][17].This description, although not unique, makes possible the investigation of gravitational perturbations, getting above the study of test-field perturbation on a fixed background.
The continuity between RBHs and horizonless objects of our models enlighten the fact that stable photon spheres are already present in the borderline cases, that is the extremal RBHs and null-throat wormholes, and its position (in Eddington-Finkelstein coordinates) coincides with that of the extremal horizon [7,18].Actually, this is a general feature of any extremal horizon: it coincides with an extreme point of the potential in the equation for null geodesics [19][20][21].
It has already been argued that extremal BHs should be asymptotically unstable [22,23], and it is natural to ask whether the instability associated to the presence of a photon sphere and that associated to extremal horizons, are actually different ways to describe the same phenomenon.It seems that trapped orbits are indeed present near extremal horizons [24,25].While a linear analysis, as the one reported here, is not sufficient to provide a definitive answer to this issue, it can nonetheless enlighten us on the possible relation between the two aforementioned phenomena.In this sense, we investigated if the long-living modes associated with the lightring instability are present also in the extremal case.However, we found that the damping times in the extremal case are orders of magnitude shorter than in the ultracompact cases, suggesting that the photon sphere instability is not triggered or possibly partially suppressed for extremal RBHs.We argue that this partial suppression could be due to the fact that extremal horizons, being indeed horizons, are not the location of a truly stable orbits but can be considered metastable photon spheres.Indeed, the presence of a horizon, even if extremal, introduce a source of dissipation, i.e. the energy that enters the horizon is completely lost.
As a general feature, the QNM spectra for the regular models that we have considered present deviations from the spectrum of a Schwarzschild BH.Assuming that the effect of rotation in more realistic models does not change the picture significantly, we find that for sufficiently large values of the regularization parameter, and for gravitational-wave events with large signal-to-noise ratio, these deviations could be detectable with next generation detectors [26][27][28].
The paper is organized as follows.In Section II we describe the two families of regular spacetimes, their main features and field sources.In Section III we illustrate the study of perturbations on these spacetimes, we report the obtained field equations and the methods used to find the corresponding QNMs.In Section IV we show and comment our results, while in Section V we discuss how the differences between the obtained spectrum for regular models and the spectrum of singular BHs can be detectable with the next generation of gravitational-waves detectors.The technical details of the derivation of the perturbative equations are given in the Appendix.

II. MODELS
In this work we consider spherically symmetric and static spacetimes with line element Most of the Bardeen-inspired models are described by the line element (1) with ϕ(r) = 0 and some given mass function m(r) which contains a regularization parameter ℓ.Besides Bardeen's original proposal, a non-comprehensive list of widely explored models includes Hayward [4], Dymnikova [5], Fan-Wang [6], and many others.A notable exception is the SV spacetime which is obtained from the Schwarzschild line element with the substitution r → √ r2 + ℓ 2 .After a change of coordinates, even its line element can be written in the form (1).As examples of the two possible families of regular geometries, in this paper we consider the Bardeen and SV models, Bardeen: Simpson-Visser: Remarkably, it can be shown that these two families of solutions pretty much cover all the possible regularized, spherically symmetric, static, BH spacetimes [1]. 2 Furthermore, depending on the value of ℓ, they can also describe ultra-compact, horizonless, objects [7] which will also be considered in this study.

A. Horizons and photon spheres
A first relevant observation for our analysis is that the Bardeen-like and SV models have very different features when interpolating from RBHs to horizonless objects.However, in both cases, there exist two special values of the regularization parameter ℓ, say ℓ ext and ℓ light with ℓ ext < ℓ light , which determine the existence and position of horizons and photon spheres.This is visually illustrated by comparing Fig. 1, where the horizon and photon-sphere structure of the spacetimes is represented according to the value of ℓ.
On the one hand, for ℓ < ℓ ext a Bardeen-like line element describes a RBH with two horizons and one unstable photon sphere; for ℓ = ℓ ext the spacetime becomes an extremal RBH, in which the two horizons and the unstable photon sphere coincide; for ℓ ext < ℓ < ℓ light the horizon disappears, the spacetime describes an ultracompact object with two photon spheres whose distance decreases with increasing values of ℓ, and one of them is stable.Finally, for ℓ > ℓ light the two photon spheres disappear and nor stable or unstable null circular orbits are anymore possible around the object.In particular, for the Bardeen  5 the two photon spheres finally merge leaving a simple compact object.For the SV metric the horizon becomes a wormhole throat for ℓ = 2M over which a stable photon sphere resides.For ℓ = 3M the two photon spheres merge and the wormhole throat becomes an unstable photon sphere.
On the other hand, the SV metric for ℓ < 2M describes a RBH geometry with a single horizon shielding a one-way spacelike throat, surrounded by an unstable photon sphere; for ℓ = 2M the spacetime represents a one-way wormhole with an extremal null throat and two photon spheres, one of which is stable and located at the throat; for 2M < ℓ < 3M the wormhole becomes traversable both ways, the throat at r = 0 is timelike and there are two accessible photon spheres; for ℓ > 3M the spacetime has only one unstable photon sphere located at the throat r = 0.

B. Field sources
As said before, the above introduced static solutions, can be considered the outcome of a transient regularization of the gravitational collapse due to quantum gravity.The implicitly assumption is that such non-classical regime gives way, at late times, to a stationary configuration that should be a solution a some gravitational theory: a low energy, effective field theory limit of quantum gravity, whatever this might be.As our solutions mimic GR ones better and better as one gets away from the objects cores, so we do expect that any such effective field theory of gravity should be encoding deviations from GR in strong gravity regimes.Also, it is well known that such theories can often be recast as GR with non-trivial, and sometimes exotic, matter sources.It is hence reasonable to explore the interpretation of our geometries as solutions of GR and check their associated matter content as this is a crucial step for considering their behaviour under perturbations.
Within GR, the effective stress-energy tensor associated with the line element (1) is given by its Einstein tensor, i.e., T µ ν = G µ ν /8π.Then, for any given RBH model, one might question a posteriori the existence of some matter distribution yielding the same stress-energy tensor.
Notice that the Einstein tensor computed from Eq. (1) has three independent components, meaning that the matter source cannot be uniquely a scalar field (for which T t t = T θ θ ), nor an electromagnetic field (for which T t t = T r r ).Nonetheless, when ϕ(r) = 0, G t t = G r r and Bardeen-like RBHs are often interpreted as solutions of GR coupled to some non-linear electrodynamics with action [13,14] where the electromagnetic Lagrangian is a non-linear function of the electromagnetic field strength being A µ the electromagnetic potential.The Maxwell field is frequently assumed purely magnetic and its magnetic charge coincides with the regularization parameter, which implies that the only non-vanishing component of the Maxwell field is F θφ = ℓ sin θ (alternatively, the only non-vanishing component of the potential is A φ = ℓ cos θ) and F = ℓ 2 /2r 4 .
The modified Maxwell field equation being L F ≡ ∂L/∂F, is trivially satisfied, while the gravitational equations imply that the electromagnetic Lagrangian is given in term of the metric functions of the spacetime as in Eq. ( 1) (with ϕ = 0) where r = r(F).
In particular, for the model considered in this paper On the other hand, when ϕ 0, to model the source it is necessary to introduce other matter fields.In particular, the SV spacetime could be sourced by a combination of non-linear electrodynamics and a self-interacting scalar field [15,17].
where ε = ±1, and the positive (negative) sign corresponds to a canonical (phantom) scalar field with positive (negative) kinetic energy.
Even in this case, we assume the Maxwell field to be purely magnetic with its magnetic charge equal to the regularization parameter, so that the modified Maxwell equation is trivially satisfied.The computation of the gravitational field equations reveals that the scalar field is phantom and satisfies the derivative of the electromagnetic Lagrangian reads which, once integrated, can be substituted in the expression for the scalar potential Finally, the Klein-Gordon equation is a consequence of the Einstein equations.In particular, using the metric functions for the SV spacetime we get Notice that we have used a different convention with respect to Ref. [15]; in particular we have chosen the scalar field to vanish at spatial infinity.

III. STUDY OF PERTURBATIONS
Assuming the gravito-scalar-magnetic interpretation given in Section II B, we can study the full effect of linear perturbations expanding the metric and the matter fields around their background values.According to their parity symmetry, even or odd, the metric and matter perturbations can be decomposed respectively in polar and axial contributions.However since the background metric and the background scalar field are even, while the background magnetic field is odd, axial electromagnetic perturbations and polar scalar perturbations are coupled to polar gravitational perturbations, while polar electromagnetic perturbations are coupled solely to axial gravitational perturbations (being impossible to have axial scalar perturbations).
If this parity coupling is not taken into account, as it commonly happened in several recent investigations [29][30][31][32][33][34][35], one obtains an incompatible system of equations for the perturbation functions, which admits only a trivial solution.This has been quite systematically overlooked in the previously mentioned literature, tacitly assuming that one of the equations can be obtained from the other two.To understand the fine details, the interested reader can follow the full derivation of the perturbative equations in Appendix A and in particular the comment in Footnote 4. On the other hand, other authors have analyzed linear perturbations carefully, but specialized to non-linear electrodynamics without scalar fields or viceversa [36][37][38][39][40][41].Our perturbative analysis extends these results to a generic spacetime described by the line element (1), interpreted as an exact solution of GR coupled to non-linear electrodynamics and scalar fields.

A. Full perturbative analysis
For each parity sector, gravitational, scalar and electromagnetic harmonic perturbations satisfy a system of coupled nonhomogeneous wave equations, which schematically read where r * is the tortoise coordinate defined as dr * /dr ≡ e ϕ / f , for I, J = {A, E} in the sector in which axial gravitational perturbations are coupled to polar electromagnetic perturbations, and I, J = {P, B, S} in the sector in which polar gravitational, axial electromagnetic and polar scalar perturbations are coupled.The variables {A, P, B, E, S} are given combinations of the metric, the electromagnetic potential, and the scalar field perturbation functions and their derivatives.The potentials V I and the coefficients c I,J are given functions of the background metric and fields, and also depend on the harmonic number l associated to the spherical-harmonics expansion.
For the sake of conciseness, we shall only summarise, in Section IV, the outcome of such full perturbative analysis, which is instead explicitly carried on in Appendix A. The latter turns out to be quite involved and dependent on the details of the matter distribution so, as a complementary analysis, we also present in what follows a test-field perturbations analysis which, albeit less accurate, has the merit to avoid assumptions on the matter distribution supporting the geometry.In Section IV we shall see that, reassuringly, the outcomes between the two kinds of analysis turn out to be qualitatively in agreement.

B. Test-field perturbations
Dealing with metrics for which a specific distribution of matter is not specified, test-field perturbations represent a simple but informative proxy.Often, the first step is to consider scalar field perturbations on top of these spacetimes.For spherically symmetric spacetimes, it is possible to extend the analysis to other spin-s fields, to include axial gravitational spin-2 perturbations [42].
The crucial point made in such analyses is that standard matter fields -such as canonical scalars and Maxwellian electric fields -couple with the polar gravitational perturbations, while in the axial sector the source stress-energy tensor is left unperturbed.If true this would imply that the axial gravitational sector could already capture some features of the QNM spectrum.However, this is not the case for a purely magnetic source, and one should be then careful in drawing conclusions.
Within this context, the perturbative equation for scalar, electromagnetic and gravitational axial perturbations for the spacetime described by Eq. ( 1) reads [42,43] where ψ s is related to the spin-s perturbed functions, the tortoise coordinate r * is still defined as dr * /dr ≡ e ϕ / f , and the potential depends on the spin-weight of the perturbation and the metric functions It is interesting to note the formal similarity between equations Eqs. ( 15) and ( 16) modulo the last term on the r.h.s. of Eq. ( 15).

C. Computation of the quasi-normal modes
We now want to solve Eqs. ( 15) and ( 16) for ω, to compute the quasi-normal mode (QNMs) frequencies, i.e. the late-time response of the compact object to an initial perturbation that is localized in space.After providing suitable boundary conditions, that depend on the physical process and on the compact object properties, we use standard and matrix-valued direct integration techniques [44,45] for the test-field case and for the full gravitational case, respectively.
For the RBH cases, the two boundaries from which we integrate are spatial infinity, where we impose the solution to be a purely outgoing wave, and the horizon, where we impose the solution to be a purely ingoing wave.
For the horizonless cases, we still impose the solution to be purely outgoing at spatial infinity, but for the other boundary condition we make a different choice for the two families.The Bardeen-like metrics with ℓ > ℓ ext describe ultracompact stars thus we impose regularity conditions at the origin.Note that in this way we are assuming that the test field perturbations can travel through the entire object with negligible interaction with matter (while in the gravitational perturbations case such interaction is self-contained in the equations of motion).Of course, this assumption may at this point seem unjustified, it is nonetheless the only one that we can do without introducing a specific, and at this stage arbitrary, coupling between the object's matter and our test field (see however our comment about absorption below).The SV metric with ℓ > 2M represents instead a traversable wormhole.Its throat, differently from a horizon, is traversable in both directions.Since the geometry on the two sides of the wormhole throat is symmetric, we assume that the perturbation will inherit the symmetry of the background.This assumption translates into perfect reflection at the throat, which we implement by demanding the perturbation to vanish there, i.e. ψ(ℓ) = 0.
Both the above assumptions can in principle be modified, e.g. for the limit ultra compact object of Bardeen-like geometries we could introduce an absorption coefficient associated to the star matter or in the wormhole case we could assume asymmetric stimulation of the wormhole mouth.We leave these extensions of the present study for future investigations.
The direct integration method we used requires an initial guess for the value of the QNM frequency.While in the RBH case we track the mode continuously starting from its "quasi-Schwarzschild" value obtained for small values of ℓ, in the ultracompact case because of the discontinuity in the boundary conditions (there is no horizon) and the large values of ℓ, we do not have any value as a reference to start from.Thus we explored carefully the (ω I , ω R ) plane in order to find the mode with smaller imaginary part, that is the fundamental one.

IV. RESULTS
In what follows we report the QNM spectra for the considered two families of spherically symmetric regular spacetimes.We focus on the quadrupolar l = 2 fundamental mode, which is the dominant in the gravitational-wave ringdown signal.Note however, that in the ultracompact horizonless cases, these QNMs could become dominant only in the late-time ringdown signal being preceded by a first part of the signal that is very similar to the Schwarzschild one [46].
For test-field perturbations we explore both the RBH and horizonless branches.For the Bardeen metric we vary the regularization parameter from ℓ = 0, that is Schwarzschild, to roughly the maximum value for which the object still possesses a photon sphere.In the SV spacetime a photon sphere is always present at the throat and thus there is no upper bound on the value of the regularization parameter, so we let it span in [0, 3.5M].We show our results in Figs. 2 and 3.
Let us note that some results in the test-field approximation were already present in literature, in a specific branch and for specific values of s.Our results are in agreement with those presented e.g. in Refs.[  3 M (extremal RBH).On the right results for values of ℓ in the horizonless branch (ℓ > ℓ ext ).Note that, in this branch, for values of the regularization parameter near (but not equal to) the extremal one, the imaginary part is extremely small and thus we have very long living modes this is not true for the extremal RBH case, indicated by the vertical line in the left panel.On the left results for values of ℓ in the RBH branch, that is from ℓ = 0 (Schwarzschild) to ℓ = 2M (one-way wormhole with an extremal null throat).On the right results for values of ℓ in the horizonless branch (ℓ > 2M).It is worth noticing the relative flatness of the real part curves which highlights weak deviations from the singular GR solution behaviour recovered for ℓ = 0. Note that, in the horizonless branch (right panel), for values of the regularization parameter near (but not equal to) the extremal one, the imaginary part is extremely small and thus we have very long living modes, this is not true for the extremal RBH case, indicated by the vertical line in the left panel.
For the full perturbative analysis the computation in the horizonless branch presents some technical difficulties and numerical instabilities, therefore we only report the more solid results for the RBH branch, shown in Figs. 4 and 5.However, in advance with the discussion in Section V, we only need the numerical values of gravitational QNMs in the RBH branch to assess the possible detectability of these deviations with the next generation of gravitational-wave detectors.

A. Summary
The results for the two families of regular models presents some differences.For what regards test-field perturbations, the SV spacetime seems to be a better mimicker since, given a certain value of the regularizing parameter ℓ, its spectrum is more similar to the Schwarzschild one (i.e |(ω S V − ω S )/ω S | < |(ω Bard − ω S )/ω S |).We must say however that, for SV, ℓ can span a bigger intervals of values and thus the spectrum can reach higher deviations from Schwarzschild in the imaginary part (some  numerical examples are reported in Table I).Furthermore the corrections to the real part of the frequency in the RBHs branch are negative (except for s = 1) while for the Bardeen spacetime are always positive.The reason for this is clear in Schutz-Will WKB approximation [51] in which ω R ∼ V(r peak ) 1/2 (where r peak is the location of the maximum of the potential).Indeed, compared to the Schwarzschild spacetime, the peak of the potential in the SV spacetime is smaller, whereas in the Bardeen spacetime it is higher.This holds for any spin s of the perturbation except for s = 1.Indeed, for this value of the spin, in the test-field approximation V SV = V Schw and the small positive corrections in the QNMs of the SV spacetime are only due to the different location of the peak in tortoise coordinates.For what regards full gravitational perturbations instead, the real part of the frequency for SV RBHs presents stronger deviations from the Schwarzschild one in the axial sector.
For both families of regular models, in the ultracompact branch we found long living modes associated to the trapping of perturbations near the stable photon sphere.The damping time grows exponentially with the harmonic number and it is longer for values of the regularization parameter near the extremal case, that is for more compact configurations.This is not surprising, a more efficient trapping is expected in these cases since there is more distance between the two photon spheres and a deeper potential well.The aforementioned conclusions stand robust within our framework; however, it is crucial to note that they may be influenced by potential interactions between the test field and matter (e.g., through an absorption coefficient).This consideration is particularly significant as the stable photon sphere seats comfortably within the region where the matter stress-energy tensor is non-negligible [7].
Finally, we also found that the isospectrality between the axial and the polar sector is broken for both families, mainly in the real part of the frequencies, with deviations that, as expected, are greater for greater values of the regularization parameter.

B. A connection between the photon sphere instability and the Aretakis instability?
We found that also in our new class of horizonless ultracompact objects there are long living perturbations modes, associated to the presence of a stable photon sphere.This is usually assumed to be the hint of a non-linear instability. 3Indeed the presence of these long living modes in the frequency decomposition is associated with a total perturbation in time domain that decay slower than 1/t and this leads to the breaking of linear approximation.
An intuitive way to see it is the following.In perturbation theory each order is the source of the next one in the linearized Einstein field equations, then if h (n) is the perturbation at the n-th order, one has thus if q ⩽ 1 then h (2) will be increasing with t, so eventually breaking the perturbative order-expansion.Furthermore, a pseudospectrum analysis [11] showed that these long living modes can be easily perturbed into unstable modes, i.e. modes with a positive imaginary part of the frequency.This means that small fluctuations in the system may trigger growing modes and thus lead to an instability.
with ω S M = 0.37367 − 0.08896i, for s = 2 test-field and linear gravitational perturbations, both in the axial and polar sectors, for selected valued of the regularization parameter.
Results are shown for the Bardeen and SV spacetimes, on the left for the RBH branch and on the right for horizonless configurations.For the Bardeen metric there are no results for ℓ/M = 1.6 and δ = 0.2, with δ ≡ ℓ/ℓ ext − 1, since for those values of compactness the spacetime not only lose the presence of the horizon but even of a photon sphere.For both spacetimes results for axial and polar gravitational perturbations are not reported for horizonless configurations because of the numerical issues present in this branch.Looking at the test field case, it is easy to see the large increment of ∆ I passing from the RBH configurations to the horizonless ones for small δ.
It should be noted that, as anticipated, the stable photon sphere responsible for the above mentioned instability is already present in the limiting case of extremal RBHs.Interestingly, this case is conjectured to be affected by another type of instability, the so-called Aretakis instability [22,23] which appears to be connected to conserved quantities of extremal horizons.
Presently, and differently from the photon sphere instability, the Aretakis instability lacks of a sound physical interpretation.In [24,25] it has been tentatively connected to the presence of null geodesics trapped near the horizon, that is geodesics that orbit arbitrarily many times around the horizon before falling in.If this connection will be confirmed then it will strongly suggest that the Aretakis instability should be interpreted as a special case of the photon sphere one.
However, we have here to notice that the former has been proven to hold also for extremal Kerr BHs [23] albeit for these BHs the photon sphere at the horizon is actually unstable.Of course, also in this case one can observe geodesics that orbit arbitrarily many times around the horizon before falling in, like it happens around any unstable photon sphere, but, usually, this is not associated to any new instability.
Furthermore, from our previous analysis, it is clear that the damping times for extremal RBHs are of the same order of magnitude of that for sub-extremal ones, while ultracompact objects with stable photon sphere presents very long living modes with damping times several order bigger.This seems to suggest that the photon sphere instability is not triggered or partially suppressed for extremal RBHs.Probably this is due to the fact that an extremal horizon, being an horizon, is not a true stable orbit but can be considered a metastable photon sphere (see Fig. 6).The presence of an horizon, even if extremal, introduces a source of dissipation: indeed the energy that enters the horizon is completely lost.
In conclusion, there seems to be no ground for a claim that the Aretakis instability is the limit of the instability associated to stable photon spheres for ultra-compact objects when an extremal trapping horizon forms (even when assuming transparent supporting matter for these solutions).

V. DETECTABILITY
At this point one may wonder if these QNMs can be distinguished from the QNMs of singular GR BHs in the observed gravitational-wave ringdown signals.In other words, will we ever be able to prove that the merging objects that produce a given ringdown signal are not singular GR BHs but RBHs? and how many observations we have to combine to do that?
A general discussion on BH spectroscopy with multiple observations is beyond the scope of this work, a complete review on the topic is [52], see also [53,54].However, as a preliminary answer to the above questions, we can here report results obtained within a particular framework for BH spectroscopy, the Parspec framework [55].One can take this initial analysis as a proof of principle of detectability of these corrections to BH QNMs.

Regular Black Hole
. Difference between a true stable photon sphere present in the spacetime of compact horizonless objects (left panel) and the "stable" photon sphere present at the horizon of extremal BHs and RBHs (right panel).The first one causes a real trapping of modes while in the second case the "trapped" modes pass through the horizon in the BH region.

A. Parspec framework
Parspec is an observable-based parametrization of the ringdown signal of rotating BHs beyond GR, it was developed for BH solutions in modified gravity but can be adapted to our phenomenological models of RBHs.We will give here a brief description of this framework.
Let us assume i = 1, . . ., N independent ringdown detections, for which q QNMs are measured.Each mode of the i-th source is parametrized as where J = 1, 2, ..., q labels the mode; M i and χ i ≪ 1 are the detector-frame mass and spin of the i-th source; D is the order of the spin expansion; ω (k) J and t (k) J are the dimensionless coefficients of the spin expansion for a Kerr BH in GR; γ i are dimensionless coupling constants that can depend on the source i but not on the specific QNM J, for γ i → 0 the GR BH case is recovered; finally δω (k)  J and δt (k) J are the "beyond Kerr" corrections, in general since all the source dependence is parametrized in γ i , these corrections are universal dimensionless numbers.
Since there is no dependence of the corrections on the source, γ i can be set to 1.We assume perturbative corrections, i.e. we assume that γ i δω (k) ≪ 1 and γ i δt (n) ≪ 1.
It should be noted that M i and χ i are extracted assuming GR BHs, i.e. computed from the full inspiral-merger-ringdown waveform within GR.One should extract mass and spin of the BH from the inspiral-merger waveform considering also GR deviations, but this can be very challenging, especially because it requires merger simulations for these RBHs.In this preliminary analysis, we shall assume the shift on the final mass and spin of the source to be negligible.
To construct the probability distribution of the beyond Kerr parameters we use a Bayesian approach: if we indicate with θ the parameters (that in our case are δω (k)  J and δt (k) J ) and with d a given set of ringdown observations, from the Bayes' theorem we have where L( d| θ) is the likelihood function and P 0 ( θ) is the prior on the parameters.Thus from the likelihood we can obtain the full posterior probability distribution P( θ| d) through a Markov chain Monte Carlo (MCMC) method based on the Metropolis-Hastings algorithm.
For each event, the likelihood is chosen to be Gaussian: where the vector ⃗ µ i is where each ⃗ µ i J is a two component vector that depends on the difference between the observed J = 1, ..., q modes and the parametrized templates in Eq. ( 20): and Σ i is the covariance matrix that includes errors and correlations between the frequencies and damping times of the i-th source.
Since the observed QNMs correspond to different values of l and m, i.e. they are "quasi-othonormal", the covariance matrix Σ i = diag(Σ (1)  i , ...., Σ (q) i ) is block-diagonal with each block corresponding to the J-th mode, and thus the likelihood function can be written as a product of Gaussian distributions: Moreover, since we consider N independent detections, the combined likelihood function of the parameters can be further factorized as:

B. Results
We considered only one mode (l = m = 2) and we stick to 0 order in the spin thus we have: The analysis can be generalized to higher order in the spin once computed the gravitational QNMs for this rotating RBHs.We considered the signal coming from the merger remnant of N binary coalescences as observed by a ground-based 3G detector (ET in the so-called ET-D configuration [26]).The 2N masses of the binary components are drawn from a log-flat distribution between [5,95] M ⊙ and the 2N spins from a uniform distribution between [−1, 1].We do not include supermassive BHs in the range of masses since ET will be poorly sensitive to them.We fix the source distance by choosing the signal-tonoise ratio (SNR) of the mode to be 10 2 .We then compute the mass and the spin of the final BH formed after merger using semianalytical relations based on numerical relativity simulations in GR [56].From the final mass of the source we compute the l = 2 frequency and damping time of a RBHs with that mass, we did the analysis for both Bardeen and SV RBHs.We compute the errors on the modes through a Fisher-matrix approach.
Like we can see in Fig. 7, O(100) observations with SNR ≈ 10 2 are enough to exclude with 90% confidence level the hypotesis of GR singular BHs that is δω = δτ = 0.This even for small, but not Planckian, values of the regularization parameter and for both families of regular models.
For the Bardeen metric, the strongest constraints come from the real part of the frequency in the polar sector: when ℓ > 0.13 its deviation from Schwarzschild real frequency allows to exclude the GR hypothesis at 90% confidence level.For the SV metric the strongest constraints come from the real part of the frequency in the axial sector and to exclude the GR hypothesis we need ℓ ⩾ 0.19.Of course these results depend on the number of sources and their SNR.Here we referred to the case of O(100) observations with SNR ∼ 10 2 which should be routinely available with third generation ground-based detectors [26,27].
Note that from the posterior probability distributions is also possible to extract a value for the observed δω and δτ (with associated errors).This should be the value at which the posterior is peaked.Thus if one knows the dependence of these corrections from the regularization parameter (δω(ℓ) and δτ(ℓ)) it is also possible to infer the value of ℓ from the posterior probability distributions.This dependence could be obtained for example fitting the numerical results for the RBHs QNMs computed in Section IV.

VI. CONCLUSIONS
In this paper we have studied test-field and gravitational perturbations on top of the two possible families of spherically symmetric black-hole mimickers, that can be modeled by the Bardeen and SV metric.Both families smoothly interpolate between RBHs and horizonless objects depending on the value of the regularization parameter ℓ that enters the metric.The results for these two families of regular models presents some differences.For what regards test-field perturbations the SV spacetime seems to be a better mimicker since, given a certain value of the regularizing parameter ℓ, its spectrum is more similar to the Schwarzschild one (i.e |(ω S V − ω S )/ω S | < |(ω Bard − ω S )/ω S |).We noticed however, that due to the larger span allowed for ℓ in the SV case, this can produce higher deviations from Schwarzschild in the imaginary part, as explicitly reported in Table I.Furthermore, the corrections to the real part of the frequency in the RBH branch are negative in the SV case while for the Bardeen spacetime are positive.For what regards full gravitational perturbations instead, the real part of the frequency for SV RBHs presents stronger deviations from the Schwarzschild one in the axial sector.We also proved that isospectrality between axial and polar QNMs is broken.
For both families of regular models, in the ultracompact branch, we found long living modes whose damping time grows exponentially with the harmonic index l and is longer for values of the regularization parameter near the extremal case, that is for more compact configurations.These modes are associated with the presence of a stable photon sphere in these spacetime and are usually considered a hint for non-linear instability.
Also the Aretakis instability is expected to affect the extremal RBH case, that is the limiting case between RBHs and horizonless objects.A linear mode analysis is insufficient to confirm it, indeed we find damping times for this case to be of the same order of magnitude of the sub-extremal case.
In general our analysis demonstrates that there are deviations of the QNM spectrum of these spacetimes from that of a Schwarzschild BH due to the non-zero value of the regularization parameter ℓ.So, we analysed the possible detectability of these deviations in the observed gravitational-wave ringdown signals.The detectability of such deviations depends on several aspects such as: the number of observations, their SNR and obviously the size of the regularization parameter.Using the Parspec framework for the analysis we showed that these deviations should be detectable in the near future for Bardeen RBHs with ℓ/M > 0.13 and SV RBHs with ℓ/M > 0.19.Indeed with a hundred of observations with SNR ∼ 100, which should be routinely available with third generation ground-based detectors [26,27], it will be possible to exclude the hypothesis of GR singular BHs with 90% confidence level or to cast constraints on the quantum gravity-induced regularization parameter ℓ.This analysis in only preliminary and we plan to extend it in several ways: using corrections at higher order in the spin, using a more realistic binary population for the sources, and treating also the final mass and spin of the remnant as unknown parameters.We also plan to extend this study on gravitational perturbations to rotating BH mimickers and to better investigate the presumed instability of the extremal case.
In conclusion, in spite of their preliminary nature we do think that the results of the investigations carried on in this work should be taken as a strong encouragement that third generation gravitational-wave experiments have the potential not only to further advance our astrophysical understanding but as well to open a whole new channel into quantum gravity phenomenology.polar µν + h axial µν , and for the matter field.In the Regge-Wheeler gauge, h µν can be written as being S lm b ≡ −Y lm ,φ / sin θ, sin θY lm ,θ with b = {θ, φ}.Likewise, we expand the electromagnetic potential as being . Finally, we decompose the scalar perturbation as In what follows, we drop the symbol l,m and the superscript lm to avoid cluttering the notation.We also assume harmonic time dependence for the perturbation functions, i.e. for any perturbative quantity δF(t, r) we write δF(t, r) = e −iωt δ F(r), but we will omit the tilde.The background metric and scalar field are even under parity transformations, while the background magnetic field is odd.Hence, to linear order, the axial gravitational perturbations couple to the polar electromagnetic perturbations (Sector I), while the polar gravitational perturbations couple to the axial electromagnetic and the polar scalar perturbations (Sector II).
It is relatively easy to check that the equations derived in the next subsections reproduce well known results in the appropriate limits, e.g. the Regge-Wheeler-Zerilli gravitational equations for the Schwarzschild spacetime for L = 0, ϕ(r) = 0 and f (r) = 1 − 2M/r, or those of Ref. [39] for ϕ(r) = 0.

Sector I: axial gravitational-polar electromagnetic
In this sector the axial gravitational perturbations couple with the polar electromagnetic perturbations.Let us begin by considering the modified Maxwell equation for a polar electromagnetic perturbation.At linear order, the field strength squared is unperturbed, F ≈ F (0) .It follows that when computing linear perturbations the Lagrangian and its derivatives are unperturbed as well, e.g.L F ≈ L F (0) , where L F (0) ≡ ∂L/∂F (0) .However, to avoid an excessive cluttering of the equations, in what follows we drop the "(0)" superscript.
The t, r and θ components of the modified Maxwell equation read Equation (A6b) can be solved for u 2 and substituting in Eq. (A6a) gives an equation for u 1 with non-homogeneous terms proportional to h 0 and h 1 .Equation (A6c) is a consequence of the first two equations.The independent components of the perturbed axial gravitational equations are the tθ, rθ and θφ Solving for h 0 in Eq. (A7c) and substituting in Eq. (A7b) we obtain a dynamical equation for h 1 , while Eq.(A7a) is automatically satisfied as a consequence of the previous equations and the modified Maxwell equations. 4o make the four independent equations more readable, it is helpful to introduce a new variable u, which corresponds to the perturbation of the tr component of the Maxwell tensor, instead of the perturbations of the potential u 1 and u 2 .
Taking the first and second derivative of Eq. (A8), solving for u ′ 1 and u ′′ 1 , and substituting into Eqs.(A6a) and (A6b) we get Solving Eqs.(A7c), (A9a) and (A9b) for h 0 , u 1 and u 2 and substituing in Eq. (A7b), we obtain the gravitational dynamical equation for h 1 while the electromagnetic dynamical equation for u is obtained from Eq. (A8) with the substitutions above and using Eq.(A10) Finally, Eqs.(A10) and (A11) can be written as wave equations by performing the substitutions h 1 = re ϕ A/ f and u = e −ϕ E/r 2 √ L F , and by introducing a tortoise-like coordinate dr * /dr = e ϕ / f , to get the coupled system where ) 2. Sector II: polar gravitational-axial electromagnetic-polar scalar In this sector the polar gravitational perturbations couple with the axial electromagnetic and polar scalar perturbations.Let us begin with the Klein-Gordon equation The field strength squared for an axial perturbation is In this case, when computing linear perturbations to Eq. ( 5) we also expand L F around F (0) , e.g.L F ≈ L F (0) + L F (0) F (0) δF, and similarly for higher derivatives.
With the further gauge choice u 3 = 0, the θ component of the modified Maxwell equations is the only non-vanishing, and reads Lastly, let us consider a polar gravitational perturbation.The θφ component of the perturbed gravitational equation requires Using the background equations, the other six independent gravitational equations, namely the tt, tr, tφ, rr, rφ and θθ components, are The off-diagonal equations are first-order differential equations in the metric perturbations and can be solved for H ′ 0 , H ′ 1 and K ′ , hence the rr component (A19d) gives an algebraic relation among the metric perturbation functions, which can be used to eliminate H 0 from the other equations.
Using these relations as well as the background equations, Eqs.(A19a) and (A19f) are automatically satisfied.Let H 1 = ωR, then the relevant equations are the tr and tφ components, which can be written as a system of two non-homogeneous coupled differential equations where Now, the procedure to obtain the equation that governs polar gravitational perturbations follows Zerilli's original derivation.The task now is to find a new couple of functions R and Ĥ to transform Eq. (A20) into where the new radial variable r is given by dr/dr = 1/n(r).To find such transformation we write Let us introduce the matricial notation then Eq. (A20) can be written as dψ/dr = Aψ + J, Eq. (A23) as ψ = F ψ, which combined with our request Eq. (A22) give the system together with the new source terms Ĵ = nF −1 J. Equation (A25) represents four equations that relate g 1 , g 2 , k 1 , k 2 , n and V in terms of α 1,2 and β 1,2 .By equating the coefficients of ω 0 and ω 2 we get eight equations, supplemented by the condition det F 0, for six unknown functions.Yet, the system is consistent and admits a solution where we have introduced λ = (l − 1)(l + 2)/2.Equation (A26a) means that the new variable r is nothing but the tortoise-like coordinate r * .The new source terms read The above system and the source terms simplify when the integral in Eq. ( A26c) is given in a closed form, and this depends strongly on the explicit form of the background metric functions.Remarkably, for ϕ = ϕ 0 + 1 2 log 1 − ℓ 2 /r 2 , with ϕ 0 being an arbitrary constant, as for the SV spacetime, and for any choice for f , we find g 2 = e ϕ .We assume it in what follows.
Finally, combining Eq. (A22) we get a master equation for the polar gravitational perturbations coupled with the axial electromagnetic and polar scalar perturbations The very last step is to use the solutions for the gravitational equations to rewrite the Klein-Gordon and modified Maxwell equations; they read Equations (A29) to (A31) can be written as wave equations by introducing new variables where with c 1 c 4 − c 3 c 6 0, so that for I, J = {P, B, S}.The potentials V I and the coefficients c I,J can be given in closed form: where 4ℓ(l − 1)l(l + 1)(l + 2) ∆ = (l − 1)(l + 2)r 2 + 4ℓ 2 L F .(A54) Once the solutions for H 0 , u 4 and δΦ are known, the other metric function perturbation is given by

Figure 1 .
Figure 1.Radii of the photon spheres (solid red lines for the inner stable one and solid purple line the outer unstable one) and horizons (dashed black line) for the Bardeen (left panel) and SV (right panel) metric.For the Bardeen metric the two horizons merge for ℓ = 4M/3 √ 3 giving way to a stable photon sphere inside the usual unstable one.For ℓ = 48M 25

Figure 3 .
Figure3.Quadrupolar l = 2 fundamental QNMs of the SV metric for test-field perturbations, s = 0 (blue), s = 1 (light purple) and s = 2 (red).On the left results for values of ℓ in the RBH branch, that is from ℓ = 0 (Schwarzschild) to ℓ = 2M (one-way wormhole with an extremal null throat).On the right results for values of ℓ in the horizonless branch (ℓ > 2M).It is worth noticing the relative flatness of the real part curves which highlights weak deviations from the singular GR solution behaviour recovered for ℓ = 0. Note that, in the horizonless branch (right panel), for values of the regularization parameter near (but not equal to) the extremal one, the imaginary part is extremely small and thus we have very long living modes, this is not true for the extremal RBH case, indicated by the vertical line in the left panel.

Figure 5 .
Figure 5. Axial (blue) and polar (light purple) l = 2 gravitational QNMs for the SV metric for values of ℓ in the RBH branch that is from ℓ = 0 (Schwarzschild) to ℓ = 2M (one-way wormhole with an extremal null throat).

Figure 7 .
Figure 7. Probability distribution functions for corrections to the Schwarzschild l = m = 2 mode when we inject as observations axial/polar QNMs of the Bardeen metric with ℓ/M = 0.3 (left upper/bottom panel) and axial/polar QNMs of the SV metric with ℓ/M = 0.3 (right upper/bottom panel).Different colors represents results obtained with different numbers N of observed sources.

Table I .
Relative deviations from the quadrupolar fundamental Schwarzschild frequency ∆ R/I