# Dark photon superradiance: Electrodynamics and multimessenger signals

Nils Siemonsen<EMAIL_ADDRESS>Perimeter Institute for
Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada  Arthur B. McDonald
Canadian Astroparticle Physics Research Institute, 64 Bader Lane, Queen’s
University, Kingston, Ontario K7L 3N6, Canada  Department of Physics &
Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada  Cristina
Mondino<EMAIL_ADDRESS>Perimeter Institute for Theoretical
Physics, Waterloo, Ontario N2L 2Y5, Canada  Daniel Egaña-Ugrinovic
<EMAIL_ADDRESS>Perimeter Institute for Theoretical Physics,
Waterloo, Ontario N2L 2Y5, Canada  Junwu Huang<EMAIL_ADDRESS>Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
Masha Baryakhtar<EMAIL_ADDRESS>Physics Department, University of
Washington, Seattle, WA 98195-1560, USA  William E. East
<EMAIL_ADDRESS>Perimeter Institute for Theoretical Physics,
Waterloo, Ontario N2L 2Y5, Canada

(August 28, 2024)

###### Abstract

We study the electrodynamics of a kinetically mixed dark photon cloud that
forms through superradiance around a spinning black hole, and design
strategies to search for the resulting multimessenger signals. A dark photon
superradiance cloud sources a rotating dark electromagnetic field which,
through kinetic mixing, induces a rotating visible electromagnetic field.
Standard model charged particles entering this field initiate a transient
phase of particle production that populates a plasma inside the cloud and
leads to a system which shares qualitative features with a pulsar
magnetosphere. We study the electrodynamics of the dark photon cloud with
resistive magnetohydrodynamics methods applicable to highly magnetized plasma,
adapting techniques from simulations of pulsar magnetospheres. We identify
turbulent magnetic field reconnection as the main source of dissipation and
electromagnetic emission, and compute the peak luminosity from clouds around
solar-mass black holes to be as large as $10^{43}\,{\rm erg/s}$ for
observationally-allowed dark photon parameter space. The emission is expected
to have a significant X-ray component and to potentially be periodic, with
period set by the dark photon mass. The luminosity is comparable to the
brightest X-ray sources in the Universe, allowing for searches at distances of
up to hundreds of Mpc with existing telescopes. We discuss observational
strategies, including targeted electromagnetic follow-ups of solar-mass black
hole mergers and targeted continuous gravitational wave searches of anomalous
pulsars.

###### Contents

  1. I Introduction
  2. II Black hole superradiance for vector fields
  3. III Kinetically-mixed superradiance clouds: an overview
  4. IV Plasma production
    1. IV.1 Synchrotron radiation
    2. IV.2 Photon assisted Schwinger pair production
    3. IV.3 Dynamics leading to a quasi-steady state
  5. V Field configurations
    1. V.1 Plasma model
    2. V.2 Field solutions: Large scale behaviors
    3. V.3 Field solutions: Small scale turbulence
    4. V.4 Summary of turbulent plasma scales
  6. VI Electromagnetic emission
    1. VI.1 Radiation and dissipation channels
    2. VI.2 Power output
    3. VI.3 Dissipation mechanism
    4. VI.4 Periodicity of emission
    5. VI.5 Emission spectra
  7. VII Multimessenger Signals
    1. VII.1 Electromagnetic follow-ups of black hole mergers
    2. VII.2 Gravitational follow-ups of anomalous pulsars
    3. VII.3 Concluding remarks
  8. VIII Discussion and Conclusions
  9. A Notation
  10. B Construction of massive vector cloud
  11. C Numerical evolution setup
  12. D Resistive force-free currents
  13. E Charge distribution and small conductivity regime
  14. F Dark photon basis
  15. G Flux discussion

##  I Introduction

Ultralight fields arise in abundance in Beyond the Standard Model (SM)
theories of particle physics. The most well-known and well-motivated such
particle is the QCD axion Weinberg (1978); Wilczek (1978), proposed to solve
the discrepancy between the observed and predicted magnitude of the neutron
electric dipole moment arising from CP violation in the strong sector of the
SM Peccei and Quinn (1997). Beyond the QCD axion, light bosonic fields have
been found to be ubiquitous in string theory Arvanitaki  _et al._ (2010);
Svrcek and Witten (2006); Abel  _et al._ (2008); Goodsell  _et al._ (2009),
and provide excellent candidates for the dark matter particle or a dark matter
mediator Preskill  _et al._ (1983); Abbott and Sikivie (1983); Dine and
Fischler (1983); Graham  _et al._ (2016); Nelson and Scholtz (2011); Arias
_et al._ (2012); Essig  _et al._ (2013); Adams  _et al._ (2022); Antypas  _et
al._ (2022), making this class of particles one of the most exciting
candidates for new physics.

Black hole (BH) superradiance Zeldovich (1971); Misner (1972); Starobinskii
(1973) is a unique mechanism that enables searches for weakly interacting
ultralight bosons Arvanitaki  _et al._ (2010); Arvanitaki and Dubovsky (2011)
that relies only on the boson’s gravitational interaction. If a new light
boson with Compton wavelength of order the BH horizon size exists in the
theory—whether or not there is an initial abundance of the particle in the
environment—the BH will spin down and source macroscopic, coherent,
gravitationally-bound states of ultralight bosons Arvanitaki  _et al._ (2010);
Arvanitaki and Dubovsky (2011); Brito  _et al._ (2015a). These bosonic
“clouds” carry up to several percent of the BH’s initial mass, and have an
energy density comparable to that of neutron star matter for stellar mass BHs
Arvanitaki and Dubovsky (2011); East and Pretorius (2017); East (2018). The
resulting large energy density of the cloud has time-dependent components,
rotating around the BH axis at a frequency fixed by the particle mass,
resulting in coherent, monochromatic gravitational wave (GW) radiation that
depletes the cloud over parametrically longer times Arvanitaki and Dubovsky
(2011); Yoshino and Kodama (2014); Arvanitaki  _et al._ (2015); Goldreich and
Julian (1969).

The signatures of GW emission and BH spindown have been proposed to constrain
and search for ultra-light bosons Arvanitaki and Dubovsky (2011); Rosa and
Dolan (2012); Pani  _et al._ (2012a, b); Yoshino and Kodama (2014); Brito  _et
al._ (2013); Arvanitaki  _et al._ (2015); Brito  _et al._ (2015b, a);
Arvanitaki  _et al._ (2017); East (2017); Brito  _et al._ (2017a, b);
Baryakhtar  _et al._ (2017); East (2018); Baumann  _et al._ (2019); Siemonsen
and East (2020); Brito  _et al._ (2020); Zhu  _et al._ (2020). Bosons in the
$10^{-13}-10^{-11}$ eV range can lead to up to thousands of GW signals
originating from our Galaxy alone Arvanitaki  _et al._ (2015, 2017); Brito
_et al._ (2017a, b). Blind continuous wave searches for monochromatic GW from
scalar boson clouds Abbott  _et al._ (2022a, b); Palomba  _et al._ (2019);
D’Antonio  _et al._ (2018); Dergachev and Papa (2019); Zhu  _et al._ (2020),
as well as stochastic searches for an excess of GW power from spin-0 Tsukada
_et al._ (2019) and spin-1 Tsukada  _et al._ (2021) boson clouds around yet
undiscovered BHs have been carried out with LIGO-Virgo-KAGRA (LVK) Aasi  _et
al._ (2015); Acernese  _et al._ (2015); Aso  _et al._ (2013) data. These
searches have produced some constraints; however, a robust underlying BH natal
spin distribution is needed to conclusively exclude particle parameter space.
Another search strategy is to follow up BHs with a measured mass and spin,
which are newly born from binary BH mergers Arvanitaki  _et al._ (2017); then
the expected signals can be precisely computed, and a conclusive search
performed. Directed searches for continuous GWs from a potential scalar boson
cloud around Cygnus X-1 have also been carried out Sun  _et al._ (2020).
Currently, the sensitivity of GW searches is not sufficient to see follow-up
signals from spin-0 bosons around binary BH remnants Isi  _et al._ (2019), but
they are promising for spin-1 bosons in upcoming observation runs, and
especially in next-generation observatories Chan and Hannuksela (2022);
Siemonsen  _et al._ (2022).

BH spin measurements have set constraints on ultralight bosons using
measurements of BH properties from X-ray binaries Arvanitaki  _et al._ (2015);
Cardoso  _et al._ (2018); Baryakhtar  _et al._ (2021); Mehta  _et al._ (2020),
as well as measurements of binary BH constituents using LVK observations
Cardoso  _et al._ (2018); Ng  _et al._ (2021a, b). The latter have produced
constraints on spin-0 bosons of a factor of $\sim 2$ in mass Ng  _et al._
(2021a, b), and we expect slightly stronger constraints for spin-1 bosons,
although such an analysis has not been carried out for the full data set. The
X-ray binary measurements depend on BH accretion disk modeling, which may
introduce additional systematics Li  _et al._ (2005); Huang  _et al._ . In
addition, while the gravitational interaction of the BH superradiance cloud
and the accretion disk does not significantly affect the constraints
Arvanitaki  _et al._ (2015), non-gravitational interactions of the cloud can
perturb the disk dynamics, invalidating the constraints on the dark photon
mass from BH spin measurements for the parameters considered in this paper.

While the gravitational aspects of superradiance have been studied
extensively, making contact with particle physics models of ultralight spin-0
or spin-1 particles can dramatically change this picture. For spin-0 axions,
the relevant interaction at the next-order after the mass term is a quartic
coupling (see e.g. Arvanitaki and Dubovsky (2011); Yoshino and Kodama (2012,
2015); Gruzinov (2016); Fukuda and Nakayama (2020)). This results in energy
exchange between levels in the cloud Gruzinov (2016); Baryakhtar  _et al._
(2021); Omiya  _et al._ (2022), slowing down the spin extraction and resulting
in lower-frequency gravitational waves from transitions and axion wave
emission Baryakhtar  _et al._ (2021). More complicated dark sectors can result
in production of new dark states Baryakhtar  _et al._ (2017); Fukuda and
Nakayama (2020); Mathur  _et al._ (2020); East (2022); East and Huang (2022);
Cannizzaro  _et al._ (2022). Beyond interactions within the dark sector
itself, interactions with SM particles can lead to additional energy loss
channels Sen (2018), although for axion-like particles these are subdominant
to the dynamics of self-interactions Baryakhtar  _et al._ (2021).

Figure 1: (left) We show the visible magnetic field strength $B^{2}$
(normalized by its maximum) in the equatorial plane of the central BH of mass
$M$, dimensionless spin $a_{*}=0.86$, and a dark photon mass $\mu=0.3/(GM)$.
The dark photon of the superradiance cloud forces the pair plasma into a
circular motion resulting in magnetic field line twisting, which is released
through magnetic field line reconnection, resulting in a turbulent plasma
state (shown here) and efficient energy dissipation into the plasma, driving
the luminous electromagnetic emissions from the system. (right) Kinetically
mixed dark photon parameter space of interest in this work. The solid (dashed)
black lines are contours of constant electromagnetic luminosity emitted from
the superradiance cloud around a BH of mass 10 (100) $\mathrm{M}_{\odot}$ and
initial spin $a_{*}=0.9$. The region above the blue contour is relevant for
electromagnetic follow-ups of compact binary mergers, discussed in Sec. VII.1
(the shaded band on top of the contour is due to uncertainties on the merger
rate). The area within the green contour is of interest for continuous
gravitational waves searches targeted on anomalous pulsars, as described in
Sec. VII.2. The gray shaded region is excluded by existing measurements of the
CMB spectrum by COBE/FIRAS Fixsen  _et al._ (1996); Caputo  _et al._ (2020).

In this work, we focus on studying the effects of the lowest-order
interactions one can write down for spin-1 dark photons: kinetic mixing with
the SM photon Okun (1982); Holdom (1986). In the presence of such a mixing,
the huge energy density of the cloud picks up a visible electromagnetic field
component that interacts directly with electrons, leading to cascade
production of charged particles and to the formation of a plasma. To study the
plasma dynamics, we analyze an isolated, relativistic superradiance dark
photon cloud, and compute the evolution of the visible electric and magnetic
fields using a resistive-magnetohydrodynamic description, valid in the limit
of a strongly magnetized, tenuous plasma, that we adapt from simulations of
pulsar magnetospheres. See Fig. 1 for an example visualization of the
resulting magnetic field strength around a rotating BH.

Our simulations show that the resulting system is a luminous multimessenger
source: a BH system which emits an enormous electromagnetic flux, up to
several orders of magnitude brighter than pulsars and magnetars. This
radiation is generated by turbulent field and plasma dynamics in the
superradiance cloud, and is expected to have a large high-energy component. We
find partial evidence for an intrinsic periodicity set by the mass of the dark
photon particle, giving rise to a novel object that we call a “new pulsar.” In
Fig. 1, we show the parameter space of dark photon particles and the expected
peak luminosity for illustrative BHs as a function of dark photon mass and
kinetic mixing parameter.

Our results motivate a variety of novel astrophysical searches to discover
these systems. These include electromagnetic follow-ups of BH mergers that
result in rotating BHs, most promising in the X-ray and radio bands. Another
target is gravitational wave follow-ups of pulsars with coincident frequencies
or positive frequency drifts, which could be superradiance cloud signals
lasting thousands of years or more; see Fig. 1. In much of the parameter
space, we find that the evolution of the cloud is still dominated by its
gravitational dynamics, making the overall evolution free of electromagnetic
modeling uncertainties. At small dark photon masses and large kinetic mixings,
the electromagnetic emission has a larger power than the GWs, where we also
have an exceptionally bright sources.

Some aspects of dark photon superradiance with a non-zero kinetic mixing have
been explored in Caputo  _et al._ (2021), and superradiance of the SM photon
itself has been treated in Pani  _et al._ (2012b); Blas and Witte (2020);
Cannizzaro  _et al._ (2021a, b). Ours is the first work to consistently take
into account the dynamics of the SM plasma that is automatically generated by
the kinetically-mixed cloud. The interactions with the plasma completely alter
the behavior of the visible electromagnetic fields in the vicinity of the BH
and the resulting signatures.

This paper is organized as follows; in Sec. II, we review gravitational spin-1
superradiance. In Sec. III, we provide an executive summary of the dynamics of
kinetically mixed superradiance, which we explore in detail in the subsequent
sections. In Sec. IV, we describe the processes by which an isolated dark
photon cloud generates its own plasma density. In Sec. V, we study the
dynamics of the coupled system of electromagnetic fields and charged currents.
In Sec. VI, we describe the key electromagnetic emission mechanisms, including
electromagnetic radiation and power dissipation in the plasma due to turbulent
dynamics. In Sec. VII, we summarize the observational signatures and propose
several detection strategies for this new class of astrophysical object,
concluding and outlining future directions in Sec. VIII.

This work spans the areas of particle physics, strong field electrodynamics,
gravity, and astrophysical systems, thus introducing much notation, some non-
standard; we collect definitions in App. A. We describe details of the
numerical simulations of the superradiance cloud and electromagnetic fields
and currents in App. B and C, respectively. We present the resistive current
prescription in App. D, the small conductivity regime in App. E, and aspects
of the dark photon basis in App. F. We use the mostly-plus metric signature
$(-,+,+,+)$ and natural units, with $\hbar=c=1$ and non-reduced Planck mass
$M_{\mathrm{pl}}=1/\sqrt{G}$.

##  II Black hole superradiance for vector fields

We begin by reviewing the BH superradiance of a massive vector (spin-1) boson
that interacts predominantly through gravity. The kinetic and mass terms for
this dark photon $A^{\prime\mu}$ are given by

$\mathcal{L^{\prime}}=-\frac{1}{4}F^{\prime}_{\mu\nu}F^{\prime\mu\nu}-\frac{1}{2}\mu^{2}A^{\prime\mu}A^{\prime}_{\mu}.$
(1)

We assume that the dark photon mass $\mu$ arises from the Stueckelberg
mechanism 200 (2009) so in what follows, we do not discuss any dynamics that
could originate from a Higgs sector East (2022); East and Huang (2022).

The superradiant instability is a purely gravitational process that can lead
to the production of an exponentially large number of massive bosons around
spinning BHs by extracting the BH’s energy and angular momentum. The bosons
occupy hydrogenic clouds characterized by a gravitational fine-structure
constant $\alpha\equiv r_{g}\mu=GM\mu$, a principal quantum number $n$, and
total, orbital, and magnetic angular momentum numbers $j$, $l$, and $m$. The
total and orbital angular momentum can differ due to the boson’s intrinsic
spin: $j\in\\{l-1,l,l+1\\}$. Amongst the different cloud levels, the fastest-
growing one for vector bosons is the $(j,n,l,m)=(1,1,0,1)$ mode. Given its
dynamical dominance, we focus for brevity exclusively on the study of this
level. This will be sufficient for exploring the features that we wish to
highlight in this work. The boson’s energy in this level at leading-order in
the gravitational coupling is given by

$\omega\simeq\mu\bigg{(}1-\frac{\alpha^{2}}{2}\bigg{)}.$ (2)

After the birth of the source BH, the number of dark photons in the cloud
grows exponentially at a leading-$\alpha$ rate that, for our dynamically
dominant mode, is set by

$\Gamma_{\mathrm{SR}}\equiv\tau_{\mathrm{SR}}^{-1}\simeq
4\alpha^{7}(\Omega_{\mathrm{BH}}-\omega)\simeq 4a_{*}\alpha^{6}\mu,$ (3)

where $a_{*}$ is the BH’s dimensionless spin and $\Omega_{\mathrm{BH}}$ its
angular velocity

$\Omega_{\mathrm{BH}}=\frac{1}{2}\left(\frac{a_{*}}{1+\sqrt{1-a_{*}^{2}}}\right)r_{g}^{-1}.$
(4)

In the last equality of (3), we approximated $\Omega_{\mathrm{BH}}\gg\omega$
and took large-spin BHs $(1-a_{*}^{2})\ll 1$. If, on the other hand, the BH’s
spin is small so that its angular velocity falls below the boson’s energy

$\Omega_{\mathrm{BH}}\leq\omega,$ (5)

then $(j,n,l,m)=(1,1,0,1)$ superradiance does not occur. Equation (5) implies
a maximum possible value for the fine-structure constant (saturated for
maximally spinning BHs)

$\alpha\lesssim 1/2.$ (6)

The above condition, together with the strong suppression of the superradiant
rate at small $\alpha$ [see (3)] indicate that superradiance is most effective
for gravitational couplings of order $\alpha\sim 10^{-1}$ or, equivalently,
boson masses $\mu\sim 0.1/r_{g}$, which for stellar BHs corresponds to
$\mu\sim 10^{-12}\,\mathrm{eV}$.

The growth of the cloud stops when sufficient spin has been extracted so that
the condition Eq. (5) is saturated. The number of dark photons in the cloud
can reach $10^{77}$ or more for a $10$ solar mass BH, with the cloud mass

$M_{c}\simeq 10^{-2}\bigg{(}\frac{\Delta
a_{*}}{0.1}\bigg{)}\bigg{(}\frac{\alpha}{0.1}\bigg{)}M,$ (7)

for $\alpha\ll 1$, where $\Delta a_{*}$ is the difference between the initial
BH spin and the final spin which saturates the superradiance condition. The
cloud mass reaches up to $10\%$ of the mass of the BH for large $\alpha$ and
high initial BH spin East and Pretorius (2017). The vector field profile
around the BH, on the other hand, is given at leading order in the fine-
structure constant by

$\displaystyle A^{\prime}_{0}$
$\displaystyle=\frac{\sqrt{M_{c}}}{\sqrt{\pi}\mu^{2}r_{c}^{5/2}}e^{-r/r_{c}}\sin\theta\sin(\omega
t-\phi),$ $\displaystyle\boldsymbol{\mathbf{A^{\prime}}}$
$\displaystyle=-\frac{\sqrt{M_{c}}}{\sqrt{\pi}\mu
r_{c}^{3/2}}e^{-r/r_{c}}\\{\cos\omega t,\sin\omega t,0\\},$ (8)

where $r_{c}={r_{g}}/{\alpha^{2}}$ is the cloud’s characteristic Bohr radius
and we have taken the BH spin direction to lie along the $z$-axis. From Eq.
(8), we see that the dark electric and magnetic fields
$\boldsymbol{\mathbf{E^{\prime}}}\equiv-\nabla
A^{\prime}_{0}-\partial_{t}\boldsymbol{\mathbf{A^{\prime}}}$ and
$\boldsymbol{\mathbf{B^{\prime}}}\equiv\nabla\times\boldsymbol{\mathbf{A^{\prime}}}$
are in a proportion
$|\boldsymbol{\mathbf{B^{\prime}}}|/|\boldsymbol{\mathbf{E^{\prime}}}|\sim\alpha$,
so the cloud is electrically dominated. At leading order in $\alpha$, the
electric field corresponding to the potential Eq. (8) is unidirectional and
equatorially oriented, and rotates on this plane at a frequency
$\omega\simeq\mu$, while the magnetic field lines form concentric tori around
the BH with a common axis perpendicular to the electric field direction and
passing through the BH. Both fields decay exponentially away from the BH. We
show these features in Fig. 2, where we present exact (in the test field
limit) solutions for both the electric and magnetic fields, obtained by
numerically solving the vector’s equations of motion in the BH’s Kerr metric.
We refer the reader to Sec. V and App. B for details on the simulations. We
note that, close to the BH, the exact field solutions differ from the ones
obtained from the approximations (8) due to corrections that arise at higher-
order in the gravitational coupling.

Figure 2: We plot representative sets of field lines of the electric (left)
and magnetic fields (right) of the superradiant cloud around the central BH in
Kerr-Schild coordinates (see App. B for details). The $m=1$ cloud is
characterized by $\alpha=0.3$, while the BH has a corresponding spin of
$a_{*}=0.86$ (further details can be found in Table 3). The BH spin-axis
points in the $z$-direction. Color indicates the field strength along each
field line, normalized by the respective maximum field strength. On the right,
we also plot the magnetic field strength inside the equatorial plane of the
BH.

Following its formation, the cloud decays via GW emission, which is the main
observable signature of superradiance clouds composed of massive bosons that
interact with the SM solely by gravitation. The GW emission power is given by

$P_{\mathrm{GW}}\simeq
17\frac{\alpha^{10}}{G}\bigg{(}\frac{M_{c}(t)}{M}\bigg{)}^{2}$ (9)

in the $\alpha\ll 1$ limit Baryakhtar  _et al._ (2017); Siemonsen and East
(2020). The quadratic dependence of the emission power on the cloud mass leads
to a power-law decay of the cloud set by

$M_{c}(t)=\frac{M_{c}(t_{0})}{1+(t-t_{0})/\tau_{\mathrm{GW}}},$ (10)

where $\tau_{\mathrm{GW}}$ is the gravitational-wave decay timescale, which is
given by

$\tau_{\mathrm{GW}}\simeq\frac{GM}{17\alpha^{11}\Delta a_{*}}\sim
30\,\mathrm{days}\bigg{(}\frac{0.1}{\Delta
a_{*}}\bigg{)}\bigg{(}\frac{0.1}{\alpha}\bigg{)}^{11}\bigg{(}\frac{M}{10M_{\odot}}\bigg{)}.$
(11)

##  III Kinetically-mixed superradiance clouds: an overview

So far, we have discussed a theory where vector bosons lack non-gravitational
interactions. Going beyond this minimal setup, dark photons may interact with
the SM at the renormalizable level via kinetic mixing with the SM U(1) gauge
boson. In an effective theory below the electroweak scale, this interaction
mixes the dark and SM photons via a Lagrangian term
$\mathcal{L}\supset\varepsilon F^{\prime}_{\mu\nu}F^{\mu\nu}/2$, where
$\varepsilon$ is a parameter that quantifies the mixing Okun (1982); Holdom
(1986). This term can be equivalently written as a mass-mixing term by
performing the field redefinition, $A^{\prime}_{\mu}\rightarrow
A^{\prime}_{\mu}+\varepsilon A_{\mu}\equiv{A}^{\prime}_{\mu}$, which results
in the Lagrangian

$\displaystyle\begin{aligned} \mathcal{L}=&\
-\frac{1}{4}{F}_{\mu\nu}{F}^{\mu\nu}-\frac{1}{4}{F}^{\prime}_{\mu\nu}{F}^{\prime\mu\nu}\\\
&\
-\frac{\mu^{2}}{2}{A}^{\prime}_{\mu}{A}^{\prime\mu}-\varepsilon\mu^{2}{A}^{\prime}_{\mu}{A}^{\mu}+I_{\mu}{A}^{\mu},\end{aligned}$
(12)

where $I^{\mu}$ is the four-dimensional spacetime current. This choice of
fields is referred to as the interaction basis; other choices of basis are
discussed in App. F. Due to the mass mixing, the dark photon field acts as a
source current for the visible fields and vice-versa, as can be seen either
from the equations of motion

$\displaystyle\nabla_{\alpha}{F}^{\alpha\beta}$ $\displaystyle=$
$\displaystyle-I^{\beta}+\varepsilon\mu^{2}{A}^{\prime\beta},$ (13)
$\displaystyle\nabla_{\alpha}{F}^{\prime\alpha\beta}$ $\displaystyle=$
$\displaystyle\mu^{2}{A}^{\prime\beta}+\varepsilon\mu^{2}{A}^{\beta},$ (14)

at leading order in the kinetic mixing parameter, or from the energy-momentum
conservation relations, which manifestly show exchange of energy between the
vector fields

$\displaystyle\nabla_{\alpha}{T}^{\alpha\beta}=$ $\displaystyle\
-{F}^{\beta\gamma}(I_{\gamma}-\varepsilon\mu^{2}{A}^{\prime}_{\gamma}),$ (15)
$\displaystyle\nabla_{\alpha}{T}^{\prime\alpha\beta}=$ $\displaystyle\
\varepsilon\mu^{2}{F}^{\prime\beta\gamma}{A}_{\gamma}\,.$

Figure 3: Schematic depiction of the evolution of a kinetically-mixed dark
photon superradiance cloud. Starting from a spinning BH (left), a vector
superradiance cloud forms on a timescale $\tau_{\mathrm{SR}}$ (center-
left).The visible electric field sourced by the cloud accelerates
environmental charged particles, leading to cascade production of electrons
and positrons on a timescale $\tau_{\mathrm{plasma}}$; note that
$\tau_{\mathrm{SR}}\gg\tau_{\mathrm{plasma}}$ and the cascade production
occurs and completes before the superradiance instability completes (center-
right). The cloud finally decays by GW emission on a timescale
$\tau_{\mathrm{GW}}$, and by transferring energy to the plasma, which loses
energy through electromagnetic emission on a timescale
$\tau_{\mathrm{EM}}$(right). See text for further details.

In the context of superradiance, the kinetic mixing term allows for the
superradiance cloud to source electromagnetic fields111The induced coupling to
SM electrons also results in a higher-dimensional self-interaction term for
the dark photons of the Euler-Heisbenberg Lagrangian. Approximately
extrapolating the results of self-interacting scalars Baryakhtar  _et al._
(2021), we estimate that the induced quartic coupling would start to affect
the growth of the cloud for $\varepsilon\gtrsim\mathcal{O}(1)$, a value far
greater than relevant for the dynamics discussed here, and that is already
excluded.. Our objective in this work is to study the corresponding
electrodynamics, and we numerically solve Maxwell’s equations with a
superradiant source term in curved spacetime. This task is technically
complex, but most of our results can be understood in simple physical terms,
so to guide the reader through the discussion presented in the following
sections, we will begin here by providing a simplified overview of our
findings.

The evolution of a kinetically-mixed dark-photon superradiance cloud can be
separated into several stages that are schematically depicted in Fig. 3.
Starting with a spinning BH (leftmost panel), these stages correspond to the
initial growth of the cloud (center-left panel), creation of a conducting
plasma via particle acceleration and pair-creation due to the visible electric
field induced by the cloud (center-right panel), and the establishment of an
electromagnetic field and plasma configuration, which decays as the plasma
radiates electromagnetically and the cloud emits GWs (rightmost panel). The
initial growth of the cloud was reviewed in Sec. II for non-interacting dark
photons. From Eq. (14), we see that the inclusion of kinetic mixing affects
the dynamics of the superradiance cloud at order $\varepsilon$, which in turn
leads to effects at order $\varepsilon^{2}$ in the visible fields via Eq.
(13). Here, we will limit ourselves to computing the visible fields at leading
(linear) order in $\varepsilon$, so in what follows we ignore the effects of
kinetic mixing on the growth of the superradiance cloud. Thus, the evolution
of the cloud is governed purely by the gravitational dynamics as in the
previous section (with a growth timescale $\tau_{\mathrm{SR}}$ given by Eq.
(3)), and we may move on to the description of the creation of the plasma by
the cloud-induced electric field. See also the text below Eq. (27) for a
discussion of the negligible sub-leading corrections to the dark photon mass.

As discussed in Sec. II, the cloud is dominated by an equatorially-oriented
dark electric field, with a direction that rotates in the equatorial plane at
frequency $\mu/(2\pi)\sim 10^{2}\,\mathrm{Hz}\,(\mu/10^{-12}\,\mathrm{eV})$.
Due to kinetic mixing, the superradiance fields act as a source term in the
visible field equations of motion, Eq. (13), and induce a visible electric
field that is equal to the dark electric field times the mixing parameter. As
the cloud grows through superradiance, the visible field grows concurrently. A
fully-formed cloud would, in the absence of charged particles, have a visible
electric field of magnitude

$\displaystyle|\varepsilon\textbf{E}^{\prime}|\simeq\frac{\varepsilon\sqrt{\Delta
a_{*}}\,\alpha^{5/2}\mu}{\sqrt{G}}$ $\displaystyle\simeq 2\cdot
10^{13}\,\,\mathrm{V/m}\,\sqrt{\Delta
a_{*}}\bigg{(}\frac{\varepsilon}{10^{-7}}\bigg{)}\bigg{(}\frac{\alpha}{0.1}\bigg{)}^{5/2}\bigg{(}\frac{\mu}{10^{-12}\,\mathrm{eV}}\bigg{)}$
(16)

at distances of order the Bohr radius from the BH ($r_{c}\sim 1/\alpha\mu$).
These large fields, however, cannot be achieved due to plasma screening.

Before the cloud reaches its full size, and when visible fields are still only
a fraction of the value (16), environmental charged particles are accelerated
to ultra-relativistic velocities. As illustrated in the middle-right panel in
Fig. 3, the rotation of the electric field with the cloud curves the charged
particle trajectories, which then emit synchrotron photons. These photons, in
turn, interact non-perturbatively with the background electric field and
produce additional electrons and positrons. The charge acceleration and pair
production processes repeat in a cascade, until a conducting plasma is
created. This mechanism is reminiscent of cascade production of electron-
positron pairs in the strong magnetic field around supermassive BHs described
by Blandford-Znajek Blandford and Znajek (1977), with the important difference
that our system is electrically instead of magnetically dominated. An extended
discussion of the plasma creation will be presented in Sec. IV. The effects of
the plasma on the dark photon interactions with the visible fields are
discussed in App. F.

Once the plasma is created, the electrodynamics can be studied by encoding the
microscopic particle physics in an effective conductivity $\sigma$, which
allows for the computation of spatial plasma currents from the electromagnetic
fields using Ohm’s law,

$\displaystyle\mathbf{J}=\sigma(\mathbf{E}+\mathbf{v}\times\mathbf{B}),$ (17)

where $\mathbf{v}$ is the plasma fluid velocity. This Ohmic prescription,
which will be studied in detail in Sec. V, is commonly used in
magnetohydrodynamics Eckart (1940), and has been proposed to treat dissipative
currents in pulsars Komissarov (2006); Gruzinov (2007, 2008); Li  _et al._
(2012); Palenzuela (2013). Pulsar magnetospheres are highly conducting,
$\sigma/\omega\gg 1$ (where $\omega$ is the plasma angular frequency), so that
in the bulk of the magnetosphere electromagnetic fields are shorted out and
combine to cancel the Lorentz forces entering Eq. (17), resulting in a so-
called “force-free” system that is mostly dissipationless Goldreich and Julian
(1969); Blandford (2002); Meier (2004). Dissipative effects, however, do arise
on specific two-dimensional planes called current sheets of thickness $\sim
1/\sigma$, which may be modelled using Eq. (17).

Our kinetically-mixed superradiance cloud shares many similarities with
pulsars, and due to screening is also expected to have force-free regions. In
fact, from numerical simulations presented in Sec. V, we will see that if the
plasma is highly conducting, $\sigma/\mu\gg 1$, plasma currents effectively
redistribute charge to screen the rotating electric field induced by
superradiance $\varepsilon\textbf{E}^{\prime}$. This leads to charge being
separated into a dipole-like distribution, with a characteristic density that
at the cloud radius is approximately given by,

$\displaystyle\rho$
$\displaystyle\simeq\varepsilon\nabla\cdot\textbf{E}^{\prime}\simeq\pm\frac{\varepsilon\sqrt{\Delta
a_{*}}\alpha^{7/2}\mu^{2}}{\sqrt{G}}$ (18) $\displaystyle\simeq\pm 5\cdot
10^{7}\textrm{cm}^{-3}\,\sqrt{\Delta
a_{*}}\bigg{(}\frac{\varepsilon}{10^{-7}}\bigg{)}\bigg{(}\frac{\alpha}{0.1}\bigg{)}^{7/2}\bigg{(}\frac{\mu}{10^{-12}\,\mathrm{eV}}\bigg{)}^{2},$

where the plus and the minus signs correspond to opposite ends of the dipole-
like pattern. A crude, non-relativistic estimate indicates that the large
magnitude of the charge density is consistent with large conductivities: for a
non-relativistic collisionless plasma the conductivity is
$\sigma\simeq\rho/\mu m_{e}$, which, using Eq. (18), gives $\sigma/\mu\sim
10^{12}$, where $m_{e}$ is the electron mass.

Despite the utility of the pulsar analogy, the resemblance with our system is
limited. First, in the absence of a plasma, the kinetically mixed
superradiance cloud is electrically, instead of magnetically, dominated.
Second, while in a pulsar the magnetic field is dipolar and decays away from
the neutron star, in Sec. V we will show that in our system the resulting
visible fields remain strong well outside of the light-cylinder
$r_{\mathrm{LC}}\equiv 1/\mu$ (the radius out to which the plasma can corotate
with the BH), up to the Bohr radius $r\gtrsim 1/\mu\alpha\gg r_{\mathrm{LC}}$.
This means that in the bulk of our system, charges cannot move fast enough to
perfectly screen the rotating source field. Electric dominance and imperfect
screening suggest that in our cloud a steady-state force-free solution does
not exist, unlike in pulsars where dissipative effects are confined to the
current sheets. Instead, our numerical simulations, presented in Sec. V, show
a dynamical interplay between resistive and force-free regions where electric
fields have been mostly screened.

Up to now, we have only discussed the electric field dynamics. Complementary
insight into the electrodynamics can be gained by studying instead the
magnetic fields induced by the plasma currents. The magnetic dynamics can be
analyzed by combining Ohm’s law Eq. (17), Faraday’s, and Ampere’s law (derived
from Eq. (13)) to obtain the magnetic induction equation

$\displaystyle\partial_{t}\boldsymbol{\mathbf{B}}=\frac{\varepsilon\mu^{2}\boldsymbol{\mathbf{B^{\prime}}}}{\sigma}+\frac{1}{\sigma}\nabla^{2}\boldsymbol{\mathbf{B}}+\nabla\times(\boldsymbol{\mathbf{v}}\times\boldsymbol{\mathbf{B}}),$
(19)

where the first term on the right hand side accounts for the background
superradiant magnetic field. The induction equation is used to study magnetic
fields in a wide variety of astrophysical plasmas, where the electric
displacement currents are smaller than plasma currents and can be neglected,
an assumption that our simulations show to be valid. The three terms on the
right-hand side of the induction equation describe different characteristic
regions of the system. Closest to the BH, the superradiant driving field is
large, and the first term on the right-hand side dominates the morphology of
the magnetic fields. Away from the neighborhood of the BH, the magnetic field
is non-trivially related to the superradiant driving fields, and the two
remaining terms become dynamically relevant in a proportion set by a magnetic
Reynolds number $R_{m}=\sigma|\textbf{v}|\ell$, where $1/\ell$ characterizes
the magnetic field gradients. In zones where $R_{m}\gg 1$, the last term on
the right-hand side of Eq. (19), which represents pure field advection, is
largest. The simulations presented in Sec. V show that large regions in the
bulk of the plasma are dominated by advection, and are characterized by
magnetic flux conservation, tight-coupling of the plasma and the magnetic
fields, mostly screened electric fields, and some emission of electromagnetic
radiation due to the time-dependent plasma charge and current densities.

Our simulations also show time-dependent regions, especially outside of the
light-cylinder, where the plasma cannot corotate with the driving fields as
advection would impose. This leads to differential rotation within the plasma
and to the twisting and shearing of magnetic field lines, as well as to
regions where the second (diffusive) term in Eq. (19) dominates due to large
field gradients and/or small plasma velocities that result in $R_{m}\ll 1$. In
these regions, we find that the interplay of advection and diffusion drives
turbulent effects, such as breaking and reconnection of field lines,
schematically shown in the rightmost panel of Fig. 3. Unscreened electric
fields along the direction of plasma currents, expected from the simple
kinematic arguments outlined above, are found at these sites. These electric
fields lead to significant Ohmic dissipation $\textbf{J}\cdot\textbf{E}$,
which in our simulation represent conversion of electromagnetic field energy
into particle acceleration and radiation.

Our simulations thus show that dissipation is associated with magnetic field
reconnection and unscreened electric fields, as in the pulsar current sheets.
In contrast to the pulsar system, however, in the kinetically mixed
superradiance cloud the resistive effects are realized in dynamically evolving
regions throughout the bulk of the plasma. As a consequence, while in pulsars
most of the power emitted is due to the time-dependent nature of the currents
in the force-free bulk, we find that in our system the comparative
preponderance of dissipative effects leads to an emission power that is
dominated by Ohmic losses.

Importantly, the simulations presented in Sec. VI.2 suggest that while
increasing the conductivity reduces the size of the dissipative regions, it
also increases their number, i.e. larger conductivities lead to
“fragmentation” of the dissipative regions without changing their volumetric
fraction. As a result, we find that the dissipative power tends to a
$\sigma$-independent value (at large $\sigma$), allowing us to provide a
prediction for the emitted power that is set entirely by the dark-photon model
parameters and the BH mass. The power emitted by our system typically exceeds
the emission power of pulsars by several orders of magnitude, and for clouds
around stellar BHs can be as large as $L\simeq 10^{43}$ erg/s (an exact
expression can be found in Eqns. (54) or (61)). We ascribe this difference to
the rapid falloff of the dipolar magnetic field of the pulsar away from the
neutron star, the large volume of the superradiance cloud when compared with
the pulsar’s emission regions
$V_{\mathrm{cloud}}/V_{\mathrm{lightcylinder}}\sim 1/\alpha^{3}$, and to our
system’s dissipative features. Given the periodic rotation of the cloud, it is
possible that the emitted power will have a pulsating component, and our
simulations indeed show some limited evidence that supports this hypothesis
(see Sec. VI.4). From our simulations we cannot compute the spectral
decomposition of the emitted power; however, we can speculate based on results
of kinetic treatments of turbulent plasmas (analogous to pulsar current sheet
simulations) that charged particles will be highly boosted by the large
electric fields resulting in a large component of high-energy radiation in the
form of X- and gamma-rays Werner  _et al._ (2016); Zhdankin  _et al._ (2017);
Cerutti  _et al._ (2015, 2016); for further discussion see Sec. VI.5.

In the final stage of our system’s evolution the cloud decays predominantly by
gravitational-wave emission, accompanied by the novel electromagnetic emission
outlined here, as depicted in the rightmost panel of Fig. 3. These emission
channels lead to concrete observational signatures that we describe in Sec.
VII, such as performing EM follow-up observations of compact binary mergers,
searching for a population of same-frequency and/or positive-frequency drift
pulsars, and targeting such anomalous pulsars with GW follow-up searches.

With this short summary in hand, we now move on to provide an in-depth
discussion of the plasma and field dynamics at each stage of their evolution,
starting with the production of the plasma.

##  IV Plasma production

In this section we describe the production of the conducting plasma within the
superradiance cloud and determine the values of the mixing parameter
$\varepsilon$ for which a plasma is plausibly created. We identify two main
processes that are crucial for the formation of the plasma, synchrotron
radiation emitted by environmental electrons that are accelerated by the
superradiance cloud, and subsequent photon-assisted Schwinger pair production
in the background electric field. Here we estimate the rates of these two
processes and show that they can effectively create the conducting plasma even
for kinetic mixing parameters that are several orders of magnitude below
current experimental bounds and of the region of interest for the
observational prospects discussed later in this work. Several other mechanisms
can produce charged particles in background fields and additionally contribute
to the formation of the plasma, but for brevity we do not discuss them here
(for a comprehensive list we refer the reader to Mészáros (1992)).

###  IV.1 Synchrotron radiation

Any stray charged particle entering the kinetically mixed superradiance cloud
will experience strong electromagnetic forces. Since the cloud’s magnetic
field is subdominant,
$|\boldsymbol{\mathbf{B}}^{\prime}|\sim\alpha|\boldsymbol{\mathbf{E}}^{\prime}|<|\boldsymbol{\mathbf{E}}^{\prime}|$,
we simply consider the motion of $e^{\pm}$ accelerated by the electric field
(equivalently we can perform a boost into a frame with vanishing magnetic
field and electric field amplitude reduced by a factor of $1-\alpha\simeq 1$).
Inside a fully grown cloud (at distance $r$ such that $r_{g}\ll r\lesssim
1/\alpha\mu$) the electric field has approximately constant amplitude given by
Eq. (16) and rotates with angular velocity $\omega\simeq\mu$. The
electrons/positrons are then approximately linearly accelerated over a time
scale of $1/\mu$, reaching a maximum boost factor of

$\displaystyle\gamma_{e}$
$\displaystyle\simeq\frac{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}{m_{e}\mu}\simeq
e\varepsilon\alpha^{5/2}\sqrt{\Delta a_{*}}\frac{M_{\rm pl}}{m_{e}}$
$\displaystyle\simeq
10^{12}\left(\frac{\varepsilon}{10^{-7}}\right)\left(\frac{\alpha}{0.1}\right)^{5/2}\left(\frac{\Delta
a_{*}}{0.1}\right)^{1/2},$ (20)

where in the first line we made use of Eq. (16). As the electric field
rotates, the electron/positrons trajectories bend with approximate radius of
curvature
$r_{c}\simeq\gamma_{e}m_{e}/(e\varepsilon|{\boldsymbol{\mathbf{E}}^{\prime}}|)\simeq
1/\mu$. During this circular motion the charged particles radiate synchrotron
photons, predominantly at frequency $\omega_{\rm
syn}=\gamma_{e}^{3}/r_{c}\simeq\gamma_{e}^{3}\mu$. We can estimate the rate
for synchrotron emission at this frequency as $P_{\rm syn}/\omega_{\rm syn}$,
which gives

$\Gamma_{\rm
syn}\simeq\frac{2}{3}\frac{e^{2}\gamma_{e}}{r_{c}}\simeq\frac{2}{3}\frac{e^{3}\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}{m_{e}}\simeq\frac{2}{3}e^{3}\varepsilon\,\alpha^{5/2}\sqrt{\Delta
a_{*}}\mu\frac{M_{\textrm{Pl}}}{m_{e}},$ (21)

where again, in the last equality we used (16). In order for the plasma to be
phenomenologically relevant it must be created before the cloud is depleted by
gravitational wave emission, on a timescale given by Eq. (11). To ensure that
this occurs, we impose the sufficient requirement that the synchrotron and
photon-assisted pair production rates (discussed in the next section) occur
before any particle can escape the cloud, i.e., that the synchrotron and pair-
production timescales are shorter than the light-crossing time of the cloud,
$1/\alpha\mu$, which is much smaller than the GW decay time (11). For the
synchrotron emission rate of Eq. (21), this leads to the requirement

$\varepsilon>\frac{1}{e^{3}\,\alpha^{5/2}\sqrt{\Delta
a_{*}}}\frac{m_{e}}{M_{\textrm{Pl}}}\simeq
10^{-18}\left(\frac{0.1}{\alpha}\right)^{5/2}\left(\frac{0.1}{\Delta
a_{*}}\right)^{1/2}$ (22)

As evident from Eq. (IV.1), the above requirement also ensures that the
accelerated electrons are highly relativistic ($\gamma_{e}\gg 1$).

###  IV.2 Photon assisted Schwinger pair production

A static electric field can decay to electron-positron pairs through quantum
tunneling, a process known as Schwinger pair production. The probability of
scalar $e^{\pm}$ pair creation was first computed in Schwinger’s seminal work
Schwinger (1951), from vacuum decay in an external, slowly varying electric
field E. The rate per unit volume $V$ is given by

$\frac{\Gamma_{e^{\pm}}}{V}=\frac{(e|\textbf{E}|)^{2}}{4\pi^{3}}\sum_{n}^{\infty}\frac{1}{n^{2}}\exp{\left(-\frac{\pi
m_{e}^{2}}{e|\textbf{E}|}n\right)},$ (23)

and is exponentially suppressed for electric fields below the critical value
$m_{e}^{2}/e\simeq 10^{18}\ \mathrm{V/m}$. Even the large electric field
generated by the dark photon superradiance cloud, given in Eq. (16), falls
short by a few orders of magnitude, making Schwinger pair production unlikely
in our setup. However, pair creation can be greatly enhanced in the presence
of highly energetic photons Schutzhold  _et al._ (2008), such as the
synchrotron photons described in the previous section. Photon assisted
Schwinger pair production is similar to magnetic pair production Erber (1966)
invoked in Blandford-Znajek processes Blandford and Znajek (1977), where
radiation with energy above the threshold $2m_{e}$ can produces electron-
positron pairs by scattering off of strong magnetic field.

Photon-assisted Schwinger pair creation can be viewed as a semiclassical
tunneling process and the production rate has been computed with methods
similar to the one used for metastable vacuum decay in Ref. Coleman (1977).
The exponential factor in the rate is given by $e^{-S_{B}}$, where $S_{B}$ is
the Euclidean action evaluated on the bounce solution (the classical
trajectory that extremizes the action). In our case Monin and Voloshin (2010);
Dunne  _et al._ (2009)

$S_{\rm
B}=-\gamma_{\theta}+\left(\frac{2m_{e}^{2}}{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}+\frac{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}{2m_{e}^{2}}\gamma_{\theta}^{2}\right)\arctan\frac{2m_{e}^{2}}{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|\gamma_{\theta}},$
(24)

where $\gamma_{\theta}=\sin\theta\
m_{e}\omega_{\gamma}/(e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|)$,
$\omega_{\gamma}$ is the photon frequency, and $\theta$ is the angle between
the direction of the photon and the background electric field. As the
synchrotron photons travel in the cloud, they will encounter electric fields
that are almost perpendicular to their direction of propagation within a time
scale of $1/\mu$, when the production rate is maximized. As a result, for a
simple estimate of the rate we can take $\sin\theta\simeq 1$.

There are two limiting cases of Eq. (24) depending on the photon frequency. If
$\omega_{\gamma}\ll 2m_{e}$, $S_{B}\simeq\pi
m_{e}^{2}/(e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|)$, which reduces to
the standard Schwinger result of Eq. (23). Therefore if the photon energy is
below the pair production threshold, the electric field still needs to be
super-critical for the process not to be exponentially suppressed. On the
other hand, if $\omega_{\gamma}\gg 2m_{e}$, as is the case for most of the
synchrotron photons described in the previous section, $S_{B}\simeq
2m_{e}^{3}/(e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|\omega_{\gamma})$.
The additional, potentially very small, factor of $2m_{e}/\omega_{\gamma}$
significantly enhances the probability of pair production. The full expression
of the rate (including the prefactor of the exponential term) is given in Ref.
Dunne  _et al._ (2009) and is larger for photons with polarization
perpendicular to the electric field. For the highly energetic synchrotron
photons with perpendicular polarization, we have

$\displaystyle\Gamma^{\gamma}_{e^{\pm}}$ $\displaystyle=\frac{\alpha_{\rm
EM}}{2\pi}\frac{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}{m_{e}}\exp\left(-\frac{2m_{e}^{2}}{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}\frac{2m_{e}}{\omega_{\rm
syn}}\right)$ $\displaystyle=\frac{\alpha_{\rm
EM}}{2\pi}\frac{e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|}{m_{e}}\exp\left[-\frac{4m_{e}^{6}\mu^{2}}{(e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}|)^{4}}\right].$
(25)

For cascade production to occur in our system the term in the exponential must
reach a magnitude of order unity when the cloud has reached its full size (or
before) so that the exponential suppression of pair-production is lifted. This
translates into a minimum value for the mixing parameter $\varepsilon\gtrsim
m_{e}^{3/2}\mu^{1/2}/e|\boldsymbol{\mathbf{E^{\prime}}}|$. More precisely, in
what follows we impose that the pair-production rate for a fully grown cloud
is faster than the light crossing time of the cloud $1/\alpha\mu$, so that the
synchrotron photons split into $e^{\pm}$ before escaping the superradiance
cloud, ensuring a cascade production of the plasma, which translates into a
minimal mixing parameter

$\displaystyle\varepsilon\gg$ $\displaystyle\left(\ln\frac{\alpha_{\rm
EM}\gamma_{e}}{2\pi\alpha}\right)^{-1/4}\frac{m_{e}^{3/2}\sqrt{2\mu}}{e|\boldsymbol{\mathbf{E^{\prime}}}|}$
$\displaystyle\approx$ $\displaystyle\left(\ln\frac{\alpha_{\rm
EM}\gamma_{e}}{2\pi\alpha}\right)^{-1/4}\frac{1}{e\alpha^{5/2}\sqrt{\Delta
a_{*}}}\frac{m_{e}}{M_{\rm pl}}\sqrt{\frac{2m_{e}}{\mu}}\simeq
10^{-10}\left(\frac{0.1}{\alpha}\right)^{5/2}\left(\frac{0.1}{\Delta
a_{*}}\right)^{1/2}\left(\frac{10^{-12}\ {\rm
eV}}{\mu}\right)^{1/2}\left(\frac{\log\frac{\alpha_{\rm
EM}\gamma_{e}}{2\pi\alpha}}{20}\right)^{-1/4},$ (26)

where in the logarithmic term we made use of (IV.1), and in going to the
second line we used Eq. (16). Notice that the above requirement also
guarantees that $\omega_{\rm syn}\gg 2m_{e}$. For kinetic mixing parameters
saturating the lower bound in Eq. (IV.2), plasma production will be triggered
when the cloud has reached a close-to maximal size, while for mixing
parameters above this lower bound the plasma will be created before the cloud
has fully grown. We show the smallest values of kinetic mixing parameters that
allow for cascade pair creation in the superradiance cloud as a function of
dark photon mass in the left panel of Fig. 4.

Figure 4:  (left) Smallest values of the kinetic mixing parameter
$\varepsilon$ that allows for efficient $e^{\pm}$ pair production in the
superradiance cloud as a function of dark photon mass $\mu$, for BH masses of
$10\,M_{\odot}$ (blue) and $100\,M_{\odot}$ (orange), with initial BH spin of
$a^{*}=0.9$. The rate for photon stimulated Schwinger pair production, given
in Eq. (25), is required to be greater than the size of the cloud,
$\alpha\mu$, when the cloud has fully grown. Smaller dark electric fields at
small $\mu$, require larger $\varepsilon$ to initiate the cascade. The sharp
cutoff corresponds to the highest dark photon mass that satisfies the
superradiance condition for the fastest growing level. (right) Ratio of the
time needed to populate the plasma over the superradiance e-folding time as a
function of the kinetic mixing parameter $\varepsilon$ and the gravitational
coupling $\alpha$ for a BH mass of $10\ M_{\odot}$ and initial BH spin of
$a^{*}=0.9$ (the ratio is independent of $M$ and only mildly dependent on
$a^{*}$). An estimate of the ratio is given in Eq. (27), while in the plot
$\tau_{\rm{plasma}}$ is evaluated using the electric field value at the time
that the cascade pair production is initiated. In the dark gray shaded region,
the electric field is always too small to produce $e^{\pm}$. When the cascade
is efficient, the plasma is filled within a small fraction of one
superradiance e-folding time ($\tau_{\rm{SR}}$ grows steeply at small
$\alpha$). In both panels, the light gray shaded region is excluded by
measurements of the CMB spectrum by COBE/FIRAS Fixsen  _et al._ (1996); Caputo
_et al._ (2020)

###  IV.3 Dynamics leading to a quasi-steady state

The plasma begins to be populated once the pair-production cascade initiates,
which as noted previously happens when the superradiance cloud has grown to a
size such that the pair-production rate becomes of the order of the cloud’s
Bohr radius, i.e. $\Gamma^{\gamma}_{e^{\pm}}\simeq\alpha\mu$. For this to
occur and up to a logarithmic correction, the superradiant field must reach a
critical value
$e\varepsilon|\boldsymbol{\mathbf{E}}^{\prime}_{\mathrm{crit}}|\simeq\sqrt{2}m_{e}^{3/2}\mu^{1/2}$.
The superradiant field grows to
$\boldsymbol{\mathbf{E}}^{\prime}_{\mathrm{crit}}$ in a few superradiance
times $\tau_{\mathrm{SR}}\simeq 1/4\alpha^{6}a^{*}\mu$. After that, the plasma
is created by cascade production on the much shorter light-crossing timescale,
over which the superradiant field and the cascade production rates are
approximately fixed to
$\boldsymbol{\mathbf{E^{\prime}}}\simeq\boldsymbol{\mathbf{E}}^{\prime}_{\mathrm{crit}}$
and $\Gamma^{\gamma}_{e^{\pm}}\simeq\alpha\mu$ respectively.

During cascade production, the charge density grows exponentially as
$n_{e}=n^{0}_{e}\exp\left(2\Gamma^{\gamma}_{e^{\pm}}t\right)$, where
$n_{e}=n^{-}_{e}+n^{+}_{e}$ is the total number density of electrons and
positrons. Pair production stops when the charged plasma effectively screens
the critical electric field due to charge separation, which happens when the
electron number density reaches
$n^{f}_{e}\simeq\varepsilon\nabla\cdot\boldsymbol{\mathbf{E}}^{\prime}_{\mathrm{crit}}/e\simeq\sqrt{2}\alpha(m_{e}\mu)^{3/2}/e^{2}$.
Assuming that when the cascade begins we start from one single electron in the
cloud, $n^{0}_{e}\simeq(\alpha\mu)^{3}$, the plasma grows
$n^{f}_{e}/n^{0}_{e}$ e-folds before the cascade stops, so that the plasma
formation time can be estimated as

$\displaystyle\tau_{\rm plasma}\simeq$
$\displaystyle\frac{1}{2\Gamma^{\gamma}_{e^{\pm}}}\ln\bigg{(}\frac{n^{f}_{e}}{n^{0}_{e}}\bigg{)}\simeq\frac{1}{2\alpha\mu}\ln\frac{m_{e}^{3/2}}{e^{2}\alpha^{2}\mu^{3/2}}\quad.$

The plasma production time is thus parametrically shorter than the
superradiance timescale by a factor

$\displaystyle\frac{\tau_{\rm plasma}}{\tau_{\rm SR}}$
$\displaystyle\simeq{2a^{*}\alpha^{5}}\ln\frac{m_{e}^{3/2}}{e^{2}\alpha^{2}\mu^{3/2}}$
(27)

In the right panel of Fig. 4 we show the above ratio of timescales for the
values of $\varepsilon$ and $\alpha$ that satisfy the plasma pair production
requirement from Eq. IV.2.

As the dark photon cloud continues to grow and the electric field increases,
more charged particles will be created and the plasma will rearrange itself in
the screening configuration, until $\boldsymbol{\mathbf{E}}^{\prime}$ has
reached its maximum value Eq. (16) after $\approx 180\ \tau_{\rm SR}$, at
which point the charge density achieves its maximal value Eq. (18). The
numerical simulation presented in the next sections will show how charge
separation in the cloud and electric screening are indeed good approximations.
Note that the formation of the plasma induces a plasma mass for the SM photon,
but does not significantly affect the dark photon mass. In fact, in the limit
$\mu\rightarrow 0$, a massless mode must remain in the theory even in the
presence of the plasma, which indicates that the leading contribution to the
dark photon mass is simply $\mu$, up to $\varepsilon^{2}$ corrections. A more
detailed discussion of plasma effects on the dark photon is presented in App.
F. Note also that the total mass of the plasma,

$\displaystyle M_{\rm
plasma}\approx\frac{m_{e}\rho}{e(\alpha\mu)^{3}}\simeq\varepsilon\alpha^{1/2}(\Delta
a_{*})^{1/2}\frac{m_{e}M_{\rm pl}}{\mu}\simeq
10^{-29}M_{\odot}\left(\frac{\epsilon}{10^{-7}}\right)\left(\frac{\alpha}{0.1}\right)^{1/2}\left(\frac{\Delta
a_{*}}{0.1}\right)^{1/2}\left(\frac{10^{-12}\ {\rm eV}}{\mu}\right)$ (28)

is much smaller than the mass of the cloud, $M_{c}=\alpha^{2}\Delta
a_{*}M_{\rm pl}^{2}/\mu$, leaving the gravitational potential unaltered. We
can then safely assume that the growth and dynamics of the superradiance cloud
is not affected by the presence of the standard model plasma.

In the left panel of Fig. 4, we show the range of $\varepsilon$ and $\mu$
where fast cascade production occurs and the plasma is populated,
$\Gamma_{e^{\pm}}^{\gamma}>\alpha\mu$. Much slower processes can also populate
this plasma in the parameter spaces where the fast cascade production is
inactive. Assuming, for example, Bondi accretion and $\mathcal{O}(1)$ sound
speeds $c_{s}$, it takes roughly

$\displaystyle\tau_{\rm acc}=\frac{M_{\rm plasma}}{\dot{M}_{\rm
Bondi}}\approx\varepsilon\alpha^{-3/2}(\Delta
a_{*})^{1/2}\frac{c_{s}^{3}M_{\rm pl}\mu}{\pi n_{M}}\simeq 10\ {\rm
years}\left(\frac{\varepsilon}{10^{-7}}\right)\left(\frac{\alpha}{0.1}\right)^{-3/2}\left(\frac{\Delta
a_{*}}{0.1}\right)^{1/2}\left(\frac{\mu}{10^{-12}\ {\rm
eV}}\right)\left(\frac{1/{\rm
cm^{3}}}{n_{M}}\right)\left(\frac{c_{s}}{1}\right)^{3}$ (29)

to populate enough charged particles inside the superradiance cloud, where
$n_{M}$ is the matter number density in the interstellar medium. Such a time
scale suggests that for parameters where the cascade production is active,
accretion from interstellar medium can be safely ignored during dark photon
superradiance and plasma generation. On the other hand, in the parameter space
where the cascade production is inactive, such processes can be fast enough to
populate a plasma inside the superradiance cloud around isolated BHs inside
our galaxy, which are most likely more than thousands if not millions of years
old.

Finally, it is worth pointing out that the transient process of cascade
particle production discussed in this section only produces a small and
unobservable amount of emission. The transient effects discussed in Caputo
_et al._ (2021), for example, occur when $|\boldsymbol{\mathbf{E}}|^{2}\sim
m_{e}^{3}\mu$, which amounts to a total energy of
$10^{32}(10^{-12}\leavevmode\nobreak\ \mathrm{eV}/\mu)^{2}(0.1/\alpha)^{3}$
ergs (independent of $\varepsilon$). As we will demonstrate in section VI.2,
this is about 20 orders of magnitude smaller than the total
electromagnetically dissipated power from the superradiance cloud. Similarly,
photon superradiance Pani  _et al._ (2012b); Blas and Witte (2020); Cannizzaro
_et al._ (2021a, b) will saturate long before the field strengths and energy
densities in the photon superradiance cloud reaches sizes relevant for
observation.

##  V Field configurations

We have shown that a pair production cascade ensues on short timescales, once
the superradiance dark photon cloud surpasses a critical electric field
strength during the exponential growth of the cloud. The generated charges
form a highly conducting plasma that is subject to the electromagnetic fields
of the dark photon cloud. In this section, we study the macroscopic state this
plasma equilibrates into. To that end, we consider, numerically, the
superradiance cloud of a kinetically mixed dark photon on a fixed Kerr
spacetime of mass $M$ and dimensionless spin parameter $a_{*}$. As outlined
throughout secs. III and IV, it is crucial to work with an Ohm’s law that
accounts for the energy dissipation into the plasma inside the cloud. Building
on the analytical discussions in the previous section, we first examine the
structure of the visible electromagnetic fields of the cloud-plasma system on
large scales. On small scales, we demonstrate that, in the high conductivity
limit, turbulent dynamics emerges accompanied by efficient magnetic field line
reconnection in the bulk of the dark photon cloud. Throughout the section, we
study the system considering all relativistic effects of the background
spacetime in the interaction basis.

We begin by introducing the plasma model in the context of a kinetically mixed
massive vector field with the SM photon in Sec. V.1. We establish that within
this model, a quasi-stationary end state of the pair production cascade is
reached and characterize its large scale behavior in Sec. V.2, and small-scale
dynamics in Sec. V.3. In Sec. V.4, we briefly summarize the main findings.

###  V.1 Plasma model

We study the kinetically mixed field equations (13) of the visible field
${A}_{\mu}$ in the interaction basis (12). A strongly magnetized, highly-
conducting pair plasma, is well-described by the force-free limit of ideal
magnetohydrodynamics. Specifically, force-free electrodynamics is applicable
when (i) resistive effects are negligible, (ii) the magnetic field strength is
larger than the electric field’s, i.e., ${B}^{2}>{E}^{2}$, and (iii) the
plasma mass density $\rho_{p}$ and pressure $P_{p}$ are much smaller than the
electromagnetic energy density, i.e., ${B}^{2}\gg\rho_{p},P_{p}$ (see e.g.,
Goldreich and Julian (1969); Blandford (2002)). In force-free electrodynamics,
one can numerically evolve only the electromagnetic fields, without keeping
track of the fluid quantities, with the current being uniquely determined from
the electromagnetic fields by the requirement that the Lorentz force vanishes.
In essence, $F_{\alpha\beta}I^{\alpha}=0$ allows
$\nabla_{\beta}F^{\alpha\beta}=I^{\alpha}$ to be rewritten as
$F_{\alpha\gamma}\nabla_{\beta}F^{\alpha\beta}=0$. (More details, include the
equations of force-free electrodynamics in terms of electric and magnetic
fields can be found in App. D.)

As discussed in the previous sections, the superradiant system considered here
satisfies condition (iii) [see (28)], but a priori violates conditions (i) and
(ii). Exactly how and why these violations persist, even at large
conductivities is the subject of the following subsections. However, one can
already anticipate that in the limit of vanishing backreaction of the pair
plasma onto the visible electromagnetic fields, the magnetic dominance is lost
by virtue of the superradiant solution being electrically dominated:
${B}^{2}=\varepsilon^{2}B^{\prime 2}<\varepsilon^{2}E^{\prime 2}={E}^{2}$.
Since electric dominance implies that there is no frame where the electric
field, and hence the acceleration on any charges, vanishes, it is furthermore
not surprising that dissipative processes also become important. Therefore, in
order to relax assumptions (i) and (ii), we modify the force-free equations to
explicitly allow for dissipation by introducing an Ohm’s law with finite
conductivity $\sigma$, while still requiring the current to be a function of
the electromagnetic fields. This is done in such a way that, under a certain
set of assumptions, the force free limit can be recovered as
$\sigma\rightarrow\infty$. Slightly abusing terminology, we shall refer to
this as resistive force-free electrodynamics. Several variations of this
approach have been applied to simulating pulsar magnetospheres Komissarov
(2006); Gruzinov (2007); Li  _et al._ (2012); Parfrey  _et al._ (2017);
Mahlmann  _et al._ (2021), and in particular the electrically dominated
current sheet Li  _et al._ (2012); Mahlmann  _et al._ (2021), and here we
generalize these to the kinetically mixed case (see also App. D for more
details). In the following, we will describe how we set up the coordinate
system, the Maxwell equations, the fluid of charged particles and the Ohm’s
law in a fully relativistic simulation of the plasma of electrons and
positrons inside the superradiance cloud.

#### Space-time decomposition:

We begin by discussing how the kinetically-mixed Maxwell’s equations for the
visible fields can be decomposed according to a given choice of time slicing
of a spacetime, which for this study will be given by the Cartesian Kerr-
Schild coordinate time $t$. Using the (future-pointing) unit normal to slices
of constant time $n^{\mu}$, the visible electric and magnetic fields defined
with respect to this slicing are

$\displaystyle E^{i}\equiv n_{\nu}F^{i\nu},$ $\displaystyle B^{i}\equiv
n_{\nu}(*F)^{i\nu}=\frac{1}{2}n_{\nu}\varepsilon^{i\nu\alpha\beta}F_{\alpha\beta}.$
(30)

where $\varepsilon^{\alpha\beta\gamma\delta}$ is the Levi-Civita tensor. These
are the quantities that we evolve on the BH spacetime. Projecting the
kinetically mixed Maxwell equations (13) into components orthogonal and
parallel to the time slice, one obtains the evolution equations in terms of
three-dimensional spatial quantities. We give the explicit form of these
equations that we use to carry out the numerical evolution in App. C.

In order to evolve the electric and magnetic fields, we also need to specify
the electromagnetic current. The four dimensional current $I^{\mu}$ can be
decomposed into a spatial component $J^{i}$, and a component perpendicular to
the slices of constant time

$\displaystyle J^{i}=I^{i}-\rho_{q}n^{i},$
$\displaystyle\rho_{q}=-n_{\mu}I^{\mu},$ (31)

where $\rho_{q}$ is the Eulerian frame (i.e., with respect to the slices of
constant time) charge density. We directly calculate the charge density from
the divergence of the electric field,

$\displaystyle\rho_{q}=$ $\displaystyle\
D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu},$ (32)

using the Gauss’s law constraint equation obtained from projecting the
kinetically mixed Maxwell equations for the visible fields (13) onto the time
slice. However, we still need to specify the spatial part of the current
$J^{i}$.

#### Ohm’s law:

As discussed above, we introduce the effect of finite conductivity using a
simple Ohm’s law in the frame of the fluid (plasma)

$\displaystyle j^{\mu}=\sigma e^{\mu},$ (33)

where $j^{\mu}$ and and $e^{\mu}$ are, respectively, the electromagnetic
current and visible electric field in the fluid frame. The latter is defined
in terms of the fluid velocity $u^{\mu}$,

$\displaystyle e^{\mu}\equiv u_{\nu}F^{\mu\nu},$ $\displaystyle j^{\mu}\equiv
I^{\mu}+(u^{\nu}I_{\nu})u^{\mu},$ (34)

in an analogous way to the Eulerian frame quantities. With this prescription,
we are neglecting anisotropic magnetic field effects. However, as the
superradiance cloud system is characterized by strong electric fields, this is
the dominant contribution away from the force-free limit. Note that the Ohm’s
law (33) instantaneously relates currents and electric fields, a prescription
that is valid for low-inertia plasmas such as ours. In eq. (33), the
conductivity $\sigma$ is a phenomenological parameter that allows for energy
dissipation via mechanisms that are set by the microphysics, which is left
unspecified. While without a microphysical description it is not possible to
determine the value of the conductivity, we expect that the conductivity of
the pair plasma considered here is large when measured in terms of the
system’s natural length scale (the inverse dark-photon mass), i.e.,
$\sigma/\mu\rightarrow\infty$ [see also a brief discussion below eq. (18)].
Note that in this limit, other charge transport mechanism such as diffusion
due to charge gradients can be safely neglected.

The resistive relation (33) allows us to compute the fluid-frame currents from
the visible electromagnetic fields. However, since we are not directly
evolving the fluid, and in particular its velocity, this is not sufficient to
give the Eulerian current entering into the evolution equations. In ideal
magnetohydrodynamics (including the limiting case of force-free), given the
electromagnetic fields, one can reconstruct one component of the velocity,
referred to as the drift velocity

$\displaystyle v^{i}_{d,{\rm
ideal}}=\frac{\varepsilon^{ijk}E_{j}B_{k}}{B^{2}},$ (35)

where here we use $v^{i}$ to refer to the Eulerian spatial velocity of the
fluid $v_{i}=-u_{i}/(n_{\mu}u^{\mu})$, though not the component of the
velocity parallel to the magnetic field. In our resistive extension, following
Gruzinov (2008) (see also Ref. Komissarov (2006)), we use a drift velocity
that is augmented with an electric component to allow for electrically
dominated regions to be treated self-consistently,

$\displaystyle v^{i}_{d}$
$\displaystyle=\frac{\varepsilon^{ijk}E_{j}B_{k}}{B^{2}+E_{0}^{2}},\quad
E_{0}^{2}=B_{0}^{2}+E^{2}-B^{2},$ $\displaystyle B_{0}^{2}$
$\displaystyle=\frac{1}{2}\left[B^{2}-E^{2}+\sqrt{(B^{2}-E^{2})^{2}+4(E_{i}B^{i})^{2}}\right].$
(36)

This ensures that the drift velocity is bounded by the speed of light, even in
electrically dominated regions, and is further reduced in regions with non-
vanishing resistive effects, i.e., $E_{i}B^{i}\neq 0$. We note that, while the
quantities above are written in terms of $E$ and $B$ fields,
$B^{2}-E^{2}=F_{\mu\nu}F^{\mu\nu}/2$ and
$E_{i}B^{i}=F_{\mu\nu}(*F)^{\mu\nu}/4$ are spacetime scalars.

In general, the full fluid velocity cannot be determined from the
electromagnetic fields alone, without extra conditions. Here, following Li
_et al._ (2012), we identify $v^{i}=v_{d}^{i}$, i.e., we set the non-drift
velocity component to zero in the BH frame defined by $n^{\mu}$. With the
fluid velocity specified, we can transform (33) into the Eulerian frame,
giving

$\displaystyle J^{i}=\rho_{q}v^{i}_{d}+\sigma
E_{0}\sqrt{\frac{B^{2}+E_{0}^{2}}{B_{0}^{2}+E_{0}^{2}}}\left(\frac{E_{0}E^{i}+B_{0}B^{i}}{B^{2}+E_{0}^{2}}\right).$
(37)

This is the kinetically mixed extension of the current considered in Li  _et
al._ (2012). In App. D, we further discuss the advantages and limitations of
this choice of current, and contrast it with other currents developed in the
literature. Notice, however, this ansatz is merely a prescription to extend
the regime of validity of the force-free paradigm to resistive regions (such
as current sheets), which has physical meaning only in the high-conductivity
limit, where $B^{2}>E^{2}$, $E_{i}B^{i}=0$, and $\sigma E_{0}<\infty$, such
that $v_{d}^{i}\rightarrow v_{d,\text{ideal}}^{i}$. Hence, while different
choices of currents may result in different physics at moderate conductivity,
for $\sigma\rightarrow\infty$ all prescripts converge towards the well-defined
force-free limit. Therefore, we primarily focus on the trend towards this
limit.

#### Numerical setup:

The visible fields $E^{i}$ and $B^{i}$ are numerically evolved on a fixed Kerr
spacetime in Cartesian Kerr-Schild coordinates. The massive vector field
solutions on this background spacetime, which enter as source terms in the
evolution equations, are constructed numerically (this is discussed in App.
B). We restrict our attention to solutions where the superradiant instability
has been saturated, which occurs when the cloud oscillation frequency is
synchronized with the horizon frequency $\omega_{R}=\Omega_{\rm BH}$. The
computational domain extends from the BH horizon to spatial infinity through
the use of compactified coordinates, and mesh refinement is used to
concentrate resolution in the central region, enabling us to resolve the BH-
cloud system sufficiently up to $r=10r_{c}$. The conductivity $\sigma$ is set
to be constant in space and time, as is serves merely as a proxy for the local
resistivity present in the cloud. In light of the lack of a microphysically
motivated conductivity, and the success of analogous choices in the case of
the pulsar magnetosphere (see e.g., Li  _et al._ (2012)), this choice, while
unphysical, is a first step towards a more complete analysis. Finally, we
evolve the system of equations forward in time using a higher-order explicit
Runge-Kutta algorithm. The limitations of this choice in the context of stiff
equations in the large conductivity limit, as well as further details of the
numerical methods, are discussed in App. C. We evolve the system, starting
from suitably chosen initial data, for a sufficiently long period of time
($\sim 200/\mu$ or longer) such that it relaxes towards an approximate
equilibrium, as measured by the Poynting flux at large radii becoming nearly
constant. Details are given in App. C.

###  V.2 Field solutions: Large scale behaviors

Figure 5: We show the visible electric field geometry, $E^{i}$, and magnitude
as a function of the plasma conductivity $\sigma$ in the equatorial plane of
the BH. The superradiance cloud’s phase is the same in each panel. Color
(red/blue) correspond to the magnitude of the component along the spin axis,
i.e., in the $z$-direction, normalized by the maximal magnitude of that
component at the given conductivity. Field lines are projections of the
electric field onto the equatorial plane, while the color of the field lines
(yellow/green) indicates the magnitude of the visible electric field
normalized by the maximal magnitude at $\sigma=0$. The BH and cloud parameters
are as in Figure 2, i.e., $\alpha=0.3$ and $a_{*}=0.86$. With increasing
conductivity, the electric field magnitude decreases compared with the vacuum
limit, and $E^{z}/E\sim\mathcal{O}(1)$ for $\sigma/\mu=20$. The field geometry
undergoes a phase-shift of $\pi/2$ between vacuum and large conductivity
limits.

In the following, we demonstrate how in the quasi-steady state solution, the
strong electric field of the superradiance cloud is mostly screened by locally
produced charges at large conductivities, and the visible magnetic field
begins to play an important role in the system. As noted above, the physical
value of the plasma conductivity $\sigma$ will be set by a microphysical scale
(due to scattering, synchrotron radiation, pair production and annihilation,
etc.) that we expect to be much smaller than the other physical length or time
scales in the system that we consider (e.g., $1/\mu$). However, due to
numerical limitations, we only consider values up to $\sigma/\mu=20$.
Therefore, in order to make qualitative and quantitative statements about the
properties of the quasi-endstate, we proceed by discussing the behavior of the
system as a function of conductivity, focusing primarily on $\sigma/\mu\geq 1$
(though we include the low-conductivity limit in App. E for completeness), and
extrapolate trends towards $\sigma\rightarrow\infty$ if possible. This is the
approach typically used to study resistive effects in pulsar magnetospheres Li
_et al._ (2012); Parfrey  _et al._ (2017); Mahlmann  _et al._ (2021); we
discuss possible shortcomings of these methods below.

We begin by considering the behavior of the visible electric field $E^{i}$ of
the quasi-stationary endstate of the pair cascade as function of conductivity
in Figure 5. In the vacuum limit, $\sigma/\mu=0$, there is no electric field
generated by a charged plasma, and the depicted field lines are just an
equatorial slice of the electric field $E^{i}=\varepsilon E^{\prime i}$ shown
in Figure 2. For $\sigma/\mu\geq 1$, the main qualitative feature, as
summarized in Sec. III, is the buildup of the dipolar screening charge density
leading to a significant reduction of the visible electric field compared to
the vacuum case. As the conductivity increases, this suppression of $E^{2}$
grows and the component of $E^{i}$ along the BH spin-axis becomes more
important, as $E^{z}/E\sim\mathcal{O}(1)$ in the $\sigma/\mu=20$ panel of
Figure 5.

Another important qualitative feature we find is that the visible electric
field exhibits a global de-phasing of $\pi/2$ with respect to the dark
electric field at large conductivities. This de-phasing can be understood
analytically in the non-relativistic limit, $\alpha\ll 1$. In this limit, a
spatial derivative, which is set by the inverse cloud size
$(\nabla\sim\alpha\mu)$, is much smaller than the time derivative
$(\partial_{t}\sim\mu)$. For electric fields with similar or larger strengths
compared to the magnetic field, as indicated by our numerical simulations at
small and moderate conductivities (see Figure 7 below), the Maxwell equations
in the interaction basis (73) reduce to

$\displaystyle\partial_{t}{E}^{i}\approx-\sigma
E^{i}+\varepsilon\mu^{2}A^{\prime i},$ (38)

where we assume the same Ohm’s law as before in (33). With the non-
relativistic superradiant field solutions $A^{\prime}_{i}$ in (8), we see that
the visible electric field is driven towards the inhomogeneous solution

$\displaystyle\textbf{E}\propto\begin{pmatrix}\sigma+i\mu\\\ \mu-i\sigma\\\
0\end{pmatrix}e^{-i\omega t}+c.c.$ (39)

over timescales of $1/\sigma$. Hence, as conductivity increases towards
$\sigma\sim\mu$, the visible electric field direction rotates with respect to
a fixed superradiance cloud phase, such that at large conductivities,
$\sigma/\mu\gg 1$, the visible and superradiant electric fields exhibit a
phase-offset of $\pi/2$. Notice, however, at very large conductivity, this
analytic approximation, in principle, is no longer valid, as it neglects the
magnetic field effects, which become important as the electric field is
screened, for $\sigma/\mu\gg 1$. These effects induce the appearance of small
scale structures visible in the last panel of Figure 5. We will elaborate on
the break down of this approximation and the emergence of small scale features
in more detail in Sec. V.3.

Figure 6: Magnetic field lines $B^{i}$ and magnitudes in a coordinate slice
spanned by the BH spin (pointing in the $z$-direction), and an arbitrarily
chosen superradiance cloud phase. The six panels show the field configurations
in the same slice for successively larger conductivities $\sigma$. The
background colors (red/blue) indicate the magnitude of the component
perpendicular to the slice $B^{\varphi}$ in the $\varphi$-direction around the
BH normalized by the magnitude of the visible magnetic field. The colors of
the field lines (yellow/green) indicate the magnitude of the visible magnetic
field along the field lines normalized by the maximal magnitude in the vacuum
case $B_{\rm max}^{\sigma=0}$. The BH and cloud parameters in all panels are
as in Figure 2, i.e., $\alpha=0.3$ and $a_{*}=0.86$. The magnitude of the
magnetic field, while exponentially decaying in the vacuum limit, is roughly
uniform at large conductivities $\sigma/\mu=20$. The small-scale features are
discussed in detail in Sec. V.3

As the conductivity grows, the electric field decreases in amplitude and the
magnetic field plays a more important role. The field line geometry and
magnitude of the visible magnetic fields are shown as a function of
conductivity in Figure 6. At vanishing $\sigma$, the solutions are identical
to the vacuum solution, and the snapshot in Figure 6 simply represents a slice
of the geometry shown in Figure 2. The magnetic null line (i.e., where
$B^{i}=0$) crosses this slice once on either side of the BH as they spiral
away from the BH. The field lines close around this null line, and the
magnetic field strength is largest close to the BH and decays exponentially
towards spatial infinity. In the vacuum limit, the dark photon field exhibits
an exact helical symmetry about the BH spin-axis222The helical Killing field
is $k^{\mu}=\xi^{\mu}-\Omega_{\rm BH}\varphi^{\mu}$ in terms of the stationary
and axisymmetric Killing vectors. The superradiant field $A^{\prime}_{\mu}$
strictly retains this symmetry, at leading order in the kinetic mixing, i.e.
$\mathcal{L}_{\textbf{k}}A^{\prime}_{\mu}\sim\mathcal{O}(\varepsilon^{2})$..
For $\sigma=\mu$, the magnetic field pattern still exhibits this helical
symmetry approximately on the spatial scales depicted in Figure 6. This
symmetry is broken for $\sigma/\mu>1$. Therefore, the last two snapshots of
Figure 6, while representing the magnetic field geometry qualitatively, are
not indicative of the full three-dimensional field geometry. Qualitatively, at
large conductivities, the plasma turns into a highly conducting pair plasma
attaching to the visible magnetic field lines. This implies a differential
rotation of the magnetic field lines at large and small distances from the BH,
which breaks the helical geometry into a more complex configuration presenting
small-scale features, that are further discussed in the next section. For both
the electric and magnetic field, the presence of the plasma leads to field
strengths that are relatively uniform in magnitude within the Bohr radius (as
can be seen in the $\sigma/\mu=20$ panel of Figure 5). The small-scale
features, crucial for the high-conductivity dynamics of the system, are
discussed in detail in the next section.

Figure 7: (top row) The ratio between visible magnetic and electric field
strengths, $B^{2}$ and $E^{2}$, respectively, as function of conductivity
$\sigma$ in the same slices as in Figure 6 (i.e., spanned by the BH spin-axis
in the z-direction and an arbitrary superradiant phase). Contour lines
indicate, where $B^{2}/E^{2}=1$. (bottom row) The magnitude of the visible
electric field component in the direction of the visible magnetic field,
$|E_{i}B^{i}|$, normalized by both magnitudes. The slices of the top and
bottom rows are identical. The BH and cloud parameters in all panels are as in
Figure 2, i.e., $\alpha=0.3$ and $a_{*}=0.86$. For $\sigma/\mu\lesssim 1$, the
electric field is dominant everywhere and the violations of $|E_{i}B^{i}|=0$
is strong, while for $\sigma/\mu\gtrsim 1$, the magnetic field begins to
dominate in some regimes and $|E_{i}B^{i}|=0$ is violated only in isolated
regions.

In order to understand the degree to which our solutions approach a force-free
solution in the $\sigma\rightarrow\infty$ limit, we consider how violations of
the force-free conditions $E_{i}B^{i}=0$ and $B^{2}>E^{2}$ (respectively,
conditions (i) and (ii) discussed in the beginning of subsection V.1), change
with increasing conductivity. In Figure 7, we show, in representative slices,
pointwise measures of the violations of these conditions, while in Figure 8 we
show how volume integrated measures of these violations decrease with
increasing conductivity.

Examining a volume integral of $E_{i}B^{i}$ as a function of conductivity,
shown in rightmost panel of Figure 8, we find that it begins to decrease like
$1/\sigma$ for $\sigma/\mu\geq 1$. From the bottom panels of Figure 7, we can
see that, in contrast to low and moderate value of $\sigma/\mu$, at high
conductivity, large values of $E_{i}B^{i}$ (relative to the magnitude of the
fields) occur only in isolated, smaller-scale regions.

We also find that the fraction of the volume that is magnetically dominated
increases with increasing conductivity, as shown in the middle panel of Figure
8, in particular, for a coordinate sphere of radius $4/(\alpha\mu)$. For
$\sigma=0$, none of this volume is magnetically dominated (as expected, since
$B^{\prime}\sim\alpha E^{\prime}$, $B=\varepsilon B^{\prime}$, and
$E=\varepsilon E^{\prime}$), while for $\sigma/\mu=20$, approximately one-
third of the volume is. The spatial extent of these magnetically dominated
regions can be seen in the top panels of Figure 7. Similarly, as shown in the
left panel of Figure 8, the global maximum visible electric and magnetic field
strengths are comparable, while the value of $\max B^{2}/E^{2}$ increases to
large values with large conductivities.

Figure 8: (left) We plot the maximal ratio of visible magnetic to electric
field magnitudes $\max B^{2}/E^{2}$, and the maximal magnetic and electric
field magnitudes, $B_{\rm max}^{2}$ and $E_{\rm max}^{2}$, normalized by their
maximal vacuum values, as a function of plasma conductivity
$\sigma/\mu\in\\{0.001,0.01,0.1,1,2,5,10,20\\}$. (middle) The fractional
coordinate volume $\mathcal{V}_{B^{2}>E^{2}}$ of magnetically dominated
regions inside a coordinate sphere of radius $4/(\alpha\mu)$ around the
central BH as function of conductivity. (right) We show the behavior of the
volume integral of $|E^{i}B_{i}|$ over a coordinate sphere of radius
$10r_{c}$, $\mathcal{I}_{EB}$, as a function of conductivity, normatlized to
its vacuum value $\mathcal{I}_{EB}^{\sigma=0}$. As above, we consider here a
BH-cloud system with $\alpha=0.3$ and spin $a_{*}=0.86$ in all panels.

Naïvely, one might expect that in the infinite conductivity limit, the visible
electric field in the fluid frame is completely shorted out by large scale
charge separation, leading to $e^{\mu}\rightarrow 0$, and hence a magnetically
dominated solution everywhere (recalling that magnetic dominance is equivalent
to the existence of a frame where the electric field vanishes, i.e.,
$B^{2}/E^{2}>1$). The results above suggest that a force-free solution might
exist for $\sigma\rightarrow\infty$, at least in a significant fraction of the
total volume taken up by the plasma and superradiance cloud. However, given
the slow increase in the magnetically dominated fraction of the volume with
increasing conductivity, it seems plausible that electrically dominated
regions with non-zero volume may persist as $\sigma\rightarrow\infty$.
Electric dominance implies an unscreened electric field in the fluid frame,
allowing for strong particle acceleration, and the dissipation of field
energy. This is consistent with our force-free simulations (without a guide
field), which always evolve towards developing electrically dominated regions
(see Figure 22 in App. D for a discussion on these force-free simulations),
and will be discussed further in the following section. In the next section,
we will also discuss how the violations of the force-free conditions are
connected to the turblent behavior the plasma.

###  V.3 Field solutions: Small scale turbulence

In the previous section, we found that a large-scale charge separation screens
the superradiant electric field, lifting the importance of the visible
electric field and leading to magnetically dominated regions. In what follows,
we focus on the magnetic field dynamics in the large-conductivity regime. We
show that, in this regime, the dark photon superradiance cloud-plasma system
is characterized by turbulent plasma dynamics in the bulk of the system. A
trend towards small-scale features can already be seen in Figure 6, and
becomes more apparent in Figure 7 in the previous section. The turbulent
regions emerge not in isolated and clearly structured lower-dimensional
regions, but rather across the bulk of the cloud. This is in contrast to the
pulsar magnetosphere, where, at least in the high conductivity limit, small
scale features and dissipation are expected to be contained mostly in a two-
dimensional current sheet.

In order to understand the turbulent behavior of the visible magnetic field
dynamics at moderate and large conductivities, it is instructive to consider
the visible magnetic induction equation in the presence of a finite kinetic
mixing with the massive vector field. For simplicity, we focus on the flat
spacetime limit only, noting that the following qualitative arguments are
unchanged on curved backgrounds. Furthermore, for clarity, we assume that the
fluid (i.e., the plasma) is mostly non-relativistic, and that at large
$\sigma$, the plasma is conduction (as opposed to advection) dominated, as
explicitly shown in App. E. Finally, we can neglect the displacement current
$\partial_{t}E^{i}\sim\mu E^{i}$, as it is suppressed compared to the
conduction term $\sim\sigma E^{i}$, if $\sigma\gg\mu$. All qualitative
arguments outlined below translate to the fully-relativistic case. Making
these assumptions, the evolution equations of the (kinetically mixed) Maxwell
equations (73) together with the current (33) reduce to the visible magnetic
induction equation

$\displaystyle\partial_{t}{B}^{i}=\frac{\varepsilon\mu^{2}B^{\prime
i}}{\sigma}+\frac{1}{\sigma}\partial_{j}\partial^{j}{B}^{i}+\varepsilon^{ijk}\varepsilon_{klm}\partial_{j}v^{l}{B}^{m},$
(40)

where $B^{\prime i}=\varepsilon^{ijk}\partial_{j}A^{\prime}_{k}$ and $v^{i}$
is the fluid velocity that we identify with the drift velocity
$v^{i}=v^{i}_{d}$, such that $v^{i}{E}_{i}=0$. For later convenience, we
define the Cartesian Kerr-Schild coordinate radius
$\hat{\rho}=(x^{2}+y^{2}+z^{2})^{1/2}$ (see App. C for details). With a
characteristic length scale $\ell=1/\mu$ of the system, we are able to define
an effective magnetic Reynolds number, $R_{m}=\ell v_{d}\sigma=\sigma
v_{d}/\mu$, of the effective plasma defined in (37)333One could instead chose
$\ell=r_{c}$. However, the qualitative arguments are unaffected by the precise
choice of magnetic Reynolds number.. For $R_{m}\gg 1$, the magnetic field
dynamics is dominated by the third term in Eq. (40), which represents field
advection. In this regime the magnetic field dynamics is entirely determined
by the plasma, as both are strongly coupled (i.e., the magnetic field is
comoving with the plasma). For $R_{m}\ll 1$, the second term in the induction
equation is most important, which accounts for magnetic field diffusion. In
this case the magnetic field decouples from the plasma motion and relaxes to a
diffusive state. The effects of the superradiant driving fields are included
explicitly in the first term of Eq. (40) and implicitly in the plasma
velocity, which depends on the driving electric fields, and in the visible
magnetic field themselves, which are sourced by electrically induced currents.

To illustrate the different domains of the magnetic field dynamics, in Figure
9 we show the magnetic field magnitude inside the equatorial plane as a
function of conductivity. In the vacuum limit, $\sigma/\mu=0$, the magnetic
field coincides with the superradiant magnetic driving field shown in Figure
2. With increasing conductivity, i.e., $\sigma/\mu\gtrsim 1$, from Figure 9 we
see that the regions where the morphology of the magnetic field resembles the
superradiant magnetic fields are confined to distances from the BH that are
smaller than a characteristic radius $r_{*}$, which we heuristically find to
be

$r_{*}\approx 80\mu GM/\sigma.$ (41)

Inside this critical radius, the superradiant driving field $B^{\prime i}$ is
exponentially large, and the first term of (40) dominates, compared with the
terms it induces in the diffusion and advection contributions444Note, the
superradiant components of the visible magnetic fields in the second and third
term in the induction equation (40), are $\alpha^{2}$ and $\alpha$ suppressed,
respectively, to the leading contribution at intermediate conductivities..
Note also that at large conductivities an overall phase-offset of $\pi/2$
between the superradiant and visible magnetic fields appears, similar to the
behavior of the electric field case (see Figure 5); at the level of the
magnetic field, this phase-offset emerges from the first term in the induction
equation (40).

For radii larger than $r_{*}$ we see that the magnetic fields are non-
trivially related to the driving superradiant electric and magnetic fields. At
large distances from the BH and on scales of order $1/\mu$, the magnetic field
strength increases with growing conductivity. On smaller scales, on the other
hand, and especially at large conductivities, we see that a series of small
scale features appear. These features, which arise on scales of order
$1/\sigma$, emerge from the interplay of the diffusion and advection
components of the magnetic induction equation. At large conductivities and due
to the large plasma velocities, most regions with $\hat{\rho}>r_{*}$ are
advection dominated. As a result, the oscillating superradiant driving fields
source visible magnetic fields that couple strongly to the plasma, while
conversely, the plasma cannot back react onto the superradiant driving fields.
Our results show that for radii $\hat{\rho}>r_{*}$, and on spatial scales of
the entire cloud (i.e., $>1/\mu$), the plasma is unable to corotate with the
driving fields; hence, differential rotation between plasma cells at different
radii occurs. This shear velocity results in twisting of the magnetic field
lines on scales of the cloud. The twist builds up magnetic energy that is
dissipated through turbulent magnetic reconnection in regions of large
magnetic diffusion, i.e., $R_{m}\ll 1$ (as we illustrate in detail below).
Thus, the small scale features in Figure 9 are a result of this turbulent
reconnection. The radius $r_{*}$ may saturate at the light cylinder of the
system for $\sigma\rightarrow\infty$ as inside it the plasma could rotate
rigidly with the superradiant driving fields. However, we do not find the
light cylinder to be a location of special importance for the largest
conductivity that we considered: $\sigma/\mu=20$.

Figure 9: We plot the visible magnetic field strength $B^{2}$ in the
equatorial plane of the system in the vacuum limit, $\sigma=0$, as well as at
moderate to high plasma conductivities, i.e., $\sigma/\mu\gtrsim 1$. The BH
and cloud parameters are as in Figure 2, i.e., $\alpha=0.3$ and $a_{*}=0.86$.
The superradiance cloud phase is identical in each of the panels. The color is
normalized by the maximal visible magnetic field strength at each
conductivity. The white dashed line indicates the critical coordinate radius
$r_{*}=80\mu GM/\sigma$, discussed in the main text. The region
$\hat{\rho}<r_{*}$ is dominated by superradiant driving, while the regions
with $\hat{\rho}>r_{*}$ are characterized by an interplay of advective and
diffusive regions. The flat spacetime light cylinder for this system is
roughly $R_{\rm LC}=GM/\alpha\approx 3.33GM$. Notice, the resolution of our
numerical methods decreases with increasing coordinate distances $|x|$ and
$|y|$, resulting in, for instance, a suppression of small-scale features in
the $\sigma/\mu>2$ cases for $|x|,|y|>50GM$.

Another qualitative feature is the disappearance of the magnetic null line of
the superradiant magnetic field for moderate conductivities, $\sigma/\mu>5$,
outside the critical radius $r_{*}$, as can be seen in Figure 9. Vanishing
magnetic field strength implies vanishing plasma bulk velocity, i.e.,
$R_{m}\approx 0$, and equivalently the presence of strong magnetic diffusion.
Hence, we find that the magnetic null line is quickly filled by magnetic field
lines diffusing into the null line from surrounding areas with finite magnetic
field strength. This process efficiently removes the null line outside the
critical radius $r_{*}$.

In addition to the large scale differential rotation about the BH, in
advection dominated regions (where the plasma and magnetic field are co-
moving) we observe localized roughly uniform oscillatory motion of the plasma
within the equatorial plane (with oscillation radius given by $1/\mu$), as
well as periodic longitudinal motion of the plasma along the BH spin-
axis555Even in the force-free simulations, this periodic oscillatory motion of
the magnetic field strength in the equatorial plane can be observed.. This
periodic motion is likely driven by the large scale superradiant electric
field throughout the plasma, in conjunction with large scale charge separation
of the pair plasma for $\sigma/\mu\gtrsim 1$. Charges are accelerated along
the large scale superradiant electric field. However, the orbital frequency
$\mu$ of the field’s direction forces the charges into a circular trajectory
with radius of $1/\mu$. This is reflected in the circular motion of features
in the magnetic field of scale $1/\mu$ inside the equatorial plane. The
circular motion of negative and positive charges is exactly out of phase by
$\pi$, resulting in out-of-phase oscillatory motion of the plasma on either
side of the BH due to the large scale charge dipole screening the superradiant
electric field. The longitudinal periodic motion along the spin axis is more
complex, and likely a result of the electric field driving within the
equatorial plane. We will discuss the observational consequence of this motion
of charge densities in Sec. VI.4.

Figure 10: We plot a selection of visible magnetic field lines of the
superradiance cloud-plasma system with conductivity $\sigma/\mu=20$,
$\alpha=0.3$, and $a_{*}=0.86$. In panels (a), (b) and (c) we also plot the
visible magnetic field strength within the equatorial plane, while in panel
(d) we show the visible magnetic field strength in the plane spanned by the BH
spin and an arbitrary superradiance cloud phase. We discuss this plot in
detail in the main text. The main macroscopic scales involved are the BH-
scale, set by the mass $M$, the superradiance cloud’s oscillation timescale
$1/\mu\approx 3.33GM$, and the cloud’s Bohr radius $r_{c}\approx 11.1GM$.

Let us illustrate some of these observations explicitly in Figure 10. In panel
(a), we show the three-dimensional geometry of the visible magnetic field
lines in the vicinity of the central BH. The field lines are mostly closed
around the BH, and only occasionally thread the event horizon. This confirms
the discussion above, as these field lines are roughly inside
$\hat{\rho}<r_{*}\approx 4GM$ for the choice of parameters in the figure, and
hence, are approximately set by the superradiant driving field (compare also
to Figure 2).

The large-scale differential rotation and twisting of magnetic field lines on
scales of the entire cloud is shown in panel (b) of Figure 10. There, we focus
on a few representative field lines crossing the equatorial plane. Both above
and below the equatorial plane, we find that the azimuthal angular velocity of
the magnetic field lines decreases with increasing distance from the BH’s spin
axis, leading to a lag of portions of the field lines far away from the BH
compared to those closer to the center. This lag causes twisting of the field
lines roughly around the spin-axis666This is only roughly true, since the
plasma motion is more complex as pointed out above. on scales of the
superradiance cloud, which, ultimately, results in the opening of the visible
magnetic field lines at large distances.

In addition to this large scale feature, the field line geometry also exhibits
features on scales of roughly $1/\mu$, as can be seen in panel (c) of Figure
10 (notice, this scale compared with the BH mass is $1/(\mu GM)\approx 3.33$).
There we show two collections of field lines exhibiting variations on spatial
scales set by the superradiance oscillation frequency. These arise likely as a
result of the periodic motion of the plasma in the large scale superradiant
electric field, twisting the field lines on scales of $1/\mu\ll r_{c}$ as
well. Both the twisting on scales of the cloud, and small scales, $1/\mu$,
builds up magnetic energy that is dissipated by means of magnetic
reconnection.

In panel (d) of Figure 10, we isolate one such reconnection sight,
representative for a class of reconnection processes active throughout the
bulk of the cloud. We leave the details to the discussion in Sec. VI.3, and
just point out here that the visible magnetic field lines entering the
reconnection region along the white arrow, diverge away into two distinct
directions along the black arrows. This indicates that the connectivity of the
field lines is discontinuously changed at the location labelled as A. Another
indication of magnetic reconnection at location A in panel (d) of Figure 10,
is the appearance of four bundles of field lines from the region around A. In
the highly conducting and highly magnetized plasma limit, (at least one of)
the dimensions of the reconnection regions are expected to scale as
$1/\sigma$. Speculating about the $\sigma\rightarrow\infty$ limit, we
hypothesize that these reconnection sites turn into one- or two-dimensional
highly fragmented localized filaments and current sheets, where potentially
large amounts of electromagnetic energy is injected into the plasma (as
typically observed in turbulent highly magnetized plasmas, e.g., Zhdankin  _et
al._ (2017); Comisso and Sironi (2018); Pontin  _et al._ (2011)). We discuss
the connection between significant energy dissipation into the plasma and
magnetic reconnection in Sec. VI.3 in detail.

###  V.4 Summary of turbulent plasma scales

Before concluding this section, we provide a summary of the main features of
the quasi-equilibrium endstate of the pair production cascade. We begin with
the largest scales first, and work towards small scales:

  * •

On spatial scales of the superradiance cloud, $\sim r_{c}=1/(\mu\alpha)$, the
superradiant electric field is efficiently screened by a roughly dipolar
charge distribution. A significant fraction of the cloud’s volume is
magnetically dominated at $\sigma/\mu=20$ with trend towards larger fractional
volumes for larger plasma conductivities. This increases the importance of the
magnetic field dynamics for the cloud-plasma system. For sufficiently high
conductivity, the magnetic field and the plasma become strongly coupled except
in isolated diffusion regions. Hence, at large distances from the central BH’s
spin axis, the plasma rotates around the BH with period much longer than the
superradiant cloud’s period, inducing differential rotation on the scale of
the entire cloud. The resulting shearing magnetic field lines reconnect inside
the bulk of the cloud.

  * •

On spatial and temporal scales set by the Compton wavelength of the dark
photon, $1/\mu$, a variety of features appear. The plasma orbits with the
superradiance cloud in circular motion with radius given by roughly $1/\mu$.
This is likely due to the large scale electric field set by the superradiance
cloud. Hence, negative and positive components of the local charge density
orbit exactly out of phase due to the large charge dipole. In the large
conductivity regime, this circular motion implies circular motion of the
visible magnetic field due to the strong coupling in advection dominated
regions. Features of size $1/\mu$ in the global visible magnetic field
geometry appear due to the built-up magnetic field twisting, which is released
in magnetic field line reconnection sites.

  * •

Besides the two macroscopic scales discussed above, the conductivity $\sigma$
sets the size of non-ideal features, which are expected to be of microscopic
size. We found that at moderately large conductivities $\sigma/\mu=20$, these
non-ideal regions begin to form filaments inside the superradiance cloud’s
plasma, setting the scale of the turbulent behavior. Speculating, for very
large conductivities, $\sigma/\mu\gg 20$, the non-ideal regions may fragment
into a large number of current sheet-like structures filling the turbulent
plasma777This is typically found in treatments of turbulent magnetized
plasmas, e.g., Comisso and Sironi (2018); Zhdankin  _et al._ (2017); Pontin
_et al._ (2011).. Below, in Sec. VI.2, we elaborate on this and identify these
regions as sites of enhanced energy dissipation.

Therefore, the superradiance cloud-plasma system is characterized by
differential rotation, as well as periodic plasma motion with period given by
the boson mass scale $1/\mu$, leading the magnetic field lines to be twisted
both on cloud size scales $1/(\mu\alpha)$, and on scales set by the dark
photon mass $1/\mu$. This twisting is relaxed through magnetic reconnection in
the bulk of the superradiance cloud in features with size set by $1/\sigma$.
These processes likely lead to strong electromagnetic transients with
periodicity set by the dark photon mass $\mu$. Characterizing the power and
periodicity of these transients is the subject of the next section.

##  VI Electromagnetic emission

The pair production cascade within the superradiance cloud saturates in a
turbulent, differentially-rotating plasma surrounding the central BH with a
partially screened electric field and magnetic field line reconnection in the
bulk. In highly-magnetized astrophysical plasmas, particles are efficiently
accelerated at reconnection sites, leading to high-energy electromagnetic
emission. Therefore, we expect strong electromagnetic signatures from the
superradiance cloud system.

In the following, we illustrate the radiation and dissipation channels in our
setup in Sec. VI.1, quantify the emitted electromagnetic luminosity in Sec.
VI.2, identify the dominant emission mechanism in Sec. VI.3, discuss the
periodicity of the emission pattern in Sec. VI.4, and comment on the possible
emission spectra in Sec. VI.5.

###  VI.1 Radiation and dissipation channels

The effective description of the pair plasma that we use, introduced in the
previous sections, includes only the leading-order resistive correction to the
force-free limit of ideal magnetohydrodynamics. In the context of this
formalism, any microphysical processes (e.g. pair production, scattering,
photon emission, or other QED effects) are averaged over, or only roughly
approximated by the macroscopic conductivity $\sigma$ that characterizes the
local dissipation in the plasma, and not included from first principles888To
some degree, this could be achieved within the context of kinetic theory and
particle-in-cell simulations. However, we leave exploring this avenue to
future work.. However, our approximation is sufficient to reliably estimate
the total electromagnetic power output of the system through the outgoing
Poynting flux, as well as through dissipation due to macroscopic spatial
currents along the visible electric field. The Poynting flux is typically
invoked, within the force-free paradigm to estimate the rotational energy
extraction rates of pulsars Goldreich and Julian (1969); Contopoulos  _et al._
(1999); Komissarov (2006); Spitkovsky (2006); Ruiz  _et al._ (2014); Petri
(2016); Carrasco  _et al._ (2018) (which have been confirmed within kinetic
theory in Cerutti  _et al._ (2016); Philippov  _et al._ (2015); Philippov and
Spitkovsky (2018)). For typical pulsars, this macroscopic coherent energy flux
is expected to be emitted from the system in the form of lower-energy radio
waves, as well as dissipate in particle acceleration processes in current
sheets resulting in X-rays and gamma-rays. While in the pulsar magnetosphere,
dissipation is mainly confined to roughly two-dimensional current sheets, the
superradiance cloud exhibits reconnection in the bulk, enabling efficient
energy transfer into high-energy emission. The dissipation of electromagnetic
energy can be interpreted as sourcing local particle acceleration, synchrotron
and curvature radiation, and plasma heating. Therefore, estimating the total
emitted Poynting flux and dissipative energy losses of the superradiant system
is crucial to understand the overall electromagnetic signatures. In the
following, we briefly outline how these quantities are computed in our setup.

#### Modified Poynting theorem:

The conservation of energy in the interaction basis, (15), can be used to
identify the macroscopic sources and types of energy flows present in the
system. The background Kerr spacetime has an asymptotically-timelike Killing
field $\xi^{\mu}$, endowing the system with a local energy conservation law.
Therefore, we define the total energy $\mathcal{E}$ of $A_{\mu}$, with respect
to $\xi^{\mu}$, within the (coordinate) domain $D$ as

$\displaystyle\mathcal{E}=\int_{D}d^{3}x\sqrt{\gamma}T^{\alpha}{}_{\mu}n_{\alpha}\xi^{\mu},$
(42)

with the volume form $d^{3}x\sqrt{\gamma}$ of a $t=$const. slice of Kerr
spacetime (see App. C for details), and the energy-momentum tensor
$T_{\mu\nu}$ of the visible fields

$\displaystyle
T_{\mu\nu}=F_{\mu}{}^{\lambda}F_{\lambda\nu}-\frac{1}{4}g_{\mu\nu}F^{\alpha\beta}F_{\alpha\beta}.$
(43)

Throughout this work, the domain $D$ of consideration is the exterior of the
BH up to a coordinate sphere $S_{\hat{\rho}}^{2}$ at coordinate radius
$\hat{\rho}=\sqrt{x^{2}+y^{2}+z^{2}}$ in Kerr-Schild coordinates (defined in
App. C). In the following, we focus entirely on the visible electromagnetic
field. Intuitively, this visible field is a superposition of the massless
(i.e., the SM photon) and the massive vector fields. The former is propagating
freely, and sourced only by the plasma, while the latter is gravitationally
bound to the BH, and non-radiative. Given this, and the energy-momentum
conservation (15), we can relate the change of the total energy of the visible
electromagnetic fields within $D$, $\partial_{t}\mathcal{E}$, to the energy
fluxes across the boundary of the domain, $\partial D$, as well as the work
done on the plasma within $D$, by the modified Poynting theorem, written as

$\displaystyle\partial_{t}\mathcal{E}=-P_{\rm EM}-\dot{\mathcal{E}}_{\rm
BH}-L_{\rm diss}+\dot{\mathcal{E}}_{\rm A^{\prime}}.$ (44)

The first two terms on the right-hand side correspond to the Poynting flux
emitted towards infinity, $P_{\rm EM}$, and the visible electromagnetic field
flux across the event horizon of the BH, $\dot{\mathcal{E}}_{\rm BH}$,
respectively. The third term describes the energy loss of the visible fields
to the pair plasma through resistive processes, $L_{\rm diss}$. Lastly, the
source of energy of the superradiant system is the energy injection of the
massive vector field $\dot{\mathcal{E}}_{A^{\prime}}$. Notice, we assumed that
the energy of the visible electromagnetic fields is much larger than the
energy contained in the pair plasma, i.e., $T_{\mu\nu}\gg T^{\rm
plasma}_{\mu\nu}$. Hence, any finite mass loss due to the accretion or
emission of fermions is not contained in the above analysis, which was shown
to be a good approximation in (28). In the following, we consider each
component on the right hand side of Eq. (44) and take the flat spacetime limit
to connect to familiar expressions.

#### Poynting fluxes:

The electromagnetic luminosity—the Poynting flux—through a coordinate sphere
at radius $\hat{\rho}$, is

$\displaystyle P_{\rm
EM}=-\oint_{S^{2}_{\hat{\rho}}}d\Omega_{\mu}T^{\mu}{}_{\nu}\xi^{\nu}\overset{\rm
flat}{=}\oint_{S^{2}_{\hat{\rho}}}d\Omega\
\boldsymbol{\hat{\rho}}\cdot(\boldsymbol{E}\times\boldsymbol{B}).$ (45)

Here, $d\Omega_{\mu}$ is the oriented area element of $S^{2}_{\hat{\rho}}$
pointing outwards, $\boldsymbol{\hat{\rho}}$ the radial unit vector, and
$d\Omega$ the solid angle. Hence, positive $P_{\rm EM}$ implies visible
electromagnetic energy leaving the domain $D$. Since the massive linear
combination of $A_{\mu}$ and $A^{\prime}_{\mu}$ is gravitationally bound to
the BH and decays as $\sim\exp(-r/r_{c})$, at large distances, the visible
Poynting flux, $P_{\rm EM}$, will receive a contribution only from the
massless linear combination (corresponding to the SM photon) for sufficiently
large $\hat{\rho}$. Analogously, the energy flux across the BH’s event horizon
is

$\displaystyle\dot{\mathcal{E}}_{\rm BH}=-\oint_{S^{2}_{\rm
BH}}d\Omega_{\mu}T^{\mu}{}_{\nu}\xi^{\nu},$ (46)

where $S^{2}_{\rm BH}$ is the event horizon, and $d\Omega_{\mu}$ the oriented
area element pointing outwards. Hence, negative $\dot{\mathcal{E}}_{\rm BH}$
implies visible electromagnetic energy accreting onto the BH. This, of course,
vanishes in the $\alpha\ll 1$ limit (i.e., in the flat spacetime limit). At
the saturation point, i.e., if $\omega=\Omega_{\rm BH}$, the massive field has
vanishing flux across the horizon, such that $\dot{\mathcal{E}}_{\rm BH}$
contains only massless fluxes. We can therefore interpret (46) as a measure of
the amount of accretion, or energy extraction from the BH (e.g. the Blandford-
Znajek process Blandford and Znajek (1977)), triggered by the plasma and
superradiance cloud.

#### Dissipative energy losses:

Besides these fluxes of energy across the boundary of the domain, the
resistive pair plasma is able to dissipate energy in the bulk of $D$. This
macroscopic dissipation is captured by the dissipative losses999Again, due to
the kinetic mixing, both the massive and massless linear combinations of
$A_{\mu}$ and $A^{\prime}_{\mu}$ dissipate energy.

$\displaystyle L_{\rm
diss}=-\int_{D}d^{3}x\sqrt{-g}F^{\alpha\beta}\xi_{\alpha}I_{\beta}\overset{\rm
flat}{=}\int_{D}d^{3}x\boldsymbol{E}\cdot\boldsymbol{J}.$ (47)

Here $g$ is the metric determinant of Kerr spacetime. In the flat spacetime
limit, this expression reduces to the Joule heating within $D$. For later
convenience, we define the local dissipation density

$\displaystyle\rho_{\rm diss}=NF^{\alpha\beta}\xi_{\alpha}I_{\beta},$ (48)

where $N=\sqrt{-g/\gamma}$ is the lapse providing a macroscopic measure of the
local rate at which energy is lost to heating, particle acceleration, etc.

#### Energy transfer from superradiance cloud:

The main energy source, driving the radiation and dissipation, is the
superradiance cloud, which extracted a non-negligible amount of rotational
energy from the BH. The rate of replenishment of $\mathcal{E}$ from the
superradiance cloud $A^{\prime}_{\mu}$ is given by the last term in (44):

$\displaystyle\dot{\mathcal{E}}_{\rm
A^{\prime}}=-\varepsilon\mu^{2}\int_{D}d^{3}x\sqrt{-g}F^{\alpha\beta}\xi_{\alpha}A^{\prime}_{\beta}\overset{\rm
flat}{=}\varepsilon\mu^{2}\int_{D}d^{3}x\boldsymbol{E}\cdot\boldsymbol{A}^{\prime}.$
(49)

This describes the energy transfer from superradiant to the visible fields.

Most important for our discussion in the following sections are the Poynting
flux $P_{\rm EM}$ and the dissipation rate $L_{\rm diss}$. These determine the
total electromagnetic power output of the system, and provide insights into
the characteristics of the emission, such as the primary emission mechanism,
the time-dependence, and the emission spectrum.

###  VI.2 Power output

We find that the cloud-plasma system settles into a driven turbulent state
with a large electric dipole screening the superradiant electric field and
bulk magnetic field reconnection. We demonstrate below that the dipole results
in coherent electromagnetic Poynting flux, while the magnetic reconnection is
associated with strong energy dissipation from the visible electromagnetic
fields into the plasma. As we show in the following, due to the turbulent
nature of the dissipation in the bulk of the cloud, the dissipative losses
dominate over the Poynting flux from the system by orders of magnitude.

In the vacuum limit, i.e., for $\sigma=0$, there is no dissipation and
electromagnetic modes propagate freely. However, both the dark and visible
electromagnetic fields fall off exponentially at large distances away from the
BH, and there is no flux to infinity. At non-zero conductivities of the
medium, on the other hand, any Poynting flux is re-absorbed by the plasma on
scales set by the skin-depth of the effective fluid, which is a complex
function of propagating mode frequencies, conductivity, background
electromagnetic field strengths, and local Ohmic losses. As we show below,
this leads the Poynting flux to go to zero in the intermediate regime
$\sigma/\mu\sim\mathcal{O}(1)$. In this regime, the energy dissipation is
expected to be largest, while any freely propagating electromagnetic modes are
re-absorbed on the skin depth length scales. In the limit where $\sigma/\mu\gg
1$, the regime in which the cloud-plasma system is expected to reside, the
Poynting flux is expected to mostly decouple from the plasma, except in
locations of large dissipation into the plasma, and propagate freely.

Figure 11: (left) We plot the total visible, time-averaged, Poynting flux
$P_{\rm EM}$, defined in (45), through a coordinate sphere at radius
$\hat{\rho}$ around the BH, for various conductivities $\sigma$. Solid lines
are positive (locally outgoing) fluxes, whereas dotted lines are negative
(locally ingoing) fluxes. The interior of the BH and the ergosphere (ES) in
the equatorial plane are indicated by shaded regions; the smallest radius
value indicates the flux through the event horizon $\dot{\mathcal{E}}_{\rm
BH}$. (right) We plot the total energy dissipation rate due to Ohmic losses
$L_{\rm diss}$, defined in (47), everywhere outside a coordinate radius
$\hat{\rho}$ for various conductivities. In both panels, we focus on an
$\alpha=0.3$ cloud with a BH of spin $a_{*}=0.86$, and Bohr radius of the
superradiance cloud of $r_{c}=1/(\mu\alpha)$; notice,
$\sqrt{2(\mu^{2}-\omega^{2})}\rightarrow\alpha\mu$ for $\alpha\ll 1$. Note,
our simulations assume a conductivity constant everywhere in space. At
intermediate conductivities, $\sigma\sim\mu$, the Poynting flux is efficiently
absorbed by the effective plasma, while towards large conductivity, the
electromagnetic modes propagate freely. The energy injection into the plasma
$L_{\rm diss}$ follows the profile of the cloud for all but the highest
conductivities considered here.

To understand the high-conductivity regime, we compute the quantities $P_{\rm
EM}$, $L_{\rm diss}$, and $\dot{\mathcal{E}}_{\rm BH}$ in our numerical
simulations for $\alpha=0.3$ and various different conductivities
$\sigma/\mu\gtrsim 1$ (note, in App. E, we discuss the small-$\sigma$ regime
for completeness). In Figure 11, we present the visible Poynting flux and
dissipative losses as functions of conductivity and coordinate radius for the
superradiant cloud-plasma system. We postpone a discussion on time-dependence
of the electromagnetic emission to Sec. VI.4, and focus here on quantities
time-averaged over one period of the superradiance cloud. As anticipated, the
energy dissipation into the plasma is largest for intermediate conductivities,
$\sigma\sim\mu$. As a result, the visible Poynting flux is efficiently re-
absorbed by the fluid and decays exponentially as it propagates away from the
BH. The sinusoidal features of $P_{\rm EM}$ for $\sigma/\mu\lesssim 1$ are
discussed in detail in App. G. The local energy dissipation follows the radial
profile of the superradiant cloud
$\sim\exp(-\sqrt{2(\mu^{2}-\omega^{2})}\hat{\rho})$, driving the energy
injection into the plasma, at large distances from the BH. Moving towards
larger conductivities, $\sigma/\mu>1$, the weakening of the dissipative losses
roughly follows a $\sim 1/\sigma$ scaling, however, with important corrections
at $\sigma/\mu\gtrsim 5$ discussed below. As can seen in the left panel of
Figure 11, energy flows into the BH for $\sigma/\mu\geq 1$ (though the horizon
is actually a source of energy for lower conductivies, see App. E) at rate
comparable to the total Poynting flux to the wavezone. The radial re-
absorption length scale of the electromagnetic flux increases significantly
with conductivity for $\sigma/\mu>1$, enabling efficient transfer of
propagating modes from the center of the superradiance cloud to the emission
zone far away from the BH. Focusing on the outgoing Poynting fluxes at large
distances, $\hat{\rho}=10r_{c}$, a trend emerges, from small fluxes at
intermediate conductivities to large power at large conductivity, saturating
at a conductivity-independent outgoing electromagnetic emission.

Figure 12:  We show the energy dissipation rate, integrated over the entire
cloud $L_{\rm diss}$, and the total (time-averaged) Poynting flux $P_{\rm EM}$
extracted at $\hat{\rho}=10r_{c}=10/(\mu\alpha)$. (left) Focusing on
$\alpha=0.3$ and $a_{*}=0.86$, we plot the quantities obtained from our
resistive force-free simulations, black circles representing $L_{\rm
diss}^{\rm num.}$, and black squares representing $P_{\rm EM}^{\rm num.}$, as
functions of conductivity. The orange and blue bands, labelled $L_{\rm
diss}^{\rm fits}$ and $P_{\rm EM}^{\rm fits}$, respectively, are a series of
fits of the form $a_{1}+a_{2}(\mu/\sigma)^{p}$ to the simulation results with
the three largest conductivities. The fits are motivated by the discussion in
the main text that the energy dissipation remains finite at infinite
conductivities. The bands are bounded by the most optimistic and pessimistic
fits to the data. The black dotted and dashed lines show the force-free
estimates for the emitted Poynting flux and total dissipation, labelled
$P_{\rm EM}^{\rm ff}$ and $L_{\rm diss}^{\rm ff}$, valid formally at
$\sigma\rightarrow\infty$ (how these are obtained is discussed in the main
text). Lastly, we show the analytical approximations (52), labelled as $L_{\rm
diss}^{\rm ana.}$, for comparison. (right) We also show $L_{\rm diss}^{\rm
num.}$ and $P_{\rm EM}^{\rm num.}$ from the simulations, but now fixing
$\sigma/\mu=20$ and varying $\alpha$, assuming $\omega=\Omega_{\rm BH}$. The
two fits in (54) to the numerical data $L_{\rm diss}^{\rm num.}$ (dash-dotted
line) and $P_{\rm EM}^{\rm num.}$ (sparse-dashed line) are labelled as $P_{\rm
EM}^{\rm fit}$ and $L_{\rm diss}^{\rm fit}$ (and use the ansatz
$a_{1}\alpha^{3}+a_{2}\alpha^{4}$ and $a_{1}\alpha^{1}+a_{2}\alpha^{2}$,
respectively). The orange and blue bands, labelled $L_{\rm diss}^{\rm
extrap.}$ and $P_{\rm EM}^{\rm extrap.}$, are the
$\sigma/\mu\rightarrow\infty$ extrapolations of the corresponding bands in the
plot on the left (there for $\alpha=0.3$) applied to the two fits $P_{\rm
EM}^{\rm fit}$ and $L_{\rm diss}^{\rm fit}$. Lastly, the analytic estimate
(52) is indicated as $L_{\rm diss}^{\rm ana.}$. A discussion of both plots can
be found in the main text.

The conductivity of the pair plasma within the superradiance cloud is expected
to be set by a micro-physical scale far smaller than any macroscopic scale of
the system, $\sigma\gg\mu$. Therefore, in the left panel of Figure 12, we
consider the trends of total dissipation $L_{\rm diss}$ and Poynting flux at
large radii $P_{\rm EM}$ towards the large-$\sigma$ limit. As pointed out
above, the coherent electromagnetic flux emitted from the system increases
rapidly from $\sigma\sim\mu$ towards a non-zero value for $\sigma\gg\mu$. The
blue band in the left panel of Figure 12, indicates possible fits with various
$\sigma$-scalings of the trend, extrapolating to physically relevant regimes,
$\sigma\gg\mu$. We discuss the behavior of $P_{\rm EM}$ for $\sigma/\mu<1$ in
App. G. Turning to the dissipation of energy into the plasma, the behavior in
the low to medium conductivity regime, $\sigma/\mu\lesssim 1$, is as expected.
The energy injection at the macroscopic level increases as $\sim\sigma$ from
the vacuum limit towards intermediate resistivity levels. Beyond the peak
dissipation power around $\sigma\sim\mu$, a decay following $L_{\rm
diss}\propto 1/\sigma$ is predicted by simple arguments of the bulk
dissipation inside the cloud (outlined below). However, instead of following
this behavior, the energy dissipation rate deviates from this trend. In order
to extrapolate to large $\sigma$, we can therefore split the dissipation into
two different components:

$\displaystyle L_{\rm diss}=L_{\rm diss}^{\rm bulk}+L_{\rm diss}^{\rm turb},$
(50)

heuristically representing the bulk and the turbulent dissipation,
respectively. Based on analytic estimates (discussed below) the bulk
dissipation $L_{\rm diss}^{\rm bulk}$ is expected to decrease as $\sim
1/\sigma$ towards small resistivity, while the turbulent dissipation component
$L_{\rm diss}^{\rm turb}$ will have a non-zero value in the infinite
conductivity limit. To capture this, in Figure 12, we fit the results from the
simulations with $L_{\rm diss}\sim a_{1}+a_{2}(\mu/\sigma)^{p}$, considering
$p\in[1,0.025]$, finding $L_{\rm diss}=L_{\rm diss}^{\rm turb}\approx 5\times
10^{-2}(\varepsilon^{2}/G)(M_{c}/M)$ for $\alpha=0.3$ and
$\sigma\rightarrow\infty$. The orange band in Figure 12 represents the range
of values for the conductivity dependence, and is bounded by the most
optimistic and pessimistic fits considered, to illustrate the uncertainty of
this extrapolation.

We compare these extrapolations with the results from force-free simulations,
valid formally at $\sigma\rightarrow\infty$, (see App. D for details on the
numerical implementation and setup; in particular the current is given by
(78)). In the force-free simulations, the continually development of
electrically dominated regions (where the evolution equations breakdown) must
be handled in an ad hoc manner, by reducing the magnitude of the electric
field by hand, which gives rise to an artificial type of dissipation in
regions where current sheets might develop in a more complete description of
the plasma dynamics. Nevertheless, we can determine the effective dissipation
rate by assuming energy conservation (44), as is common in the literature
(e.g. Most and Philippov (2020)), and also compute the Poynting flux given by
the force-free simulations at large distances and across the BH horizon.
Encouragingly, these infinite conductivity results are in good agreement with
the extrapolation of the resistive plasma simulations towards large
conductivity, shown as the orange bands in the left panel of Figure 12. As we
discuss in below, this turbulent dissipation component is associated with
magnetic reconnection and other small scale features of the solution, and
hence, is expected to remain finite even at very large bulk conductivities,
$\sigma\rightarrow\infty$. The Poynting flux at large distances from the
central BH in the force-free setting, shown in Figure 12, is consistent with
the large-$\sigma$ extrapolations of $P_{\rm EM}$ (i.e., is within the blue
band in the left panel of Figure 12).

The bulk dissipation is linked to the large-scale visible electric field
induced by the superradiance cloud, while the turbulent dissipation emerges
from higher-order magnetic field corrections. The former can be understood
analytically in the non-relativistic limit, i.e., for $\alpha\ll 1$, by means
of the solutions derived in (39). Given this solution and neglecting magnetic
field effects, we can determine the electromagnetic current density of the
plasma to be

$\displaystyle\begin{aligned} \textbf{J}&\ =\sigma\textbf{E}\\\ &\
=-\frac{e^{-r/r_{c}}\sqrt{M_{c}\mu}\alpha^{3/2}\varepsilon\sigma\omega^{2}}{2\sqrt{\pi}(\sigma^{2}+\omega^{2})}\begin{pmatrix}\sigma+i\omega\\\
-\omega+i\sigma\\\ 0\end{pmatrix}e^{-i\omega t}+c.c.\end{aligned}$ (51)

In the non-relativistic limit, the spatial extend $r_{c}$ of the superradiance
cloud is large compared with the oscillation frequency $1/\mu$, implying local
charge neutrality $\rho_{q}=0$. Hence, the spatial dependence of $J^{i}$ is a
large scale modulation of the locally neutral plasma for $\alpha\ll 1$. Using
the flat spacetime limit of (47), together with the above current density and
the visible electric field solution (39), the bulk dissipation rate of the
cloud in the $\alpha\ll 1$ limit is

$\displaystyle L_{\text{diss}}^{\rm
bulk}=\frac{\sigma\alpha\varepsilon^{2}}{\mu(1+(\sigma/\mu)^{2})}\frac{M_{c}}{GM}\
.$ (52)

As expected, this quantity goes to zero in both the insulating and highly
conducting limits, leading to free propagation of Poynting flux away from the
system. From Figure 12, we can see that this expression provides a good
approximation for the total dissipation for $\sigma\lesssim 2\mu$ for
$\alpha=0.3$, while for higher conductivities, the turbulent contributions,
i.e., the magnetic field driven component, to $J^{i}$ are more important.

Most relevant for determining the potentially observable electromagnetic
signatures of superradiant systems are the $\sigma\rightarrow\infty$ estimates
for $L_{\rm diss}$ and $P_{\rm EM}$ as functions of $\alpha$. In the right
panel of Figure 12, we show our results for $L_{\rm diss}$ and $P_{\rm EM}$ in
simulations with $\sigma/\mu=20$ as functions of $\alpha$. Focusing on the
numerical results first, it is evident from the right panel of Figure 12 that
the Poynting flux and the total dissipation power have different scalings with
$\alpha$. The analytic estimate (52) for $L_{\rm diss}^{\rm bulk}\sim\alpha
M_{c}/M$ suggests a leading order $\alpha$-scaling of $L_{\rm diss}^{\rm
num.}\sim\alpha M_{c}/M$, which we find in Figure 12 to provide the best fit
to the data. For the Poynting flux, we find a leading-order scaling of $P_{\rm
EM}\sim\alpha^{3}$ to fit best. The two fits shown in the right panel of
Figure 12, are101010Note, fits of the form $L^{\rm fit}_{\rm
diss}\sim\alpha^{p}+\alpha^{2p}$ with $p\in(0.7,1)$ are plausible based on the