# Plasma-photon interaction in curved spacetime. I:
Formalism and quasibound states around nonspinning black holes

Enrico Cannizzaro  Dipartimento di Fisica, “Sapienza” Università di Roma &
Sezione INFN Roma1, Piazzale Aldo Moro 5, 00185, Roma, Italy  Andrea Caputo
School of Physics and Astronomy, Tel-Aviv University, Tel-Aviv 69978, Israel
Department of Particle Physics and Astrophysics, Weizmann Institute of
Science, Rehovot 7610001, Israel  Max-Planck-Institut für Physik (Werner-
Heisenberg-Institut), Föhringer Ring 6, 80805 München, Germany  Laura Sberna
Max Planck Institute for Gravitational Physics (Albert Einstein Institute) Am
Muühlenberg 1, 14476 Potsdam, Germany  Perimeter Institute, 31 Caroline St N,
Waterloo, ON N2L 2Y5, Canada  Paolo Pani  Dipartimento di Fisica, “Sapienza”
Università di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 5, 00185, Roma,
Italy

###### Abstract

We investigate the linear dynamics of an electromagnetic field propagating in
curved spacetime in the presence of plasma. The dynamical equations are
generically more involved and richer than the effective Proca equation adopted
as a model in previous work. We discuss the general equations and focus on the
case of a cold plasma in the background of a spherically symmetric black hole,
showing that the system admits plasma-driven, quasibound electromagnetic
states that are prone to become superradiantly unstable when the black hole
rotates. The quasibound states are different from those of the Proca equation
and have some similarities with the case of a massive scalar field, suggesting
that the linear instability can be strongly suppressed compared to previous
estimates. Our framework provides the first step toward a full understanding
of the plasma-photon interactions around astrophysical black holes.

##  I Introduction

###  I.1 Motivation

The propagation of fundamental fields on black hole (BH) spacetimes is a
decades-long fascinating subject, which has received considerable attention in
the last years for its relevance for tests of fundamental physics with BH
observations Barack:2018yly . A particularly relevant case is the one of an
electromagnetic (EM) field propagating in a BH spacetime in the presence of
plasma. The study of plasmas is crucial in many branches of physics and
astrophysics, including the study of BHs, for at least two reasons: (i)
astrophysical BHs can be surrounded by accretion disks due to the outward
transfer of angular momentum of accreting matter 1973A&A….24..337S ;
Novikov:1973kta ; (ii) photons emitted by accreting BHs can interact with the
plasma of the interstellar medium, potentially affecting their propagation at
large distances.

It is common lore to assume that a photon propagating in a plasma is simply
dressed with an effective mass, induced by the interaction with the background
medium. A photon in a plasma is therefore commonly described with a Proca-like
equation $\nabla_{\sigma}F^{\sigma\nu}=\omega_{\rm pl}^{2}A^{\nu}$, where
$F_{\sigma\nu}=\partial_{\sigma}A_{\nu}-\partial_{\nu}A_{\sigma}$ is the EM
field strength and111Throughout this paper, we use $G=c=1$ and rationalized
Heaviside units for the Maxwell equations. The metric signature is mostly
positive.

$\omega_{\rm pl}=(ne^{2}/m_{e})^{1/2}$ (1)

is the plasma frequency of the background medium Sitenko:1967 , with $n$ the
particle density of electrons with charge $e$ and mass $m_{e}$. The Proca
equation has been studied in detail on both Schwarzschild and Kerr metrics
Rosa:2011my ; Pani:2012vp ; Pani:2012bp ; Baryakhtar:2017ngi ; Cardoso:2018tly
; Frolov:2018ezx ; Dolan:2018dqv ; Baumann:2019eav . In the case of a rotating
BH, the system can develop a superradiant instability, similar to the one
exhibited by massive scalar fields (see Ref. Brito:2015oca  for a review).
This sows the seeds for striking observational signatures, which are
particularly significant for primordial BHs dark matter Pani:2013hpa  or in
the presence of extra degrees of freedom such as axions Blas:2020nbs . This
mechanism has been also advocated as a possible explanation for the origin of
fast radio bursts Conlon:2017hhi .

Given the interesting phenomenological implications of plasma-photon
interactions around BHs, it is of crucial importance to understand first of
all whether or not the Proca equation represents a good description of the
dynamics of the system. The fact that the Proca equation can at most be an
approximation is already clear from the counting of the degrees of freedom: a
massive spin-1 field propagates three polarizations (two transverse modes and
a longitudinal one), whereas photons in a cold plasma propagate only two
transverse modes, since the putative longitudinal mode is electrostatic, see
details below.

In this work we take the first step toward a detailed investigation of the
propagation of EM waves in a plasma in curved spacetime, and show that even in
the case of a Schwarzschild BH – where of course superradiance cannot develop
– the spectrum of the quasibound states deviates from the one of a simple
Proca equation. This will impact the plasma-driven superradiant instability,
since the quasibound states around a static BH are those prone to become
superradiantly unstable when the BH rotates Pani:2012bp ; Pani:2012vp ;
Brito:2015oca . Our results call for further investigation, and are likely
pointing toward a strong suppression of the expected superradiant instability
in rotating BH spacetimes.

Other plasma effects have been pointed out recently and indicate that plasma
may generically hinder our ability to perform strong-field gravity tests, for
example by quenching the superradiance instability Cardoso:2020nst ;
Blas:2020kaa . The present work shows that the suppression of the instability
will be even more severe (and already present at the linear level) when one
accounts precisely for the presence of the plasma in the dynamical equations.

In principle, the above problems could be studied in the framework of general-
relativity-magneto-hydrodynamics (GRMHD) (see, e.g., Refs. Abramowicz:2011xu ;
Yuan:2014gma  for a review). However, GRMHD simulations typically focus on
highly-magnetized plasmas in the absence of an electric field, and are also
computationally expensive, which makes them unsuitable to follow the
(potentially long) evolution of the microscopic EM field on different
backgrounds. Indeed, to the best of our knowledge, a detailed understanding of
the propagation of EM waves in a plasma in curved spacetime problem was not
considered in previous work.

###  I.2 Cold and hot plasma

A plasma is a gas composed of electrons and ions (ionized by heating or
photoionization) which obey the classical Maxwell-Boltzmann rather than Fermi-
Dirac or Bose-Einstein statistics Sitenko:1967 . A plasma can be further
classified in terms of its temperature. In a cold plasma the thermal velocity
of the electrons is larger than that of ions, but still much smaller than the
speed of an EM wave propagating in it, namely

$v_{\rm thermal}\equiv\sqrt{\frac{2T_{e}}{m_{e}}}\ll\frac{\omega}{k}\,,$ (2)

where $\omega$ and $k$ are the frequency and wave number of the EM field, and
$T_{e}$ is the temperature of the electrons. In this case the thermal pressure
is negligible. In a cold plasma, the longitudinal modes are plasma
oscillations which do not propagate and do not transport energy. The presence
of this electrostatic mode already shows that the Proca equation cannot fully
describe the interaction of photons with a cold plasma.

On the other hand, in a hot plasma, the electron thermal velocity cannot be
neglected and the electrons cannot be considered at rest with respect to
propagating waves. In this case the electrostatic modes are converted into
propagating, energy-transporting, longitudinal modes called Langmuir waves
2017mcp..book…..T . These modes propagate at the speed of sound in the plasma
and the (nonrelativistic) equations resemble the Proca one in the case in
which the electrons are ultrarelativistic PhysRevD.45.525 ; PhysRevE.62.2989
(although this approximation is beyond the regime of validity of the original
equations).

In the following we will limit ourselves to the study of cold plasma, since
the latter provides a good description of accretion disks around BHs
1974ApJ…191..499P ; 1973A&A….24..337S . Indeed, in the inner region of a
typical accretion disk the temperature can be estimated as Abramowicz:2011xu

$T_{\rm disk}\simeq 4\times
10^{3}\alpha^{-1/4}\Big{(}\frac{M_{\odot}}{M}\Big{)}^{1/4}\Big{(}\frac{r}{M}\Big{)}^{-3/8}\text{eV}\ll
m_{e},$ (3)

where $\alpha\sim{\cal O}(1)$ is a dimensionless coefficient relating the
kinematic viscosity of the fluid with the velocity of turbulent elements.
Therefore, $v_{\rm thermal}\approx 0.06$ or smaller. On the other hand, for a
quasibound state around a BH of mass $M$, $\omega\sim\omega_{\rm pl}$ and
$k=2\pi/\lambda$, where $\lambda\sim M/(M\omega_{\rm pl})^{2}$ is the typical
length scale of the mode Arvanitaki:2010sy ; Brito:2014wla . This gives
$\omega/k\sim 0.3$ or larger for the most interesting case $M\omega_{\rm
pl}\lesssim 0.5$. We also stress that in the present work we treat the plasma
as spherically symmetric, and static. The first is a simplifying assumption
and should be relaxed to accommodate more realistic accretion disk geometries.
Once spherical symmetry is assumed, the plasma is static to very good
approximation: the time scales of interest to this work are much shorter than
the time scales of other important astrophysical phenomena. In particular, the
BH accretion timescale can be conservatively estimated to be given by a
fraction of the Eddington accretion timescale $\tau_{\rm accr}=M/\dot{M}\sim
f_{\rm Edd}^{-1}\,10^{15}\,{\rm s}$, where $f_{\rm Edd}$ is an accretion
efficiency factor. For $f_{\rm Edd}\sim\mathcal{O}(1)$, which is conservative
as accretion may be much less efficient, one sees that the time scale of the
plasma radial motion is much longer than the time scale of the quasi-trapped
long-lived perturbations studied below, $\tau_{\rm accr}\gg\omega_{\rm
pl}^{-1}$.

The remaining of this manuscript is organized as follows. In Sec. II we define
the general setup of our work, the main equations in spherical symmetry (Sec.
II.2) and the plasma profiles we considered (Sec. II.4). We show our numerical
results for the quasibound states in Sec. III and comment on the implications
for spinning BHs in Sec. IV. Finally, we conclude in Sec. V.

##  II Setup

###  II.1 General equations

The general study of the propagation of EM waves through a cold plasma in
curved spacetime was pionereed in 1981A&A….96..293B , where the authors
derived the system of nonlinear equations governing the plasma and the EM
field. In this work we start from those equations and specialize to a
background Schwarzschild metric, for which we study the quasibound states for
different plasma configurations.

Consider a two-component plasma made of electrons and ions. Let us denote the
number density and four- velocity of the electrons as $n$ and $u^{\mu}$, while
$J^{\mu}$ stands for the ion current density. The system of differential
equations for the plasma quantities reads 1981A&A….96..293B

$\displaystyle\nabla_{\nu}F^{\mu\nu}=enu^{\mu}+J^{\mu},$ (4) $\displaystyle
u^{\mu}\nabla_{\mu}u^{\nu}=e/m_{e}F^{\nu}{}_{\mu}u^{\mu},$ (5) $\displaystyle
u^{\mu}u_{\mu}=-1,$ (6) $\displaystyle\nabla_{\mu}(nu^{\mu})=0.$ (7)

These are Maxwell’s equations together with the momentum and particle
conservation equations, in covariant form.

We study the propagation of a perturbation through the plasma by introducing
the small perturbations $\tilde{n},\tilde{u}^{\mu},\tilde{F}_{\mu\nu}$, e.g.
${F}_{\mu\nu}={F}_{\mu\nu}^{\rm background}+\tilde{F}_{\mu\nu}$ (and likewise
for other quantities). Here we neglect second-order perturbations of the
plasma and EM field, as well as any perturbation of the background metric
$g_{\mu\nu}$ (since the gravitational backreactions of these fields is small).
We also neglect perturbations of the ions, since they will be suppressed with
respect to those of the electrons by a factor $\propto m_{e}/m_{\rm ion}\ll
1$. In the hot plasma case an extra term would appear in the momentum equation
(5), due to the pressure of the fluid.

The presence of a plasma implies the existence of a preferred rest frame.
Locally, the plasma defines surfaces of simultaneity for the observer, whose
effective metric tensor is

$h_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}.$ (8)

The tensor $h_{\mu\nu}$ projects orthogonally onto the tangent rest plane of
the (electron) plasma. Then, the kinematic of an electron fluid is described
by two matrices: the rate of rotation (the vorticity)
$\omega^{\mu\nu}\equiv-\omega^{\nu\mu}$, and the rate of deformation
$\theta^{\mu\nu}\equiv\theta^{\nu\mu}$. This follows from the general
decomposition Ellis:1971pg

$\nabla_{\mu}u_{\nu}=\nabla{}_{(\mu}u_{\nu)}+\nabla{}_{[\mu}u_{\nu]}={\omega_{\nu}}_{\mu}+\theta_{\nu\mu}-u_{\mu}u^{\alpha}\nabla_{\alpha}u_{\nu}\,,$
(9)

from which we get

$\displaystyle\omega_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(v_{\mu\nu}-v_{\nu\mu})\,,$ (10)
$\displaystyle\theta_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(v_{\mu\nu}+v_{\nu\mu})\,,$ (11)

where we defined the tensor
$v^{\mu\nu}=h^{\mu\alpha}h^{\nu\beta}u_{\alpha;\beta}$. We can also define the
plasma frequency as in Eq. (1), the electric component $E^{\mu}\equiv
F^{\mu}{}_{\nu}u^{\nu}$, the magnetic component $B_{\mu\nu}\equiv
h_{\mu}{}^{\alpha}h_{\nu}{}^{\beta}F_{\alpha\beta}$, and the Larmor tensor
${\omega_{\rm L}}^{\mu\nu}=-\frac{e}{m_{e}}B^{\mu\nu}$. With these
definitions, by differentiating Maxwell’s equation (4) and using the momentum
equation (5), Ref. 1981A&A….96..293B  obtained the perturbed equation for the
vector potential perturbation $\tilde{A}^{\mu}$ in the Landau gauge,
$u_{\mu}\tilde{A}^{\mu}=0$, containing both the influence of the gravitational
potential and that of the moving plasma:

$\displaystyle\Big{(}h^{\alpha\beta}u^{\nu}\nabla_{\nu}(\nabla^{\mu}{}_{\beta}-\delta^{\mu}{}_{\beta}\nabla^{\gamma}{}_{\gamma})+(\omega^{\alpha\beta}+\omega_{\rm
L}^{\alpha\beta}+\theta^{\alpha\beta}$ (12) $\displaystyle+\theta
h^{\alpha\beta}+\frac{e}{m_{e}}E^{\alpha}u^{\beta})\times(\nabla^{\mu}{}_{\beta}-\delta^{\mu}{}_{\beta}\nabla^{\gamma}{}_{\gamma})$
$\displaystyle+\omega_{\rm
pl}^{2}h^{\alpha\mu}u^{\gamma}\nabla_{\gamma}+\omega_{\rm
pl}^{2}(\theta^{\alpha\mu}-\omega^{\alpha\mu})$
$\displaystyle\Big{)}\tilde{A}_{\mu}=0,$

where $\nabla_{\mu\nu}\equiv\nabla_{\mu}\nabla_{\nu}$ and
$\theta=\theta^{\mu}{}_{\mu}$. The above equation is the starting point for a
rigorous analysis of the linearized photon dynamics in a cold plasma in curved
spacetime. It is clearly very different from an effective Proca equation which
would have the form $\nabla_{\alpha}\nabla^{\alpha}\tilde{A}_{\mu}=\omega_{\rm
pl}^{2}\tilde{A}_{\mu}$. In particular, note that Eq. (12) contains third-
order derivatives.

In the flat spacetime limit and in the Fourier domain, the spatial part of Eq.
(12) reads

$\bm{k}(\bm{k}\cdot\tilde{\bm{E}})-k^{2}\tilde{\bm{E}}+\Big{(}1-\frac{\omega_{\rm
pl}^{2}}{\omega^{2}}\Big{)}\omega^{2}\tilde{\bm{E}}=0\,,$ (13)

which is indeed the equation that regulates the propagation of an EM wave in
an isotropic plasma with dielectric tensor $\epsilon=1-\omega_{\rm
pl}^{2}/\omega^{2}$. Taking $\bm{k}$ along the $z$ direction one can write the
standard dispersion tensor $\bm{D}$ and find its determinant, which reads
$\text{det}|\bm{D}|=(\omega^{2}\epsilon-k^{2})^{2}\omega^{2}\epsilon=0$. There
are two types of solutions to this dispersion relation: the first one
corresponds to the longitudinal mode with $\epsilon=0$, i.e.
$\omega=\omega_{\rm pl}$; while the second one corresponds to the two
transverse modes $(\omega^{2}\epsilon-k^{2})=0$, i.e.
$\omega^{2}=k^{2}+\omega_{\rm pl}^{2}$ Raffelt:1996wa .

###  II.2 Spherically symmetric spacetimes

We now specialize to the symmetries of the Schwarzschild background. We work
in the coordinates $(t,r,\theta,\phi)$, in which the line element reads

$ds^{2}=-fdt^{2}+f^{-1}dr^{2}+r^{2}d\Omega_{2}^{2}$ (14)

with $f(r)=1-2M/r$, where $M$ is the BH mass. In this case both the vorticity
and the deformation tensors are zero, as can be easily checked from Eqs. (10)
and (11). The four velocity of a static plasma is
$u^{\alpha}=(u^{0},\vec{0})$, with $u^{0}=f^{-1/2}$ satisfying the
normalization condition $u_{\mu}u^{\mu}=-1$. From Eq. (5), the electric field
has then only one nonvanishing radial component
$E^{\alpha}=(0,m_{e}/e\,\Gamma^{r}_{00}(u^{0})^{2},0,0)$, where
$\Gamma^{\mu}_{\alpha\beta}$ are the standard Christoffel’s symbols. We assume
an unmagnetized plasma $B_{\mu\nu}=0$ (and therefore also ${\omega_{\rm
L}}^{\mu\nu}=0$).

Moreover, in any spherically symmetric spacetime it is possible to separate
the angular part of the fields from the radial one by performing a multipolar
expansion. Following Ref. Rosa:2011my , we introduce a basis of four vector
spherical harmonics:

$\displaystyle Z_{\mu}^{(1)lm}$ $\displaystyle=[1,0,0,0]Y^{lm},$ (15)
$\displaystyle Z_{\mu}^{(2)lm}$ $\displaystyle=[0,f^{-1},0,0]Y^{lm},$ (16)
$\displaystyle Z_{\mu}^{(3)lm}$
$\displaystyle=\frac{r}{\sqrt{l(l+1)}}[0,0,\partial_{\theta},\partial_{\phi}]Y^{lm},$
(17) $\displaystyle Z_{\mu}^{(4)lm}$
$\displaystyle=\frac{r}{\sqrt{l(l+1)}}[0,0,\frac{\partial_{\phi}}{\sin\theta},-\sin\theta\partial_{\theta}]Y^{lm},$
(18)

where $Y^{lm}(\theta,\phi)$ are the standard scalar spherical harmonics. These
vector spherical harmonics satisfy the orthogonality condition

$\int d\Omega
Z^{(i)lm}_{\mu}\hat{\eta}^{\mu\nu}Z^{(i^{\prime})l^{\prime}m^{\prime}}_{\nu}=\delta^{ii^{\prime}}\delta^{ll^{\prime}}\delta^{mm^{\prime}},$
(19)

where $d\Omega=\sin\theta d\theta d\phi$ and $\hat{\eta}^{\mu\nu}={\rm
diag}[1,f^{2},1/r^{2},1/(r^{2}\sin^{2}\theta)]$. The perturbation of the
vector potential can be decomposed in this basis as

$\tilde{A}_{\mu}(r,t,\theta,\phi)=\frac{1}{r}\sum_{i=1}^{4}\sum_{l,m}c_{i}u_{(i)}^{lm}(t,r)Z_{\mu}^{(i)lm}(\theta,\phi),$
(20)

where $c_{1}=c_{2}=1$ and $c_{3}=c_{4}=1/\sqrt{l(l+1)}$. Using this
decomposition and a frequency-domain representation
$u_{(i)}^{lm}(t,r)=u_{(i)}^{lm}(r)e^{-i\omega t}$, the field equations become

$\displaystyle u_{(1)}$ $\displaystyle=0\,,$ (21)
$\displaystyle\left(rf(l+l^{2}+r^{2}\omega_{\rm
pl}^{2})-r^{3}\omega^{2}\right)u_{(2)}-r^{2}f^{2}u_{(3)}^{\prime}$
$\displaystyle=0\,,$ (22) $\displaystyle
l(1+l)rfu_{(2)}+r^{3}(\omega^{2}-f\omega_{\rm pl}^{2})u_{(3)}$
$\displaystyle-l(1+l)r^{2}fu_{(2)}^{\prime}+2Mrfu_{(3)}^{\prime}+r^{3}f^{2}u_{(3)}^{\prime\prime}$
$\displaystyle=0\,,$ (23) $\displaystyle\left(rf(l+l^{2}+r^{2}\omega_{\rm
pl}^{2})+r^{3}\omega^{2}\right)u_{(4)}$
$\displaystyle-2Mrfu_{(4)}^{\prime}-r^{3}f^{2}u_{(4)}^{\prime\prime}$
$\displaystyle=0\,,$ (24)

where $u_{(i)}^{\prime}=\partial_{r}u_{(i)}$, we have suppressed the $l$
superscript, and the radial dependence of $\omega_{\rm pl}=\omega_{\rm
pl}(r)$. Owing to the spherical symmetry of the background, the equations do
not depend on the angular number $m$. Note that, despite the fact that Eq.
(12) contains third-order derivatives, the final system of equations in the
frequency domain is of second differential order.

From Eqs. (21)–(24) we can immediately notice that the polar (i.e., even-
parity) sector, described by the functions $u_{(1)},u_{(2)}$ and $u_{(3)}$ is
completely decoupled from the axial (i.e., odd-parity) sector, described by
the function $u_{(4)}$. This resembles the case of a massive vector field
discussed in Ref. Rosa:2011my  and is in fact a consequence of the spherical
symmetry of the background. We can therefore treat the two sectors
independently, as in the following.

####  II.2.1 Axial sector

Interestingly, the axial equation is identical to the axial Proca equation
found in Rosa:2011my . Thus, the axial sector can be easily reduced to a
Schrödinger-like equation analogous to the massive vector case Rosa:2011my

$\mathcal{D}_{2}u_{(4)}(r)=0,$ (25)

with
$\mathcal{D}_{2}\equiv\frac{d^{2}}{dr_{*}^{2}}+\omega^{2}-f\Big{(}\frac{l(l+1)}{r^{2}}+\omega_{\rm
pl}^{2}\Big{)}$, and in terms of the tortoise coordinate defined by
$\frac{dr_{*}}{dr}=f^{-1}$. In this case the plasma frequency plays indeed the
role of an effective mass for the component $u_{(4)}$.

####  II.2.2 Polar sector

On the other hand, the polar sector is different from the Proca case.

In the monopole case ($l=0$), only the first two spherical harmonics are
defined and the field equations reduce to $u_{(1)}=0$ and

$(\omega^{2}-\omega_{\rm pl}^{2}f)u_{(2)}=0.$ (26)

In general, this equation admits only the trivial solution $u_{(2)}=0$,
whereas in the flat spacetime limit ($M\rightarrow 0$) we recover the presence
of electrostatic degrees of freedom, $\omega^{2}=\omega_{\rm pl}^{2}$, which
are typical of photons propagating in cold plasma.

For $l\neq 0$, the polar equations for $u_{(2)}$ and $u_{(3)}$ can be reduced
to a single Schrödinger-like equation,

$\frac{d^{2}}{dr_{*}^{2}}\psi-V(r)\psi=0,$ (27)

where the complicated form of the effective potential $V(r)$ is given in
Appendix A. As expected, $V(r\to 2M)\to-\omega^{2}$ and
$V(r\to\infty)\to\omega_{\rm pl}^{2}(r\to\infty)-\omega^{2}$. The fact that
the polar sector reduces to a single second-order differential equations
implies the presence of only one dynamical degree of freedom. Note that the
potential $V$ depends on the plasma frequency $\omega_{\rm pl}(r)$ and also on
its radial derivatives.

###  II.3 The hydrogenic spectrum for Proca modes

In the next section we shall compute the quasibound states of our problem by
solving Eqs. (25) and (27) numerically with suitable boundary conditions. The
latter select an infinite set of complex eigenfrequencies:
$\omega=\omega_{R}+i\omega_{I}$. It will be useful to compare the results with
the case of a Proca field with mass $\hbar\omega_{\rm pl}={\rm const}$ around
a Schwarzschild BH. In the latter case, the spectrum has a clear physical
interpretation in the Newtonian limit, which corresponds to the Compton
wavelength of the Proca field, $\sim 1/\omega_{\rm pl}$, being much larger
than the horizon size. This requires $M\omega_{\rm pl}\ll 1$. To leading order
in this limit the spectrum of quasibound states has a hydrogenic form
Galtsov:1984ixy ; Rosa:2011my ; Pani:2012vp ; Pani:2012bp ; Baryakhtar:2017ngi

$\displaystyle\omega_{R}$ $\displaystyle\sim$ $\displaystyle\omega_{\rm
pl}\left(1-\frac{(M\omega_{\rm pl})^{2}}{2(l+S+1+n)^{2}}\right)\,,$ (28)
$\displaystyle\omega_{I}$ $\displaystyle\sim$
$\displaystyle-\frac{1}{2}C^{(1)}_{lS}(M\omega_{\rm
pl})^{4l+2S+5}\omega_{R}\,,$ (29)

where $l$ is the total angular momentum222This can be seen by applying the
angular-momentum operator to the state Baryakhtar:2017ngi . Note that $l$ was
erroneously identified with the orbital angular momentum in Ref. Rosa:2011my .
See Ref. Brito:2015oca  for a discussion of different spherical bases used in
the literature. of the state with spin projections $S=-1$, 0, 1 (with $S=0$
for axial modes and $S=\pm 1$ for the two polarizations of polar modes), $n$
is the overtone number ($n=0$ for the longest-lived, fundamental mode), and
$C^{(1)}_{lS}$ are constants (given, e.g., in Ref. Brito:2015oca ). The most
unstable mode is the polar dipole with $S=-1$, $l=1$, which has
$C^{(1)}_{1-1}=16$ Baryakhtar:2017ngi . The dominant slope of the imaginary
part of the mode is $\omega_{I}\propto(M\omega_{\rm pl})^{10}$ for the axial
dipole, and $\omega_{I}\propto(M\omega_{\rm pl})^{8}$ and
$\omega_{I}\propto(M\omega_{\rm pl})^{12}$ for the two polar dipole modes. The
above analytical approximation in the Newtonian limit is in excellent
agreement with the exact numerical results Detweiler:1980uk ; Dolan:2007mj ;
Baryakhtar:2017ngi ; Cardoso:2018tly ; Dolan:2018dqv ; Baumann:2019eav ;
Brito:2013wya ; Brito:2020lup .

In the spinning case, the imaginary part acquires a factor
$\omega_{I}\propto(\omega_{R}-m\Omega)$, which depends on the BH angular
velocity $\Omega$. Therefore, in the superradiant regime,
$\omega_{R}<m\Omega$, the modes with the smallest slope in the static case
(namely, the polar dipole with $S=-1$) become the ones with the shortest
instability timescale, $\tau=1/\omega_{I}$ Brito:2015oca  (see Sec. IV for
more details).

In the next section we shall compare our numerical results for the full
plasma-photon system with the above hydrogenic behavior of a Proca field.

###  II.4 Plasma profiles

We consider two different plasma profiles: a homogeneous density profile and a
Bondi-like accretion disk model. We will examine first the case of homogeneous
density, and therefore homogeneous plasma frequency, $\omega_{\rm pl}={\rm
const}$. This approximation is obviously not realistic, especially close to
the BH. However, it will serve as a warm-up to elucidate the structure of the
equations and identify the correct limit asymptotically far from the BH.

Then, we will consider a radial dependence $\omega_{\rm pl}=\omega_{\rm
pl}(r)$, as expected in the surroundings of a spherically symmetric BH. In
particular, we will consider the Bondi-like accretion disk model, which
describes the accretion dynamics of a non-self-interacting gas around a
spherically symmetric compact object 1952MNRAS.112..195B . This model predicts
a power-law density profile for the gas and consequently a plasma frequency of
the type

$\omega_{\rm pl}^{2}(r)=\omega_{\rm
B}^{2}\Big{(}\frac{2M}{r}\Big{)}^{\lambda}+\omega_{\rm\infty}^{2},$ (30)

where $\sqrt{\omega_{\rm B}^{2}+\omega_{\infty}^{2}}$ is the plasma frequency
at the horizon (since in the relevant regime $\omega_{\infty}\ll\omega_{\rm
B}$, with a little abuse of notation we shall refer to $\omega_{\rm B}$ as the
horizon plasma density). The slope $\lambda$ depends on the adiabatic index of
the gas (e.g., $\lambda=3/2$ for monoatomic species). The constant term
$\omega_{\infty}$ is the asymptotic plasma frequency at infinity, i.e the
interstellar medium plasma frequency far away from the central BH.

##  III Plasma-driven quasibound EM modes of a Schwarzschild BH

###  III.1 Numerical method

We compute the characteristic frequencies of our system using a direct
integration shooting method Ferrari:2007rc ; Pani:2012bp ; Rosa:2011my . The
main idea is to integrate the system from the horizon outwards to infinity,
imposing suitable boundary conditions. Close to the horizon, the solution must
be a purely ingoing wave, as the horizon behaves as a one-way membrane,

$u_{(i)}\sim e^{-i\omega r_{*}}\sum_{n}b_{(i)\,n}(r-2M)^{n},$ (31)

where the coefficients $b_{(i)\,n}$ can be computed in terms of the arbitrary
coefficient $b_{(i)\,0}$ by expanding the relevant equations near the horizon.
At infinity, the solution at the leading order can be written in its most
generic form as

$u_{(i)}\sim B_{(i)}e^{-k_{\infty}r_{*}}+C_{(i)}e^{+k_{\infty}r_{*}},$ (32)

where $k_{\infty}=\sqrt{\omega_{\rm pl}^{2}(r\to\infty)-\omega^{2}}$. We are
interested in strongly localized, quasibound state solutions for
$\omega<\omega_{\rm pl}$, and we will therefore impose the condition
$C_{(i)}=0$, implying exponentially damped solutions at infinity, and solve
the associated eigenvalue problem. This method does not rely on a specific
shape of the effective potential and is thus extremely flexible. Moreover, the
direct integration method works particularly well for the computation of
quasibound states, while it is expected to be less precise for the computation
of normal modes (corresponding to the boundary condition $B_{(i)}=0$),
especially for the highly damped modes.

Since the axial and polar sectors are decoupled and each one is described by a
single second-order differential equation, for a given value of $l$ we expect
to find two families of modes (axial and polar), each one defined by an
overtone number $n$. In the following we shall mainly focus on the fundamental
($n=0$) modes, although tracking a fixed overtone number is difficult in the
small $M\omega_{\rm pl}$ regime. As we shall see, in the same regime the
imaginary part of the mode can become extremely small and sensitive to small
numerical errors. The latter can be reduced by increasing the numerical
accuracy and the truncation order of the series expansions at the horizon and
at infinity.

This procedure is also suitable to find unstable modes, i.e. those with
$\omega_{I}>0$ and whose time dependence is exponentially growing as
$e^{\omega_{I}t}$. We have searched for such modes in different configurations
and found none, confirming the reasonable expectation that the system is
stable in the static case.

###  III.2 Constant-density plasma

Figure 1: Imaginary (upper panel) and real (lower panel) part of the plasma-
driven quasibound EM mode of a Schwarzschild BH in the axial sector for the
case $\omega_{\rm pl}={\rm const}$. The imaginary part of the frequency is
always negative, while the real part is always smaller than the plasma
frequency. We show the fundamental mode ($n=0$, orange) and the first overtone
($n=1$, blue) for $l=1$. The dashed curves correspond to the hydrogenic
spectrum, which is recovered as expected in the limit $M\omega_{\rm pl}\ll 1$.
Note that the axial sector is equivalent to that of a Proca field Rosa:2011my
.

Figure 1 shows the real part and the imaginary part of the axial eigenvalues
$\omega=\omega_{R}+i\omega_{I}$, normalized to the plasma frequency. We
consider $l=1$ and two different overtone numbers, respectively $n=0$ (orange)
and $n=1$ (blue). As expected, the imaginary part of the frequency is always
negative and represents a mode which is exponentially decaying in time, while
the real part is always smaller than $\omega_{\rm pl}$, a necessary condition
to obtain solutions that are confined in the vicinity of the BH. Since the
axial sector is the same as in the Proca case, this plot is equivalent to what
found in Ref. Rosa:2011my  for the axial quasibound states of a Proca field
around a Schwarzschild BH. The dashed curves shown in Fig. 1 correspond to the
hydrogenic spectrum expected in the $M\omega_{\rm pl}\ll 1$ limit.

While in the axial sector the potential is simple and equivalent to that of a
Proca equation, in the polar sector the situation is different. Indeed, not
only are the equations more involved than in the Proca case, but the polar
sector only contains a single propagating mode, in contrast with the two polar
modes of a Proca field. Thus, as expected, the structure of the plasma-driven
polar quasibound states is different and cannot be mapped into a hydrogenic
spectrum as in the axial case. This is shown in Fig. 2.

Figure 2: Same as in Fig. 1 but for the polar sector. In this case the modes
do not reduce to the hydrogenic Proca spectrum in the $M\omega_{\rm pl}\ll 1$
limit.

In this case the numerical computation of the eigenfrequencies becomes
increasingly more challenging for values of the coupling $M\omega_{\rm
pl}\lesssim 0.4$ and we therefore do not show the lower part of the spectrum
in Fig. 2. However, in order to investigate the small-mass coupling regime,
$M\omega_{\rm pl}\ll 1$, we separately studied smaller intervals around
$M\omega_{\rm pl}\simeq 0.15$, as shown in Fig. 3. Although the numerical
results are noisy in this limit, we find that the scaling of $\omega_{I}$ is
significantly different from the hydrogenic behavior predicted in the Proca
case for $S=-1$. Indeed, our best fit yields $\omega_{I}/\omega_{\rm
pl}\sim(M\omega_{\rm pl})^{10.84}$, to be compared with
$\omega_{I}/\omega_{\rm pl}\sim(M\omega_{\rm pl})^{7}$ for the $S=-1$ Proca
mode. In fact, we find a scaling that is quite close to the (subleading) polar
mode with polarization $S=1$ for the massive vector case Rosa:2011my .

Figure 3: Imaginary part of the plasma-driven EM quasibound mode in the polar
sector in $M\omega_{\rm pl}\ll 1$ regime. The blue dots indicate the numerical
data; the black solid line is the best-fit power law, while the red line
denotes the hydrogenic scaling in the Proca case.

In some sense, the plasma-driven EM quasibound states are phenomenologically
closer to a massive scalar field than to a Proca field. In both plasma-driven
EM fields and massive scalar fields, the minimum scaling of the imaginary part
(which sets the shortest, and therefore most interesting, mode lifetime) is
$\omega_{I}/\omega_{\rm pl}\sim(M\omega_{\rm p})^{9}$. This scaling is set by
the axial dipole mode in the plasma-photon system. Thus, the lifetime of the
fastest quasibound states in the plasma-driven case is parametrically longer
than for a Proca field – for which the shortest lifetime is set by the polar
modes – by a factor $(M\omega_{\rm pl})^{-2}\gg 1$.

###  III.3 Plasma profile from Bondi accretion

We now analyze the quasibound states for a more realistic plasma density
profile, as predicted by Bondi accretion. In this case the axial sector is
also described by Eq. (25), except for a radial dependence in the plasma
frequency $\omega_{\rm pl}\rightarrow\omega_{\rm pl}(r)=\sqrt{\omega_{\rm
B}^{2}(2M/r)^{\lambda}+\omega_{\infty}^{2}}$. The polar sector is again more
involved and its potential depends nontrivially on the radial derivatives of
the plasma frequency. The equations for the polar case are presented in
Appendix A.

In Fig. 4 we show the results for the imaginary and real part of the frequency
in the axial sector as a function of $\omega_{\infty}$, the quantity that
would asymptotically corresponds to a mass term. In all cases we fix
$\lambda=3/2$ for the Bondi model (corresponding to monoatomic gas) and we
show the results for different values of the horizon plasma frequency
$\omega_{\rm B}$. We notice that, as in the homogeneous plasma case, the real
(imaginary) part of the frequency decreases (increases) monotonically with
$M\omega_{\infty}$. The effect of $\omega_{\rm B}$ is the opposite: larger
values of the density in the vicinity of the BHs lead to a smaller imaginary
part and a larger real part. Indeed, for larger values of $\omega_{\rm B}$ it
is more difficult to find minima in the effective potentials, i.e. it is more
difficult to support quasibound states. In fact, we were not able to find
quasibound states solutions for $\omega_{\rm B}M\gtrsim 1.7$, in good
agreement with the numerical study in Ref. Dima:2020rzg  for a scalar toy
model. Furthermore, this behavior is consistent with the
$\omega_{B}\rightarrow 0$ limit, as in this case the spectrum shown in Fig. 1,
with smaller values of the real part of the eigenfrequencies, must be
recovered.

Figure 4: Imaginary (upper panel) and real (lower panel) part of the
quasibound state frequency in the axial sector for a photon in a Bondi-
accretion plasma model, see Eq. (30). The imaginary part is always negative,
while the real part is always smaller than the plasma frequency, as expected.
The imaginary part of the frequency decreases (the modes become short-lived)
as the horizon plasma $\omega_{\rm B}$ density increases.

Note that, in general, the imaginary part of the quasibound states for a Bondi
plasma profile is smaller than in the homogeneous-plasma case. This makes the
numerical computation harder, especially for polar modes for which the
imaginary part is even smaller than in the axial case. Although not shown, we
were also able to compute the plasma-driven polar mode in the case of a Bondi
profile, at least when $M\omega_{\infty}={\cal O}(1)$. As in the homogeneous
case, the axial modes are found to be the shortest lived for small values of
the plasma frequency.

##  IV The role of the BH spin

In this work we focused on finding the plasma-driven EM quasibound states
around a nonspinning (Schwarzschild) BH. A natural and paramount extension is
to consider the spinning case, namely solving the equations on a Kerr
background metric.

This case is significantly more involved, but also more interesting in view of
the fact that the quasibound states can become (superradiantly) unstable in
the spinning case Brito:2015oca . A detailed analysis of the photon-plasma
interactions around a Kerr BH is left for future work. Here we anticipate the
main differences relative to the static case considered here:

  * •

Vorticity. In a Kerr BH the vorticity tensor does not vanish, i.e. (from Eq.
(10)) $\omega_{ab}\sim{\cal O}(\chi)$, where $\chi=J/M^{2}$ is the
dimensionless angular momentum of a Kerr BH. On the other hand, the rate of
expansion is still zero, $\theta_{ab}=0$. Therefore, the linearized field
equations (12) contain extra terms which vanish in the static case;

  * •

Plasma profile. A spinning BH drags plasma in its vicinity. In the stationary
and axisymmetric case, the plasma density can generically depend both on the
radius and on the angle $\theta$, i.e. $\omega_{\rm pl}=\omega_{\rm
pl}(r,\theta)$;

  * •

Mode coupling. In the spinning case, axial and polar modes with different
harmonic index $l$ are coupled to each other Kojima:1992ie ; Pani:2012bp ;
Pani:2012vp ; Pani:2013pma . Thus, even if the axial sector is the same as in
the Proca case when the BH is static, in the spinning case it would get
corrections arising from the different polar sector. For a Kerr BH we
therefore expect that the entire quasibound spectrum will differ from the
Proca one.

Overall, for the above reasons do not expect the linearized equations to be
separable, at variance with the Proca case Frolov:2018ezx  (unless possibly
for specific profiles, see e.g. Refs. Cardoso:2013opa ; Cardoso:2013fwa ). The
linear dynamics can still be computed perturbatively in the slow-rotation
approximation Pani:2012bp ; Pani:2012vp ; Pani:2013pma , or fully numerically
(see, e.g., Ref. Cardoso:2018tly  for the Proca case), or using the methods
adopted in Refs. Dolan:2012yt ; Baumann:2019eav ; Dima:2020rzg .

It is reasonable to expect that, for a Kerr BH, the quasibound spectrum is
modified by corrections proportional to the superradiant Brito:2015oca  factor
$\omega-m\Omega$, where $m$ is the azimuthal number of the perturbation and
$\Omega$ is the BH angular velocity. In particular, the imaginary part of the
quasibound state might change sign when $\omega<m\Omega$, turning the mode
from stable to unstable. This correction is evident in the hydrogenic spectrum
of massive bosonic perturbations of a Kerr BH Detweiler:1980uk ; Pani:2012bp ;
Pani:2012vp ; Baryakhtar:2017ngi ; Cardoso:2018tly ; Brito:2020lup  in the
(Newtonian) small-mass limit. Nonetheless, for the reasons listed above one
cannot exclude that the quasibound spectrum of a Kerr BH can be qualitatively
different from the Proca case. A detailed analysis is required and will appear
elsewhere.

It is important to stress that the plasma-driven superradiant instability can
be strongly quenched by nonlinear and relativistic effects Cardoso:2020nst ;
Blas:2020nbs . Thus, although we expect that the linear spectrum of plasma-
driven EM modes around a Kerr BH is significantly different than in the Proca
case, the nonlinear and relativistic quenching might be efficient also for
this more realistic configuration.

For the Schwarzschild background considered in this work, expanding Eqs.
(4)–(7) to second order in the perturbations, we find that the linearized
analysis is valid provided

$|\partial_{i}E^{i}|\ll en_{e}\,.$ (33)

This equation implies $e\omega E/v\ll\omega_{\rm pl}^{2}m_{e}$, where $E$ is
the amplitude of the electric field, and $v=\omega/k$ is the velocity of the
EM wave in the plasma. For $v\sim 1$ and $\omega\sim\omega_{\rm pl}$ we
recover the condition $e|E|\ll m_{e}\omega_{\rm pl}$, which implies that the
plasma backreaction on the dispersion relation is negligible. For a stable
quasibound state that decays in time, this condition is always satisfied if it
holds initially. However, for a superradiantly growing state the condition is
eventually violated before the EM field can extract enough energy from the BH
Cardoso:2020nst .

##  V Conclusion and extensions

We initiated a detailed investigation of the equations governing the
linearized dynamics of an EM field propagating in a plasma in curved
spacetime. Our main goal was to explore the differences with the case of Proca
theory which has been so far used as a proxy for the plasma-photon interaction
around a BH. We focused on the case of cold plasma around a background
Schwarzschild BH and computed the quasibound spectrum of EM modes. We showed
that, while the axial sector is identical to the Proca one, polar
perturbations propagate a single degree of freedom (at variance with the two
polar modes of a Proca field) which is described by a new field equation. This
drastically affects the polar quasibound-mode spectrum and the scaling of
these modes in the small-frequency limit.

We also studied the stability of the system, by looking for unstable modes
numerically. We explored a large region of the parameter space and found no
unstable modes, as expected in the static case.

As discussed in the previous section a natural and urgent extension is to
consider the case of plasma-driven EM modes of a spinning (Kerr) BH, which is
likely to be drastically different from the Proca case. We expect
quantitatively (and perhaps qualitatively) different results for what concerns
the superradiant instability. In this context, our analysis suggests that the
linearized analyses performed in the past (e.g., Pani:2013hpa ; Conlon:2017hhi
; Dima:2020rzg ) in the context of plasma-driven superradiant instabilities
around spinning BHs should be revisited.

Although a detailed analysis in the spinning case is required, our results
already suggest that the dominant unstable modes could be more akin to the
scalar massive modes of a Kerr BH rather than to the vector massive modes. If
confirmed, this would imply that the toy model adopted in Ref. Dima:2020rzg
could be more accurate than what one would in principle expect.

Another interesting extension of our work is to study nonlinear plasma-photon
effects in a coherent and covariant framework. It was recently argued that
nonlinear and relativistic effects can strongly quench the superradiant
instability Cardoso:2020nst ; Blas:2020nbs . The framework adopted here could
be used to study this interesting problem quantitatively. For example, one
could expand the exact field equations to second order in the perturbations
and study the backreaction of the plasma onto the EM field.

We also neglected magnetic fields and focused on the case of cold plasma.
Considering a magnetized/hot plasma is another natural extension.

Finally, we assumed that the plasma was static, but this is not the case for
an accreting BH. While this approximation is valid for modes that are short-
lived with respect to the typically accretion timescale of the BH, this might
not always be the case, especially for problems (e.g. the superradiant
instability) which could be characterized by extremely long timescales. Thus,
another interesting extension would be to solve the linearized photon-plasma
dynamics in the time domain, also taking into account the dynamics of the
plasma density, or perform a multiscale adiabatic approximation Brito:2014wla
.

We hope to report on these interesting problems in the near future.

###### Acknowledgements.

We thank Richard Brito and Vitor Cardoso for useful conversations and
comments. P.P. acknowledges financial support provided under the European
Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480. We also
acknowledge support under the MIUR PRIN and FARE programmes (GW-NEXT, CUP:
B84I20000100001), and from the Amaldi Research Center funded by the MIUR
program “Dipartimento di Eccellenza” (CUP: B81I18001170001). A.C. acknowledges
support from the the Israel Science Foundation (Grant No. 1302/19), the US-
Israeli BSF (grant 2018236) and the German Israeli GIF (grant I-2524-303.7).
A.C. acknowledges hospitality from the MPP of Munich. L.S. acknowledges that
research at Perimeter Institute is supported in part by the Government of
Canada through the Department of Innovation, Science and Economic Development
Canada, and by the Province of Ontario through the Ministry of Colleges and
Universities.

##  Appendix A Polar sector

Here we derive the full potential for the polar sector. We first consider a
homogeneous plasma ($\omega_{\rm pl}={\rm const}$) and then a generic radially
dependent density profile, i.e. $\omega_{\rm pl}=\omega_{\rm pl}(r)$.

###  A.1 Homogeneous plasma case

From Eq. (22), it is possible to find an expression for $u_{(2)}$ in terms of
$u_{(3)}^{\prime}$. This expression can be inserted in Eq. (23) to get a
single second order equation for the decoupled variable $u_{(3)}$:

$\displaystyle
u_{(3)}^{\prime\prime}-u^{\prime}_{(3)}\frac{2f^{2}l(l+1)\omega^{2}+(f-1)f^{2}r^{2}\omega_{\text{pl}}^{4}-f\omega_{\text{pl}}^{2}\left(f(f+1)l(l+1)+2(f-1)r^{2}\omega^{2}\right)+(f-1)r^{2}\omega^{4}}{fr\left(f\omega_{\text{pl}}^{2}-\omega^{2}\right)(f\left(l^{2}+l+r^{2}\omega_{\text{pl}}^{2}\right)-r^{2}\omega^{2})}$
$\displaystyle+u_{(3)}\left(\omega^{2}-\frac{f\left(l^{2}+l+r^{2}\omega_{\text{pl}}^{2}\right)}{r^{2}}\right)f^{-2}$
$\displaystyle=0.$ (34)

Equation (A.1) is still not in the form of a Schroedinger-like equation. We
therefore perform the following substitution: $u_{(3)}=G(r)\psi(r)$, where $G$
is an unknown function such that with this transformation equation Eq. (A.1)
assumes the form $d^{2}\psi/dr_{*}^{2}-V(r)\psi=0$. We find
$G(r)=r^{-1}\sqrt{fr(l^{2}+l+r^{2}\omega_{\rm
pl}^{2})-r^{3}\omega^{2}}\sqrt{r(f\omega_{\rm pl}^{2}-\omega^{2})}^{-1}$ and
we finally arrive at a Schrodinger-like equation in terms of the tortoise
coordinate where the potential has the following cumbersome form

$\displaystyle
V(r)=\frac{d(r)}{4\left(\omega^{2}-f\omega_{\text{pl}}{}^{2}\right)\left(fl(l+1)r+fr^{3}\omega_{\text{pl}}{}^{2}-r^{3}\omega^{2}\right)}\,,$
(35)

with

$d(r)=-\omega^{4}\Big{[}f^{4}\left(3l^{2}(l+1)^{2}+72l(l+1)r^{2}\omega_{\text{pl}}^{2}\right)+2f^{3}\Big{(}-18l(l+1)\left(l^{2}+l+1\right)r^{2}\omega_{\text{pl}}^{2}\\\
-l^{2}(l+1)^{2}(2l(l+1)+1)-36l(l+1)r^{4}\omega_{\text{pl}}^{4}-20r^{6}\omega_{\text{pl}}^{6}\Big{)}-f^{2}l^{2}(l+1)^{2}\Big{]}\\\
-\omega^{8}\left(-12fl(l+1)r^{4}-20fr^{6}\omega_{\text{pl}}^{2}\right)-4r^{6}\omega^{10}-\omega^{2}\Big{[}2f^{5}\omega_{\text{pl}}^{2}\left(-3l^{2}(l+1)^{2}-29l(l+1)r^{2}\omega_{\text{pl}}^{2}\right)\\\
+f^{4}\left(4l^{2}(l+1)^{2}(2l(l+1)+1)\omega_{\text{pl}}^{2}+48l(l+1)r^{4}\omega_{\text{pl}}^{6}+4l(l+1)(9l(l+1)+5)r^{2}\omega_{\text{pl}}^{4}+20r^{6}\omega_{\text{pl}}^{8}\right)\\\
+2f^{3}\left(l^{2}(l+1)^{2}\omega_{\text{pl}}^{2}+l(l+1)r^{2}\omega_{\text{pl}}^{4}\right)\Big{]}-4f^{6}l(l+1)\omega_{\text{pl}}^{4}\left(l^{2}+l+4r^{2}\omega_{\text{pl}}^{2}\right)\\\
-2f^{5}\omega_{\text{pl}}^{2}\left(-2l^{2}(l+1)^{2}\left(l^{2}+l+1\right)\omega_{\text{pl}}^{2}-6l(l+1)r^{4}\omega_{\text{pl}}^{6}-2l(l+1)(3l(l+1)+1)r^{2}\omega_{\text{pl}}^{4}-2r^{6}\omega_{\text{pl}}^{8}\right)\\\
-\omega^{6}\left(f^{2}\left(48l(l+1)r^{4}\omega_{\text{pl}}^{2}+4l(l+1)(3l(l+1)+5)r^{2}+40r^{6}\omega_{\text{pl}}^{4}\right)-30f^{3}l(l+1)r^{2}-2fl(l+1)r^{2}\right)\,.$
(36)

Note that $V\to-\omega^{2}$ as $r\to 2M$ and $V\to\omega_{\rm
pl}^{2}-\omega^{2}$ as $r\to\infty$.

###  A.2 Nonhomogeneous plasma case

The case of a generic radially dependent plasma density proceeds in the same
fashion as the homogeneous case. We find the same expression for the function
$G$ with $\omega_{\rm pl}$ replaced by $\omega_{\rm pl}(r)$. As in the
previous case we then arrive to a Schrodinger-like equation with an effective
potential

$\displaystyle
V(r)=\frac{d(r)}{4\left(\omega^{2}-f\omega_{\text{pl}}(r){}^{2}\right)\left(fl(l+1)r+fr^{3}\omega_{\text{pl}}(r){}^{2}-r^{3}\omega^{2}\right)}\,,$
(37)

where

$d(r)=\omega^{4}\Big{[}36f^{3}l(l+1)r^{2}\left(-2f+l^{2}+l+1\right)\omega_{\text{pl}}(r){}^{2}+f^{2}l^{2}(l+1)^{2}(f(-3f+4l(l+1)+2)+1)\\\
-4f^{3}l(l+1)r\omega_{\text{pl}}(r)\left((9f-5)r^{2}\omega_{\text{pl}}^{\prime}(r)-fr^{3}\omega_{\text{pl}}^{\prime\prime}(r)\right)\\\
+4f^{4}l(l+1)r^{4}\omega_{\text{pl}}^{\prime}(r){}^{2}+72f^{3}l(l+1)r^{4}\omega_{\text{pl}}(r){}^{4}+40f^{3}r^{6}\omega_{\text{pl}}(r){}^{6}\Big{]}+\\\
\omega^{2}\Big{[}2f\omega_{\text{pl}}(r){}^{2}\left(f^{2}l^{2}(l+1)^{2}(f(3f-4l(l+1)-2)-1)+4f^{4}l(l+1)r^{4}\omega_{\text{pl}}^{\prime}(r){}^{2}\right)\\\
-4f^{5}l^{2}(l+1)^{2}r^{2}\omega_{\text{pl}}^{\prime}(r){}^{2}+4f^{4}l(l+1)r\omega_{\text{pl}}(r){}^{3}\left(2(7f-3)r^{2}\omega_{\text{pl}}^{\prime}(r)-2fr^{3}\omega_{\text{pl}}^{\prime\prime}(r)\right)\\\
-48f^{4}l(l+1)r^{4}\omega_{\text{pl}}(r){}^{6}-2f^{3}l(l+1)r^{2}(f(-29f+18l(l+1)+10)+1)\omega_{\text{pl}}(r){}^{4}\\\
-4f^{3}l(l+1)r\omega_{\text{pl}}(r)\left(f^{2}l(l+1)r\omega_{\text{pl}}^{\prime\prime}(r)+2fl(l+1)\omega_{\text{pl}}^{\prime}(r)\right)-20f^{4}r^{6}\omega_{\text{pl}}(r){}^{8}\Big{]}\\\
-2f^{2}\omega_{\text{pl}}(r){}^{4}\left(6f^{4}l(l+1)r^{4}\omega_{\text{pl}}^{\prime}(r){}^{2}-2f^{3}l^{2}(l+1)^{2}\left(-f+l^{2}+l+1\right)\right)\\\
+\omega^{6}\left(-48f^{2}l(l+1)r^{4}\omega_{\text{pl}}(r){}^{2}-40f^{2}r^{6}\omega_{\text{pl}}(r){}^{4}+2fl(l+1)r^{2}(f(15f-6l(l+1)-10)+1)\right)\\\
+4f^{5}l(l+1)r^{3}\omega_{\text{pl}}(r){}^{5}\left((1-5f)\omega_{\text{pl}}^{\prime}(r)+fr\omega_{\text{pl}}^{\prime\prime}(r)\right)+12f^{5}l(l+1)r^{4}\omega_{\text{pl}}(r){}^{8}\\\
+4f^{4}l(l+1)r\omega_{\text{pl}}(r){}^{3}\left(f^{2}l(l+1)r\omega_{\text{pl}}^{\prime\prime}(r)+(f+1)fl(l+1)\omega_{\text{pl}}^{\prime}(r)\right)-4r^{6}\omega^{10}\\\
-4f^{5}l(l+1)r^{2}(4f-3l(l+1)-1)\omega_{\text{pl}}(r){}^{6}+4f^{5}r^{6}\omega_{\text{pl}}(r){}^{10}+\omega^{8}\left(12fl(l+1)r^{4}+20fr^{6}\omega_{\text{pl}}(r){}^{2}\right)$
(38)

In this case the asymptotic behavior reads $V\to-\omega^{2}$ as $r\to 2M$ and
$V\to\omega_{\rm pl}^{2}(r\to\infty)-\omega^{2}$ as $r\to\infty$. Obviously,
the effective potential in Eq. (37) reduces to Eq. (35) when $\omega_{\rm
pl}={\rm const}$.

## References

  * (1) L. Barack et al., “Black holes, gravitational waves and fundamental physics: a roadmap,” Class. Quant. Grav. 36 no. 14, (2019) 143001, arXiv:1806.05195 [gr-qc].
  * (2) N. I. Shakura and R. A. Sunyaev, “Reprint of 1973A&A….24..337S. Black holes in binary systems. Observational appearance.,” A&A 500 (June, 1973) 33–51.
  * (3) I. D. Novikov and K. S. Thorne, “Astrophysics and black holes,” in Proceedings, Ecole d’Eté de Physique Théorique: Les Astres Occlus, pp. 343–550.  1973\.
  * (4) A. G. Sitenko, Electromagnetic Fluctuations in Plasma.  Academic, New York, 1976.
  * (5) 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].
  * (6) P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, “Black hole bombs and photon mass bounds,” Phys. Rev. Lett. 109 (2012) 131102, arXiv:1209.0465 [gr-qc].
  * (7) P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, “Perturbations of slowly rotating black holes: massive vector fields in the Kerr metric,” Phys. Rev. D 86 (2012) 104017, arXiv:1209.0773 [gr-qc].
  * (8) M. Baryakhtar, R. Lasenby, and M. Teo, “Black Hole Superradiance Signatures of Ultralight Vectors,” Phys. Rev. D96 no. 3, (2017) 035019, arXiv:1704.05081 [hep-ph].
  * (9) 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,” JCAP 03 (2018) 043, arXiv:1801.01420 [gr-qc].
  * (10) V. P. Frolov, P. Krtous, D. Kubiznak, and J. E. Santos, “Massive Vector Fields in Rotating Black-Hole Spacetimes: Separability and Quasinormal Modes,” Phys. Rev. Lett. 120 (2018) 231103, arXiv:1804.00030 [hep-th].
  * (11) S. R. Dolan, “Instability of the Proca field on Kerr spacetime,” Phys. Rev. D98 no. 10, (2018) 104006, arXiv:1806.01604 [gr-qc].
  * (12) D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, “The Spectra of Gravitational Atoms,” JCAP 1912 no. 12, (2019) 006, arXiv:1908.10370 [gr-qc].
  * (13) R. Brito, V. Cardoso, and P. Pani, Superradiance: New Frontiers in Black Hole Physics, vol. 906.  Springer, 2015.  arXiv:1501.06570 [gr-qc].
  * (14) P. Pani and A. Loeb, “Constraining Primordial Black-Hole Bombs through Spectral Distortions of the Cosmic Microwave Background,” Phys. Rev. D 88 (2013) 041301, arXiv:1307.5176 [astro-ph.CO].
  * (15) D. Blas and S. J. Witte, “Imprints of Axion Superradiance in the CMB,” arXiv:2009.10074 [astro-ph.CO].
  * (16) J. P. Conlon and C. A. Herdeiro, “Can black hole superradiance be induced by galactic plasmas?,” Phys. Lett. B 780 (2018) 169–173, arXiv:1701.02034 [astro-ph.HE].
  * (17) V. Cardoso, W.-d. Guo, C. F. Macedo, and P. Pani, “The tune of the universe: the role of plasma in tests of strong-field gravity,” arXiv:2009.07287 [gr-qc].
  * (18) D. Blas and S. J. Witte, “Quenching Mechanisms of Photon Superradiance,” arXiv:2009.10075 [hep-ph].
  * (19) M. A. Abramowicz and P. Fragile, “Foundations of Black Hole Accretion Disk Theory,” Living Rev. Rel. 16 (2013) 1, arXiv:1104.5499 [astro-ph.HE].
  * (20) F. Yuan and R. Narayan, “Hot Accretion Flows Around Black Holes,” Ann. Rev. Astron. Astrophys. 52 (2014) 529–588, arXiv:1401.0586 [astro-ph.HE].
  * (21) K. S. Thorne and R. D. Blandford, Modern Classical Physics: Optics, Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics.  2017\.
  * (22) R. Kulsrud and A. Loeb, “Dynamics and gravitational interaction of waves in nonuniform media,” Phys. Rev. D 45 (Jan, 1992) 525–531. https://link.aps.org/doi/10.1103/PhysRevD.45.525.
  * (23) J. T. Mendonça, A. M. Martins, and A. Guerreiro, “Field quantization in a plasma: Photon mass and charge,” Phys. Rev. E 62 (Aug, 2000) 2989–2991. https://link.aps.org/doi/10.1103/PhysRevE.62.2989.
  * (24) D. N. Page and K. S. Thorne, “Disk-Accretion onto a Black Hole. Time-Averaged Structure of Accretion Disk,” ApJ 191 (July, 1974) 499–506.
  * (25) A. Arvanitaki and S. Dubovsky, “Exploring the String Axiverse with Precision Black Hole Physics,” Phys.Rev. D83 (2011) 044026, arXiv:1004.3558 [hep-th].
  * (26) R. Brito, V. Cardoso, and P. Pani, “Black holes as particle detectors: evolution of superradiant instabilities,” Class. Quant. Grav. 32 no. 13, (2015) 134001, arXiv:1411.0686 [gr-qc].
  * (27) R. A. Breuer and J. Ehlers, “Propagation of electromagnetic waves through magnetized plasmas in arbitrary gravitational fields,” A&A 96 no. 1-2, (Mar., 1981) 293–295.
  * (28) G. Ellis, “Relativistic cosmology,” Gen. Rel. Grav. 41 (2009) 581–660.
  * (29) G. Raffelt, Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles.  5, 1996.
  * (30) D. Gal’tsov, G. Pomerantseva, and G. Chizhov, “Behavior of massive vector particles in a Schwarzschild field,” Sov. Phys. J. 27 (1984) 697–700.
  * (31) S. L. Detweiler, “Klein-Gordon equation and rotating black holes,” Phys. Rev. D 22 (1980) 2323–2326.
  * (32) S. R. Dolan, “Instability of the massive Klein-Gordon field on the Kerr spacetime,” Phys.Rev. D76 (2007) 084001, arXiv:0705.2880 [gr-qc].
  * (33) R. Brito, V. Cardoso, and P. Pani, “Massive spin-2 fields on black hole spacetimes: Instability of the Schwarzschild and Kerr solutions and bounds on the graviton mass,” Phys. Rev. D88 (2013) 023514, arXiv:1304.6725 [gr-qc].
  * (34) R. Brito, S. Grillo, and P. Pani, “Black Hole Superradiant Instability from Ultralight Spin-2 Fields,” Phys. Rev. Lett. 124 no. 21, (2020) 211101, arXiv:2002.04055 [gr-qc].
  * (35) H. Bondi, “On spherically symmetrical accretion,” MNRAS 112 (Jan., 1952) 195.
  * (36) V. Ferrari, L. Gualtieri, and S. Marassi, “A New approach to the study of quasi-normal modes of rotating stars,” Phys. Rev. D 76 (2007) 104033, arXiv:0709.2925 [gr-qc].
  * (37) A. Dima and E. Barausse, “Numerical investigation of plasma-driven superradiant instabilities,” Class. Quant. Grav. 37 no. 17, (2020) 175006, arXiv:2001.11484 [gr-qc].
  * (38) Y. Kojima, “Equations governing the nonradial oscillations of a slowly rotating relativistic star,” Phys. Rev. D 46 (1992) 4289–4303.
  * (39) P. Pani, “Advanced Methods in Black-Hole Perturbation Theory,” Int. J. Mod. Phys. A28 (2013) 1340018, arXiv:1305.6759 [gr-qc].
  * (40) V. Cardoso, I. P. Carucci, P. Pani, and T. P. Sotiriou, “Matter around Kerr black holes in scalar-tensor theories: scalarization and superradiant instability,” Phys. Rev. D 88 (2013) 044056, arXiv:1305.6936 [gr-qc].
  * (41) V. Cardoso, I. P. Carucci, P. Pani, and T. P. Sotiriou, “Black holes with surrounding matter in scalar-tensor theories,” Phys. Rev. Lett. 111 (2013) 111101, arXiv:1308.6587 [gr-qc].
  * (42) S. R. Dolan, “Superradiant instabilities of rotating black holes in the time domain,” Phys.Rev. D87 (2013) 124026, arXiv:1212.1477 [gr-qc].