sectionSection #1 equationEq. (#1) figureFig. #1 subfigureFig. 0#1 tableTable
#1

# Accreting Schwarzschild-like compact object: Plasma-photon interaction and
stability

Avijit Chowdhury<EMAIL_ADDRESS>Department of Physics,
Indian Institute of Technology Guwahati, Guwahati 781039, India  School of
Physical Sciences, Indian Association for the Cultivation of Science, Kolkata
700032, India  Shauvik Biswas<EMAIL_ADDRESS>School of Physical
Sciences, Indian Association for the Cultivation of Science, Kolkata 700032,
India  Sumanta Chakraborty<EMAIL_ADDRESS>School of Physical Sciences,
Indian Association for the Cultivation of Science, Kolkata 700032, India

###### Abstract

Accretion is a common phenomenon associated with any astrophysical compact
object, which is best described by plasma, a state of matter composed of
electrons and heavy ions. In this paper, we analyze the linear dynamics of
electromagnetic (EM) fields propagating through the accreting plasma around a
static and spherically symmetric horizon-less, exotic compact object (ECO).
The general equations governing the propagation of EM waves in such a
background exhibits quasi-bound states, with a characteristic oscillation
around the BH values, for both the axial and the polar modes, as well as for
homogeneous and inhomogeneous plasma distributions. The amplitude of these
oscillations depend on the non-zero reflectivity of the surface of the compact
object, while the oscillation length depends on its compactness. This results
into slower decay of the quasi-bound states with time for certain parameter
space of the plasma frequency, compared to BHs, making these ECOs more prone
to instabilities.

##  I Introduction

Black holes (BHs) are one of the most remarkable constructs of General
Relativity (GR), unmatched in their simplicity and elegance. The
characteristic feature of a BH is the existence of an event horizon, acting as
a one-way membrane, which causally shields the exterior region from the
interior, and hides the singularity within. The observation of gravitational
waves (GWs), from the merger of two compact objects, by the LIGO-Virgo-Kagra
(LVK) collaboration [1, 2, 3, 4, 5] and the shadow measurement of compact
objects at the centre of M87 and our own galaxy, by the Event Horizon
Telescope (EHT) collaboration [6, 7], are consistent with these compact
objects being BHs. Even though, both of these observations have proved beyond
doubt that whatever be the nature of the compact object, it must be more
compact than the photon sphere, but the existence of an event horizon has not
been directly verified [8, 9, 10]. Therefore, the question remains whether BHs
are the only possible objects which are more compact than a photon
sphere111There can actually be stars, made out of normal matter, which are
more compact than the photon sphere, known as the Buchdahl stars [11, 12, 13,
14, 15]. However, in the presence of accretion, these stars would describe a
significant outflow, which are absent in the EHT observation [16, 17]
(however, see [18]). This suggests that the objects must be more compact than
the Buchdahl stars and hence must be made out of exotic matter fields., or can
there exist other exotic compact objects (ECOs) [19, 20, 21, 22, 23, 24]. The
exotic-ness of these ECOs stem from the absence of the event horizon in these
spacetimes. The ECOs are rather characterised by a surface, close to the event
horizon, that is partially/fully reflecting. The reflectivity of the ECO
surface, along with its position relative to the BH horizon, encodes
information about the object’s interior. The ECOs behave exactly as a BH for
scales larger than the size of the photon sphere, but the physics at the
horizon is completely different. As a consequence there are important effects
in the GW sector arising from these ECOs — (a) existence of echoes in the
ringdown waveform of GWs [25, 26, 27, 28, 29, 30, 31, 27, 32, 33, 34, 35, 36],
(b) non-zero tidal Love numbers [37, 38, 39, 40], (c) modifications in the GW
phasing due to tidal heating [41, 42, 43], among many others. Interestingly,
most of the above analysis have been performed assuming the ECOs to be
isolated, while in reality these objects are in astrophysical environments, be
it a dark matter halo, or accreting plasma. The consequences of ECOs evolving
in a dark matter environment have been analysed recently [44], while in the
present work we focus on the effect of accretion on static and spherically
symmetric ECOs.

Just as galactic BHs are supposed to be immersed in a dark matter halo, any
astrophysical BHs are also surrounded by accreting matter. Since the accreting
matter is essentially made out of plasma, it is expected that electromagnetic
(EM) emission from the accreting matter will interact with the plasma, thereby
affecting its propagation at large distances. Hence, the most non-trivial
behaviour is expected from the EM fields moving in the spacetime. This has
already been observed in [45], regarding the propagation of EM wave around an
accreting Schwarzschild BH, leading to quasi-bound states. This happens
because, the photons propagating through the plasma effectively gains a mass,
dependent on the plasma frequency $\omega_{\rm pl}$ of the background
accreting medium [46],

$\displaystyle\omega_{\text{pl}}=\sqrt{\frac{n_{\text{e}}e^{2}}{m_{\text{e}}}}~{}.$
(1)

Here $n_{\text{e}}$ is the number density of the electrons, having charge $e$
and mass $m_{\text{e}}$. In case of inhomogeneous distribution of plasma, the
electron number density and the plamsa frequency become dependent on the
spacetime coordinates. In the case of static and spherically symmetric
background, they become functions of the radial coordinate $r$. The existence
of effective mass for photons drastically changes the implications regarding
the stability of the compact object [47, 48, 49, 50, 51, 52, 53, 54, 55, 56].
Since it is well-known that non-zero effective mass of the photons may lead to
superradiant instabilities for astrophysical BHs [57], and more so for the
ECOs, owing to their reflective nature [58]. Therefore, it is important to ask
what happens to the photons propagating around an accreting ECO, in
particular, whether they lead to any instabilities. As we will demonstrate,
the dynamics of the axial part of an EM field is exactly like a Proca field,
while the dynamics of the polar sector is very different, which is consistent
with the earlier findings in the literature [59, 45]. Another striking
difference between the Proca field and EM waves moving in a plasma is that,
the motion of EM wave through plasma is still characterised by only two
polarization modes (both transverse), whereas a Proca field is always
characterised by three polarization modes (one longitudinal and two
transverse). Therefore, one must revisit the plasma-photon interaction around
ECOs, which have been well understood for BHs [60, 61, 45, 62, 63]. This is
also motivated from the quasi-normal mode (QNM) spectrum for ECOs, where the
decay time of the QNMs associated with ECOs exceeded the BH values by several
orders of magnitude, making the ECOs prone to instabilities. However, QNMs are
eigenvalues of scattering states associated with the compact object, while the
EM fields around accreting ECOs lead to quasi-bound states. Therefore,
implications for the spectrum of these quasi-bound states associated with an
ECO need to be carefully assessed to understand the stability of these
objects. As a first step in this direction, we have determined the relevant
equations and the frequency spectrum of the quasi-bound states associated with
EM waves moving in the background of accreting Schwarzschild-like ECO in this
work. We hope to extend our analysis for a rotating ECO and study its
stability.

For this purpose, we start by generalizing the dynamical equations governing
the propagation of the EM field through unmagnetised, cold and collision-less
plasma on an arbitrary static and spherically symmetric background.
Subsequently, we model the spacetime exterior to the ECO surface by the
Schwarzschild metric and investigate how the quasi-bound frequencies respond
to the nonzero reflectivity of the ECO surface and its location relative to
the horizon. We have considered both constant and frequency-dependent
reflectivity, as well as, both the axial and polar part of the EM field. In
addition, the case of Bondi accretion, involving inhomogeneous plasma
distribution has also been discussed. We have also corroborated our numerical
findings with analytic calculations.

The structure of the paper is as follows: we start by defining the general
equations governing the motion of an EM wave in a cold and collision-less
plasma in II. These equations have been separated into axial and polar
sectors, and spelled out for a generic static spherically symmetric spacetime,
along with the plasma profiles in III. Subsequently, in IV we discuss the
numerical setup, methodology and the results obtained. Finally, we conclude
with a summary of our results in V, highlighting their importance and
prospects. Some detailed analytical computations have been delegated to A, B
and C, respectively.

_Notations and conventions:_ Throughout the paper, we will use units in which
$G=1=c$ and the signature of the flat spacetime Minkowski metric will be taken
to be mostly positive, i.e., $\eta_{\mu\nu}=\textrm{diag.}(-1,1,1,1)$.
Moreover, all the Greek indices $\mu,\nu,\alpha,\cdots$ denote spacetime
indices, and we will exclusively work with four spacetime dimensions.

##  II Plasma-photon interaction in curved spacetime: Basic equations

In this section we will present the basic equations governing the interaction
between plasma and photon in a curved background. We will keep this section
brief and the equations general, since they have appeared in various other
contexts before [64, 65, 66, 67], and shall apply them subsequently for a
general static and spherically symmetric spacetime. As emphasized before, the
interaction between plasma and photon is typical of any accreting system,
where the accreting material is best described by a plasma, comprising of
electrons and ions, which emits EM waves, interacting with the plasma itself.
If within the plasma, the thermal velocity of the electrons is small compared
to the propagation speed (phase velocity) of the EM waves, one can classify it
as cold plasma. The plasma will be referred to as hot in the reversed
situation. Therefore in the cold plasma the following inequality holds true

$\displaystyle v^{\text{thermal}}_{\rm e}=\sqrt{\frac{2k_{\rm B}T_{\rm
e}}{m_{\rm e}}}\ll\frac{\omega}{k}~{},$ (2)

where $\omega$ and $k$ are the frequency and wave number of the EM field,
respectively, and $T_{\rm e}$ is the temperature of the electrons within the
plasma. In what follows we will concentrate on the case of cold plasma.

The propagation of EM waves through an unmagnetized cold plasma consisting of
electrons and ions (since the ions are much heavier than electrons they barely
move compared to the EM waves) in curved spacetime is described by three four
vectors — (a) the vector potential $A_{\mu}$, describing the EM waves, (b) the
four velocity $u^{\mu}$ of the electrons, (c) the current density $J^{\mu}$ of
the ions, and a scalar — the number density $n_{\rm e}$ of electrons. The
dynamics of the system is thus described by the Maxwell equation together with
the momentum and particle conservation equation, [68],

$\displaystyle\nabla_{\nu}F^{\mu\nu}$ $\displaystyle=$ $\displaystyle
en_{\text{e}}u^{\mu}+J^{\mu}~{},$ (3) $\displaystyle
u^{\mu}\nabla_{\mu}u^{\nu}$ $\displaystyle=$ $\displaystyle\frac{e}{m_{\rm
e}}F^{\nu}_{\ \mu}u^{\mu}~{},$ (4)
$\displaystyle\nabla_{\mu}\big{(}n_{\text{e}}u^{\mu}\big{)}$ $\displaystyle=$
$\displaystyle 0~{},$ (5)

where $F_{\mu\nu}=\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}$ is the EM field
strength tensor and $u^{\mu}u_{\nu}=-1$ with $m_{\text{e}}$ being the electron
mass.

The background spacetime metric will get perturbed by the presence of the
accreting plasma, but only at the second order, since the energy-momentum
tensor is quadratic in the matter degrees of freedom. For the EM field, this
can be easily understood through the following analysis: the vector potential
behaves as $A_{\mu}=\epsilon A_{\mu}^{\rm acc}$, where $\epsilon$ is a small
parameter and $A_{\mu}^{\rm acc}$ is the vector potential generated by the
equilibrium distribution of the accreting plasma. Therefore, the energy
momentum tensor of the accreting EM field behaves as
$T_{\mu\nu}\sim\epsilon^{2}T_{\mu\nu}^{\rm acc}$, and hence the metric gets
perturbed as $g_{\mu\nu}=g_{\mu\nu}^{\rm bg}+\mathcal{O}(\epsilon^{2})$222An
explicit example would be the case of the Reissner-Nordström BH, where the
vector potential scales as the charge $Q$, while the metric is affected only
at $\mathcal{O}(Q^{2})$., where $g_{\mu\nu}^{\rm bg}$ is the background static
and spherically symmetric metric.

To assess the stability of the system, we introduce small perturbations over
and above the accreting plasma in equilibrium, in particular, we let
$A_{\mu}=\epsilon A_{\mu}^{\rm acc}+\epsilon^{2}A_{\mu}^{\rm pert}$. Similar
perturbing effects are present to the plasma quantities, namely the electron
density and the electron velocity, respectively. Therefore, the dynamical
equations for the perturbed EM field $A_{\mu}^{\rm pert}$ can be obtained in
terms of $g_{\mu\nu}^{\rm bg}$, $A_{\mu}^{\rm acc}$ and unperturbed quantities
associated with plasma in equilibrium. Thus any effect of backreaction on the
spacetime metric due to accretion does not affect the dynamics of the photon
field $A_{\mu}^{\rm pert}$. As emphasized before, the ions in the plasma are
assumed to be heavier than the electrons, and hence any perturbation of the
ion quantities is suppressed compared to that of the electron quantities. For
this purpose, we project the background metric on the plane orthogonal to the
four-velocity of electrons, by introducing $h_{\mu\nu}=g^{\rm
bg}_{\mu\nu}+u_{\mu}u_{\nu}$. Thus, any generic tensor can be decomposed into
components along and orthogonal to $u_{\mu}$, e.g. [69],

$\nabla_{\mu}u_{\nu}=\omega_{\mu\nu}+\theta_{\mu\nu}-u_{\mu}u^{\alpha}\nabla_{\alpha}u_{\nu}$
(6)

where

$\displaystyle\omega_{\mu\nu}=\frac{1}{2}\left(v_{\mu\nu}-v_{\nu\mu}\right)~{},$
(7)
$\displaystyle\theta_{\mu\nu}=\frac{1}{2}\left(v_{\mu\nu}+v_{\nu\mu}\right)~{}.$
(8)

Here, $v^{\mu\nu}\equiv h^{\mu\alpha}h^{\nu\beta}\nabla_{\alpha}u_{\beta}$ is
the projected four-velocity gradient, with $\omega_{\mu\nu}$ and
$\theta_{\mu\nu}$ being the symmetric and the anti-symmetric parts of
$v_{\mu\nu}$ (referred to as the vorticity and the deformation tensors,
respectively). The electric part of the Maxwell’s tensor is given by
$E^{\mu}\equiv F^{\mu}_{\ \nu}u^{\nu}$, while in the present context involving
un-magnetized cold plasma, the magnetic field identically vanishes. Given 4,
it immediately implies that the acceleration of the electrons are generated by
the electric field, as it should. With the above definitions, and following
Ref. [68], we obtain the dynamical equations for the perturbed EM field
$A^{\rm pert}_{\mu}$ in an unmagnetized cold plasma, within the Landau gauge
($u^{\mu}A_{\mu}^{\rm pert}=0$) as,

$\displaystyle\Big{[}h^{\alpha\beta}u^{\nu}\nabla_{\nu}\big{(}\nabla^{\mu}\nabla_{\beta}-\delta^{\mu}_{\beta}\square\big{)}+\big{(}\omega^{\alpha\beta}+\theta^{\alpha\beta}+\theta
h^{\alpha\beta}$
$\displaystyle\qquad+\frac{e}{m_{\text{e}}}E^{\alpha}u^{\beta}\big{)}\big{(}\nabla^{\mu}\nabla_{\beta}-\delta^{\nu}_{\beta}\square\big{)}+\omega_{\text{pl}}^{2}h^{\alpha\mu}u^{\gamma}\nabla_{\gamma}$
$\displaystyle\qquad+\omega_{\text{pl}}^{2}\big{(}\theta^{\alpha\mu}-\omega^{\alpha\mu}\big{)}\Big{]}A_{\mu}^{\rm
pert}=0~{}.$ (9)

Here, all the raising and lowering of indices have been done by the background
metric, e.g., $\theta=g^{\mu\nu}_{\rm bg}\theta_{\mu\nu}$. Note that the
Larmour tensor does not appear in the above equation, since we are assuming
the plasma to have zero magnetic field. Moreover, II is in contrast to the
effective Proca-like equation, traditionally used to describe the dynamics of
photon in plasma, as already pointed out in [64].

##  III EM waves in plasma: Dynamical equations in a generic static and
spherically symmetric spacetime

Following the dynamical equations satisfied by an EM field moving in a plasma
surrounding an arbitrary curved background, in this section, we apply the same
to a generic static and spherically symmetric spacetime. We take the following
line element for the static and spherically symmetric background,

$ds^{2}=-f(r)dt^{2}+\frac{1}{g(r)}dr^{3}+r^{2}d\Omega^{2}~{},$ (10)

where $d\Omega^{2}$ denotes the line element on a unit two-sphere. The four
velocity of an electron within the plasma in equilibrium is given by
$u^{\alpha}=(u^{0},0,0,0)$, with $u^{0}=(1/\sqrt{f})$. In this case, the
vorticity and the deformation tensor vanish identically. The spherical
symmetry of the background spacetime allows us to decompose the vector
potential $A_{\mu}^{\rm pert}$ into two independent sectors — (a) the axial
sector, which transforms as $(-1)^{l+1}$ and (b) the polar sector, which
transforms as $(-1)^{l}$, under parity. In addition, every component of the
vector potential can be expressed as the multiplication of two parts — one
depending on $(t,r)$ coordinates and another on $(\theta,\phi)$ coordinates,
respectively. Therefore, after decomposing the vector potential in the
spherical harmonic basis, we obtain,

$\displaystyle A_{\mu}^{\rm pert}$ $\displaystyle=\sum_{lm}A_{\mu,~{}lm}^{\rm
polar}+A_{\mu,~{}lm}^{\rm axial}~{};$ (11) $\displaystyle A_{\mu,~{}lm}^{\rm
polar}$
$\displaystyle=u^{lm}_{1}(t,r)Y_{lm}\delta_{\mu}^{t}+u^{lm}_{2}(t,r)Y_{lm}\delta_{\mu}^{r}$
$\displaystyle\qquad+\frac{u^{lm}_{3}(t,r)}{l(l+1)}\left(0,0,\partial_{\theta}Y_{lm},\partial_{\phi}Y_{lm}\right)~{},$
(12) $\displaystyle A_{\mu,~{}lm}^{\rm axial}$
$\displaystyle=\frac{u^{lm}_{4}(t,r)}{l(l+1)}\left(0,0,\frac{1}{\sin\theta}\partial_{\phi}Y_{lm},-\sin\theta\partial_{\theta}Y_{lm}\right)~{}.$
(13)

The above angular functions ensure that under parity transformation the
functions $u^{lm}_{i}$ belong to the polar sector for $i=1,2,3$, and to the
axial sector for $i=4$. Thus, the EM field has been neatly separated into
axial and polar sectors.

###  III.1 Axial Sector

In the axial sector, we have a single equation for $u^{lm}_{4}$ to solve. This
equation can be obtained by substituting $A_{\mu}^{\rm axial}$ from 13 to II,
resulting into the following dynamical equation for $u^{lm}_{4}(t,r)$,

$\displaystyle
f(r)\Big{[}-2l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Big{]}\frac{\partial
u_{4}^{lm}(t,r)}{\partial t}$
$\displaystyle\quad+r^{2}\bigg{[}g^{\prime}\frac{\partial^{2}u_{4}^{lm}(t,r)}{\partial
t\partial r}+2g\frac{\partial^{3}u_{4}^{lm}(t,r)}{\partial t\partial
r^{2}}\bigg{]}$
$\displaystyle\quad+r^{2}\left[gf^{\prime}\frac{\partial^{2}u_{4}^{lm}(t,r)}{\partial
t\partial r}-2\frac{\partial^{3}u_{4}^{lm}(t,r)}{\partial t^{3}}\right]=0~{}.$
(14)

Here, ‘prime’ denotes derivative with respect to $r$. Note that, the plasma
frequency $\omega_{\rm pl}$ depends on the radial coordinate $r$ through the
number density of electrons. Defining,
$\partial_{t}u^{lm}_{4}\equiv\Psi^{lm}_{\text{ax}}(t,r_{*})$, we can rewrite
the III.1 as,

$\frac{d^{2}\Psi^{lm}_{\text{ax}}}{dr_{*}^{2}}-\frac{d^{2}\Psi^{lm}_{\text{ax}}}{dt^{2}}-\frac{f(r)}{r^{2}}\Big{[}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Big{]}\Psi^{lm}_{\text{ax}}=0~{},$
(15)

where $r_{*}$ is the tortoise coordinate, defined as

$\frac{dr_{*}}{dr}=\frac{1}{\sqrt{f(r)g(r)}}~{}.$ (16)

The above equation can be expressed as a Schrödinger-like equation by
transforming to the frequency domain. This is achieved by the expressing the
axial master variable as, $\Psi^{lm}_{\rm ax}(t,r)=\int d\omega e^{-i\omega
t}\Psi^{lm}_{\rm ax}(\omega,r)$, such that we obtain the following equation,

$\frac{d^{2}\Psi^{lm}_{\text{ax}}}{dr_{*}^{2}}+\left(\omega^{2}-V^{lm}_{\text{ax}}\right)\Psi_{\text{ax}}=0~{}.$
(17)

where, the effective potential takes the following form,

$V^{lm}_{\text{ax}}=\frac{f(r)}{r^{2}}\Big{[}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Big{]}~{}.$
(18)

Intriguingly, the axial potential does not depend on $m$, and in the limit
$\omega_{\rm pl}\to 0$, we get back the effective potential experienced by a
photon in a static and spherically symmetric background. Even though the
differential equation satisfied by $A_{\mu}^{\rm pert}$ is complicated, at
least for the axial sector, it resembles the effective potential of a Proca
field with mass $m\sim\omega_{\rm pl}$ [47, 45].

###  III.2 Polar Sector

In the polar sector, we have three unknown functions $u^{lm}_{1}$,
$u^{lm}_{2}$ and $u^{lm}_{3}$, whose evolution equations are to be determined
from the general differential equation for $A_{\mu}^{\rm pert}$. Since we have
chosen the four velocity of electrons in equilibrium to have only the zeroth
component, the Landau gauge condition ensures that,

$u^{lm}_{1}=0~{}.$ (19)

Therefore, instead of three we now have two EM degrees of freedom. In order to
determine the differential equations associated with the evolution of
$u^{lm}_{2}$ and $u^{lm}_{3}$, we substitute 12 in II, which along with the
condition presented in 19, yields the following two coupled equations,

$\displaystyle\Big{[}l(l+1)$
$\displaystyle+r^{2}\omega_{\text{pl}}^{2}\Big{]}\frac{\partial
u_{2}^{lm}}{\partial t}-\frac{\partial^{2}u_{3}^{lm}}{\partial t\partial
r}+\frac{r^{2}}{f(r)}\frac{\partial^{3}u_{2}^{lm}}{\partial t^{3}}=0~{},$ (20)
$\displaystyle f(r)\Bigg{[}2\omega_{\text{pl}}^{2}$
$\displaystyle\frac{\partial u_{3}^{lm}}{\partial
t}+g^{\prime}\Big{\\{}l(l+1)\frac{\partial u_{2}^{lm}}{\partial
t}-\frac{\partial^{2}u_{3}^{lm}}{\partial t\partial r}\Big{\\}}\Bigg{]}$
$\displaystyle+g\Bigg{[}f^{\prime}\Big{\\{}$ $\displaystyle
l(l+1)\frac{\partial u_{2}^{lm}}{\partial
t}-\frac{\partial^{2}u_{3}^{lm}}{\partial t\partial
r}\Big{\\}}+2f\Big{\\{}-\frac{\partial^{3}u_{3}^{lm}}{\partial t\partial
r^{2}}$ $\displaystyle+l(l+1)\frac{\partial^{2}u_{2}^{lm}}{\partial t\partial
r}\Big{\\}}\Bigg{]}+2\frac{\partial^{3}u_{3}^{lm}}{\partial t^{3}}=0~{}.$ (21)

Note that the above equations are third order differential equations, which is
consistent with the order of the differential equation for $A_{\mu}^{\rm
pert}$, as presented in II. To proceed further, and to obtain a second order
differential equation, we switch to the frequency domain, i.e., we express the
two variables $u_{2}^{lm}$ and $u_{3}^{lm}$ as, $u_{2}^{lm}(t,r)=\int d\omega
U_{2}^{lm}(\omega,r)e^{-i\omega t}$ and $u_{3}^{lm}(t,r)=\int d\omega
U_{3}^{lm}(\omega,r)e^{-i\omega t}$. Therefore, 20 can be used to express
$U_{2}^{lm}$ in terms of the radial derivative of $U_{3}^{lm}$ as,

$U^{lm}_{2}=\frac{f{U^{lm}_{3}}^{\prime}}{-r^{2}\omega^{2}+f\\{l(l+1)+r^{2}\omega_{\text{pl}}^{2}\\}}~{}.$
(22)

Since in the limit, $f(r)\to 0$, the EM degree of freedom $U_{2}^{lm}$
identically vanishes, we should find out the master differential equation in
terms of $U_{3}^{lm}$. Substituting the above expression for $U^{lm}_{2}$ in
the frequency domain representation of 21, the following decoupled equation
for $U_{3}^{lm}$ can be obtained,

$\displaystyle{U_{3}^{lm}}^{\prime\prime}$
$\displaystyle+\frac{1}{2rfg\Bigl{(}-\omega^{2}+f\omega_{\text{pl}}^{2}\Bigr{)}\biggl{[}-r^{2}\omega^{2}+f\Bigl{\\{}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}\biggr{]}}\Bigg{[}r^{3}\omega^{4}gf^{\prime}+f\biggl{(}r\omega^{2}g\Bigl{\\{}l(l+1)-2r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}f^{\prime}+r^{3}\omega^{4}g^{\prime}\biggr{)}$
$\displaystyle+f^{2}\biggl{(}g\Bigl{\\{}-4l(l+1)\omega^{2}+l(l+1)rf^{\prime}\omega_{\text{pl}}^{2}+r^{3}f^{\prime}\omega_{\text{pl}}^{4}\Bigr{\\}}-r\omega^{2}\Bigl{\\{}l(l+1)+2r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}g^{\prime}\biggr{)}+f^{3}\omega_{\text{pl}}\biggl{(}r\omega_{\text{pl}}\Bigl{\\{}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}g^{\prime}$
$\displaystyle+4l(l+1)g\bigl{\\{}\omega_{\text{pl}}+r\omega_{\text{pl}}^{\prime}\bigr{\\}}\biggr{)}\bigg{]}{U_{3}^{lm}}^{\prime}+\frac{r^{2}\omega^{2}-f\left\\{l(l+1)+r^{2}\omega_{\text{pl}}^{2}\right\\}}{r^{2}gf}U_{3}^{lm}=0~{}.$
(23)

Note that the above equation involves derivatives of the plasma frequency,
since we have assumed the plasma properties to be a function of the radial
coordinate. In the case of homogeneous plasma, satisfying
$\omega_{\text{pl}}=\textrm{constant}$, in a Schwarzschild background, with
$g=f=1-(2M/r)$, one can easily verify that III.2 reduces to Eq. (A1) of Ref.
[45]. In the general case as well, i.e., with $f\neq g$, it is possible to
recast III.2 in a Schrodinger-like form, which reads

$\frac{d^{2}\Psi^{lm}_{\text{pol}}}{dr_{*}^{2}}-V^{lm}_{\text{pol}}\Psi^{lm}_{\text{pol}}=0~{},$
(24)

where, $r_{*}$ is the Tortoise coordinate, defined in 16. In order to express
the polar sector differential equation to the above form, one must rescale
$U_{3}^{lm}$, so that the first derivative term in III.2 goes away. This is
achieved by introducing the master function $\Psi^{lm}_{\text{pol}}$, for the
polar sector, which is defined as, $\Psi_{\text{pol}}\equiv G(r)U_{3}^{lm}$
with333The function $G(r)$ is defined as the inverse of the auxiliary function
defined in [45]; the homogeneous plasma profile can be obtained from the above
equation by considering $\omega_{\text{pl}}$ to be a constant.

$G(r)=r\sqrt{\frac{f\omega_{\text{pl}}^{2}-\omega^{2}}{f\big{\\{}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\big{\\}}-r^{2}\omega^{2}}}~{}.$
(25)

Therefore, alike the odd sector, in the even sector as well, we have only one
propagating degree of freedom. Therefore, as in the massless case, here also
we have two propagating degrees of freedom. This is in striking contrast with
the Proca case, where one has three propagating degrees of freedom. The
complete expression for the effective potential $V_{\rm pol}$ in the polar
sector can be found in A. Therefore, 17 and 24, together with 18 and A in A
govern the photon dynamics through cold and unmagnetised plasma distribution
in the background of a generic static and spherically symmetric spacetime
geometry.

###  III.3 Plasma Profile

After describing the evolution equations for the axial and polar parts of EM
waves, we present here the details of the plasma profile considered in this
work. We will work with two distinct plasma profiles:

  * •

A homogeneous density profile, characterised by a constant electron/ion
density, yielding a constant plasma frequency $\omega_{\rm pl}$. Therefore,
the derivative $\omega_{\rm pl}^{\prime}$, along with all the higher
derivatives identically vanishes.

  * •

An inhomogeneous distribution of plasma, dependent on the radial coordinate
$r$, which may arise in a Bondi-like accretion disk model, describing the
accretion dynamics of a non-self-interacting gas around a spherically
symmetric compact object. The Bondi-like accretion model predicts a power-law
density profile of the electron/ion and consequently a plasma frequency of the
form,

$\omega^{2}_{\text{pl}}(r)=\omega^{2}_{B}\left(\frac{r_{0}}{r}\right)^{\lambda}+\omega^{2}_{\infty}~{},$
(26)

where, $r_{0}$ is the radius of the compact object, such that
$\sqrt{\omega^{2}_{B}+\omega^{2}_{\infty}}$ is the plasma frequency at the
surface of the compact object. For example, in the case of a BH as well as
exotic-compact object, one may consider $r_{0}$ to be replaced by $r_{+}$ in
26. Moreover, $\omega_{\infty}$ is the plasma frequency at infinity, such that
$\omega_{B}\ll\omega_{\infty}$. Finally, the parameter $\lambda$ depends on
the adiabatic index of the gas, e.g., $\lambda=(3/2)$ for a monoatomic gas.

###  III.4 Nature of the central compact object

In this section, we will describe the nature of the central compact object
around which the plasma is accreting, and whose exterior geometry is described
by the static and spherically symmetric line element as in 10. The analysis
assuming the compact object to be a BH has been performed in [45], while here
we focus on the case of ECO. These objects are ultra compact and have radii
close to the BH horizon. The difference between the surface radius of an ECO
from that of a BH is one of the characteristic properties of the ECO and is
traditionally denoted by $\epsilon$, with the following definition:

$\epsilon\equiv\frac{r_{0}-r_{+}}{r_{+}}~{}.$ (27)

Here, $r_{0}$ is the surface radius of the ECO and $r_{+}$ denotes the
location of the would be horizon444This is obtained by solving the equation
$f(r)=0=g(r)$. If the zeroes of $f(r)$ and $g(r)$ do not coincide, then rather
than a BH it will denote a wormhole spacetime, another class of ECO.. We will
consider cases with $\epsilon\ll 1$, and this requires either some exotic
matter field or, some quantum effects at play, since normal matter fields
cannot give rise to such compact structures. Furthermore, since the horizon is
replaced by a surface, it is expected that the EM waves will not be fully
absorbed, as it would have been for a BH, rather will also have a reflected
component. The reflectivity denotes another characteristic feature of an ECO.

To quantify the reflectivity, we note that the effective potential $V_{\rm
ax}^{lm}$ in the axial sector, presented in 18, vanishes in the limit $f(r)\to
0$, i.e., $r\to r_{+}$. Therefore, 17 has plane waves $e^{\pm i\omega r_{*}}$
as solutions. The same holds true for the polar sector as well, where $V_{\rm
pol}^{lm}$ will go to $-\omega^{2}$ in the limit $r\to r_{+}$. Therefore, 24
will also have plane wave modes as solutions. For a BH, only ingoing waves
will exist near the horizon, and hence the solution will be given by,
$e^{-i\omega r_{*}}$. While for an ECO there will be a reflected part and
hence the general solution will be,

$\Psi_{\rm ax/pol}^{lm}=Ae^{-i\omega r_{*}}+Be^{i\omega r_{*}}~{}.$ (28)

Following which, one may define the reflectivity of an ECO in the present
context of static and spherically symmetric background geometry as,

$R(\omega)\equiv\left[\frac{1-\dfrac{i}{\omega}\left(\dfrac{d\ln\Psi^{lm}_{\rm
ax/pol}}{dr_{*}}\right)}{1+\dfrac{i}{\omega}\left(\dfrac{d\ln\Psi^{lm}_{\rm
ax/pol}}{dr_{*}}\right)}\right]_{r_{0}}=\frac{B}{A}e^{2i\omega r^{0}_{*}}~{},$
(29)

where, $r^{0}_{*}=r_{*}(r_{0})$. Therefore, the ratio of the outgoing and
ingoing amplitudes of EM waves near the surface of the ECO is proportional to
the reflectivity, modulo an overall phase factor. Therefore any model of an
ECO will be characterized by the two parameter family $(\epsilon,R)$. Note
that, in general, the reflectivity is a function of the frequency $\omega$ of
EM waves. For example, in the case of BHs with quantum corrections, it
acquires a reflectivity: $R=\exp(-\hbar\omega/k_{\rm B}T_{\rm H})$, where
$T_{\rm H}$ is the Hawking temperature. On the other hand, in the context of
the membrane paradigm, the dominant contribution to the reflectivity is given
by, $|R|=(1-16\pi\eta)/(1+16\pi\eta)$. In the present context as well, we will
present our results for constant reflectivity, as well as for the Boltzman
reflectivity.

##  IV Plasma-driven quasi-bound modes for Schwarzschild-like ECO

Having described all the necessary ingredients, we now present the analysis
involving the stability of quasi-bound states for the EM waves propagating
within a plasma, which is acccreting around an ECO. For the QNMs, which are
scattering states, it is well-known that ECOs lead to a very small value for
the imaginary part of QNMs, responsible for the echoes in the GW ringdown. On
the other hand, for the quasi-bound states, the difference between BHs and
ECOs have not been studied extensively (however, see [10, 70]). In particular,
if any compact object is accreting, which is common for most of the compact
objects in an astrophysical environment, such quasi-bound states will
naturally form. Thus it is important to understand these quasi-bound states
and their implication for stability, if the central compact object is not a
BH, but an ECO.

To model the spherically symmetric ECO, we assume that the spacetime outside
the surface of the ECO is described by a Schwarzschild metric, with the metric
elements

$f(r)=g(r)=1-\frac{2M}{r}~{},$ (30)

where $r\geq r_{0}$, with $r_{0}=2M(1+\epsilon)$ being the location of the
surface of the ECO. Note that our previous analysis involving the evolution
equations for the EM field in the axial and the polar sector are completely
general and can be applied to any generic static and spherically symmetric
background. However, to see the basic features of the ECO on the quasi-bound
states we have kept the background spacetime simplest, namely Schwarzschild.
It is possible to generalize our analysis to any other static and spherically
symmetric geometry, including wormholes.

###  IV.1 Numerical Method

In this section we provide a brief analysis of the numerical method we have
adopted for determining the characteristic frequencies of plasma-driven quasi-
bound states for EM field in the background of an ECO. Throughout this paper,
we use the direct integration (shooting) method [47]. The idea is to integrate
the system of differential equations forward from the surface of the compact
object and backwards from a large radius (treated as infinity) to a common
matching point. The matching between the two solutions is achieved by
demanding that the Wronskian of the two solutions vanish at the matching
point. This yields the frequencies of the quasi-bound modes. The procedure
requires suitable boundary conditions to be specified close to the surface of
the compact object and near asymptotic infinity, however, it is not specific
to the nature of the effective potential and thus has versatile applicability.

Near the surface of the compact object, which has both ingoing and outgoing
modes present in the analysis, the solution of the differential equations
governing the axial and polar modes, may be written as the sum of a purely
ingoing wave and a purely outgoing wave,

$\displaystyle\Psi^{lm}_{\text{ax}/\text{pol}}$ $\displaystyle\sim
e^{-i\omega(r_{*}-r^{0}_{*})}\sum_{n}h^{\text{ax}/\text{pol}}_{n}(r-2M)^{n}$
$\displaystyle+R(\omega)e^{i\omega(r_{*}-r^{0}_{*})}\sum_{n}k^{\text{ax}/\text{pol}}_{n}(r-2M)^{n}~{}.$
(31)

Here, $R(\omega)$ is the reflectivity of the surface of the ECO, defined in
29. Note that $r=2M$ is a regular singular point of the evolution equations
for the axial, 17, and the polar, 24, sectors. Hence the series around $r=2M$
has a radius of convergence $\gg\epsilon$, so that the above series applies to
the ECO as well.

Just as the above series describes the master functions for polar and axial
sectors near the surface of the ECO, near asymptotic infinity as well, similar
power series expansion of the solutions are possible and they are generically
given by,

$\Psi^{lm}_{\rm
ax/pol}\sim\sum_{n}\frac{B_{n}^{\text{ax}/\text{pol}}}{r^{n}}e^{-k_{\infty}r_{*}}+\sum_{n}\frac{C_{n}^{\text{ax}/\text{pol}}}{r^{n}}e^{k_{\infty}r_{*}}~{},$
(32)

where $k_{\infty}=\sqrt{\omega_{\infty}^{2}-\omega^{2}}$, with
$\omega_{\infty}$ being the plasma frequency at asymptotic infinity. Since we
are interested in bound states, which are strongly localized near the ECO,
their amplitudes at large distance must be exponentially small. Therefore, for
the existence of quasi-bound states we assume that — (a)
$\omega<\omega_{\infty}$, and (b) $C_{n}^{\text{ax}/\text{pol}}=0$. This
corresponds to the exponentially decaying solution near the spatial infinity,
which are highly localized near the compact object.

Using the approach and the boundary conditions specified above, in the
following section, we will discuss the fundamental quasi-bound frequencies of
the EM waves propagating on the background of a Schwarzschild-like ECO in the
presence of accreting plasma. We would like to point out that tracking the
fundamental mode (along with overtones) using the direct integration method
becomes extremely challenging for small $M\omega_{\text{pl}}$, as the
imaginary part of the quasi-bound frequency becomes very small, and hence
prone to numerical errors. This applies to both BHs and ECOs. However, in the
case of ECOs, there are additional numerical challenges, to be described
later. The numerical accuracy of the results in general is also dependent on
the truncation order of the series at both the boundaries. To ensure the
robustness of our numerical results: (a) we verified that the frequencies
obtained are independent of the location of the matching point, (b) for each
case, we ensured that the change induced in the frequencies due to
(reasonable) change in the truncation order is below the precision of the
obtained results and (c) in each case, we scanned the parameter space in small
steps to ensure that same overtone mode is tracked through the parameter
space.

###  IV.2 Quasi-bound states: Axial sector

With the numerical scheme and appropriate boundary conditions at hand, let us
determine the quasi-bound states associated with the axial sector of EM waves.
We start by assuming a homogeneous plasma profile, whose effective potential
for a Schwarzschild-like ECO with mass $M$ is given by

$V^{lm}_{\text{ax}}=\left(1-\frac{2M}{r}\right)\left(\omega_{\text{pl}}^{2}+\frac{l(l+1)}{r^{2}}\right)~{}.$
(33)

For quasi-bound state to exist for both BHs and ECOs, the above potential must
have two extrema — one maxima and one minima, and the maxima must be closer to
the compact object than the minima (see 1 for an example). Taking derivatives
of the effective axial potential $V^{lm}_{\text{ax}}$, it follows that, the
above conditions are satisfied for $M\omega_{\rm pl}<\sqrt{l(l+1)/12}$ (for a
derivation, see B). Thus, there exists an upper bound for the plasma
frequency.

Figure 1: Plot of the effective potential $V_{\rm ax}$ in the axial sector for
a Schwarzschild-like ECO (with $M=1$) as a function of the Tortoise coordinate
$r_{*}$, with $M\omega_{\text{pl}}=0.2$, $\epsilon=10^{-5}$ and $l=1$. The
dashed line represents the surface of the ECO, while the inset figure shows
the minima of the effective potential resulting in bound states. Note that the
potential vanishes near the surface of the ECO, but approaches to a constant
value $\omega_{\rm pl}^{2}=0.04$ asymptotically.

Figure 2: Plots of the imaginary and real parts of the frequencies of
fundamental quasi-bound states associated with EM waves propagating in the
background of Schwarzschild-like ECO in the axial sector, with the plasma
frequency $M\omega_{\text{pl}}$ have been presented. For these plots we have
chosen $\epsilon=10^{-5}$ and three possible choices of constant reflectivity,
$R=0.2,0.5$ and $-0.5$, respectively. The corresponding frequencies associated
with a Schwarzschild BH have also been depicted (the black line in each plot).
The left panel shows the frequencies of quasi-bound states for positive
reflectivity, whereas the right panel shows quasi-bound state frequencies for
both positive and negative reflectivities. The inset in each panel shows the
deviation of the real part of the quasi-bound state frequencies from that of a
Schwarzschild black hole.

Given all these inputs and the numerical method depicted above, the
fundamental quasi-bound frequencies for the Schwarzschild black hole (which
corresponds to the case $R=0$) and for Schwarzschild-like ECOs have been
plotted in 2. We observe that for a fixed value of the reflectivity, both the
decay rate (inverse of the imaginary part of the quasi-bound frequency) and
the frequency of the quasi-bound state oscillate about the corresponding black
hole values. This oscillation is reminiscent of the “two-well” structure of
the effective potential, one between the maxima of the potential and the
surface of the ECO and the other originating from the minima of the effective
potential $V_{\rm ax}^{lm}$, which arises due to the effective mass of the
photon as it propagates through the plasma. Akin to the QNMs, the relative
deviation of the real part of the quasi-bound frequencies for ECOs, from the
black hole values, is smaller compared to the corresponding deviation of the
decay rate (see 2). Thus, depending on the plasma frequency,
$M\omega_{\text{pl}}$, which in turn depends on the electron number density,
the quasi-bound states for EM waves around an ECO can be either long-lived or
short-lived, compared to a BH. In particular, the deviation of the imaginary
part of the quasi-bound states for ECOs from BHs is at most by a factor of 10,
as evident from 2, while for QNMs such a deviation is $\mathcal{O}(10^{6})$
and higher depending on the model of ECO. Therefore, in striking contrast to
QNMs, the frequencies of quasi-bound states of ECOs are more closer to the
corresponding values for BHs. These results are also supported by analytical
computations, which have been presented in detail in C.

Dependence of the quasi-bound frequencies on the reflectivity can also be
assessed from 2. As evident, with increasing reflectivity, the deviation of
the decay rate from the black hole value increases. For higher values of
$M\omega_{\text{pl}}$, the real part of the frequency also shows significant
deviation from the black hole results, while the imaginary part shows stronger
decay (compared to black holes). Negative values of reflectivity, presented in
the second panel of 2 imply a phase reversal of the reflected wave from the
ECO surface. So far we have limited our analysis for $M\omega_{\rm pl}>0.2$.
This is because, extending the results for $M\omega_{\text{pl}}<0.2$ is
extremely challenging, particularly for larger reflectivity due to numerical
inaccuracies and it becomes exceedingly difficult to track the fundamental
mode, mostly because of the oscillatory nature. We would like to emphasize
that the above numerical instability arises solely in the case of ECOs, as
argued in the previous section.

Figure 3: Left: plots of the imaginary and real parts of the plasma-driven
quasi-bound states frequencies for EM fields associated with a Schwarzschild-
like ECO in the axial sector have been presented against the plasma frequency
$M\omega_{\text{pl}}$. Here we have fixed the reflectivity to be $R=0.5$, but
consider three different positions of the reflective boundary. The inset shows
the deviation of the real part of the quasi-bound state frequency from that of
a Schwarzschild black hole. Right: plots of the imaginary and real parts of
the quasi-bound state frequencies with $\epsilon$ have been presented for
fixed reflectivity, $R=0.5$ and three possible choices of the plasma frequency
$M\omega_{\rm pl}$. The inset in the lower panel shows the zoomed-in view of
the real part of quasi-bound frequencies, for $M\omega_{\text{pl}}=0.22$.

.

The quasi-bound state frequencies for a given reflectivity also strongly
depend on the position of the ECO surface, as can be observed from 3. As the
ECO surface moves closer to the maxima of the potential, i.e., as $\epsilon$
increases, the trapping region of the potential in the interior of the maxima
becomes smaller (see 1), and the oscillation length of the imaginary part of
the frequency extends further, covering more regions of the
$M\omega_{\text{pl}}$ parameter space. From the right hand side plot of 3 it
is evident that for larger values of $M\omega_{\rm pl}$, the variation of the
real and imaginary parts of the frequencies of quasi-bound states with
$\epsilon$ is the strongest, while for smaller values of $M\omega_{\rm pl}$,
the variations with $\epsilon$ are weaker. Also the oscillation length of the
imaginary part of quasi-bound state frequencies depend on $M\omega_{\rm pl}$,
increasing with decreasing values of $M\omega_{\rm pl}$. All of these results
are consistent with analytical estimation of the bound state frequencies,
presented in C, showing oscillatory behaviour, as observed in the numerical
analysis.

Figure 4: Plots of the imaginary and the real parts of the fundamental quasi-
bound frequencies for axial EM field around an accreting Schwarzschild-like
ECO have been presented, against the plasma frequency $\omega_{\text{pl}}$. We
have taken $\epsilon=10^{-5}$, and the reflectivity of the ECO has been taken
to be Boltzmann reflectivity. The insets in each plots show the deviation from
the black hole results.

So far our analysis assumes constant reflectivity, however as emphasized
before, the reflectivity of an ECO can in general be frequency dependent. To
illustrate implications for such a frequency dependent reflectivity, we
consider the case of Boltzmann reflectivity. The Boltzmann reflectivity is
akin to any quantum corrected BH, which emits Hawking radiation at temperature
$T_{\rm H}$, and since one can assume the BHs to be not perfectly black,
rather reflecting as a black body at temperature $T_{\rm H}$, it is natural to
assume $R=\exp(-\hbar\omega/k_{\rm B}T_{\rm H})$. The quasi-bound state
frequencies for an ECO with Boltzmann reflectivity have been presented in 4.
As expected, in this context as well, we observe that both the frequency and
the decay rate of the plasma-driven EM quasi-bound states oscillate about the
black hole values. The deviation from the black hole values increases with
increasing $M\omega_{\text{pl}}$ coupling, and hence with increasing the
electron number density in the plasma. This is because, increasing
$\omega_{\rm pl}$ increases the effective mass of the photon, leading to
deeper potential well. This leads to significant deviation from the BH quasi-
bound state frequencies.

Figure 5: Plots of the imaginary part of the plasma-driven quasi-bound
frequencies of an axial EM field, in the background of a Schwarzschild-like
ECO, has been presented against the asymptotic frequency $\omega_{\infty}$ for
Bondi accretion. We have assumed $\epsilon=10^{-5}$ with $\lambda=(3/2)$,
i.e., monoatomic gas for the accreting matter. As evident, increasing values
for the near-horizon frequency Here we have compared two cases for
$\omega_{B}=\frac{1}{M}$ and $\omega_{B}=\frac{0.9}{M}$.

####  IV.2.1 Bondi Accretion and quasi-bound state frequencies

We have described a homogeneous plasma profile and the quasi-bound state
frequencies associated with axial EM fields. In this section, we would like to
discuss the more realistic Bondi-like accretion, which fits well the spherical
symmetry of the background spacetime. The plasma frequency, in the present
context is given by 26. In 5, we show the behaviour of the imaginary part of
the axial quasi-bound state frequency with $M\omega_{\infty}$ for Bondi
accretion. We fix $\lambda=3/2$, corresponding to monoatomic gas. From 5 it is
clear that as in the homogeneous case, the imaginary part of the quasi-bound
state frequency for a given reflectivity oscillates about the monotonically
increasing BH results. The real part of the quasi-bound state frequency
behaves similarly as in the homogeneous distribution, however, the deviation
from the BH result is extremely small and hence has not been depicted here. As
in the BH case [45], for ECOs as well, the decay rate of the quasi-bound state
around a Schwarzschild-like ECO decreases with increasing the surface plasma
frequency $\omega_{\rm B}$. The corresponding deviation from the BH results,
on the other hand, decreases with increasing $\omega_{\rm B}$. Thus, unlike
BHs, for certain choices of $\omega_{\rm B}$, these quasi-bound states
associated with ECOs decay sufficiently slowly, making these plasma
frequencies susceptible to instabilities. Again this is a feature unique to
the ECOs. In general, for a given surface plasma frequency $\omega_{\rm B}$,
the imaginary part of the quasi-bound state frequency is smaller compared to
the homogeneous plasma distribution. This makes it more susceptible to
numerical instability.

###  IV.3 Quasi-bound states: Polar sector

In this final section we discuss the frequencies of quasi-bound states
associated with the polar sector. The main difficulty with the polar sector is
the complicated structure of the effective potential, presented in 24, and in
particular, the nontrivial dependence of the effective potential on $\omega$.
Therefore, we work exclusively on the homogeneous plasma distribution and
depict the variation of the plasma-driven EM quasi-bound frequencies with
$M\omega_{\text{pl}}$ in 6.

Figure 6: The imaginary and the real parts of the fundamental plasma-driven
quasi-bound frequencies have been depicted, against the plasma frequency, for
the polar sector of an EM wave moving in the background of a Schwarzschild-
like ECO. We have taken the compactness of the ECO, such that
$\epsilon=10^{-5}$, and the reflectivity is taken to be $R=\pm 0.0001$. The
black curve corresponds to the frequencies associated with Schwarzschild BH.
The inset in each plot shows the corresponding deviation from the BH results.

We note that, as in the axial case, both the frequency and the decay rate
oscillate about the BH result, and the amplitude of oscillation increases with
$M\omega_{\text{pl}}$. This is due to deepening of the potential well with
increasing plasma frequency. Due to the nontrivial dependence of the effective
potential in the polar sector on $\omega$, the numerical evolution becomes
increasingly unstable as the reflectivity is increased, and hence we had to
work with small reflectivities. The general behaviour of the quasi-bound state
frequencies and the decay rates remain similar to that in the axial sector.
For Boltzmann-reflectivity, on the other hand, the quasi-bound state
frequencies shows similar behaviour as in the axial sector, however, since the
amplitude of deviation from the BH results is even smaller it has not been
depicted here.

For the Bondi accretion, with inhomogeneous plasma frequency, the effective
potential in the polar sector becomes even more complicated and the numerical
accuracy becomes worse. The small values of the imaginary part of the quasi-
bound state frequencies add up to the difficulty.

As evident from 6, even in the polar sector, depending on the value of the
plasma frequency, the quasi-bound states associated with ECO either decay with
a longer or, shorter time scale compared to those around a BH. In particular,
the oscillating pattern around the BH value is seen in the polar sector as
well. It is important to comment that throughout our analysis we have not
observed any unstable modes, i.e., modes with positive imaginary parts. Hence
within the purview of our analysis, the plasma-driven EM quasi-bound states
around static and spherically symmetric ECOs are stable, but the degree of
stability is weaker or stronger compared to the corresponding BH.

##  V Conclusion

Proving the existence of BHs and exploring its properties has been one of the
pinnacle for research on gravitational physics. With the detection of GWs by
the LVK collaboration, we are merely one step away from unveiling the true
nature of BHs, which are abundant in our universe, from the center of our
galaxies to various astrophysical scenarios. In order to really test the BH
hypothesis, we need alternatives to BHs, and these objects are collectively
dubbed as ECOs, which are horizonless but possess a light ring similar to that
of a BH. These ECOs are characterized by a real surface with a radius greater
than the event horizon of a corresponding BH. As far as an external observer
is concerned, the interior properties of an ECO are represented by two
physical quantities — (a) the reflectivity $R$ of the ECO surface, and (b) the
location of its surface, characterized by $\epsilon$. There have been plethora
of works, discussing various properties of these ECOs, but mostly in the
absence of any environment. In the present work, we have investigated the
effect of accretion onto the stability of such an ECO. In detail, we have
studied the propagation of EM waves through a cold and collisionless plasma
accreting onto a Schwarzschild-like ECO. The symmetry of the background
geometry allowed us to decompose the evolution equations for the EM field in
axial and polar sectors, which decouples from one another. We have presented,
for the first time, the equations governing the dynamics of the EM wave
propagating in a plasma within a general static and spherically symmetric
background, both for the polar and the axial sectors. Though, in the axial
sector the governing dynamical equations are identical to that of a Proca
field, in the polar sector, the dynamical equations for the EM field are
completely different. This observation matches with [45].

Subsequently, we have applied our general equations to study the quasi-bound
states of the EM waves and their stability for Schwarzschild-like ECO. In
particular, our interest lies in understanding the effect of reflectivity and
the location of the ECO surface on the frequencies of the plasma-driven quasi-
bound states. These states are particularly interesting in a rotating compact
object as they are prone to superradiant instability. Of course, a spherically
symmetric compact object is not subjected to superradiance but a study of the
quasi-bound states even in such a system may provide valuable insight into the
stability of such compact objects.

To study the frequency spectrum of quasi-bound states, we employ the direct
integration method with reflective boundary conditions near the surface of the
ECO and exponentially decaying behaviour asymptotically. In all the cases
studied, we observed the general trend that the decay rate of the quasi-bound
states increases with the plasma frequency $M\omega_{\text{pl}}$, whereas the
real frequency of the quasi-bound state decreases with the same. For a given
reflectivity both the real and the imaginary parts of the quasi-bound state
frequencies oscillates about the BH results. The amplitude of oscillation
increases with the reflectivity. On the other hand, as the compactness of the
ECO increases (smaller values of $\epsilon$), the oscillation length of the
frequencies for quasi-bound states decreases. Therefore, the
$M\omega_{\text{pl}}$ parameter space gets populated with more short and long
lived modes compared to BHs. Thus, for a Schwarzschild-like ECO with a given
mass, depending on the electron number density, the quasi-bound states may be
short or long-lived compared to a BH. This is true for both homogeneous and
Bondi-like plasma distribution, with constant as well as frequency-dependent
reflectivity, in both the polar and axial sectors. It is important to note
that the relative deviation of the real part of the quasi-bound frequencies of
the ECO from the BH quantities are much smaller than the corresponding
deviation in the decay rate. Moreover, the axial modes are short-lived
compared to the polar modes, in general. All of these numerical results are
consistent with our analytical findings presented in C. Finally, we have also
observed that, if the reflectivity of the ECO surface is frequency-dependent,
as in the case of Boltzmann reflectivity, the deviation from the BH result
decreases.

Through our analysis we were able to show that, unlike the QNM frequencies,
the frequencies of quasi-bound states of EM waves associated with ECOs differs
from the BH values by an order or two. However, the frequencies of quasi-bound
states around ECOs describe a unique oscillation pattern around BH values,
whose amplitude is related to the reflectivity of the ECO, and the oscillation
length is dependent on $\epsilon$ — a new feature for ECOs. Our analyses
enabled us to explore a considerable portion of the parameter space within
which, we did not find any unstable modes indicating the stability of the ECOs
under plasma driven EM perturbation. The results also strongly indicate that
in the case of rotating ECOs, even in the linear regime, the presence of
plasma will result in stronger quenching (compared to black holes) of
superradiant instability for a certain part of the parameter space, where the
imaginary part of quasi-bound state frequency is larger compared to a BH,
while for some other parameters the superradiant instability will be stronger
for ECOs. Thus, it will be interesting to study how the presence of plasma
affects the instability time scale in the case of rotating ECOs [70], which we
hope to address in a separate work. Along these lines, it will also be
interesting to study the effect of nonlinearity and magnetic field in the
photon-plasma interaction on the background of ECOs. Similar studies are
already underway for black holes [60, 62, 71] and we hope to report the
consequences for ECOs elsewhere.

###### Acknowledgements.

SB thanks IACS for financial support and Somsubhra Ghosh and Pratick Sarkar
for some helpful discussion. AC and SC acknowledges hospitality from IIT
Gandhinagar and IUCAA, where a substantial part of the work was done. AC
acknowledges financial support from SERB, Government of India (File No.
PDF/2023/000550). Research of SC is supported by the MATRICS and the Core
research grants from SERB, Government of India (Reg. Nos. MTR/2023/000049 and
CRG/2023/000934).

sectionAppendix #1

##  Appendix A Effective potential in the Polar sector

In this appendix we will present the full expression for the effective
potential in the polar sector. The expression involves the metric function, as
well as the plasma frequency and their derivatives, which is given by

$\displaystyle V^{lm}_{\text{pol}}$
$\displaystyle=\frac{1}{4r^{2}\Bigl{(}\omega^{2}-f\omega_{\text{pl}}^{2}\Bigr{)}^{2}\Bigl{[}r^{2}\omega^{2}-f\Bigl{\\{}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}\Big{]}^{2}}\Bigg{[}-r^{4}\omega^{6}\Bigl{\\{}4r^{2}\omega^{4}-l(l+1)gf^{\prime
2}\Bigr{\\}}+r^{2}\omega^{4}f\Biggl{\\{}2r^{2}\omega_{\text{pl}}^{2}\Bigl{(}10r^{2}\omega^{4}$
$\displaystyle\qquad+l(l+1)gf^{\prime
2}\Bigr{)}+l(l+1)\biggl{[}r^{2}\omega^{2}\bigl{(}12\omega^{2}+f^{\prime}g^{\prime}\bigr{)}+2g\Bigl{(}-5r\omega^{2}f^{\prime}+l(l+1)f^{\prime
2}+r^{2}\omega^{2}f^{\prime\prime}\Bigr{)}\biggr{]}\Biggr{\\}}$
$\displaystyle\qquad-r\omega^{2}f^{2}\Bigg{\\{}r\omega^{2}\biggl{(}40r^{4}\omega^{2}\omega_{\text{pl}}^{4}+2l(l+1)r^{2}\omega_{\text{pl}}^{2}\bigl{(}24\omega^{2}+f^{\prime}g^{\prime}\bigr{)}+l(l+1)\Bigl{[}12l(l+1)\omega^{2}+\bigl{(}2r\omega^{2}+l(l+1)f^{\prime}\bigr{)}g^{\prime}\Bigr{]}\biggr{)}$
$\displaystyle\qquad+l(l+1)g\Biggl{(}3r\omega_{\text{pl}}^{2}\Bigl{[}l(l+1)+r^{2}\omega_{\text{pl}}^{2}\Bigr{]}f^{\prime
2}+2\omega^{2}f^{\prime}\Bigl{[}l(l+1)-11r^{2}\omega_{\text{pl}}^{2}-9r^{3}\omega_{\text{pl}}\omega^{\prime}_{\text{pl}}\Bigr{]}+2r\omega^{2}\biggl{[}-6\omega^{2}$
$\displaystyle\qquad+\Bigl{\\{}l(1+l)+2r^{2}\omega_{\text{pl}}^{2}\Bigr{\\}}f^{\prime\prime}\biggr{]}\Biggr{)}\Bigg{\\}}+2f^{5}\omega_{\text{pl}}^{3}\Biggl{\\{}6l(l+1)r^{4}\omega_{\text{pl}}^{5}+2r^{6}\omega_{\text{pl}}^{7}+l(l+1)r^{2}\omega_{\text{pl}}^{3}\Bigl{(}6l(l+1)-6g+rg^{\prime}\Bigr{)}$
$\displaystyle\qquad+l(l+1)\omega_{\text{pl}}\biggl{(}l(l+1)rg^{\prime}+2\Bigl{\\{}l^{2}(l+1)^{2}-3r^{4}gf^{\prime
2}\Bigr{\\}}\biggr{)}+l^{2}(l+1)^{2}r\Bigl{\\{}rg^{\prime}\omega^{\prime}_{\text{pl}}+2g\bigl{(}2\omega^{\prime}_{\text{pl}}+r\omega^{\prime\prime}_{\text{pl}}\bigr{)}\Bigr{\\}}$
$\displaystyle\qquad+l(l+1)r^{3}\omega_{\text{pl}}^{2}\Bigl{[}rg^{\prime}\omega^{\prime}_{\text{pl}}+g\bigl{(}-8\omega^{\prime}_{\text{pl}}+2r\omega^{\prime\prime}_{\text{pl}}\bigr{)}\Bigr{]}\Biggr{\\}}$
$\displaystyle\qquad-2f^{4}\Biggl{\\{}10r^{6}\omega^{2}\omega_{\text{pl}}^{8}+l(l+1)r^{3}\omega_{\text{pl}}^{6}\bigl{(}24r\omega^{2}-gf^{\prime}\bigr{)}+l(l+1)r\omega_{\text{pl}}^{4}\Bigl{[}-g\bigl{(}18r\omega^{2}+l(l+1)f^{\prime}\bigr{)}+3r\omega^{2}\bigl{(}6l(l+1)+rg^{\prime}\bigr{)}\Bigr{]}$
$\displaystyle\qquad-l(l+1)r^{4}g\omega_{\text{pl}}^{5}f^{\prime}\omega^{\prime}_{\text{pl}}+2l^{2}(l+1)^{2}r^{2}\omega^{2}gf^{\prime
2}+2l(l+1)\omega^{2}\omega_{\text{pl}}^{2}\biggl{(}l(l+1)rg^{\prime}+2\Bigl{\\{}l^{2}(l+1)^{2}-r^{4}gf^{\prime
2}\Bigr{\\}}\biggr{)}$
$\displaystyle\qquad+l^{2}(l+1)^{2}r\omega^{2}\omega_{\text{pl}}\Bigl{\\{}rg^{\prime}\omega^{\prime}_{\text{pl}}+2g\bigl{(}2\omega^{\prime}_{\text{pl}}+r\omega^{\prime\prime}_{\text{pl}}\bigr{)}\Bigr{\\}}-l(l+1)r^{2}\omega_{\text{pl}}^{3}\biggl{(}-2r^{2}\omega^{2}g^{\prime}\omega^{\prime}_{\text{pl}}+g\Bigl{[}\bigl{\\{}16r\omega^{2}+l(l+1)f^{\prime}\bigr{\\}}\omega^{\prime}_{\text{pl}}$
$\displaystyle\qquad-4r^{2}\omega^{2}\omega^{\prime\prime}_{\text{pl}}\Bigr{]}\biggr{)}\Biggr{\\}}+\omega^{2}f^{3}\Biggl{\\{}40r^{6}\omega^{2}\omega_{\text{pl}}^{6}-20l(l+1)r^{4}g\omega_{\text{pl}}^{3}f^{\prime}\omega^{\prime}_{\text{pl}}+2l(l+1)\omega^{2}\biggl{(}l(l+1)rg^{\prime}+2\Bigl{[}l^{2}(l+l)^{2}+r^{4}gf^{\prime
2}\Bigr{]}\biggr{)}$
$\displaystyle\qquad+l(l+1)r^{3}\omega_{\text{pl}}^{4}\bigl{(}72r\omega^{2}-14gf^{\prime}+rf^{\prime}g^{\prime}+2rgf^{\prime\prime}\bigr{)}+l(l+1)r^{2}\omega_{\text{pl}}^{2}\Bigl{[}36l(l+1)\omega^{2}+\bigl{(}6r\omega^{2}+l(l+1)f^{\prime}\bigr{)}g^{\prime}$
$\displaystyle\qquad+g\bigl{(}-36\omega^{2}+2l(l+1)f^{\prime\prime}\bigr{)}\Bigr{]}-2l(l+1)r^{2}\omega_{\text{pl}}\biggl{(}-r^{2}\omega^{2}g^{\prime}\omega^{\prime}_{\text{pl}}+g\Bigl{[}\bigl{\\{}8r\omega^{2}+3l(l+1)f^{\prime}\bigr{\\}}\omega^{\prime}_{\text{pl}}-2r^{2}\omega^{2}\omega^{\prime\prime}_{\text{pl}}\Bigr{]}\biggr{)}\Biggr{\\}}\Bigg{]}~{}.$
(34)

This is the expression we have used in our numerical analysis to determine the
frequencies of the quasi-bound states.

##  Appendix B Relative location of maxima and minima for the axial effective
potential

The axial effective potential $V_{\rm ax}^{lm}$ is given by 18, which admits a
minima and a maxima. The maxima in the potential is akin to the light ring,
while the minima specifically originates from the effective mass associated
with the axial sector of the EM field. For existence of quasi-bound orbits it
is necessary that the minima lies at a larger radial coordinate than the
maxima, which gies a condition on the parameters of the problem. Setting
${V_{\rm ax}^{lm}}^{\prime}=0$, we obtain the following algebraic equation,

$\displaystyle
r^{2}-\frac{l(l+1)}{M\omega^{2}_{\text{pl}}}r+\frac{3l(l+1)}{\omega^{2}_{\text{pl}}}=0~{},$
(35)

whose solutions depict the locations of the extrema of the potential ${V_{\rm
ax}^{lm}}$. As evident from the above equation, both the roots will exist iff
the following condition holds,

$\displaystyle M\omega_{\text{pl}}\leq\sqrt{\frac{l(l+1)}{12}}~{}.$ (36)

The above inequality suggests that for $l=1$, for the existence of double
root, we must have $M\omega_{\text{pl}}\leq(1/\sqrt{6})$. When this condition
is satisfied, the two roots of the above algebraic equation read,

$\displaystyle r_{>}=$
$\displaystyle\frac{1}{2}\left[\frac{l(l+1)}{M\omega_{\text{pl}}^{2}}+\sqrt{\frac{l(l+1)}{\omega_{\text{pl}}^{2}}\left[\frac{l(l+1)}{M^{2}\omega_{\text{pl}}^{2}}-12\right]}\right]~{},$
(37) $\displaystyle r_{<}=$
$\displaystyle\frac{1}{2}\left[\frac{l(l+1)}{M\omega_{\text{pl}}^{2}}-\sqrt{\frac{l(l+1)}{\omega_{\text{pl}}^{2}}\left[\frac{l(l+1)}{M^{2}\omega_{\text{pl}}^{2}}-12\right]}\right]~{}.$
(38)

One can easily check that $r_{>}>r_{<}$, along with the fact that $r_{>}$ is
the location of the minima, while $r_{<}$ is the location of the maxima.
Implying that the minima is always at a larger distance than the maxima. It
can also be verified by computing the potential at the location of extrema,
which reads,

$\displaystyle
V(r)=\omega_{\text{pl}}^{2}+\left[\frac{4Ml(l+1)}{r^{3}}-\frac{l(l+1)}{r^{2}}\right]~{}.$
(39)

such that,

$\displaystyle
V(r_{<})-V(r_{>})=\frac{(r_{>}-r_{<})l(l+1)}{r_{>}^{2}r_{<}^{2}}\left[\frac{l(l+1)}{3M\omega_{\text{pl}}^{2}}-4M\right]\geq
0~{}.$ (40)

Therefore, considering the monopole mode, as long as
$M\omega_{\text{pl}}\leq(1/\sqrt{6}))\simeq 0.4$, there exists a trapping well
outside the light ring, resulting into the formation of quai-bound states.
This result has been used in the main text.

##  Appendix C Analytical estimate of the axial quasi-bound frequencies for
Schwarzschild-like ECO

In this appendix, we will provide analytical results demystifying the
oscillatory behaviour of the axial quasi-bound frequencies associated with
ECOs, derived through numerical methods. For this purpose, we start with the
following equation governing the propagation of axial EM wave,

$\displaystyle
x^{2}(1+x)^{2}{\Psi_{\text{ax}}^{lm}}^{\prime\prime}+x(1+x){\Psi_{\text{ax}}^{lm}}^{\prime}+\left[\omega^{2}r^{2}_{+}(1+x)^{4}-x(1+x)\left\\{\omega^{2}_{\text{pl}}r_{+}^{2}(1+x)^{2}+l(l+1)\right\\}\right]\Psi_{\text{ax}}^{lm}=0~{},$
(41)

where, we have introduced the dimensionless radial coordinate
$x\equiv(r/r_{+})-1$, and ‘prime’ denotes derivative with respect to $x$. In
what follows we will divide the spacetime into three regions — (i) the near
horizon region (or, near ECO region) defined as $(M\omega_{\text{pl}})x\ll 1$
and $(M\omega)x\ll 1$, (ii) the far region, which is defined as $x\gg 1$, and
most importantly, (iii) the intermediate region, defined as $1\ll
x\ll(l/M\omega)$. We will solve the above equation analytically in the near
horizon and far region, and shall match the solutions in the intermediate
region to determine the frequencies.

Near the horizon, using the approximations discussed above, 41 takes the
following form,

$\displaystyle
x^{2}(1+x)^{2}{\Psi^{lm}_{\text{ax}}}^{\prime\prime}+x(x+1){\Psi^{lm}_{\text{ax}}}^{\prime}+\left[\omega^{2}r_{+}^{2}-x(x+1)l(l+1)\right]\Psi_{\text{ax}}^{lm}=0~{}.$
(42)

For simplicity, we will henceforth set $M=1$, and the solution to the above
equation takes the following form,

$\displaystyle\Psi^{lm~{}\text{(near)}}_{\text{ax}}$
$\displaystyle=A\,x^{-2i\omega}(1+x)^{1+\delta}\,{}_{2}F_{1}(-l-2i\omega+\delta,l+1-2i\omega+\delta,1-4i\omega,-x)$
$\displaystyle\qquad+B\,x^{2i\omega}(1+x)^{1+\delta}\,{}_{2}F_{1}(-l+2i\omega+\delta,l+1+2i\omega+\delta,1+4i\omega,-x)~{},$
(43)

where $\delta\equiv\sqrt{1-4\omega^{2}}$, with $A$ and $B$ being arbitrary
constants to be fixed later by imposing appropriate boundary conditions. One
such boundary condition is imposed by taking $x\to 0$ limit of the above near
horizon solution, which takes the following form,

$\displaystyle\Psi_{\text{ax}}^{lm~{}\text{(near)}}(x\to 0)\simeq Ae^{-i\omega
r_{*}}+Be^{+i\omega r_{*}}~{},$ (44)

and then realizing that the ECO surface is reflective in nature. This allows
us to relate the ratio $(B/A)$ to the reflectivity of the ECO surface as in
29.

Further, using the asymptotic behaviour of the hypergeometric function, we can
show that in the intermediate region the near horizon solution takes the
following form,

$\displaystyle\\!\\!\\!\\!\\!\\!\Psi_{\text{ax}}^{lm~{}\text{(near)}}\left(x\rightarrow
x_{\text{int}}\right)$
$\displaystyle\simeq\Bigg{[}A\left(\frac{\Gamma(1-4i\omega)\Gamma(2l+1)}{\Gamma(l+1-2i\omega+\delta)\Gamma(l+1-2i\omega-\delta)}\right)+B\left(\frac{\Gamma(1+4i\omega)\Gamma(2l+1)}{\Gamma(l+1+2i\omega+\delta)\Gamma(l+1+2i\omega-\delta)}\right)\Bigg{]}x^{l+1}$
$\displaystyle+\Bigg{[}A\left(\frac{\Gamma(1-4i\omega)\Gamma(-2l-1)}{\Gamma(-l-2i\omega+\delta)\Gamma(-l-\delta-2i\omega)}\right)+B\left(\frac{\Gamma(1+4i\omega)\Gamma(-2l-1)}{\Gamma(-l+2i\omega+\delta)\Gamma(-l+2i\omega-\delta)}\right)\Bigg{]}x^{-l}~{}.$
(45)

Here, $1\ll x_{\rm int}\ll(l/\omega)$. As evident, the part of the solution in
the first line denotes the growing mode, while the part in the second line
depicts decaying mode.

Along identical lines, in the far region, the dynamical equation for the axial
sector of the plasma-driven EM modes becomes,

$\displaystyle
x^{2}{\Psi_{\text{ax}}^{lm}}^{\prime\prime}+\left[-q^{2}x^{2}+\nu
qx-l(l+1)\right]\Psi_{\text{ax}}^{lm}=0~{};\quad
q\equiv-\sqrt{\omega_{\text{pl}}^{2}-\omega^{2}}~{};\quad\nu\equiv\frac{2\omega^{2}-\omega_{\text{pl}}^{2}}{q}~{}.$
(46)

The general solution of the above differential equation is given by a linear
combination of confluent hypergeometric functions, which read,

$\displaystyle\Psi_{\text{ax}}^{lm~{}\text{(far)}}=e^{-2qx}\left[C\,(4qx)^{l+1}\,U(l+1-\nu,2l+2,4qx)+D\,(4qx)^{-l}\,M(-l-\nu,-2l,4qx)\right]~{}.$
(47)

For depicting a bound state, we demand the above solution to be an
exponentially decaying one near infinity, and using the properties of
confluent hypergeometric function one can show that the above condition will
hold true if and only if, we choose $D=0$. With this choice, in the
intermediate regime, we have the following behaviour for the far region
solution,

$\displaystyle\Psi_{\text{ax}}^{lm~{}\text{(far)}}\left(x\rightarrow
x_{\text{int}}\right)\simeq
C\frac{(4q)^{l}\pi}{\sin[2(l+1)\pi]}\left[\frac{x^{l+1}}{\Gamma(-l-\nu)\Gamma(2l+2)}-\frac{(4q)^{-(2l+1)}x^{-l}}{\Gamma(l+1-\nu)\Gamma(-2l)}\right]~{}.$
(48)

This also has a growing and a decaying mode. Therefore comparing the growing
and the decaying modes in C with the corresponding modes in 48, and
eliminating the arbitrary constant $C$, we arrive at the following relation
between $A$ and $B$,

$\displaystyle(4q)^{2l+1}\frac{\Gamma(-2l-1)\Gamma(-2l)\Gamma(l+1-\nu)}{\Gamma(2l+1)\Gamma(2l+2)\Gamma(-l-\nu)}=-\frac{\frac{\Gamma(1-4i\omega)}{\Gamma(l+1-2i\omega+\delta)\Gamma(l+1-2i\omega-\delta)}+\left(\frac{B}{A}\right)\frac{\Gamma(1+4i\omega)}{\Gamma(l+1+2i\omega+\delta)\Gamma(l+1+2i\omega-\delta)}}{\frac{\Gamma(1-4i\omega)}{\Gamma(-l-2i\omega+\delta)\Gamma(-l-\delta-2i\omega)}+\left(\frac{B}{A}\right)\frac{\Gamma(1+4i\omega)}{\Gamma(-l+2i\omega+\delta)\Gamma(-l+2i\omega-\delta)}}~{}.$
(49)

The above condition can be written in a more compact form as,

$\displaystyle\frac{\Gamma(-l-\nu)\Gamma(2l+2)}{\Gamma(-2l)\Gamma(l+1-\nu)}=-(4q)^{2l+1}\frac{\Gamma(-2l-1)}{\Gamma(2l+1)}\frac{\Gamma(l+1-2i\omega+\delta)\Gamma(l+1-2i\omega-\delta)}{\Gamma(-l-2i\omega+\delta)\Gamma(-l-\delta-2i\omega)}\frac{1}{\mathcal{G}}~{},$
(50)

where $\mathcal{G}$ has the following form,

$\displaystyle\mathcal{G}=\frac{1+\frac{B}{A}\frac{\Gamma(1+4i\omega)}{\Gamma(1-4i\omega)}\frac{\Gamma(l+1-2i\omega+\delta)}{\Gamma(l+1+2i\omega+\delta)}\frac{\Gamma(l+1-2i\omega-\delta)}{\Gamma(l+1+2i\omega-\delta)}}{1+\frac{B}{A}\frac{\Gamma(1+4i\omega)}{\Gamma(1-4i\omega)}\frac{\Gamma(-l-2i\omega+\delta)}{\Gamma(-l+2i\omega+\delta)}\frac{\Gamma(-l-2i\omega-\delta)}{\Gamma(-l+2i\omega-\delta)}}~{}.$
(51)

It is clear that the contribution of the reflective nature of the ECO surface
comes from the $\mathcal{G}$ term. If we had a BH instead, then $B=0$ and
hence we have $\mathcal{G}=1$. In this case the above expression matches with
the one of [72]. In order to simplify it further, we rewrite 50 as follows,

$\displaystyle\frac{\Gamma(-l-\nu)\Gamma(2l+2)}{\Gamma(-2l)\Gamma(l+1-\nu)}=-(4q)^{2l+1}\frac{F_{2}(\omega)}{F_{1}(\omega)}\underbrace{\left(\frac{1+\frac{B}{A}\frac{F_{2}(-\omega)}{F_{2}(\omega)}}{1+\frac{B}{A}\frac{F_{1}(-\omega)}{F_{1}(\omega)}}\right)}_{\mathcal{G}^{-1}}~{},$
(52)

where the functions $F_{1}(\omega)$ and $F_{2}(\omega)$ have the following
expressions,

$\displaystyle F_{1}$
$\displaystyle=\frac{\Gamma(1-4i\omega)\Gamma(2l+1)}{\Gamma(l+1-2i\omega+\delta)\Gamma(l+1-2i\omega-\delta)}~{},$
(53) $\displaystyle F_{2}$
$\displaystyle=\frac{\Gamma(1-4i\omega)\Gamma(-2l-1)}{\Gamma(-l-2i\omega+\delta)\Gamma(-l-2i\omega-\delta)}~{}.$
(54)

In order to derive the bound state frequencies analytically, we have to work
in the regime $\omega\sim\omega_{\text{pl}}\ll 1$. So we will keep terms up to
linear order in $\omega$, which implies $\delta\simeq 1$. In this case it
follows that,

$\displaystyle\frac{F_{1}(-\omega)}{F_{1}(\omega)}\simeq
1+4i\omega\big{\\{}2\Psi(1)-\Psi(l+2)-\Psi(l)\big{\\}}~{},$ (55)
$\displaystyle\frac{F_{2}(-\omega)}{F_{2}(\omega)}\simeq-1-4i\omega\big{\\{}2\Psi(1)-\Psi(l+2)-\Psi(l)\big{\\}}~{},$
(56)

where, $\Psi$ denotes the digamma function. Therefore, the function
$\mathcal{G}$ becomes,

$\displaystyle\mathcal{G}=\frac{1+\frac{B}{A}\left(1+4i\omega\big{\\{}2\Psi(1)-\Psi(l+2)-\Psi(l)\big{\\}}\right)}{1-\frac{B}{A}\left(1+4i\omega\big{\\{}2\Psi(1)-\Psi(l+2)-\Psi(l)\big{\\}}\right)}~{}.$
(57)

Moreover using the reflection formula for the gamma functions, one can show
that,

$\displaystyle\frac{F_{2}(\omega)}{F_{1}(\omega)}=\frac{\Gamma(l+2-2i\omega)\Gamma(l+2+2i\omega)\Gamma(l+2i\omega)\Gamma(l-2i\omega)}{\pi(2l)!(2l+1)!\sin(\pi+2\pi
l)\,\textrm{cosec}^{2}(l\pi+2i\omega\pi)}~{}.$ (58)

Note that, we have assumed $l\in\mathbb{C}$, and will take the
$l\to\textrm{integer}$ limit at the very end, akin to the analytic
continuation method. Similarly, using the reflection formula the left hand
side of 52 can be simplified, which yields,

$\displaystyle\frac{\Gamma(-l-\nu)\Gamma(2l+2)}{\Gamma(-2l)\Gamma(l+1-\nu)}=-\frac{(2l)!(2l+1)!}{\Gamma(l+1+\nu)\Gamma(l+1-\nu)}\frac{\sin[\pi(l+\nu)]}{\sin[2l\pi]}~{}.$
(59)

Substituting all of these results in 52, and also incorporating the fact that
$\sin(\pi+x)=-\sin x$, we arrive at the following expression,

$\begin{split}\frac{(2l)!(2l+1)!}{\Gamma(l+1+\nu)\Gamma(l+1-\nu)}\frac{\sin[\pi(l+\nu)]}{\sin(2l\pi)}=-(4q)^{2l+1}\frac{\Gamma(l+2-2i\omega)\Gamma(l+2+2i\omega)\Gamma(l+2i\omega)\Gamma(l-2i\omega)}{\pi(2l)!(2l+1)!\sin(2l\pi)~{}\text{cosec}^{2}(l\pi+2i\omega\pi)}\frac{1}{\mathcal{G}}~{}.\end{split}$
(60)

Further, using the identity
$|\Gamma(1+n+il)|^{2}=\frac{l\pi}{\text{sinh}(l\pi)}\prod_{k=1}^{n}(k^{2}+l^{2})$,
we finally arrive at the following expression,

$\displaystyle\frac{(-1)^{l}\sin(\pi\nu)}{\Gamma(l+1-\nu)\Gamma(l+1+\nu)}=\frac{(4q)^{2l+1}4\pi^{2}\omega^{2}}{\pi[(2l)!]^{2}[(2l+1)!]^{2}\mathcal{G}}\prod_{k=1}^{l+1}(k^{2}+4\omega^{2})\prod_{k=1}^{l-1}(k^{2}+4\omega^{2})~{}.$
(61)

From this expression it is clear that at leading order in the frequency
$\omega$, the left hand side of the above expression is finite, while the
right hand side is very small. To go around this we let, $\nu=\nu_{0}=(l+n+1)$
which yields,

$\displaystyle\omega\simeq\omega_{\text{pl}}\left(1-\frac{\omega_{\text{pl}}^{2}}{2(n+l+1)^{2}}\right)\equiv\omega_{0}$
(62)

This is the real part of the quasi-bound frequency and is independent of the
reflectivity of the ECO. This is why the real part of the quasi-bound state
frequencies effectively followed the BH trend. To compute the imaginary part
of the frequency, we consider $\omega=\omega_{0}+\delta\omega$. Then to the
leading order in $\delta\omega$ one can show that,
$\delta\nu=-(n+l+1)^{2}(\omega_{\text{pl}})^{-3}\delta\omega$. Under this
consideration, the left hand side of 60 becomes,

$\displaystyle\frac{(2l)!(2l+1)!}{\Gamma(l+1+\nu)\Gamma(l+1-\nu)}\frac{\sin[\pi(l+\nu)]}{\sin(2l\pi)}=\frac{(2l)!(2l+1)!n!}{(2l+1+n)!}\frac{(n+1+l)^{3}}{\mu^{3}}\delta\omega~{}.$
(63)

The right hand side of 60, on the other hand, yields, with all $\omega$ set to
$\omega_{0}$, the following result,

$\displaystyle\text{RHS}\simeq\frac{4^{2l+1}\omega_{\text{pl}}^{2l+2}}{(l+n+1)^{2l+1}}2i\omega_{0}\frac{\prod_{k=1}^{l-1}k^{2}\prod_{k=1}^{l+1}k^{2}}{(2l)!(2l+1)!}\left[\frac{1-\frac{B}{A}\left(1+4i\omega\\{2\Psi(1)-\Psi(l+2)-\Psi(l)\\}\right)}{1+\frac{B}{A}\left(1+4i\omega\\{2\Psi(1)-\Psi(l+2)-\Psi(l)\\}\right)}\right]~{},$
(64)

where, we have used $q\simeq-\\{\omega^{2}_{\text{pl}}/(n+1+l)\\}$. Using 29
and the location of the ECO surface, we can replace the ratio $(B/A)$ in terms
of reflectivity as follows,
$(B/A)=R(\omega)e^{-4i\omega\ln\epsilon}(1-4i\omega)$. With these results,
equating 63 and 64 we finally obtain,

$\displaystyle\delta\omega=2i\omega_{0}\frac{4^{2l+1\omega_{\text{pl}}^{4l+5}}}{(n+l+1)^{2l+4}}\frac{(2l+n+1)!}{n!}\left[\frac{(l+1)!(l-1)!}{(2l)!(2l+1)!}\right]^{2}\left[\frac{1-R(\omega_{0})e^{-4i\omega_{0}\ln\epsilon}\left(1+4i\omega_{0}\\{2\Psi(1)-\Psi(l+2)-\Psi(l)\\}\right)}{1+R(\omega_{0})e^{-4i\omega_{0}\ln\epsilon}\left(1+4i\omega_{0}\\{2\Psi(1)-\Psi(l+2)-\Psi(l)\\}\right)}\right]~{}.$
(65)

The oscillatory behaviour is due to the $e^{-4i\omega_{0}\ln\epsilon}$ term,
with the oscillation scale dependent on $\ln\epsilon$, and the amplitude on
$R(\omega_{0})$. Therefore, the features in the numerical spectrum is
recovered from our analytical expressions as well. The same can also be seen
from the plot of the imaginary part of $\delta\omega$ as well, presented in 7,
given the above analytic expression. It is to be emphasized that at large
values of $\omega_{\rm pl}$, the numerical and analytical results differ
considerably, which has to do with the failure of the assumption
$\omega\sim\omega_{\text{pl}}$.

Figure 7: The imaginary part of the fundamental plasma-driven axial EM quasi-
bound state frequencies has been plotted against the plasma frequency
$\omega_{\rm pl}$. We assume homogeneous plasma distribution as well as
constant reflectivity R=$\pm 0.5$. The oscillatory behaviour about the BH
value is clearly evident.

## References

  * [1] 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].
  * [2] Virgo, LIGO Scientific Collaboration, B. P. Abbott et al., “Properties of the binary black hole merger GW150914,” arXiv:1602.03840 [gr-qc].
  * [3] Virgo, LIGO Scientific Collaboration, B. P. Abbott et al., “Tests of general relativity with GW150914,” Phys. Rev. Lett. 116 no. 22, (2016) 221101, arXiv:1602.03841 [gr-qc].
  * [4] KAGRA, VIRGO, LIGO Scientific Collaboration, R. Abbott et al., “Open Data from the Third Observing Run of LIGO, Virgo, KAGRA, and GEO,” Astrophys. J. Suppl. 267 no. 2, (2023) 29, arXiv:2302.03676 [gr-qc].
  * [5] C. Cahillane and G. Mansell, “Review of the Advanced LIGO Gravitational Wave Observatories Leading to Observing Run Four,” Galaxies 10 no. 1, (2022) 36, arXiv:2202.00847 [gr-qc].
  * [6] Event Horizon Telescope Collaboration, K. Akiyama et al., “First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole,” Astrophys. J. 875 no. 1, (2019) L1, arXiv:1906.11238 [astro-ph.GA].
  * [7] Event Horizon Telescope Collaboration, “First sagittarius a* event horizon telescope results. i. the shadow of the supermassive black hole in the center of the milky way,” Astrophys. J. Lett. 930 no. 2, (May, 2022) L12.
  * [8] V. Cardoso, E. Franzin, and P. Pani, “Is the gravitational-wave ringdown a probe of the event horizon?,” Phys. Rev. Lett. 116 no. 17, (2016) 171101, arXiv:1602.07309 [gr-qc]. [Erratum: Phys.Rev.Lett. 117, 089902 (2016)].
  * [9] L. Barack et al., “Black holes, gravitational waves and fundamental physics: a roadmap,” Classical and Quantum Gravity 36 no. 14, (June, 2019) 143001.
  * [10] V. Cardoso and P. Pani, “Testing the nature of dark compact objects: a status report,” Living Reviews in Relativity 22 no. 1, (July, 2019) . http://dx.doi.org/10.1007/s41114-019-0020-4.
  * [11] H. A. Buchdahl, “General Relativistic Fluid Spheres,” Phys. Rev. 116 (1959) 1027.
  * [12] A. Alho, J. Natário, P. Pani, and G. Raposo, “Compactness bounds in General Relativity,” arXiv:2202.00043 [gr-qc].
  * [13] R. Goswami, S. D. Maharaj, and A. M. Nzioki, “Buchdahl-Bondi limit in modified gravity: Packing extra effective mass in relativistic compact stars,” Phys. Rev. D 92 (2015) 064002, arXiv:1506.04043 [gr-qc].
  * [14] S. Chakraborty and N. Dadhich, “Limits on stellar structures in Lovelock theories of gravity,” Phys. Dark Univ. 30 (2020) 100658, arXiv:2005.07504 [gr-qc].
  * [15] N. Dadhich and S. Chakraborty, “Buchdahl compactness limit for a pure Lovelock static fluid star,” Phys. Rev. D 95 no. 6, (2017) 064059, arXiv:1606.01330 [gr-qc].
  * [16] Event Horizon Telescope Collaboration, K. Akiyama et al., “First Sagittarius A* Event Horizon Telescope Results. V. Testing Astrophysical Models of the Galactic Center Black Hole,” Astrophys. J. Lett. 930 no. 2, (2022) L16, arXiv:2311.09478 [astro-ph.HE].
  * [17] Event Horizon Telescope Collaboration, K. Akiyama et al., “First Sagittarius A* Event Horizon Telescope Results. III. Imaging of the Galactic Center Supermassive Black Hole,” Astrophys. J. Lett. 930 no. 2, (2022) L14, arXiv:2311.09479 [astro-ph.HE].
  * [18] R. Carballo-Rubio, F. Di Filippo, S. Liberati, and M. Visser, “Constraints on horizonless objects after the EHT observation of Sagittarius A*,” JCAP 08 no. 08, (2022) 055, arXiv:2205.13555 [astro-ph.HE].
  * [19] V. Cardoso and P. Pani, “The observational evidence for horizons: from echoes to precision gravitational-wave physics,” arXiv:1707.03021 [gr-qc].
  * [20] S. H. Volkel and K. D. Kokkotas, “Ultra Compact Stars: Reconstructing the Perturbation Potential,” Class. Quant. Grav. 34 no. 17, (2017) 175015, arXiv:1704.07517 [gr-qc].
  * [21] N. Oshita, Q. Wang, and N. Afshordi, “On Reflectivity of Quantum Black Hole Horizons,” JCAP 04 (2020) 016, arXiv:1905.00464 [hep-th].
  * [22] P. Bueno, P. A. Cano, F. Goelen, T. Hertog, and B. Vercnocke, “Echoes of Kerr-like wormholes,” Phys. Rev. D97 no. 2, (2018) 024040, arXiv:1711.00391 [gr-qc].
  * [23] P. V. P. Cunha, C. A. R. Herdeiro, and M. J. Rodriguez, “Does the black hole shadow probe the event horizon geometry?,” Phys. Rev. D97 no. 8, (2018) 084020, arXiv:1802.02675 [gr-qc].
  * [24] R. Dey, S. Chakraborty, and N. Afshordi, “Echoes from braneworld black holes,” Phys. Rev. D 101 no. 10, (2020) 104014, arXiv:2001.01301 [gr-qc].
  * [25] E. Maggio, V. Cardoso, S. R. Dolan, and P. Pani, “Ergoregion instability of exotic compact objects: electromagnetic and gravitational perturbations and the role of absorption,” Phys. Rev. D 99 no. 6, (2019) 064007, arXiv:1807.08840 [gr-qc].
  * [26] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, “Gravitational-wave signatures of exotic compact objects and of quantum corrections at the horizon scale,” Phys. Rev. D 94 no. 8, (2016) 084031, arXiv:1608.08637 [gr-qc].
  * [27] Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, “A recipe for echoes from exotic compact objects,” Phys. Rev. D 96 no. 8, (2017) 084002, arXiv:1706.06155 [gr-qc].
  * [28] K. W. Tsang, M. Rollier, A. Ghosh, A. Samajdar, M. Agathos, K. Chatziioannou, V. Cardoso, G. Khanna, and C. Van Den Broeck, “A morphology-independent data analysis method for detecting and characterizing gravitational wave echoes,” Phys. Rev. D98 no. 2, (2018) 024023, arXiv:1804.04877 [gr-qc].
  * [29] N. Oshita and N. Afshordi, “Probing microstructure of black hole spacetimes with gravitational wave echoes,” Phys. Rev. D99 no. 4, (2019) 044002, arXiv:1807.10287 [gr-qc].
  * [30] R. A. Konoplya, Z. Stuchlik, and A. Zhidenko, “Echoes of compact objects: new physics near the surface and matter at a distance,” Phys. Rev. D99 no. 2, (2019) 024007, arXiv:1810.01295 [gr-qc].
  * [31] R. Dey, S. Biswas, and S. Chakraborty, “Ergoregion instability and echoes for braneworld black holes: Scalar, electromagnetic, and gravitational perturbations,” Phys. Rev. D 103 no. 8, (2021) 084019, arXiv:2010.07966 [gr-qc].
  * [32] E. Maggio, A. Testa, S. Bhagwat, and P. Pani, “Analytical model for gravitational-wave echoes from spinning remnants,” Phys. Rev. D 100 no. 6, (2019) 064056, arXiv:1907.03091 [gr-qc].
  * [33] S. Biswas, M. Rahman, and S. Chakraborty, “Echoes from braneworld wormholes,” Phys. Rev. D 106 no. 12, (2022) 124003, arXiv:2205.14743 [gr-qc].
  * [34] S. Chakraborty, E. Maggio, A. Mazumdar, and P. Pani, “Implications of the quantum nature of the black hole horizon on the gravitational-wave ringdown,” Phys. Rev. D 106 no. 2, (2022) 024041, arXiv:2202.09111 [gr-qc].
  * [35] S. Biswas, C. Singha, and S. Chakraborty, “Galactic wormholes: Geometry, stability, and echoes,” Phys. Rev. D 109 no. 6, (2024) 064043, arXiv:2307.04836 [gr-qc].
  * [36] A. Chowdhury and N. Banerjee, “Echoes from a singularity,” Phys. Rev. D 102 no. 12, (2020) 124051, arXiv:2006.16522 [gr-qc].
  * [37] P. Pani, “I-Love-Q relations for gravastars and the approach to the black-hole limit,” Phys. Rev. D92 no. 12, (2015) 124030, arXiv:1506.06050 [gr-qc].
  * [38] V. Cardoso, E. Franzin, A. Maselli, P. Pani, and G. Raposo, “Testing strong-field gravity with tidal Love numbers,” Phys. Rev. D95 no. 8, (2017) 084014, arXiv:1701.01116 [gr-qc].
  * [39] S. Nair, S. Chakraborty, and S. Sarkar, “Dynamical Love numbers for area quantized black holes,” Phys. Rev. D 107 no. 12, (2023) 124041, arXiv:2208.06235 [gr-qc].
  * [40] S. Chakraborty, E. Maggio, M. Silvestrini, and P. Pani, “Dynamical tidal Love numbers of Kerr-like compact objects,” arXiv:2310.06023 [gr-qc].
  * [41] S. Datta, R. Brito, S. Bose, P. Pani, and S. A. Hughes, “Tidal heating as a discriminator for horizons in extreme mass ratio inspirals,” Phys. Rev. D 101 no. 4, (2020) 044004, arXiv:1910.07841 [gr-qc].
  * [42] S. Chakraborty, S. Datta, and S. Sau, “Tidal heating of black holes and exotic compact objects on the brane,” Phys. Rev. D 104 no. 10, (2021) 104001, arXiv:2103.12430 [gr-qc].
  * [43] S. Datta, “Tidal heating of Quantum Black Holes and their imprints on gravitational waves,” Phys. Rev. D 102 no. 6, (2020) 064040, arXiv:2002.04480 [gr-qc].
  * [44] S. Mitra, S. Chakraborty, R. Vicente, and J. C. Feng, “Probing the quantum nature of black holes with ultra-light boson environments,” arXiv:2312.06783 [gr-qc].
  * [45] E. Cannizzaro, A. Caputo, L. Sberna, and P. Pani, “Plasma-photon interaction in curved spacetime I: formalism and quasibound states around nonspinning black holes,” Phys. Rev. D 103 (2021) 124018, arXiv:2012.05114 [gr-qc].
  * [46] A. Sitenko, Electromagnetic Fluctuations in Plasma.  Academic Press, 1967.
  * [47] J. G. Rosa and S. R. Dolan, “Massive vector fields on the schwarzschild spacetime: Quasinormal modes and bound states,” Physical Review D 85 no. 4, (Feb., 2012) . http://dx.doi.org/10.1103/PhysRevD.85.044043.
  * [48] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, “Black-hole bombs and photon-mass bounds,” Physical Review Letters 109 no. 13, (Sept., 2012) . http://dx.doi.org/10.1103/PhysRevLett.109.131102.
  * [49] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, “Perturbations of slowly rotating black holes: Massive vector fields in the kerr metric,” Physical Review D 86 no. 10, (Nov., 2012) . http://dx.doi.org/10.1103/PhysRevD.86.104017.
  * [50] P. Pani and A. Loeb, “Constraining primordial black-hole bombs through spectral distortions of the cosmic microwave background,” Physical Review D 88 no. 4, (Aug., 2013) . http://dx.doi.org/10.1103/PhysRevD.88.041301.
  * [51] M. Baryakhtar, R. Lasenby, and M. Teo, “Black hole superradiance signatures of ultralight vectors,” Physical Review D 96 no. 3, (Aug., 2017) . http://dx.doi.org/10.1103/PhysRevD.96.035019.
  * [52] J. P. Conlon and C. A. R. Herdeiro, “Can black hole superradiance be induced by galactic plasmas?,” Phys. Lett. B 780 (2018) 169–173, arXiv:1701.02034 [astro-ph.HE].
  * [53] V. Cardoso, O. J. Dias, G. S. Hartnett, M. Middleton, P. Pani, and J. E. Santos, “Constraining the mass of dark photons and axion-like particles through black-hole superradiance,” Journal of Cosmology and Astroparticle Physics 2018 no. 03, (Mar., 2018) 043–043. http://dx.doi.org/10.1088/1475-7516/2018/03/043.
  * [54] V. P. Frolov, P. Krtouš, D. Kubizňák, and J. E. Santos, “Massive vector fields in rotating black-hole spacetimes: Separability and quasinormal modes,” Physical Review Letters 120 no. 23, (June, 2018) . http://dx.doi.org/10.1103/PhysRevLett.120.231103.
  * [55] S. R. Dolan, “Instability of the proca field on kerr spacetime,” Physical Review D 98 no. 10, (Nov., 2018) . http://dx.doi.org/10.1103/PhysRevD.98.104006.
  * [56] D. Baumann, H. S. Chia, J. Stout, and L. t. Haar, “The spectra of gravitational atoms,” Journal of Cosmology and Astroparticle Physics 2019 no. 12, (Dec., 2019) 006–006. http://dx.doi.org/10.1088/1475-7516/2019/12/006.
  * [57] R. Brito, V. Cardoso, and P. Pani, Superradiance: Energy Extraction, Black-Hole Bombs and Implications for Astrophysics and Particle Physics.  Springer International Publishing, 2015.  http://dx.doi.org/10.1007/978-3-319-19000-6.
  * [58] E. Maggio, P. Pani, and V. Ferrari, “Exotic Compact Objects and How to Quench their Ergoregion Instability,” Phys. Rev. D96 no. 10, (2017) 104047, arXiv:1703.03696 [gr-qc].
  * [59] V. Cardoso, W.-D. Guo, C. F. B. Macedo, and P. Pani, “The tune of the Universe: the role of plasma in tests of strong-field gravity,” Mon. Not. Roy. Astron. Soc. 503 no. 1, (2021) 563–573, arXiv:2009.07287 [gr-qc].
  * [60] E. Cannizzaro, A. Caputo, L. Sberna, and P. Pani, “Plasma-photon interaction in curved spacetime. ii. collisions, thermal corrections, and superradiant instabilities,” Phys. Rev. D 104 (Nov, 2021) 104048. https://link.aps.org/doi/10.1103/PhysRevD.104.104048.
  * [61] E. Cannizzaro, Searching for new physics in the neighborhood of a black hole: fundamental interactions, plasmas and tests of gravity.  PhD thesis, Rome U., U. Rome La Sapienza (main), 2024.
  * [62] E. Cannizzaro, F. Corelli, and P. Pani, “Nonlinear photon-plasma interaction and the black hole superradiant instability,” Phys. Rev. D 109 (Jan, 2024) 023007. https://link.aps.org/doi/10.1103/PhysRevD.109.023007.
  * [63] T. F. M. Spieksma, E. Cannizzaro, T. Ikeda, V. Cardoso, and Y. Chen, “Superradiance: Axionic couplings and plasma effects,” Phys. Rev. D 108 (Sep, 2023) 063013. https://link.aps.org/doi/10.1103/PhysRevD.108.063013.
  * [64] E. Cannizzaro, A. Caputo, L. Sberna, and P. Pani, “Plasma-photon interaction in curved spacetime I: formalism and quasibound states around nonspinning black holes,” Phys. Rev. D 103 (2021) 124018, arXiv:2012.05114 [gr-qc].
  * [65] E. Cannizzaro, A. Caputo, L. Sberna, and P. Pani, “Plasma-photon interaction in curved spacetime. II. Collisions, thermal corrections, and superradiant instabilities,” Phys. Rev. D 104 no. 10, (2021) 104048, arXiv:2107.01174 [gr-qc].
  * [66] A. Caputo, S. J. Witte, D. Blas, and P. Pani, “Electromagnetic signatures of dark photon superradiance,” Phys. Rev. D 104 no. 4, (2021) 043006, arXiv:2102.11280 [hep-ph].
  * [67] G. Alonso-Álvarez and J. M. Cline, “Sterile neutrino dark matter catalyzed by a very light dark photon,” JCAP 10 (2021) 041, arXiv:2107.07524 [hep-ph].
  * [68] R. A. Breuer and J. Ehlers, “Propagation of electromagnetic waves through magnetized plasmas in arbitrary gravitational fields,” Astron. Astrophys. 96 no. 1-2, (Mar., 1981) 293–295.
  * [69] G. F. R. Ellis, “Relativistic cosmology,” Proc. Int. Sch. Phys. Fermi 47 (1971) 104–182.
  * [70] L. Zhou, R. Brito, Z.-F. Mai, and L. Shao, “Superradiant instabilities of massive bosons around exotic compact objects,” Phys. Rev. D 108 (Nov, 2023) 103025. https://link.aps.org/doi/10.1103/PhysRevD.108.103025.
  * [71] T. F. M. Spieksma, E. Cannizzaro, T. Ikeda, V. Cardoso, and Y. Chen, “Superradiance: Axionic couplings and plasma effects,” Phys. Rev. D 108 no. 6, (2023) 063013, arXiv:2306.16447 [gr-qc].
  * [72] J. G. Rosa and S. R. Dolan, “Massive vector fields on the Schwarzschild spacetime: quasi-normal modes and bound states,” Phys. Rev. D 85 (2012) 044043, arXiv:1110.4494 [hep-th].