numerical data and result in louder signals for $\alpha<0.1$. Fits with
$p\in(1,1.2)$ are equally plausible, however, are entirely consistent with the
large-$\sigma$ extrapolation uncertainty of (53) down to $\alpha\sim 10^{-4}$.

$\displaystyle L_{\rm diss}^{\rm fit}$ $\displaystyle=\
\varepsilon^{2}F(\alpha)\frac{M_{c}}{GM},$ (53) $\displaystyle P_{\rm EM}^{\rm
fit}$ $\displaystyle=\ \varepsilon^{2}G(\alpha)\frac{M_{c}}{GM},$ (54)

with

$\displaystyle\begin{aligned} F(\alpha)&=1.31\times 10^{-1}\alpha-1.88\times
10^{-1}\alpha^{2},\\\ G(\alpha)&=6.86\times 10^{-4}\alpha^{3}-1.36\times
10^{-3}\alpha^{4}.\end{aligned}$ (55)

Determining whether these scalings are also valid in the
$\sigma\rightarrow\infty$ limit requires simulations with larger
conductivities across a larger range of values for $\alpha$. Hence, we
estimate the theoretical uncertainties of the fits (53) and (54), indicated as
orange and blue bands in the right panel of Figure 12, as follows. For
$\alpha=0.3$ (i.e., the left panel of Figure 12), we obtain a series of
different $\sigma\rightarrow\infty$ extrapolations for both the total
dissipation and the Poynting flux (blue and orange bands in the left panel of
Figure 12). The spread of these $\sigma\rightarrow\infty$ extrapolating fits
in the left panel of Figure 12, corresponds to the width of the orange/blue
bands at $\alpha=0.3$ in the right panel of Figure 12. We then use this
relative uncertainty of the large-$\sigma$ extrapolation at $\alpha=0.3$ and
apply it to the fits for total dissipation and Poynting flux, i.e., (53) and
(54), for all $\alpha$, hence, obtaining the orange and blue bands in the
right panel of Figure 12. Therefore, the fit in (53) for $L_{\rm diss}$ is
likely an over-estimate of the $\sigma\rightarrow\infty$ result (the lower
bound of this extrapolation uncertainty is given by $L_{\rm diss}^{\rm
fit}/4$), while the fit in (54) for $P_{\rm EM}$, is likely under-estimating
the true flux at $\sigma\rightarrow\infty$. It is clear from Figure 12, that
the turbulent dissipation power dominates over the total Poynting flux across
the entire parameter space of the $m=1$ superradiant state. Furthermore, the
relatively flat $\alpha$-scaling of the total power output is in stark
contrast to the dependence of the total emitted gravitational wave energy flux
$P_{\rm GW}\propto\alpha^{10}{M_{c}^{2}}/{M^{2}}$ from the oscillating dark
photon cloud in the $\alpha\ll 1$ limit. Lastly, for comparison, the
superradiant instability growth rate of the $m=1$ state scales as $\Gamma\sim
a_{*}\alpha^{6}\mu$.

###  VI.3 Dissipation mechanism

Figure 13: We plot the dissipation density $\rho_{\rm diss}$ (defined in (48))
normalized by the maximal density $\rho_{\rm diss}^{\rm max}=\max\rho_{\rm
diss}$, the ratio of visible electromagnetic fields $B^{2}/E^{2}$, the
violations of the force-free condition $E_{i}B^{i}=0$ normalized by the
magnitude of the visible electric and magnetic fields, and the magntiude of
the plasma velocity $v_{d}=|\boldsymbol{v}_{d}|$ in the equatorial plane of
the central BH. All panels correspond to the same coordinate time and a BH of
spin $a_{*}=0.86$, cloud with $\alpha=0.3$, and plasma conductivity of
$\sigma/\mu=20$. We indicate where $B^{2}/E^{2}=1$ by a contour line. Regions
of small plasma velocities, i.e., large magnetic diffusion, are also sites of
large $E_{i}B^{i}\neq 0$ and locally enhanced energy injection density
$\rho_{\rm diss}$. This implies that magnetic reconnection sites are locations
of enhanced energy injection into the plasma.

Given the importance of the turbulent dissipation power, even at large
conductivity, we next discuss the spatial dependence of the energy dissipation
density $\rho_{\rm diss}$, defined in (48), and demonstrate that regions of
high dissipation are associated with magnetic reconnection. This density
captures the macroscopic energy injection of the electromagnetic fields into
the plasma, driving dissipative processes at the microscopic level. Recall,
magnetic reconnection sites are regions, where the connectivity of the
otherwise frozen-in magnetic field is changed, resulting in large spatial
current along the visible electric fields dissipating energy. In the
following, we focus on the $\alpha=0.3$, $a_{*}=0.86$ and $\sigma/\mu=20$
case, while commenting on how these results extrapolate to the physically
relevant limit of high conductivity.

In the left panel of Figure 13, we show the local dissipation rate per volume
$\rho_{\rm diss}$ [defined in (48)]. On large scales, this quantity follows
the same exponential fall-off in the radial direction as the superradiance
cloud (at $\sigma/\mu=20$), while on smaller scales, strong variation
associated with turbulent features is apparent111111These two spatial
components are naturally associated with the bulk and turbulent dissipation
powers $L_{\rm diss}^{\rm bulk}$ and $L_{\rm diss}^{\rm turb}$, driven by
electric and magnetic fields, respectively (as discussed in the previous
section).. We focus on the latter since, as argued in the previous subsection,
we expect these to persist (though develop smaller scales) in the
$\sigma\rightarrow\infty$ limit. From Figure 13, it is clear that the regions
of local enhancement in the dissipation are associated with magnetic dominance
(second panel) combined with $|E_{i}B^{i}|\sim EB$ (third panel), or with
electric dominance. From Eq. (37), the plasma either allows for a significant
component of the current parallel to the electric field. Focusing on the
insets showing the neighborhood of the BH, locations of the locally enhanced
dissipation density, in addition to having large $|E_{i}B^{i}|$ and being
magnetically dominated, are also bordered by zones of low drift velocity
$v_{d}$ (or equivalently, small magnetic Reynolds number $R_{m}$; fourth
panel), which are associated with efficient magnetic reconnection. Farther
away from the BH, large regions (at $\sigma/\mu=20$) are strongly electrically
dominated, and associated with enhanced dissipation and low drift velocity (in
contrast to the magnetically dominated regions which have $v_{d}\sim 1$),
again indicating strong magnetic diffusion and reconnection.

The magnetic field geometry within the plasma is set by the three different
spatial scales discussed in Sec. V.4. The differential rotation on scales of
the cloud induces a shearing of the magnetic field lines on scales of the
cloud $r_{c}$ around the spin-axis of the BH, while the intermediate scale
oscillations of the plasma (both in the equatorial plane and longitudinally
along the BH spin axis) drive twisting of the field lines on scales of
$1/\mu$. This combined macroscopic build-up of magnetic field twisting is
released in magnetic diffusion regions associated with currents along the
visible electric fields of thickness $1/\sigma$. As discussed in more detail
in Sec. V.3, regions of small magnetic Reynolds number $R_{m}=\sigma
v_{d}/\mu$ are characterized by efficient magnetic diffusion. In these
diffusive regions, the connectivity of the magnetic field lines changes, i.e.
reconnection occurs, which drives enhanced dissipation and accounts for the
dominant channel for the loss of electromagnetic energy in the superradiance
cloud-plasma system.

Figure 14: We show the visible magnetic field lines (in dark/yellow, color
indicating the visible magnetic field strength normalized by the global
maximum $B/B_{\rm max}$) and electric field lines (in white), as well as the
local dissipation density $\rho_{\rm diss}$ (all colors) in two different
contexts. Both panels show a close-up of the plasma roughly $15GM$ away from
the central BH of spin $a_{*}=0.86$, as well as $\alpha=0.3$ and
$\sigma/\mu=20$. In (a), the dissipation density is shown as semi-transparent
isosurfaces. The BH is located towards the bottom of the plot. In (b), the
dissipation density is plotted on a semi-transparent plane spanned by the $z$
and $y$ directions. Here the BH is located towards the top left of the plot.
The numbers in both panels indicate regions of large energy injection into the
plasma. The arrows show the divergence of magnetic field lines away from the
reconnection site. A detailed discussion can be found in the main text.

To illustrate this connection explicitly, we show two example magnetic
reconnection events in Figure 14. In panel (b) of Figure 14, magnetic field
lines enter the reconnection region (labelled “3.”) from the top along the
white arrow (and from behind the semi-transparent plane on which $\rho_{\rm
diss}$ is plotted). These same field lines exit the region in two directions
along the two black arrows (similarly for the lines entering from behind the
semi-transparent plane). The point where the field lines diverge is associated
with the magnetic field magnitude dropping to near zero (indicated by the
color of the magnetic field lines) and locally enhanced dissipation density
$\rho_{\rm diss}$ (indicated by the color in the semi-transparent plane). This
is characteristic of discontinuous reconnection, where the magnetic field
lines change connectivity discontinuously along a line or plane where the
magnetic field goes to zero. In two dimensions, X-point reconnection is the
canoncial example of discontinuous reconnection, and most prominent in current
sheets of the pulsar magnetosphere powering the high-energy component of the
electromagnetic emissions. In particular, the field line geometry shown in
panel (b) of Figure 14 resembles spine-fan type magnetic reconnection Li  _et
al._ (2021); Pontin (2012).

We show a second example of reconnection in panel (a) of Figure 14, where we
isolate two sites of large magnetic field gradients. There, we show a set of
visible magnetic field lines that are strongly twisted as they connect two
regions of locally enhanced $\rho_{\rm diss}$ (labelled as “1.” and “2.”),
separated by a distance of $\sim 1/\mu$. Within each dissipation region, the
field lines are curved on smaller scales (plausibly set by $1/\sigma$). The
strong field gradients in these regions, as well as the fact that the magnetic
field magnitude does not go to zero, make this example more consistent with
continuous reconnection, where magnetic field lines pass through each other in
a diffusion dominated region of small plasma velocity121212Continuous type
reconnection typically occurs at quasi-separatrix layers with large, but
bounded, squashing degree Li  _et al._ (2021); Pontin (2012); Démoulin (2006).
We do not attempt to identify quasi-separatrix layers, instead resort to
identifying reconnection zones based on field diffusivity, magnetic field
curvature and dissipation density. (see also Figure 13).

Therefore, we find that the visible magnetic field line connectivity changes
discontinuously (and we find evidence for continuous reconnection) at various
places in the bulk of the cloud. Both types are accompanied by enhanced energy
injection into the plasma. The reconnection is fundamentally driven by the
fixed orbital frequency of the superradiant magnetic field, suggesting that,
for even larger conductivities, $\sigma/\mu\gg 20$, the qualitative picture is
unchanged. We expect that, in this limit, though the size of the diffusion
regions may decrease as $\sim 1/\sigma$ down to micro-physical scales, the
rate of energy dissipation, which is driven by reconnection, remains roughly
constant. Identifying the changing connectivity of the magnetic field as the
driver of dissipation allows us to make a connection with existing kinetic
analyses to roughly estimate the particle and emission spectrum of the system,
which is the subject of Sec. VI.5.

###  VI.4 Periodicity of emission

Figure 15: Here we consider a system with a BH of spin $a_{*}=0.86$, a
superradiant cloud with $\alpha=0.3$, and a plasma conductivity of
$\sigma/\mu=20$. (left) We show the time dependence of the total visible
Poynting flux entering the BH $\dot{\hat{\mathcal{E}}}_{\rm EH}$ and the
outward flux $\hat{P}_{\rm EM}$ through spheres of coordinate radii of
$6r_{c}$ and $10r_{c}$. Hats indicates the rescaling
$\hat{P}=P(G/\varepsilon^{2})(M/M_{c})$. Time is normalized by the period of
the superradiance cloud $T_{\alpha}=2\pi/\omega$. (right) We show a snapshot
of the visible Poynting flux per solid-angle (centered on the BH) through a
coordinate sphere at $6r_{c}$, normalized by the maximum value. Due to the
differential rotation of the turbulent plasma, this pattern only rotates
slowly along the azimuthal direction, i.e., with period $T\gg T_{\alpha}$. At
a coordinate radii of $10r_{c}$, the periodic modulation of the amplitude of
the Poynting flux is mostly gone, indicating that the dissipation in the
interior region is periodic. The small-scale features in the angular
distribution $dP_{\rm EM}/d\Omega$ is a result of the formation of current
sheets and turbulence in the plasma.

So far, we have discussed the total time-averaged Poynting flux and turbulent
energy injection into the plasma, ignoring the time-dependence of the
emission. In the case of a pulsar, the beamed radio emission, as well as the
pulsed high-energy component of the spectrum are characteristics that emerge
from the oscillation of the magnetic dipole field around the spin-axis of the
star. Since even in the well-studied pulsar case, the radio emission mechanism
is a topic of debate, we ignore it in the following, returning to a brief
discussion of this low-frequency component in Sec. VI.5. The pulsed X-ray
component of the pulsar spectrum requires a large-scale, coherent magnetic
field geometry deep inside the light cylinder, along which charges are
accelerated and radiate, or an oscillating current sheet outside the light
cylinder Bai and Spitkovsky (2010a, b); Cerutti  _et al._ (2016); Philippov
and Spitkovsky (2018); Kalapotharakos  _et al._ (2018). In the context of the
superradiant system considered here, we do not find such coherent and
persistent field or current sheet structures, at least at the conductivities
we consider in this study. Instead, we find the plasma surrounding the BH to
be in a turbulent state without persistent large-scale magnetic or electric
fields. However, since this is driven periodically by the superradiant fields
at a frequency $\omega$ one may naïvely expect the electromagnetic emission to
still be periodic as well. Because modeling the light curve, as is done in the
pulsar case (see e.g., Refs. Bai and Spitkovsky (2010b); Cerutti  _et al._
(2016)), is challenging for the superradiant system, we consider the time-
dependence of the macroscopic Poynting flux at large distances, as well as the
energy dissipation density throughout the bulk of the cloud in order to
understand the temporal evolution of the electromagnetic emission. We find
that the (sub-dominant) Poynting flux shows signs of periodicity, and the
dissipation density locally exhibits weak evidence of periodicity. We close
this section by discussing techniques which could improve our understanding of
the temporal and viewing angle dependence of the electromagnetic emission.

We begin by discussing the time-dependence of the amplitude and angular
distribution of the total Poynting flux. In the left panel of Figure 15, we
demonstrate that both the total Poynting flux entering the BH, and the flux
passing through a coordinate sphere at distance $6r_{c}$ from the central BH,
vary periodically on timescales set by the cloud’s frequency131313Notice, the
visible Poynting flux at finite distances from the BH contains propagating
massive dark photon states that are bound to the BH. This component is
exponentially suppressed at large distances. Hence, we have made sure that the
contribution to $P_{\rm EM}^{6r_{c}}$ in Figure 15 from the massive states is
negligible.. The periodic absorption of electromagnetic energy by the BH is
driven by the longitudinal periodic plasma motion in the vicinity of the event
horizon elaborated on below. The plasma density is expected to roughly follow
the profile of the superradiance cloud (compare (18)), and thus the
conductivity should decrease exponentially away from the central BH. This,
paired with the decreasing contribution of the periodic component of $P_{\rm
EM}$ with increasing distance141414Recall, the conductivity is spatially
constant in our simulations. At $\hat{\rho}=10r_{c}$, the conductivity is
several orders of magnitude smaller compared with (18); hence, the physical
relevance of the Poynting flux at those large distances should be interpreted
with caution. (as shown in Figure 15), suggest that the Poynting flux
periodically injects energy into the plasma $\sim\mathcal{O}(r_{c})$ away from
the BH. Hence, the periodicity of the total emitted Poynting flux is
suggestive of periodic electromagnetic emission.

Figure 16: (top row) We show the dissipation density $\rho_{\rm diss}$ in the
plane spanned by the BH spin and a fixed direction in the equatorial plane at
three instances during a single superradiance cloud period. Here we focus on a
BH with spin $a_{*}=0.86$, $\alpha=0.3$, and plasma conductivity
$\sigma/\mu=20$. The arrows indicate the direction of the motion of the
features, and are discussed in detail in the main text. (bottom) We plot the
maximum of the dissipation density $\max\rho_{\rm diss}$ as a function to time
(normalized by the cloud’s period $T_{\alpha}=2\pi/\omega$) for different
values of $\alpha\in\\{0.1,0.2,0.3,0.4\\}$ (and associated saturated BH spins
satisfying $\omega=\Omega_{\rm BH}$, see Tab. 3). In the bottom panel, we also
indicate the times of the snapshots in the top panels by their corresponding
number labels. Over the course of a single period of the superradiance cloud,
the plasma undergoes periodic motion along the BH’s spin axis (as indicated by
the arrows in 1. and 3. in the top row), leading to peaks in the local
dissipation density, when the plasma from below and above the BH collide in
the equatorial plane (corresponding to snapshot 2.).

We now turn to the temporal variation of the energy injection density
$\rho_{\rm diss}$ into the plasma. Generally, the dissipation density follows
the motion of the plasma on scales of $1/\mu$ and $r_{c}$. Close to the BH,
the most relevant periodic motion of the plasma is the longitudinal motion
(discussed in Sec. V.3) along the spin-axis of the BH. In the first panel of
Figure 16, the dissipation density is peaked in pockets above and below the
BH, moving along the spin-axis towards the equatorial plane, as indicated by
the arrows. This corresponds to the time when the overall maximum of the
dissipation density is at its lowest value per period, as shown in the bottom
plot in Figure 16. Subsequently, two regions of enhanced dissipation density
(and the associated plasma) collide within the equatorial plane, as show in
the middle panel of the top row, leading to locally and temporally large
amplitudes of the energy injection rate $\rho_{\rm diss}$. Finally, the
regions of enhanced dissipation density begin to move away from the equatorial
plane along the spin-axis in the last panel in the top row of Figure 16, and
the associated maximum of the dissipation density decreases. This process
repeats on timescales of the cloud’s period. It is non-trivial to translate
this behavior directly into observable variations of the electromagnetic
signature. We may speculate, however, that this periodic enhancement of the
dissipation density could lead to a periodic flaring of the superradiance
cloud (analogous to e.g., Most and Philippov (2020); Nathanail  _et al._
(2020); Ripperda  _et al._ (2020)). It should be noted though, that the total
integrated turbulent energy dissipation does not show significant temporal
modulations.

In summary, both the Poynting flux and the local dissipation rate exhibit weak
evidence of periodicity on timescales set by the dark photon mass $1/\mu$.
Ultimately, our large-scale macroscopic description of the system is
insufficient to resolve and understand the microscopic particle acceleration
processes active in the turbulent plasma. As well, with our current treatment,
we are unable to determine how the local dissipation rate translates into
observing-angle-dependent radiation, though naïvely one expects this to be
strongly modulated by the oscillation of the superradiance cloud. This could
be improved by consider a small domain near the BH and using higher resolution
to resolve higher values of conductivity than considered in this study, e.g.,
using local resistive force-free techniques. This could determine whether
coherent magnetic and electric field geometries remain inside the light
cylinder, even at very large conductivities. These sufficiently large scale
field geometries may then be used to perform light curve modeling paralleling
the advances made in understanding high-energy pulsar light curves. Another
avenue could be to utilize particle-in-cell (PIC) simulations of the turbulent
regions of the plasma. This would recover the particle acceleration and non-
thermal heating within the plasma, and could therefore be utilized to
understand the time- and angular-dependence of the X-ray/gamma-ray sky map of
the superradiant system.

###  VI.5 Emission spectra

The plasma is characterized by differential rotation on scales of the entire
cloud, superradiant driving on scales of $1/\mu$, and turbulence down to
microscopic scales accounted for in our setup by the inverse conductivity
$1/\sigma$. Visible electromagnetic energy is dissipated into the pair plasma
by resistive processes. Likely this dissipation occurs primarily through
particle acceleration, with the subsequent synchrotron emission of highly
boosted particles leading to high energy photons that escape the system.
Besides this non-thermal component at the high-energy end of the emission
spectrum, various coherent low-energy radio emission processes may be active
in regions of the superradiance cloud. In the following, we briefly review
existing kinetic theory results for the spectra of turbulent pair plasma and
possible low-frequency radio emission mechanisms that may be relevant to the
emission spectrum for the system considered here.

In the pulsar magnetosphere, resistive processes occur mainly in current
sheets outside the light cylinder. There, magnetic dominance is lost, and
electromagnetic energy is efficiently dissipated by accelerating and heating
the plasma. In order to gain insight into the microphysical processes in these
accelerating regions, kinetic approaches based on numerical PIC methods are
typically utilized Nishikawa  _et al._ (2021). Within this framework, the
distributions of charged particles are evolved in time according to the
Lorentz force of the local electromagnetic field, while back-reacting on the
ambient fields through the charge and current they source. Local simulations
resolve microphysical scales such as the Lamour radius
$r_{L}=m_{e}\gamma/(eB)$ of an electron with mass $m_{e}$ and boost factor
$\gamma=E_{e}/m_{e}$ in an ambient magnetic field of strength $B$. Therefore,
these methods are powerful tools to determine the classical particle kinetic
spectrum self-consistently from first principles. On the other hand, radiative
corrections to the particle motion, pair production, are neglected, or added
in an ad-hoc fashion, and the expensive nature of these simulations make it
difficult to apply in a global setting, while still achieving sufficient
resolution to accurately approximating the microphysics. Nonetheless, PIC
methods have successfully recovered the global pulsar magnetosphere, the
expected non-thermal particle spectrum, and have played a central role in
studies of the radio emission mechanism of pulsars Philippov  _et al._ (2015);
Cerutti  _et al._ (2016); Zhang (2020).

In the case of magnetic reconnection, PIC approaches have found that the local
electron acceleration results in a particle distribution $N_{e}(\gamma)$ with
a high-energy, power-law tail below a cutoff $\gamma_{c}$, of the form
$dN_{e}/d\gamma\propto\gamma^{-p}e^{-\gamma/\gamma_{c}}$ Werner  _et al._
(2016); Sironi and Spitkovsky (2014); Guo  _et al._ (2014), for $\gamma\gtrsim
1$ (see also Refs. Kagan  _et al._ (2015); Werner  _et al._ (2018); Cerutti
_et al._ (2015, 2016); Philippov and Spitkovsky (2018)). The size of the
resistive region $\ell$ sets the high-energy cutoff $\gamma_{c}$ as the boost
factor where $\ell$ equals the Lamour radius. Most applicable to the
superradiant system at hand are studies focusing on three-dimensional
turbulent pair plasmas Zhdankin  _et al._ (2017); Comisso and Sironi (2018),
determining the power-law to be roughly $p=2.8$, for large plasma
magnetizations. We leave a detailed investigation of the kinetic spectrum in
the context of a kinetically mixed superradiance cloud to future work, and in
the following make a crude estimate of the high-energy component of the
emission associated with this electron kinetic spectrum based on the
characteristic length scales and field strengths. The high-energy cutoff
$\gamma_{c}=e\ell\langle B^{2}\rangle^{1/2}/m_{e}$ is set by the average
ambient magnetic field strength $\langle B^{2}\rangle^{1/2}$, defined by
$\langle B^{2}\rangle=1/S_{2r_{c}}\int_{S_{2r_{c}}}dVB^{2}$, where
$S_{2r_{c}}$ is a coordinate volume of a sphere of radius $2r_{c}$ centered on
the BH. From the resistive force-free simulations with $\sigma/\mu=20$, we
extract this root-mean-square magnetic energy for each value of $\alpha$ we
consider. Fitting the $\alpha$-dependence by151515Functions of the form
$\sim\alpha^{2}$ or $\sim\alpha^{3}$ provide worse fits, but are plausible
given the numerical and theoretical uncertainty. $\sim\alpha^{5/2}$, the
average magnetic field is

$\displaystyle\langle B^{2}\rangle^{1/2}=2.5\times 10^{8}\
\text{Gauss}\left(\frac{\varepsilon}{10^{-7}}\right)\left(\frac{M_{\odot}}{M}\right)\left(\frac{\alpha}{0.1}\right)^{5/2}.$
(56)

With this, the cutoff electron and positron boost factor, with $\ell=1/\mu$,
is given by

$\displaystyle\gamma_{c}\approx 2.2\times
10^{7}\left(\frac{\varepsilon}{10^{-7}}\right)\left(\frac{\alpha}{0.1}\right)^{3/2}.$
(57)

The size of the resistive region could be set by smaller length scales than
the $1/\mu$ value used above. However, for most the parameter space of
interest, the effects of the radiation reaction will become important for much
lower boost factors than in (57). Synchrotron backreaction becomes
significant, when the radiation reaction timescale $\tau_{R}=E_{e}/P^{\rm
sync}_{e}$, where $P_{e}^{\rm sync}$ is the total single electron synchrotron
power, is comparable to, or smaller than, the Lamour timescale $\tau_{L}=2\pi
r_{L}$. Hence, this radiation reaction becomes important for
$\gamma>\gamma_{r}=(3m_{e}^{2}/(e^{3}\langle B^{2}\rangle^{1/2}))^{1/2}$ with

$\displaystyle\gamma_{r}=3\times
10^{3}\left(\frac{10^{-7}}{\varepsilon}\right)^{1/2}\left(\frac{M}{M_{\odot}}\right)^{1/2}\left(\frac{0.1}{\alpha}\right)^{5/4}.$
(58)

Therefore, in the regime $1\lesssim\gamma<\min(\gamma_{r},\gamma_{c})$, the
power-law electron and positron kinetic spectrum results in a synchrotron
photon spectral power-law $P(\nu)\propto\nu^{-s}$ with spectral index $s=0.9$
Rybicki and Lightman (1986) (assuming $p=2.8$ Comisso and Sironi (2018)),
making up the non-thermal tail of the high-energy component of the emitted
photon spectrum, while above this range the synchrotron spectral index is
modified. The synchrotron spectrum from a single electron or positron in this
non-thermal distribution with $\gamma\leq\gamma_{r}$ peaks at emission
frequencies

$\displaystyle\nu_{\rm peak}=12\
\text{keV}\left(\frac{\gamma}{10^{2}}\right)^{2}\left(\frac{\varepsilon}{10^{-7}}\right)\left(\frac{M_{\odot}}{M}\right)\left(\frac{\alpha}{0.1}\right)^{5/2},$
(59)

where the value of $\gamma$ is chosen inspired by simulations presented in
Cerutti  _et al._ (2016) and a dedicated PIC simulation will be helpful to
determine the exact spectrum, and an electron with kinetic energy of
$m_{e}\gamma_{r}$ radiates mostly at $\nu_{r}=6.4\ \text{MeV}$.

In summary, in the superradiance cloud, the electromagnetic fields lose energy
predominantly through magnetic reconnection in a strong ambient magnetic field
with strength on the order of (56). We can expect that this efficiently
accelerates electrons and positron to large boost factors,
$\gamma\sim\mathcal{O}(10^{3})$, as given by the minimum of the values in (57)
and (58), and that these high-energy particles then radiate synchrotron
photons in the process, with spectrum ranging from a few keV up to MeV [see
(59)]. Therefore, it is likely that there will be strong non-thermal high-
energy component of the emission spectrum from the kinetically mixed
superradiance clouds.

We turning now to the low-frequency, i.e., radio, end of the spectrum, where
the emission mechanisms are far less well-understood. Even in the well-studied
pulsar case, this a topic of debate. Low-frequency electromagnetic phenomena
such as pulsar radio emissions and fast-radio bursts are thought to be sourced
through a shock induced synchrotron maser emission mechanism in the pulsar
wind, reconnection driven radio emission, or near field processes Zhang
(2020). In, for instance, Refs. Mahlmann  _et al._ (2022); Philippov  _et al._
(2019), plasmoids forming from the discontinuous reconnection of the pulsar
current sheet was demonstrated to result in the emission of fast magnetosonic
waves, plausibly escaping as radio emission to infinity. Therefore, the
efficient magnetic reconnection of the superradiant plasma, some of which
occurs through the discontinuous reconnection channel, suggests that the cloud
is a source of continuous radio flux as part of the total power output161616We
note that, near the BH, the plasma frequency is on the order of a GHz, but it
is expected to decrease exponentially away from the BH with the superradiance
cloud density (see Sec. IV.3).. Furthermore, as discussed above in Sec. VI.4,
the plasma performs periodic longitudinal motion along the BH spin axis with
frequency $1/\mu$, resulting in collisions of regions with enhanced
dissipation density within the equatorial plane close to the central BH. If
these collisions, at the microphysical level, manifest as colliding shock
waves and trigger a synchrotron maser mechanism in the process (see, e.g.
Refs. Metzger  _et al._ (2019); Plotnikov and Sironi (2019)), then one would
expect periodically enhanced radio flux from these shocks. Hence, the naïve
expectation is that the plasma-filled superradiance cloud is a source of
continuous radio flux, together with periodic peaks in the radio power with
pulse period set by the dark photon mass $1/\mu$.

##  VII Multimessenger Signals

The system studied thus far motivates a novel target for multimessenger
searches: a new bright, possibly periodic, source with specific, unusual
properties. In this section, we summarize the relevant dynamics and
observational signatures of the kinetically-mixed dark photon superradiance
cloud. The numerical simulations performed in this work give us an estimate of
the total electromagnetic power emitted, but do not directly provide the
spectrum of the emitted radiation, nor conclusively establish its periodicity.

Nevertheless, the unique properties of the system and the analogy with the
well-studied neutron star pulsars allow us to identify promising search
strategies based on our system’s combination of electromagnetic and GW
emission. Given reasonable assumptions, outlined below, we expect current and
planned telescopes and GW observatories to reveal dark photons in the
$10^{-14}-10^{-11}\ \rm{eV}$ mass range, with kinetic mixing below the current
cosmological bound $\varepsilon\lesssim 3\times 10^{-7}$.

The evolution of our new pulsar begins with the birth of a new, rotating, BH.
Around this BH, the superradiance instability populates a cloud of dark
photons in $\mathcal{O}(100)$ superradiance times171717 The results in this
section are obtained using the gravitational waveform model SuperRad Siemonsen
_et al._ (2022) for the superradiant vector cloud, with approximate
expressions for the timescales, etc. given to guide the reader.,

$t_{\rm growth}\sim\ln(M_{c}/\mu)\tau_{\rm SR}\approx 10^{4}\,{\rm
s}\left(\frac{M}{10\
M_{\odot}}\right)\left(\frac{0.7}{a_{*}}\right)\left(\frac{0.1}{\alpha}\right)^{7}.$
(60)

The resulting large electromagnetic fields will, for large enough values of
the kinetic mixing parameter, $\varepsilon\gtrsim 10^{-10}$, build up a dense
plasma of charged particles in the last few $e$-folds before saturation (see
Sec. IV). The rotation of the cloud and resulting turbulent electromagnetic
processes in the plasma lead to a large flux of electromagnetic emission from
the system, as described in Secs. V and VI.2. The energy output is dominated
by the dissipative losses in the turbulent regions, which we infer from our
numerical simulations (see Eq. (53)) to be

$\displaystyle L_{\rm{EM}}$
$\displaystyle=\varepsilon^{2}F(\alpha)\frac{M_{c}}{GM}\simeq\varepsilon^{2}\frac{\alpha^{2}\Delta
a_{*}}{G}$ $\displaystyle\simeq 4\times
10^{41}{\rm{erg}/\rm{s}}\left(\frac{\varepsilon}{10^{-7}}\right)^{2}\left(\frac{\alpha}{0.1}\right)^{2}\left(\frac{\Delta
a_{*}}{0.1}\right),$ (61)

where $F(\alpha)=0.13\alpha-0.19\alpha^{2}$ is a polynomial fit to the
simulations, and we used the mass of the superradiance cloud at its maximum
and the small $\alpha$ limit. The luminosity can be up to five orders of
magnitude brighter than the Crab pulsar’s bolometric luminosity Hester (2008),
and up to ten orders of magnitude brighter than the solar luminosity. In the
following, we will assume the luminosity to be given by Eq. (61) also for
$\alpha<0.1$, below the smallest simulated value.

The cloud slowly decays through emission of GWs and electromagnetic radiation
on a timescale generally dominated by the gravitational dissipation of the
cloud, except at small dark photon masses and large mixing, as shown in the
left panel of Fig. 17. The observable electromagnetic signal lasts for
$\sim\rm{min}\left\\{\tau_{\rm GW},\tau_{\rm EM}\right\\}$,

$\displaystyle\tau_{\rm GW}$ $\displaystyle\approx$
$\displaystyle\frac{GM}{17\alpha^{11}\Delta a_{*}}\approx 10^{6}\,{\rm
s}\left(\frac{M}{10M_{\odot}}\right)\left(\frac{0.1}{\alpha}\right)^{11}\left(\frac{0.1}{\Delta
a_{*}}\right),$ $\displaystyle\tau_{\rm EM}$ $\displaystyle\approx$
$\displaystyle\frac{GM\ln 2}{\varepsilon^{2}F(\alpha)}\approx
10^{11}{\rm{s}}\left(\frac{M}{10M_{\odot}}\right)\left(\frac{10^{-7}}{\varepsilon}\right)^{2}\left(\frac{10^{-2}}{{\it{F}}(0.1)}\right).$

The decay is power-law in time when more energy is released in GWs than
electromagnetic radiation and exponential otherwise,

$\frac{M_{c}(t)}{M_{c}(t_{0})}=\begin{cases}[1+(t-t_{0})/\tau_{\rm{GW}}]^{-1}&\
\ \tau_{\rm GW}\ll\tau_{\rm EM}\\\ e^{-(t-t_{0})/\tau_{\rm EM}\ln 2}&\ \
\tau_{\rm GW}\gg\tau_{\rm EM}.\end{cases}$ (63)

At small $\alpha$ and large enough $\varepsilon$, the superradiance growth
time $\tau_{\rm SR}$ can become slower than $\tau_{\rm EM}$, possibly
preventing the cloud from reaching its full size. A detailed study of the
cloud saturation in this case is beyond the scope of this work and we always
require $\tau_{\rm{SR}}<\tau_{\rm EM}$, which applies to most of the open
parameter space (see Fig. 17).

To summarize, the evolution of a kinetically mixed cloud is fully fixed by the
dark photon mass $\mu$, which sets the overall fundamental scale, and the
dimensionless couplings $\varepsilon$ and $\alpha$. Over the lifetime of the
cloud, the electromagnetic and GW signals grow exponentially—at a rate fixed
uniquely by the dark photon mass and $\alpha$—until the cloud reaches its
maximum size, and then decrease on longer timescales as the cloud
disappears.181818The electromagnetic signal is only present for large enough
cloud sizes $M_{c}$ and kinetic mixings $\varepsilon$ which are needed to
generate the plasma, Sec. IV. If the GW radiation dominates, there is a unique
relation between the growth and decay timescales, $\tau_{\mathrm{SR}}$ and
$\tau_{\mathrm{GW}}$. For $\alpha\ll 1$ (for which the GW timescale is given
by Eq. (11)) it takes the approximate form

$\tau_{\mathrm{SR}}\simeq
0.98\,\tau_{\mathrm{GW}}^{7/11}r_{g}^{4/11}\bigg{[}\,a^{*}\bigg{(}\frac{\Delta
a^{*}}{0.5}\bigg{)}^{7/11}\bigg{]}.$ (64)

We show the characteristic electromagnetic luminosities and their time
evolution due to GW emission in the right panel of Fig. 17.

Figure 17: (left) Lifetime of the superradiance cloud as a function of the
kinetic mixing parameter $\varepsilon$ and the gravitational coupling $\alpha$
for a BH with an initial mass of $10\ M_{\odot}$ and spin of $a_{*}=0.9$. At
large $\alpha$, the cloud decays through GW emission, and the lifetime is
independent of $\varepsilon$. When $\alpha$ is too small, the power emitted in
electromagnetic radiation overcomes the GW power and the cloud depletes faster
for larger $\varepsilon$ (see Eq. LABEL:eq:signaldecay). In both regimes, the
lifetime is proportional to the BH mass, $\tau\propto M$, so the transition is
independent of the value chosen. The initial BH spin determines the largest
value of $\alpha$ that satisfies the superradiance condition, but otherwise
has a mildly effect on the lifetime of the cloud. In the dark gray shaded
region $\tau_{\rm{SR}}>\tau_{\rm{EM}}$, while the light shading corresponds to
parameters excluded by measurements of the CMB spectrum by COBE/FIRAS Fixsen
_et al._ (1996); Caputo  _et al._ (2020). (right) Time evolution of the
superradiance cloud’s electromagnetic luminosity (see Eq. 61) for two
different values of $\varepsilon$, for $M=10\ M_{\odot}$ and $a_{*}=0.9$ (the
luminosity is independent of $M$ and only mildly dependent on $a_{*}$, while
the decay time will increase for heavier BHs). After the spinning BH is
formed, the energy emitted in radiation quickly grows exponentially with the
superradiance cloud, and later slowly decreases due to the cloud mass decay
through GW emission.

Clearly, the observational prospects of this system depend heavily on the
spectral shape of the electromagnetic radiation; unfortunately, our
simulations do not give us this information. However, given the similarities
of our system to pulsars, which have been observed across electromagnetic
bands in many systems, and have PIC simulations in agreement with aspects of
the observations, we can make educated guesses as to the expected emission in
different bands.

In analogy with neutron star pulsars, one dominant emission mechanism could be
synchrotron radiation. As discussed in Sec. VI.5, the boost factor of
electrons and positrons in the cloud peaks in the range
$\gamma\sim\mathcal{O}(10-100)$, giving typical electron radiation frequencies
ranging from keV to a few MeV. Depending on the value of the kinetic mixing
and BH mass, most of the spectrum would fall within the range of X-ray
telescopes, such as Chandra cha  and Swift Burrows  _et al._ (2005) of
$\mathcal{O}(0.1-10\,{\rm keV})$ or Fermi-GBM ($8\,{\rm keV}$ to $1\,{\rm
MeV}$) Atwood  _et al._ (2009).

GW emission from the system is monochromatic with frequency
$f_{\rm{GW}}=\omega/\pi$, and the electromagnetic emission is also expected to
have periodicity on timescales of $1/\omega$ (Sec VI.4). Here, $\omega$ is the
energy per dark photon, given by its rest mass with
$\Delta\omega/\mu\approx-\alpha^{2}/2-\mathcal{O}(\alpha^{4})$ corrections due
to the gravitational potential energy of the BH Baumann  _et al._ (2019) and
$\Delta\omega_{c}/\mu\approx-(5/8)\alpha^{2}M_{c}/M$ the gravitational self-
energy of the cloud Siemonsen  _et al._ (2022). The combination of the
decrease of the BH mass as the cloud grows, and the decrease of the cloud mass
as it decays, result in a monotonically increasing frequency correction. In
other words, the cloud period decreases during the whole evolution of the
cloud Siemonsen  _et al._ (2022), in stark contrast with conventional pulsars,
for which the period increases in time. Thus, rotating BHs can host an
anomalously bright “pulsar” which spins up over time Arvanitaki  _et al._
(2015).

This new type of pulsar has several unique features and a peculiar evolution
history. While there are a variety of signatures that can be looked for, here
we highlight two distinct observational prospects: searching for the emergence
of a bright electromagnetic source from a known rotating BH remnant (with or
without a periodic component), or searching for continuous GWs emitted by
anomalous pulsars. These two observational strategies are best-suited to
regimes of large and small $\alpha$, respectively.

Given the signal uncertainties, we propose discovery oriented searches rather
than exclusions. It is possible that the absence of a large number of
ultraluminous X-ray sources already places limits on dark photon parameter
space; however, given uncertainties in natal BH spin distributions, as well as
in the emission spectrum of the dark photon cloud, such a constraint would not
be robust. Similarly, we can speculate that such kinetically mixed
superradiance clouds may account for some of the ultraluminous X-ray sources
already observed Swartz  _et al._ (2004); Pintore  _et al._ (2018).

###  VII.1 Electromagnetic follow-ups of black hole mergers

Figure 18: Range of dark photon kinetic mixing parameter $\varepsilon$ and mass $\mu$ producing a visible signal for electromagnetic follow-up observations of LVK compact binary merger events with a BH remnant. As an example, we choose the three best target events: GW170814 (blue), GW190814 (yellow), GW200202_154313 (red) (measured parameters are given in Table 1). The thin solid lines show the regions where a signal could exist, which are bounded from below by the requirement that the visible electric field is large enough to produce the plasma (see Sec. IV and Eq. (IV.2)), to the left by the requirement that the cloud grows within 10 years (see Eq. (60)), and to the right by the superradiant condition for the fastest-growing bound state (see Sec. II). The reach is further limited on the right by the signal duration falling below a minimum observational time. The gray shaded region is excluded by measurements of the CMB spectrum by COBE/FIRAS Fixsen  _et al._ (1996); Caputo  _et al._ (2020). (left) Current prospects for an X-ray Burrows  _et al._ (2005) (solid) and radio transient Dobie  _et al._ (2022) (dashed) search. (right) Future prospects for an X-ray search (solid) and a radio search for a pulsating source (dashed). See the text for more details. The cyan contour corresponds to one merger event per year visible in the X-ray, with shaded band indicating the error due to the uncertainty in the BH merger rate. Name | Final BH mass $[M_{\odot}]$ | Final BH spin | Distance [Mpc] | $\Delta\Omega\ [\rm{deg}^{2}]$
---|---|---|---|---
GW170814 Abbott  _et al._ (2019) | 53.2 | 0.72 | 600 | 87
GW190814 Abbott  _et al._ (2021a) | 25.7 | 0.28 | 230 | 19
GW200202_154313 Abbott  _et al._ (2021b) | 16.76 | 0.69 | 410 | 170
Table 1: List of example compact binary merger events observed by LIGO-Virgo-
KAGRA that are promising candidates for a dark photon superradiance search
through electromagnetic follow up observations and corresponding central
values of their parameters. See Fig. 18 for the observational prospects.

Most of the compact binary mergers detected by the LVK observatories result in
a BH remnant with a mass between 10 to a 100 $M_{\odot}$ and high spin, due to
the capture of a significant component of the binary’s orbital angular
momentum Buonanno  _et al._ (2008). Electromagnetic follow-up observations of
the mergers could reveal the existence of a dark photon superradiance cloud
around the BH remnant. We select three illustrative events from the LSC O1-3
catalogs Abbott  _et al._ (2019, 2021a, 2021b) and forecast the sensitivity of
radio and X-ray searches in Fig. 18. The best targets and their inferred
parameters are listed in Table 1; they correspond to the closest mergers with
the best sky localization and one high mass and one low mass event, to cover
the widest range of dark photon masses. We also include the possible neutron
star–BH merger GW190814 Abbott  _et al._ (2020), for which an electromagnetic
counterpart was not found in the radio Dobie  _et al._ (2019, 2022); Gourdji
_et al._ (2021), optical, and near-infrared Kilpatrick  _et al._ (2021); de
Wet  _et al._ (2021); Thakur  _et al._ (2020); Ackley  _et al._ (2020); Vieira
_et al._ (2020), and X-ray Page  _et al._ (2020). These null results could
exclude new dark photon parameter space, provided a more reliable prediction
for the emission spectra from the superradiance system.

For our projections, we assume an $\mathcal{O}(1)$ of the superradiance cloud
luminosity (61) is emitted into X-rays and a $\mathcal{O}(10^{-4})$ fraction
into radio frequencies, in analogy with the spectrum of standard pulsars Kaspi
and Beloborodov (2017); Manchester  _et al._ (2005). For the follow-up to
detect the EM emission, the cloud needs to grow within a reasonable
observational timescale; for concreteness we impose the requirement
$t_{\mathrm{growth}}<10$ years (c.f. Eq. (60)), which translates into a
large-$\alpha$ requirement given by $\alpha\gtrsim 0.036,0.04$, and $0.03$ for
our three selected candidates. The left panel of Fig. 18 shows how current
instruments could already measure a signal, taking a radio flux sensitivity of
$40\ \mu\rm{Jy}$ at $944\ \rm{MHz}$ for an observation time of 10 h,
achievable by ASKAP searching for transient events Dobie  _et al._ (2022), and
an X-ray detection sensitivity of $2\times 10^{-14}\ \rm{erg}\
\rm{cm}^{-2}s^{-1}$ in $10^{4}\ \rm{s}$, attainable by Swift-XRT Burrows  _et
al._ (2005) as it scans through the sky (eROSITA has similar performance
Predehl (2017)). Better reach can be achieved with Chandra, which however
requires prior angular localization. For instance, for sources with angular
localization from radio observations, Chandra could probe almost an order of
magnitude smaller X-ray fluxes Jaodand  _et al._ (2019).

The sky localization is expected to improve during O4 with the Advanced LIGO,
Virgo, and KAGRA network, with about $10\%$ of the events localized within $5\
\rm{deg}^{2}$ Abbott  _et al._ (2018), improving the prospects of detecting an
electromagnetic counterpart. In the right panel of Fig. 18, we assume smaller
positional errors that would allow a Chandra-like X-ray search sensitive to
$10^{-15}\ \rm{erg}\ \rm{cm}^{-2}s^{-1}$ fluxes in $10^{5}\ \rm{s}$. In the
same panel we also indicate with a cyan contour the region of parameter space
above which more than one X-ray event per year from BH mergers could be
observed with the same X-ray sensitivity. To obtain this contour we made use
of the BH merger rate as a function of primary BH mass measured by LVK Abbott
_et al._ (2021c), assuming a final BH spin of 0.7, and a final mass equal to
twice the primary mass. The shaded band around the cyan contour indicates the
uncertainty in the merger rate.

The searches discussed so far in this subsection are aimed at a steady source
that shines for as long as the required observational time, without assuming
any periodicity. If we further assume partial or total periodicity in the
electromagnetic emission power, a search for a long lasting pulsating radio
signal could be performed. We estimate that the prospect of such a search with
a sensitivity of $10\ \mu\rm{Jy}$ at $1\ \rm{GHz}$ with $500\ \rm{MHz}$
bandwidth for 15 minute observations daily over the lifetime of the
cloud,191919Private communications with Kendrick Smith. which is comparable to
the performances of FAST Nan  _et al._ (2011) and slightly better than CHIME
CHIME/Pulsar Collaboration (2021).

The superradiance cloud could also have higher energy emission, up to
$\gamma$-rays. Telescopes with nearly all sky coverage, such as Fermi-LAT
Atwood  _et al._ (2009), are well suited to perform follow-up observations of
compact binary coalescences (which have been done and are planned during O4
Ren and Wang (2022)). Current flux sensitivities result in reach in kinetic
mixing comparable to current constraints, but a signal could be detected in
the event of an exceptionally close merger.

In the event of a positive detection of a new luminous source following a
binary BH merger, there are non-trivial cross checks that can be used to
confirm the superradiance origin of the signal. First, by measuring the peak
luminosity of the different sources, the luminosity’s unique dependence on the
parameter $\alpha$—Eq. (61)—can be verified. For a dark photon with a given
mass and kinetic mixing, the luminosity will only depend on the BH mass and
spin, which are measured in the merger. Secondly, the measured light curve
should be consistent with an exponential growth as dictated by superradiance,
and power-law decay, as expected from GW emission, with growth and decay
timescales satisfying the non-trivial relation Eq. (64). Thirdly, if a
pulsating electromagnetic signal is observed, the measurement of the period,
given by $2\pi/\mu$, can be used to extract the dark photon mass and verify
the signal growth and decay times dependence on the parameter $\alpha$—see
Eqs. (60) and (LABEL:eq:signaldecay). Multiple sources will have the same
period up to Doppler shifts of $\mathcal{O}(10^{-3})$, and gravitational
potential corrections of up to a few percent. Finally, the superradiance
system also emits continuous GWs with frequency $\omega/\pi$, which could be
detected with a search targeted on the luminous source Arvanitaki  _et al._
(2017); Isi  _et al._ (2019); Chan and Hannuksela (2022); Siemonsen  _et al._
(2022). The GW has the same signal growth and decay times as the
electromagnetic emission, which, together with the period coincidence,
constitute unmistakable signatures of dark photon superradiance.

The discovery of such an ultraluminous source depends crucially on our ability
to localize the newly formed BH. In particular, for sources with short
duration (large $\alpha$), it is important that the source is located hours
before the merger to allow telescope observation coincident with the merger.
Space-based, mid-band detectors can locate stellar mass BH merger events
similar to GW150914 to an angular area of $\sim 0.01\ \rm{deg}^{2}$ an hour
before the merger Graham and Jung (2018); Seymour  _et al._ (2022). Similar
angular localization can be reached with the Laser Interferometer Space
Antenna (LISA) and similar space missions for intermediate mass and
supermassive BH mergers Cutler (1998); Ruan  _et al._ (2019). These missions
could significantly improve the chances of finding the X-ray and radio signals
from the dark photon superradiance cloud.

###  VII.2 Gravitational follow-ups of anomalous pulsars

If the superradiance electromagnetic emission is periodic, old galactic BHs
dressed with a dark photon cloud could look like a neutron star pulsar, and be
detected by ongoing surveys of pulsating sources. Such “fake” pulsars appear
to be rotating at a single frequency $f$ set by the dark photon mass, spin up
over time, and emit continuous GWs with frequency $f_{\rm{GW}}=2f$. Their GW
strain is up to several orders of magnitude larger than the neutron star spin
down limit and can be searched for with targeted continuous GW analyses.

We select two types of potential candidates in the ATNF pulsar catalogue
Manchester  _et al._ (2005); pul . The first type is a set of frequency
multiplets, i.e. sources $i$ with at least one other pulsar $j$ that satisfies
$|f_{i}-f_{j}|/f_{i}<10^{-3}$, which takes into account the spread in
frequency due to Doppler shifts of $\mathcal{O}(10^{-3})$. The second type is
a source with $\dot{f}_{\rm{obs}}>0$ that is not known to be in a binary
system, to avoid spurious positive spin frequency derivatives due to the
source’s acceleration in a binary orbit.

Since all these pulsars are within our Galaxy, the BH formation event must
have occurred long enough ago so that it was not observed. We therefore impose
$\tau_{\rm{GW}}>10^{3}\ \rm{yr}$, which translates into a maximum value for
$\alpha\lesssim 0.03-0.05$ for the selected sources. For each system, the dark
photon mass is fixed to be $\mu=2\pi f$ (gravitational corrections to the dark
photon energy are negligible in the small $\alpha$ limit). A lower bound on
$\alpha$ is set by the smallest possible BH mass, which we take to be 3
$\rm{M}_{\odot}$, resulting in $\alpha\gtrsim 0.05-0.005$ across the range of
frequencies considered. We take the initial BH spin to be 0.5. A better choice
for the BH spin would be to sample it from the spin distribution of the
Galactic BH population; however, this is not well known.

We only consider sources with $f>50\ \rm{Hz}$, since for a stellar mass BH,
smaller frequencies would correspond to a very small $\alpha$, and
correspondingly small signal. Above this threshold, the ATNF catalogue
contains 229 sources with measured $f$, distance, and with positive or unknown
$\dot{f}_{\rm{obs}}$. Among these, there are 3 frequency triplets and 26
frequency doublets (19 of the doublets have a non-zero allowed range of
$\alpha$ between the upper and lower bound), statistically compatible with the
number of accidental multiplets expected for a uniform frequency distribution
between $50$ and $500$ Hz. For each fake pulsar in a multiplet, the emitted GW
strain cannot be uniquely predicted, since the BH mass and age (that sets the
remaining superradiance cloud mass, and thus affects the GW emission power)
are unknown. In the left panel of Fig. 19, we show the maximum possible strain
for the allowed range of $\alpha$ values, for two example BH ages of $10^{3}$
and $10^{6}$ years. As pointed out previously, we do not consider younger
systems. These young systems are excluded by all-sky searches for continuous
GWs Abbott  _et al._ (2022a, 2021d), but strains below the current bounds are
possible for older systems, and within the reach of a targeted search Abbott
_et al._ (2022c).

With the same selection criteria as above, but without requiring frequency
multiplets, we find 20 sources with $\dot{f}_{\rm{obs}}>0$, with values
between $4\times 10^{-17}$ Hz/s and $5\times 10^{-14}$ Hz/s. Only 16 of these
sources have a non-zero range of $\alpha$ values between the upper and lower
bound. Interestingly, four of these sources (J0024-7204Z Freire  _et al._
(2017), B0021-72G Robinson  _et al._ (1995); Freire  _et al._ (2017),
J1801-0857C Lynch  _et al._ (2011); Pan  _et al._ (2021), and B0021-72M
Robinson  _et al._ (1995); Freire  _et al._ (2017)) also belong to a frequency
doublet. The frequencies of the four doublets are approximately 219.6, 247.5,
267.4, and 272.0 Hz. These candidates could be further strengthened (or
disfavored) by performing spin derivative measurements of the other component
of each of the doublets, which at present do not have $\dot{f}$ measurements,
and that are J0514-4002D Ridolfi  _et al._ (2022), J1824-2452J Freire (2008),
J0024-7204ad Ridolfi  _et al._ (2021), and J0125-2327 Morello  _et al._
(2019), respectively.

If the spin-up is due to the superradiance cloud decay through GW emission, we
expect an intrinsic spin frequency derivative Siemonsen  _et al._ (2022)

$\dot{f}_{\rm{int}}\simeq\frac{5}{8\pi}\alpha\mu^{2}GP_{\rm{GW}}.$ (65)

The observed frequency derivative could differ from the intrinsic one due to
additional positive contributions from acceleration along the line of sight.
Thus, one should interpret the measured spin derivative $\dot{f}_{\rm{obs}}$
only as an upper limit on the intrinsic spin derivative. Here, however, our
intention is to provide a first example of how spin derivative measurements
could be used to discover anomalous pulsars, so in what follows we simply
neglect accelerations along the line of sight and assume
$\dot{f}_{\rm{obs}}\simeq\dot{f}_{\rm{int}}$. With this assumption, the
measured value of $\dot{f}_{\rm{obs}}$ fixes the power emitted in GWs
$P_{\rm{GW}}$ for any given $\alpha$. Since the power emitted decreases as the
cloud’s mass decays, there is a minimum $\alpha$ that allows for a spin up
rate as large as the one observed when the cloud is as young as possible,
i.e., at least 1000 yr. This lower bound on $\alpha$ is stronger than the one
described above from the minimum BH mass for $f\lesssim 270\ \rm{Hz}$ .
Therefore, the observable strain is predicted to be within a small range given
by the small spread of $\alpha$ values allowed, as shown on the left of Fig.
19 with the short segments bounded by the downward and upward pointing blue
triangles.

We find that most sources should have already been seen by all-sky searches
for continuous waves, but one candidate pulsar with a frequency around 600 Hz
remains unconstrained, which could be an interesting candidate for a targeted
search. We note that if the assumption
$\dot{f}_{\rm{obs}}\simeq\dot{f}_{\rm{int}}$ is correct, the four previously
mentioned events that are frequency multiplets and have
$\dot{f}_{\mathrm{obs}}>0$ would be excluded, since their
$\dot{f}_{\mathrm{obs}}>0$ measurement leads to a strain prediction that is
already ruled out by existing GW searches. These sources, however, could
certainly be compatible with data if $\dot{f}_{\rm{int}}<\dot{f}_{\rm{obs}}$,
due to the source accelerations along the line of sight that we have
neglected.

For each anomalous pulsar, the kinetic mixing parameter could lie within a
range of allowed values, which are shown in the right panel of Fig. 19. In the
prediction of the strain described above, we assumed that the cloud decays
through GW emission, and that the power emitted in electromagnetic radiation
is subdominant at all times, giving an upper bound on $\varepsilon$ of around
a few times $10^{-7}$ for a system that is 1000 years old (see Fig. 17)—a
different time evolution for the cloud is, in principle, allowed and would
give a different observable strain. The lower bound on $\varepsilon$ shown in
Fig. 19 comes from requiring the superradiance cloud to pair produce the
plasma, as described in Sec. IV. This bound, however, can be relaxed since the
BH was formed long ago, and charged particles could be slowly accreted and
build up the plasma over time (see Sec. IV.3). In this case, $\varepsilon$ as
small as $10^{-12}$ could produce a pulsating source that is luminous enough
to be observed at galactic distances. Finally, we notice that a few of the
pulsating sources used here have a measured luminosity, which could be used to
fix the value of $\varepsilon$ if the emission spectra were known (see Sec.
VI.5). Taking the fraction of total luminosity that goes in the radio band to
be between $10^{-4}$ and $10^{-6}$ gives $\varepsilon$ between approximately
$10^{-12}$ and $10^{-11}$ Manchester  _et al._ (2005); Kaspi and Beloborodov
(2017).

Figure 19: (left) We show the projected gravitational wave strain of observed
pulsating sources whose luminosity could be powered by the kinetically mixed
dark photon superradiance cloud. The potential candidates are selected from
the ATNF pulsar catalogue Manchester  _et al._ (2005), as described in Sec.
VII.2, and are frequency doublets (dark red), frequency triplets (light red),
and pulsars with positive measured frequency derivative (blue). The filled
(empty) triangles correspond to the largest possible strain that a source in a
frequency multiplet could produce if the cloud were created $10^{3}$
($10^{6}$) yr ago. The strain could take any value below that upper bound
(thin solid lines) down to a minimum strain outside of the range shown here.
The blue down-pointing and up-pointing triangles denote the range of strains
allowed assuming that the spin-up is due to gravitational wave emissions from
the cloud. The $95\%$ C.L. upper limits on the signal strain amplitude from
Refs. Abbott  _et al._ (2022a) and Abbott  _et al._ (2021d) are shown in dark
and light gray shading, respectively. The solid gray line corresponds to the
expected sensitivity of a targeted search with LIGO Livingston derived in
Abbott  _et al._ (2022c). (right) Range of kinetic mixing parameters allowed
for each pulsar, with frequency doublets (dark red), frequency triplets (light
red), and pulsars with positive measured frequency derivative (blue). For the
frequency multiplets, filled (empty) down-pointing triangles correspond to the
largest $\varepsilon$ that allows the cloud to decay through GW emission for
at least $10^{3}$ ($10^{6}$) yr. For the sources that are spinning up, down-
pointing triangles give the largest $\varepsilon$ that allows the cloud to
decay through GW emission up to their current age (which is fixed by
$\dot{f}_{\rm{obs}}$). For all the sources, up-pointing triangles denote the
smallest mixing parameter that allows for plasma pair production in the cloud
(see Sec. IV and Eq. (IV.2)).

###  VII.3 Concluding remarks

In the last two subsections, we discussed two detection strategies which take
advantage of the multimessenger signals from a dark photon superradiance cloud
around a spinning BH. In both cases, a combination of the electromagnetic and
GW observations allow us to perform highly non-trivial cross checks to
uniquely identify our system and measure the dark photon mass and kinetic
mixing parameter. Our analysis can be strengthened by dedicated numerical
studies which

  * •

extend our analysis to smaller gravitational coupling $\alpha$ and larger
conductivity $\sigma/\mu$,

  * •

provide robust information about the electromagnetic emission spectrum,

  * •

show robust evidence for or against periodicity in the electromagnetic
emission in different electromagnetic bands.

In addition to the signals mentioned above, our system may host a plethora of
phenomena, including ultraluminous X-ray sources (ULX) Swartz  _et al._
(2004), as well as transient processes such as fast radio bursts (FRB) Amiri
_et al._ (2019); Petroff  _et al._ (2019) and X-ray flares and giant flares
Cline  _et al._ (1980). Follow-up GW observations of these events could
identify their origin as a dark photon superradiance cloud. With
electromagnetic observations alone, more dedicated numerical studies are
needed to determine the similarities and differences in the spectral
properties and transient dynamics between our system and a neutron star pulsar
or magnetar Philippov  _et al._ (2015), or other astrophysical sources. If BH
superradiance is discovered with continuous GW searches at LVK and future GW
observatories, follow-up electromagnetic observations can discover the
kinetically mixed phenomena presented in this paper or put robust constraints
on the dark photon parameter space.

Finally, note that in this work we have exclusively focused on signals from
stellar BHs, but the computations and the results for the luminosity presented
in Eq. (61) also apply to the case of isolated supermassive BHs, which would
be of relevance for follow-up signatures of mergers detected by LISA.

##  VIII Discussion and Conclusions

In this paper, we performed a detailed analytical and numerical study of the
dynamics of a kinetically-mixed, dark photon superradiance cloud and the
resulting multimessenger signatures. If a dark photon with Compton wavelength
on the order of the radius of a stellar mass black hole exists, isolated,
rotating black holes turn into very bright electromagnetic sources for kinetic
mixings on the order of $\varepsilon\sim 10^{-10}-10^{-6}$. We summarize our
approach, discuss the unique dynamics and observational signatures of the
system, and outline future directions for investigation below.

A dark photon superradiance cloud is a dark electromagnetic field rotating at
a frequency fixed by the dark photon energy. Charged particles in the
interstellar medium enter the rotating cloud and, in the presence of a kinetic
mixing, initiate a transient phase of cascade particle production, resulting
in a dense plasma. The stable rotation of the cloud, large strength of the
visible electromagnetic fields, and the appearance of the charged plasma
resemble pulsar magnetospheres. We use resistive electrodynamic methods that
interpolate between the vacuum and the force-free limit. Our methods are
adapted from those originally developed to model resistive effects in the
pulsar magnetosphere, here applied to a fixed Kerr spacetime and the
kinetically mixed case.

Due to the system’s large conductivity, plasma currents redistribute charge in
an approximately dipolar form on large scales in order to screen the coherent
visible electromagnetic field set up by the rotating superradiance cloud. We
find that, due to the electric dominance of the oscillating superradiance
cloud, the charged plasma cannot completely screen the electromagnetic field.
Our numerical simulations of the field electrodynamics indicate that the
differential rotation between the background fields and the plasma leads to
the emergence of small-scale turbulence in the form of magnetic field
reconnection and unscreened electric fields (Sec. V.3). We establish that the
electromagnetic emission from the system is dominated by such small-scale
turbulent dissipation into the standard model plasma, with a peak luminosity
of up to 10 orders of magnitude larger than the solar luminosity, as described
in Sec. VI.2 (see Eqs. (54) and (61)). In addition, we find a Poynting flux
component to the emission. Though subdominant, $\sim 10^{2}-10^{4}$ times
smaller than the local dissipation, this may also be significant for the
observational signatures of the system. These are the main results of our
paper.

Given the rotation of the superradiance cloud, we expect that at least a
fraction of our system’s emission is periodic (with period $T\approx
2\pi/\mu$), a property that would extend the analogy with pulsars into the
observational domain. Our simulations show some evidence for such a
periodicity, as discussed in Sec. VI.4, but cannot conclusively verify this
expectation. By analogy with the results for pulsar magnetospheres and
findings from PIC simulations of turbulent plasmas, we expect that the
emission spectrum contains a non-thermal X-ray component for kinetic mixing
parameters $\varepsilon\gtrsim 10^{-7}$, softer spectra for smaller mixing
parameters, and likely a radio component (see Sec. VI.5). Our simulations
cannot directly determine the spectrum of the emitted luminosity, so given the
differences between our system and pulsars, both in field profile morphology
and strength, more investigations are essential to conclusively establish the
spectral shape.

Based on the system’s large luminosity, characteristic time evolution, and
expected periodicity, we propose several search strategies for dark photon
superradiance clouds. The first strategy relies on the extreme brightness of
our system: electromagnetic follow-up searches targeting binary BH mergers
observed by the LVK Collaboration (Sec. VII.1). Based on our expectations of
the spectrum, the most promising reach is achieved by X-ray and radio
observations. In addition to the spatial and temporal correlation of a merger
and the resulting electromagnetic emission, one could discriminate between
superradiance clouds and standard astrophysical sources by requiring a fast
exponential rise, and $1/t$ power-law fall-off of the light curve, consistent
with superradiant growth and subsequent decay by GW emission. Given the
measured remnant BH mass and spin from the gravitational waveform, these
timescales are fully determined by the dark photon mass; for the parameters
necessary to see these signals over cosmological distances, the
electromagnetic power is subdominant to the GW power and does not affect the
time evolution. This results in a non-trivial relation between the rise and
decay times (see Eq. (64)), which if experimentally confirmed, would provide a
smoking gun signature of superradiance.

Further assuming periodicity, a variety of additional signatures can be
explored. For the proposed electromagnetic follow-ups, the periodicity alone
could be used to measure the dark photon mass $\mu$, as the periodicity is set
by the dark photon mass up to a few-percent binding energy corrections from
the BH potential. This would further strengthen the evidence for the new
physics hypothesis by requiring consistency between the period-based
measurement of $\mu$ and the light curve rise and decay times.

Another observational strategy targets known pulsars, either by selecting
those with positive frequency derivative, or by selecting those with the same
measured period up to $\mathcal{O}(10^{-3})$ due to Doppler shifts.
Interestingly, by surveying existing pulsar catalogues we find four candidate
sources for which a partner with the same frequency exist, and that have a
positive frequency derivative measurement. Our computation of the GW power
with the measured “pulsar” period and spin-up rate suggest that many of these
objects would emit continuous gravitational wave with strain above the LVK
threshold if their origin is dark photon superradiance. Many potential sources
could be excluded with blind continuous wave searches while others could be
further probed with targeted continuous wave searches (section VII.2).

Apart from these multi-messenger signatures of correlated EM and GW emission,
at small $\alpha$ and large $\varepsilon$ the dark photon cloud depletion is
dominated by EM emission, leading to persistent ultra-bright sources of
X-rays. Observation of several periodic sources with periods within a few
percent of each other would be strong evidence for the origin of these objects
being a dark photon superradiance cloud. The absence of these sources in the
Universe may already imply constraints on the dark photon parameter space
assuming a natal black hole mass and spin distribution; given the
uncertainties in BH properties and the EM signal spectral shape, we leave this
study to future work.

We close with a discussion of future directions. In this paper, we studied the
case of an isolated BH and a single new particle with a minimal interaction,
the kinetically-mixed dark photon. In the presence of a more complex dark
sector, such as the existence of a dark Higgs, the dynamics may be altered
further. Recently, some of the authors showed that in the higgsed dark photon
scenario, vortex lines can form and deplete the cloud before it reaches its
maximum size East (2022); East and Huang (2022). Further study is needed to
understand how the vortex production and evolution is affected by the presence
of plasma (see, e.g., Refs. Adelberger  _et al._ (2007); East and Huang
(2022)).

Furthermore, many BHs are not isolated: they are surrounded by an accretion
disk, leading to additional dynamics. While it is a negligible perturbation to
the background dark photon field, a dense accretion disk affects the visible
electromagnetic field and resulting electromagnetic dynamics and emission.
Furthermore, preliminary studies show that even a small superradiance cloud
(or a small $\varepsilon\sim 10^{-12}$ for a maximal cloud) affects the disk
properties, and would invalidate the ISCO-based spin measurements Li  _et al._
(2005); McClintock  _et al._ (2014); Reid  _et al._ (2014); we leave further
details to a future publication Huang  _et al._ .

This work sets the stage for the exploration of kinetically-mixed dark photon
superradiance. We have established numerical simulation techniques to
understand the electrodynamics of the cloud and its electromagnetic emission.
The parameter space of possible signatures is vast, from weak, long-lasting
signals to bright, short signals, across the electromagnetic spectrum. We have
detailed several possible observations, and a more wide-ranging study of
search strategies is warranted. Utilizing different simulation approaches to
better understand the electromagnetic emission spectrum would be invaluable to
pin down observables. Finally, given the similarities between our system and a
magnetar, it is conceivable that the dark photon superradiance cloud could
host a wide range of astrophysical phenomena, such as X-ray flares and fast
radio bursts. Due to the differences in the magnetic field structures, further
investigations are clearly needed, and would be of great interest.

###### Acknowledgements.

We are extremely grateful to Andrei Gruzinov for many invaluable conversations
and suggestions. We thank Asimina Arvanitaki, Liang Dai, Savas Dimopoulos,
Anson Hook, Matthew Johnson, Robert Lasenby, Luis Lehner, Elias Most, Peter
Michelson, Michael Müller, Zhen Pan, Sasha Philippov, Geoffrey Ryan, Kendrick
Smith, and Ling Sun for many useful discussions. MB is supported by the U.S.
Department of Energy, Office of Science, Office of Basic Energy Sciences
Energy Frontier Research Centers program under Award Number DE-SC0022348 and
through the Department of Physics and College of Arts and Science at the
University of Washington. MB, DEU, and JH thank the Institute for Nuclear
Theory at the University of Washington for its kind hospitality and
stimulating research environment. The INT is supported in part by the U.S.
Department of Energy grant No. DE-FG02-00ER41132. JH and CM thank the Center
for Computational Astrophysics and New York University for hospitality. The
Center for Computational Astrophysics at the Flatiron Institute is supported
by the Simons Foundation. NS and WE acknowledge financial support by the
Natural Sciences and Engineering Research Council of Canada (NSERC). DEU is
supported by Perimeter Institute for Theoretical Physics and by the Simons
Foundation. Research at Perimeter Institute is supported in part by the
Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Colleges and Universities. This research was undertaken thanks in
part to funding from the Canada First Research Excellence Fund through the
Arthur B. McDonald Canadian Astroparticle Physics Research Institute. This
research was enabled in part by support provided by SciNet (www.scinethpc.ca),
Compute Canada (www.computecanada.ca), and the Digital Research Alliance of
Canada (www.alliancecan.ca). Simulations were performed on the Symmetry
cluster at Perimeter Institute, the Niagara cluster at the University of
Toronto, and the Narval cluster at the École de technologie supérieure in
Montreal.

##  Appendix A Notation

We list the variable definitions used throughout the text in Table 2.

Variable | Description
---|---
$\varepsilon$ | Kinetic mixing
$\mu$ | Massive vector field mass
$r_{c}$ | Cloud’s Bohr radius
$a_{*}$ | Black hole dimensionless spin
$\Delta a_{*}$ | Loss of black hole dimensionless spin
$M$ | Black hole mass
$r_{g}$ | Half black hole Schwarzschild radius
$\omega$ | Superradiance cloud angular frequency
$f$ | Superradiance cloud frequency
$\Omega_{\rm BH}$ | Black hole horizon frequency
$M_{c}$ | Superradiance cloud mass
$\alpha$ | Gravitational fine structure constant
$t,x,y,z$ | Cartesian Kerr-Schild coordinates
$\sigma$ | Plasma conductivity
$\gamma_{e}$ | Electron/positron Lorentz factor
$m_{e}$ | Electron/positron mass
$e$ | Positron charge
$\tau_{\rm plasma}$ | Timescale to populate the $e^{\pm}$ plasma
$A^{\prime}_{\mu}$ | Dark vector potential (interaction basis)
$E^{\prime i},\textbf{E}^{\prime}$ | Dark electric field (interaction basis)
$B^{\prime i},\textbf{B}^{\prime}$ | Dark magnetic field (interaction basis)
$A_{\mu}$ | Visible vector potential (interaction basis)
$E^{i},\textbf{E}$ | Visible electric field (interaction basis)
$B^{i},\textbf{B}$ | Visible magnetic field (interaction basis)
$I^{\mu}$ | Electromagnetic 4-current
$J^{i}$ | Electromagnetic spatial current
$\rho_{q}$ | Charge density
$v^{i}_{d}$ | Plasma drift velocity in Eulerian frame
$f_{\rm GW}$ | Gravitational wave frequency
$P_{\rm GW}$ | Gravitational wave luminosity
$\tau_{\rm GW}$ | Gravitational radiation timescale
$\tau_{\rm SR}$ | Superradiance instability timescale
$P_{\rm EM}$ | Visible Poynting flux
$L_{\rm diss}$ | Visible energy dissipation power
$\tau_{\rm EM}$ | Electromagnetic radiation timescale
$\nu$ | Spectral frequency of electromagnetic emissions
Table 2: List of the variables used most commonly throughout the main text, as
well as a brief description.

##  Appendix B Construction of massive vector cloud

To construct the superradiantly unstable massive vector field modes on a fixed
Kerr spacetime of mass $M$ and dimensionless spin $a_{*}$, we follow Refs.
Lunin (2017); Frolov  _et al._ (2018); Krtouš  _et al._ (2018), as well as
Dolan (2018); Siemonsen and East (2020). Neglecting all non-linear effects
both on the massive vector field $A^{\prime}_{\mu}$, as well as in the
gravitational sector, we treat $A^{\prime}_{\mu}$ as a test field on a fixed
Kerr background spacetime. In this Ricci-flat geometry, the vector field
satisfies the massive vector wave equation

$\displaystyle
g^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}A^{\prime\gamma}=\mu^{2}A^{\prime\gamma},$
(66)

with $g^{\alpha\beta}=g^{\alpha\beta}_{\text{Kerr}}$ and mass parameter $\mu$.
The Kerr family of BH spacetimes, a special case of the Kerr-NUT-(A)dS class
of spacetime, admits additional symmetries beyond stationarity and axisymmetry
generated by the timelike and axial Killing fields, $\boldsymbol{\xi}$ and
$\boldsymbol{\eta}$, respectively. These “hidden” symmetries were utilized in
Ref. Frolov  _et al._ (2018) to construct an ansatz for the field
$A^{\prime}_{\mu}$, satisfying the massive vector wave equation, that
separates radial and angular dependencies. The ansatz makes use of the
Killing-Yano symplectic 2-form $\boldsymbol{h}$, whose tensor components
satisfy $\nabla_{\mu}h_{\nu\gamma}=2g_{\mu[\nu}\xi_{\gamma]}$. The vector
field ansatz of frequency $\omega$ and azimuthal index $m$ reads Frolov  _et
al._ (2018)

$\displaystyle A^{\prime\mu}=B^{\mu\nu}\nabla_{\nu}Z,$ $\displaystyle
Z=R(r)S(\theta)e^{-i(\omega T-m\varphi)}$ (67)

in Boyer-Lindquist (BL) coordinates $(T,r,\theta,\varphi)$. The polarization
tensor $B^{\mu\nu}$ is implicitly defined by
$B^{\alpha\beta}(g_{\beta\gamma}+ih_{\beta\gamma})=\delta^{\alpha}_{\gamma}$
utilizing the Killing-Yano 2-form. Due to the presence of $\boldsymbol{\xi}$
and $\boldsymbol{\eta}$ in the spacetime, the temporal and azimuthal
dependencies are trivially satisfied. The radial and polar dependencies are
determined by solving a nonlinear ordinary differential eigenvalue problem
with eigenfunctions $R(r)$ and $S(\theta)$ and complex eigenvalues $\omega$
and $\nu$ of the form, $\mathcal{D}^{r}_{\omega,\nu}R(r)=0$ and
$\mathcal{D}^{\theta}_{\omega,\nu}S(\theta)=0$, where
$D^{r,\theta}_{\omega,\nu}$ are second order ordinary differential operators
depending on $r$ and $\theta$ only. Bound state solutions are obtained by
imposing ingoing radiation boundary conditions at the horizon, and
asymptotically flat boundary conditions at spatial infinity. For details on
how the polar equation is solved, see Ref. Dolan (2018). The radial solution
is expanded around the outer horizon, $r\geq r_{+}$, with
$r_{\pm}=M\pm\sqrt{M^{2}-a^{2}}$, in a Frobenius series of the form Dolan
(2018)

$\displaystyle
R_{\text{near}}(r)=\hat{r}^{i\kappa}(1+\hat{a}_{1}\hat{r}+\hat{a}_{2}\hat{r}^{2}+\dots),$
(68)

with $\kappa=2mr_{+}(\omega-m\Omega_{\rm BH})/(r_{+}-r_{-})$, coefficients
$\hat{a}_{i}$, and $\hat{r}=(r-r_{+})(r_{+}-r_{-})^{-1}$. The coefficients
$\hat{a}_{i}$ can be solved for by plugging $R_{\text{near}}(r)$ into the
radial equation $\mathcal{D}^{r}_{\omega,\nu}R(r)=0$. The near-horizon
solution (68) is then used to numerically integrate the radial second order
ordinary differential equation outwards from $r_{s}=r_{+}+\epsilon$ towards
large $r\gg 10/(\alpha\mu)$, with $\epsilon=10^{-4}M$ typically, in the spirit
of the shooting method. Integration cannot start at $r=r_{+}$, as the BL
coordinates are singular on the event horizon. For further details, see Refs.
Dolan (2018); Siemonsen and East (2020). The superradiantly unstable vector
cloud can then be reconstructed with Eq. (67).

Our numerical setup, outlined in App. C, utilizes the Kerr spacetime in
Cartesian Kerr-Schild (KS) coordinates $(t,x,y,z)$. Therefore, the relevant
spacelike hypersurface is the surface of constant KS coordinate time $t$, not
BL time $T$; these are two different slices of the Kerr spacetime. Therefore,
we transform the above constructed vector field $A^{\prime\mu}$ from BL
coordinates to KS coordinates. The two gauges are related by

$\displaystyle\begin{aligned} t=&\
T+\frac{M^{2}\log\frac{r-r_{+}}{r-r_{-}}}{\sqrt{M^{2}-a^{2}}}+M\log\Delta,\\\
x=&\ \sin\theta(r\cos\bar{\phi}-a\sin\bar{\phi}),\\\ y=&\
\sin\theta(a\cos\bar{\phi}+r\sin\bar{\phi}),\\\ z=&\
r\cos\theta,\end{aligned}$ (69)

where $\Delta=a^{2}-2Mr+r^{2}$. The inverse of these relations is

$\displaystyle\begin{aligned} T=&\
t-\frac{M^{2}\log\frac{r-r_{+}}{r-r_{-}}}{\sqrt{M^{2}-a^{2}}}-M\log\Delta,\\\
r=&\ 2^{-1/2}\Big{[}-a^{2}+x^{2}+y^{2}+z^{2}\\\ &\
+\sqrt{4a^{2}z^{2}+(a^{2}-x^{2}-y^{2}-z^{2})^{2}}\Big{]}^{1/2},\\\ \phi=&\
\arctan\left[\frac{rx+ay}{-ax+ry}\right]-\frac{a\log\frac{r-r_{+}}{r-r_{-}}}{2\sqrt{M^{2}-a^{2}}},\
\ \ x>0\\\ \theta=&\ \arccos\frac{z}{r},\end{aligned}$ (70)

where
$\bar{\phi}=\phi+(2\sqrt{M^{2}-a^{2}})^{-1}a\log[(r-r_{+})(r-r_{-})^{-1}]$,
valid outside the event horizon $r>r_{+}$. Using these coordinate
transformations, we transform the massive vector field to KS coordinates for
$r>r_{+}+\epsilon$. Subtleties arise for $r<r_{+}+\epsilon$. However, it is
necessary for the successful evolution of the system of equations in the
interaction basis that the source term is defined (at the very least)
everywhere outside the event horizon $r\geq r_{+}$, and while all modes on the
horizon are marginally trapped, the finite difference length scale
(potentially) allows for values at points inside the horizon to numerically
affect points just outside the horizon. Furthermore, the exponential blue-
shift captured in the transformation rules (69) and (70) of the type
$\cos\log(r-r_{+})$ for $r>r_{+}$ leads to an exponential amplification of any
truncation error in the numerical solution of the radial and angular equations
$\mathcal{D}^{r}_{\omega,\nu}R(r)=0$ and
$\mathcal{D}^{\theta}_{\omega,\nu}S(\theta)=0$. Hence, any power-like
converging numerical truncation error is exponentially enhanced for
$r\rightarrow r_{+}$, and dominates the solution for some $r_{d}$ with
$r_{d}>r>r_{+}$. To address this subtlety, we employ a $C^{4}$-transition
function $f_{t}(r)$ that matches the Frobenius solution $R_{\text{near}}(r)$,
valid for $r\gtrsim r_{+}$ and the numerical solution $R_{\text{num.}}(r)$,
valid for $r>r_{+}+\epsilon$, in the overlap region
$r_{+}+\epsilon<r<r_{+}+10^{2}\epsilon$:
$R_{\text{matched}}(r)=f_{t}R_{\text{near}}(r)+(1-f_{t})R_{\text{num.}}(r)$.
In addition to these manipulations at small radii, we also need to address the
exponential fall-off as $r\rightarrow\infty$. Due to finite floating point
precision, the shooting method will inevitably switch from the exponentially
decaying solution into the exponentially diverging solution at some large
$r=r_{\text{max}}$. Therefore, in order to provide sensible estimates also for
$r>r_{\text{max}}$ (which is necessary since we are working with a
compactified setup that includes spacelike infinity, as outlined below in App.
C), we fit an exponential of the form $ae^{-br}$ for $b>0$ to the solution
$R_{\rm num.}(r)$ in the range $r\in(0.9r_{\text{max}},r_{\text{max}})$. With
this, we obtain a $R_{\rm matched}(r)$ that is valid for $r\in(r_{+},\infty)$.

Finally, the 3+1 superradiant vector variables are projected with respect to
the $t$=const. spacelike hypersurface using the hypersurface normal $n^{\mu}$
and the projector
$\gamma^{\alpha}{}_{\beta}=\delta^{\alpha}_{\beta}+n^{\alpha}n_{\beta}$
(further details can be found in App. C). In this framework, the vector field
decomposes into $\chi_{\phi}\equiv-n_{\mu}A^{\prime\mu}$ and
$\chi^{i}\equiv\gamma^{i}{}_{\mu}A^{\prime\mu}$, which are reconstructed from
the matched and extrapolated radial solution, in conjunction with the polar
solution, the transformation rules (70) and (69), and ansatz (67), everywhere
in the $t=$ const. spacelike hypersurface with $r>r_{+}$. Lastly, we find
better convergence properties of the $D_{i}B^{i}=0$ constraint close to the BH
event horizon if a buffer region between the event horizon and the excision
surface in the BH interior is used. To that end, we utilize second order
extrapolation of all 3+1 variables along lines of constant $\theta,\varphi,t$
from $r>r_{+}$ to $r_{-}<r<r_{+}$ to ensure well-defined gradients at $r_{+}$,
we utilize zeroth order extrapolation from a distance of $\delta=10^{-4}M$
away from the spin-axis to set the cloud values on the axis in KS coordinates,
since the BL coordinates also exhibit a coordinate singularities at the poles.
The time-dependence of the cloud in KS coordinates is then simply
$\chi_{\phi},\chi_{i}\sim e^{i\omega t}$.

$\alpha$ | $\omega M$ | $a_{*}$
---|---|---
0.1 | $0.099485$ | $0.382787$
0.2 | $0.195543$ | $0.678411$
0.3 | $0.283390$ | $0.857953$
0.4 | $0.357498$ | $0.946250$
Table 3: The properties of the four clouds used in the main text. We consider
only the $m=1,\hat{n}=0$ and $S=-1$ superradiant vector boson clouds in the
language of Ref. Siemonsen and East (2020) [corresponding to
$(j,n,l,m)=(1,1,0,1)$ used in Sec. II] around spinning BHs of mass $M$ and
dimensionless spin $a_{*}$. The saturation condition, $\omega=\Omega_{\rm
BH}$, fixes the spin $a_{*}$ for each $\alpha$.

In this work, we consider only those superradiance clouds that arose from the
fast growing modes, and subsequently saturated the superradiance condition,
i.e., satisfying $\omega=m\Omega_{\rm BH}$. That is, we focus on
$m=1,\hat{n}=0$, and $S=-1$ clouds, in the language of Siemonsen and East
(2020) (corresponding to $(j,n,l,m)=(1,1,0,1)$, introduced in Sec. II).
Properties of the clouds considered in this work are given in Table 3. These
states are expected to be the endstate of the instability to a good
approximation East (2017); East and Pretorius (2017), and hence, exhibit
vanishing growth rates.

Figure 20: The norm $L_{D_{d}}(\tau)$ of the residual of the massive vector
wave equation, defined in (71), as a function of the number of grid points
$N_{p}$ considered with respect to the base resolution $N_{B}$. The default
resolution used for each of the configurations given in table 3 and presented
throughout the main text is $N_{p}/N_{B}=2$.

We test our numerical implementation of the above described reconstruction of
the cloud in KS coordinates by considering the numerical truncation error
$\tau$ of the massive vector wave equation (66) in KS coordinates:
$\tau=\sum_{\alpha}|(\square_{\text{Kerr,KS}}-\mu^{2})(\chi_{\phi}n_{\alpha}+\chi_{\alpha})|$.
Inside the event horizon, the truncation error $\tau$ is divergent with
decreasing grid spacing $h$, which is likely due to inaccuracies in the
extrapolation procedure described above. On the event horizon, the truncation
error shows marginal point-wise convergence
$\tau|_{r=r_{+}}\sim\mathcal{O}(h^{1/2})$ or better. To quantify the
convergence behavior outside the event horizon, in the coordinate domain $D$,
we introduce the norm $L_{D}(f)$ of function $f$ as

$\displaystyle L_{D}(f)\equiv\int_{D}d^{3}x\sqrt{\gamma}|f|,$ (71)

where $d^{3}x\sqrt{\gamma}=d^{3}x\sqrt{\det\gamma_{ij}}$ is the volume form of
the spacelike hypersurface. Here, and in the following, the domain $D=D_{d}$
is the coordinate shell defined by the radii $10r_{c}>r>r_{+}$ (recall,
$r_{c}=1/(\mu\alpha)$ is the cloud’s Bohr radius). The norm of the truncation
error converges approximately as $L_{D_{d}}(\tau)\sim\mathcal{O}(h^{2})$ for
$\alpha\in\\{0.2,0.3\\}$ and as $L_{D_{d}}(\tau)\sim\mathcal{O}(h^{1})$ for
$\alpha\in\\{0.1,0.4\\}$, see Figure 20. We use second order accurate methods
to compute the residual $\tau$. This convergence behavior can be explained by
considering the shooting method underlying the reconstructed solution. As
discussed above, any remaining truncation error close to the horizon is
exponentially amplified, while at large radii, the shooting method inevitable
switches from the exponentially decaying solution to the exponentially
increasing solution due to finite floating point accuracy (see also a
discussion of this in Siemonsen and East (2020)). While our shooting method
implementation makes use of higher-than double-precision floating point
arithmetic, the convergence is ultimately limited, especially for small
$\alpha$. Further, the extrapolation of the cloud to the spin-axis in KS
coordinates, and the exponential extrapolation at large radii, can further
reduce the convergence rate, especially at large $\alpha$. Therefore, even
slow convergence is a good sign of correct implementation of the superradiant
source terms in our compactified KS evolution implemented discussed next.

##  Appendix C Numerical evolution setup

Our calculations are carried out on a fixed rotating BH background of mass $M$
and dimensionless spin $a_{*}$. We adopt the Kerr metric $g_{\alpha\beta}$ in
Cartesian Kerr-Schild coordinates $(t,x,y,z)$, with line element

$\displaystyle\begin{aligned} ds^{2}=&\ -dt^{2}+dx^{2}+dy^{2}+dz^{2}\\\ &\
+\frac{2Mr^{3}}{r^{4}+a^{2}M^{2}z^{2}}\Big{[}dt+\frac{z}{r}dz\\\ &\
+\frac{r(xdx+ydy)}{r^{2}+a^{2}M^{2}}-\frac{aM(xdy-
ydx)}{r^{2}+a^{2}M^{2}}\Big{]}^{2},\end{aligned}$ (72)

where $r$ satisfies $(x^{2}+y^{2})/(r^{2}+a^{2})+z^{2}/r^{2}=1$. However, we
evolve the fields using a $3+1$ space–time decomposition, making using of
several such geometric quantities. This includes the unit normal to slices of
constant time, which can be decomposed as $n^{\mu}=(1,-\beta^{i})/N$, where
$N$ and $\beta^{i}$ are the lapse and shift vector, respectively. There is
also the spatial metric/projection operator
$\gamma_{\mu\nu}=n_{\mu}n_{\nu}+g_{\mu\nu}$, the extrinsic curvature tensor
$K_{ij}=-(D_{t}\gamma_{ij})/(2N)$, where
$D_{t}=\partial_{t}-\mathcal{L}_{\beta}$ with $\mathcal{L}_{\beta}$ the Lie
derivative along the shift. In a Minkowski spacetime in Cartesian coordinates
these 3+1 variables would be $N=1$, $\beta^{i}=0$, $\gamma_{ij}=\delta_{ij}$,
and $K_{ij}=0$. For convenience, we also define also the covariant derivative
$D_{i}$ defined with respect to $\gamma_{ij}$, as well as the trace of the
extrinsic curvature $K=K_{ij}\gamma^{ij}$.

The kinetically mixed Maxwell’s equations, presented in covariant from in
(13), in terms of the visible electric and magnetic fields are given by

$\displaystyle\begin{aligned} D_{t}E^{i}=&\
NKE^{i}+\varepsilon^{ijk}D_{j}(NB_{k})-NJ^{i}+N\varepsilon\mu^{2}\gamma^{i}{}_{\mu}A^{\prime\mu},\\\
D_{t}B^{i}=&\ NKB^{i}-\varepsilon^{ijk}D_{j}(NE_{k}),\\\ D_{i}E^{i}=&\
\rho_{q}+\varepsilon\mu^{2}n_{\mu}A^{\prime\mu},\\\ D_{i}B^{i}=&\
0,\end{aligned}$ (73)

where $\varepsilon^{ijk}=\varepsilon^{ijk\alpha}n_{\alpha}$ is the
3-dimensional Levi-Civita tensor. These evolution equations (73) are
discretized using fourth-order accurate spatial finite-difference stencils, in
conjunction with fourth-order Runge-Kutta integration in time East  _et al._
(2015); East and Yang (2018). All simulations are performed on a 3D Cartesian
grid that includes spatial infinity through the use of compactified
coordinates (details can be found in Pretorius (2005)). Sixth-order, Kreiss-
Oliger-type numerical dissipation is applied to the evolution variables for
numerical stability. This further dissipates shorter wavelength features at
large distances beyond the compactification scale, minimizing reflection off
spatial infinity. We use between ten and seven mesh-refinement levels centered
on the BH, with refinement ratio $2:1$, for $\alpha=0.1$ through $\alpha=0.4$,
respectively. We choose the finest level to have length roughly twice the
diameter of the BH in each linear dimension. The system is rescaled so that
$20\mu^{-1}\alpha^{-1}=20r_{c}$ is roughly equal to the compactification scale
(recall $r_{c}$ is the cloud’s Bohr radius). This allows us to resolve both
the scale set by the BH, as well as that set by the superradiance cloud, for
sufficiently long times as to ensure relaxation into a quasi-equilibrium
state. For all cases, we use a grid spacing of $\Delta x\approx 0.03M$ on the
finest mesh refinement level. Due to this scaling, radiation extraction can be
done up to a distance of $r\leq 10r_{c}$; beyond this distance, high-frequency
radiation modes are no longer sufficiently resolved in the wave extraction
zone due to the compactification of the domain. The Maxwell equations with an
Ohm’s law become stiff in the high conductivity limit (discussed further in
App. D). Hence, at conductivies of $\sigma/\mu>2$, we adjust the time-step
$\Delta t$ to account for this behavior. For $\alpha=0.3$, we decrease it
gradually with increasing conductivity from $\Delta t/\Delta x=0.5$ down to
$\Delta t/\Delta x=0.075$ to achieve a robust numerical evolution. For
$\alpha\in\\{0.1,0.2,0.4\\}$ and $\sigma/\mu=20$, we scale $\Delta t$, such
that $\sigma\Delta t$ remains as small, or smaller than, the value of
$\sigma\Delta t$ for $\alpha=0.3$, everywhere in the relevant computational
domain. The construction of the superradiance cloud is described in detail in
App. B. As we are neglecting backreaction of the presence of the plasma and
massless electromagnetic fields, the superradiance cloud is not evolved
numerically, rather it is a pre-prescribed function of time.

The set of equations (73) is comprised of the two constraints, the Gauss law
for electric and magnetic fields, and the Faraday equation and Ampere’s law as
evolution equations. Numerically, we damp possible violations of the
constraint equations by means of two constraint-damping fields $\Phi$ and
$\Psi$ Dedner  _et al._ (2002); Palenzuela  _et al._ (2009). To that end, we
perform the replacements $D_{t}B^{i}\rightarrow D_{t}B^{i}-ND^{i}\Phi$ and
$D_{t}E^{i}\rightarrow D_{t}E^{i}-ND^{i}\Psi$ at the level of the evolution
equations in (73). Furthermore, we promote the constraint equations to
evolution equations for these auxiliary fields, following Refs. Dedner  _et
al._ (2002); Palenzuela  _et al._ (2009):

$\displaystyle\begin{aligned} D_{t}\Psi=&\
-N(D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}-\rho_{q})-N\kappa\Psi,\\\
D_{t}\Phi=&\ -ND_{i}B^{i}-N\kappa\Phi.\end{aligned}$ (74)

This ensures that any numerical violation of the constraints
$D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}-\rho_{q}=0$ and
$D_{i}B^{i}=0$ are damped exponentially over timescales $1/\kappa$. For all
(resistive) force-free simulations, the constraint
$D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}-\rho_{q}=0$ is trivially
satisfies since the charge density is defined to be
$\rho_{q}=D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}$. Hence, unless we
explicitly assume vacuum (and in particular, set $\sigma=0$), $\Psi$ is
identically zero. In all cases, the initial conditions for these auxiliary
fields is $\Phi=\Psi=0$. In addition, we perform ideal force-free simulations
by means of two ad-hoc field modifications applied at each grid point after an
evolution step Palenzuela  _et al._ (2009, 2010) (see Ref. Alic  _et al._
(2012) for a discussion):

$\displaystyle E^{i}\rightarrow$ $\displaystyle\
E^{j}\left(\delta^{i}_{j}-\frac{B_{j}B^{i}}{B^{2}}\right),$ (75)
$\displaystyle E^{i}\rightarrow$ $\displaystyle\
E^{i}\left\\{1-\hat{\theta}(\lambda)+\hat{\theta}(\lambda)\frac{B}{E}\right\\},$
(76)

where $\lambda=E^{2}-B^{2}$ and $\hat{\theta}$ is the Heaviside function. This
prescription enforces the two force-free conditions, $E_{i}B^{i}=0$ and
$B^{2}>E^{2}$, by explicitly rescaling the electric field at each grid point.
The rescaling is a form of ad-hoc numerical dissipation that is not physically
motivated and reproduces physical dissipation behavior only in special cases.
Therefore, as pointed out in the main text, the dissipation estimates provided
by these force-free evolution schemes should be interpreted with caution.

The evolution of the system proceeds as follows. We evolve the system towards
its equilibrium state in several steps. Initially, we set the fields to
visible fields to zero $E^{i}=B^{i}=0$ and evolve until time $t=t_{s}$
assuming vacuum $I^{\alpha}=0$. With $t_{s}\approx 5/\mu$, this allows the
system to equilibrate at roughly $E^{i}=\varepsilon E^{\prime i}$ and
$B^{i}=\varepsilon B^{\prime i}$, which is purely the superradiance cloud’s
contribution to the visible fields. During this time, we utilize both the
electric and magnetic field’s Gauss constraint cleaning potentials $\Phi$ and
$\Psi$. These ensure that constraint violations in the magnetic field are kept
small, as well as efficiently remove constraint violations of the $E^{i}=0$
initial data on timescales $1/\kappa\ll t_{s}$. At $t=t_{s}$, the resistive
current (37) is discontinuously turned on, and the system is evolved until the
total Poynting flux at the largest radii where we extract it is relaxed to a
quasi-constant value. For $\alpha=0.4$, we found that this required the system
to be evolved for $\sim 8t_{\rm LC}$ light crossing times of the entire cloud,
defined as $t_{\rm LC}=10r_{c}$, whereas for $\alpha=0.1$, we evolved the
system for $\sim 3t_{\rm LC}$. During the evolution of the system, we monitor
the behavior of the Gauss constraint $D_{i}B^{i}=0$ throughout the entire
domain. This provides a measure for the rate of convergence of the numerical
solution, and the self-consistency of the numerical implementation.

Figure 21: (right and top left) We present the convergence of the norm (71) of
the Gauss constraint $D_{i}B^{i}$ within domains $D_{s}$ (solid lines) and
$D_{d}$ (dashed lines) for the $\alpha=0.3$ system at three different
conductivities $\sigma$. The former is the entire domain outside the outer
horizon with $r>1/\mu$ (i.e., neglecting the near-horizon behavior), whereas
the latter is the entire domain outside the BH $r>r_{+}$. $N_{p}/N_{B}$ is
defined as in Figure 20, and we set $\varepsilon=10^{-6}$. (bottom left) We
plot the Poynting fluxes across the ergosurface, $P_{\rm ES}$, the event
horizon, $\dot{\mathcal{E}}_{\rm BH}$, as well as the flux $P_{\rm EM}$
coordinate radii $r=8r_{c}$ and $r=10r_{c}$ in the $\alpha=0.3$ cloud with
$I^{\mu}=0$ and $\varepsilon=10^{-6}$.

To test the numerical implementation of the kinetically mixed force-free
Maxwell equations, together with the reconstruction of the massive vector
field modes, we begin by considering the $\alpha=0.3$ vacuum case $\sigma=0$.
To that end, we consider (73) with $I^{\mu}=0$, and perform a series of
simulations with increasing resolution starting from vacuum initial data:
$B_{i}=E_{i}=0$. To demonstrate the correct implementation of the equations,
we monitor the evolution of the constraints $D_{i}B^{i}=0$. In Figure 21, we
show the convergence of this Gauss constraint with grid spacing $h$ utilizing
the norm defined in (71). As can be seen there, roughly
$L_{D_{s}}(D_{i}B^{i})\sim\mathcal{O}(h^{3.5})$, whereas
$L_{D_{d}}(D_{i}B^{i})\sim\mathcal{O}(h^{2})$. Therefore, in the bulk of the
vector cloud, the convergence is roughly fourth order, as expected. Close to
the event horizon, convergence is slower. This may be attributed to the
presence of the excision surface close to the event horizon, as well as the
lower convergence order of the massive vector cloud residual $\tau$ just
outside the horizon. The constraint
$D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}=0$ (not shown here) is
violated by the initial data, but quickly becomes dominated by converging
truncation error after a few periods of the cloud. This could be improved
upon, by choosing constraint satisfying initial data. However, the goal of
this work is to investigate the system with non-vanishing charge density, and
at late times, the latter constraint is satisfied to floating point accuracy
if $I^{\mu}\neq 0$. Moving to cases with non-vanishing conductivity, we show
in Figure 21 the convergence properties of $\alpha=0.3$ systems for
$\sigma/\mu=1$ and $\sigma/\mu=20$. The former exhibits the same convergence
behavior as the vacuum case, whereas the latter follows roughly
$L_{D_{d}}(D_{i}B^{i})\sim\mathcal{O}(h^{0.6})$. This can be explained by
considering the scales of the features that need to be resolved. As we showed
in the main text, the scale of the features is roughly given by $1/\sigma$,
which translates into a length scale of $\ell\approx 0.17M$ for $\alpha=0.3$.
On the finest level, the grid spacing, $\Delta x=0.03GM$, is sufficient to
resolve these features, while on coarser levels, numerical dissipation likely
dampens these scales efficiently. This damping is numerical, rather than
physical, in nature, and does not obey the Maxwell equations, and therefore,
leads to a larger violation of the Gauss constraint and worse convergence
properties.

Apart from the convergence of the constraints, the vacuum quasi-stationary
state, after several $t_{\rm LC}$, should exhibit no energy flux across the
event horizon, as it is, by construction, synchronized with the BH angular
velocity. In practice, there are various sources of numerical error that can
spoil this property. The synchronization condition $\omega=\Omega_{\rm BH}$
can be achieved only up to finite precision, when solving for the
superradiance cloud. Finite resolution both in the evolution scheme, as well
as in the cloud construction scheme, may also leave room for the solution to
develop a small, but finite, energy flux across the horizon and towards
spatial infinity. To quantify this, and to obtain a rough estimate for the
time at which the system is truly settled, $t_{\text{settle}}$, we monitor the
energy fluxes across the horizon, the BH ergosurface, and coordinate spheres
of radius $\hat{\rho}=8r_{c}$, and $\hat{\rho}=10r_{c}$ in Figure 21. In the
continuum limit, with $\omega=\Omega_{\rm BH}$ exactly, we expect all these
fluxes to be zero. Therefore, the flux evolution presented in Figure 21 can be
used to establish $t_{\rm settle}$, i.e., when the system has reduced the
superradiance cloud’s emission powers to the degree necessary.

We briefly comment on issues related to performing simulations of the
superradiance cloud system in the small-$\alpha$ limit. In the Newtonian
limit, the massive vector wave equation on a Kerr background is obtained by
expanding in small $\alpha$ to leading order. All spin-effects are subleading
in this expansion, and the leading contribution is solely given by the far-
zone Newtonian potential of the BH $\sim GM/r$. In this limit, the vector wave
equation reduces to a radial Schrödinger-type equation with solution (8).
Within a numerical time-domain evolution setup, the singular behavior of the
Newtonian potential at the location of the BH poses challenges. However, there
are subtleties associated with replacing the far-zone weak-field metric by
Minkowski both in the interaction and mass eigenbases ((12) and (86),
respectively). Within the mass eigenbasis, the force-free condition (or
resistive generalization thereof) $F^{\alpha\beta}I_{\beta}=0$, depends on the
electric and magnetic field components of (8). The non-relativistic field (8)
and its electric and magnetic field components are multivalued, i.e.,
discontinuous, at the origin, leading to a breakdown in the validity of
numerical schemes around the origin. Additionally, the usual force-free
current would require modification, as it requires the field (8) to satisfy
the corresponding Maxwell equations on a weak-field background spacetime. As
noted above, a weak-field metric is numerically challenging to implement, such
that the choice of a flat background introduces inconsistencies when using the
(resistive) force-free current. These could be remedied, however, by including
terms involving higher order derivatives of the massive vector field (8) in
the equations, but this would add further complications at the origin. By
contrast, the evolution equations in the interaction basis depend only on
$A^{\prime}_{\mu}$, and not on its spatial derivatives; hence, the interaction
basis evolution approach allows one to evolve the system self-consistently on
a flat background. On the other hand, this choice is accompanied by subtleties
associated with the photon-dark photon interaction term in (12). The massless
state can mix into the massive state as it radiates towards infinity. In a
weak-field metric, this mixing prevents the massive component of the visible
field $A_{\mu}$ from radiating to infinity, as it is bound to the central
gravitational potential. In the flat spacetime limit, leakage of the massive
state into radiation emitted to infinity is not prevented (an illustration of
this behavior is presented in Figure 23). Therefore, in this context, any
radiated Poynting flux is to be understood as an upper bound for the total
emitted power. All these subtleties are absent in the fully-relativistic
calculations we use as our main results, where the relativistic clouds
constructed in App. B is considered on a Kerr BH background spacetime as
described above.

##  Appendix D Resistive force-free currents

In this appendix, we discuss different resistive generalizations of force-free
electrodynamics used in the literature to identify the approach most
applicable in the kinetically-mixed scenario at hand. In the main text, we
demonstrated that the system is characterized by turbulence and magnetic
reconnection with efficient energy dissipation into the plasma. In principle,
there are two feasible approaches for capturing these effects: resistive
magnetohydrodynamics and kinetic PIC methods. PIC methods, which capture the
macro- and micro-physics, are ideally suited to tackle the magnetosphere of
the kinetically mixed superradiance cloud. However, as we are interested in
the overall electromagnetic power-output and large scale features of the
system in three dimensions, PIC simulations are prohibitively computationally
expensive, especially on a curved background BH spacetime. Full resistive
magnetohydrodynamics, on the other hand, is notoriously difficult to apply to
regimes in which the plasma mass density is far below the energy density of
the electromagnetic field, which is the case for the superradiant system
considered here. Hence, we choose to use a resistive approach where the plasma
dynamics is not directly tracked, and rely only on the electromagnetic field’s
evolution. A few approaches have been developed in the literature,
particularly to model the resistive regions of pulsar magnetospheres
Komissarov (2006); Gruzinov (2007); Li  _et al._ (2012); Parfrey  _et al._
(2017); Kalapotharakos  _et al._ (2012) (see also Refs. Palenzuela  _et al._
(2009); Palenzuela (2013)). All are based on an electromagnetic current
$J^{i}$ that aims to capture the physics of a highly conducting plasma in
strong electromagnetic fields, while being specified solely in terms of the
electromagnetic fields. Generally, this current can be decomposed into a piece
describing the drift velocity of the charges, $v_{d}^{i}$, and a contribution
orthogonal to the drift velocity

$\displaystyle J^{a}=\rho_{q}v_{d}^{a}+J^{a}_{\perp}.$ (77)

In all (resistive) force-free approaches, the charge density is defined using
the (kinetically-mixed) electric Gauss law:
$\rho_{q}=D_{i}E^{i}-\varepsilon\mu^{2}n_{\mu}A^{\prime\mu}$. In the
following, we briefly review the currents considered in the literature, apply
these to kinetically-mixed superradiance clouds in the non-relativistic limit,
and compare our findings with the vacuum and force-free limits in order to
evaluate their applicability.

Figure 22: We show the visible electric and magnetic fields, charge density,
and the ratio of electromagnetic field magnitudes obtained from a force-free
simulation using the current (78), with prescription (76) for a superradiance
cloud system with $\alpha=0.3$ and a BH spin of $a_{*}=0.86$. In the case of
the visible electric field, the field line color code is the same as in Figure
5.

We begin with a discussion of the commonly invoked force-free paradigm, which
assumes vanishing Lorentz force $F^{\mu\nu}I_{\nu}$, as well as
$F_{\mu\nu}F^{\mu\nu}>0$ and $F_{\mu\nu}(*F)^{\mu\nu}=0$. In terms of electric
and magnetic fields, the last two conditions are, respectively, equivalent to
magnetic dominance $B^{2}>E^{2}$ and $E_{i}B^{i}=0$. With these assumptions,
the corresponding force-free current is given by (the relativistic
contributions in the form of the extrinsic curvature $K$ and covariant
derivative $D_{i}$ are defined in App. C)

$\displaystyle\begin{aligned} v_{d}^{a}\
&=\frac{\varepsilon^{ajk}E_{j}B_{k}}{B^{2}},\\\ J^{a}_{\perp,\text{FF}}\
&=\frac{B^{a}}{B^{2}}\Big{[}2KB_{i}E^{i}-2K_{ij}E^{i}B^{j}+B_{i}\epsilon^{ijk}D_{j}B_{k}\\\
&\ \
-E_{i}\epsilon^{ijk}D_{j}E_{k}+\varepsilon\mu^{2}B_{i}\gamma^{i}_{\mu}A^{\prime\mu}\Big{]}.\end{aligned}$
(78)

The current is perpendicular to the electric field, and there is no
dissipation of electromagnetic energy. Therefore, the force-free limit of
ideal magnetohydrodynamics excludes resistive processes or transfer of
electromagnetic energy to the plasma (e.g., particle acceleration, magnetic
reconnection, plasma heating, etc.), which however, are active throughout the
superradiance cloud due to the electric dominance of the fields
$A^{\prime}_{\mu}$ in vacuum. Notice, the plasma drift velocity $v_{d}^{a}$ is
entirely determined by the dynamics of the electromagnetic fields. While the
force-free approximation is, in principle, ideal, it can break down at current
sheets and places where magnetic dominance is lost. Numerically, this is
handled with numerical dissipation which is particularly large in turbulence
driven regimes (due to the cascade to short, unresolved wavelengths), as well
as by (as noted in App. C) enforcing the force-free conditions by rescaling
the visible electric field, as shown in (76). We perform a set of force-free
simulations of the $\alpha=0.3$ and $a_{*}=0.86$ cloud-plasma system. As the
turbulent features reach scales much below the grid scale of our simulations,
numerical dissipation and the prescription (76) efficiently remove energy that
was sourced by the superradiance cloud. Therefore, it serves as an artificial
source of dissipation, that nonetheless agrees well with the
$\sigma\rightarrow\infty$ extrapolations shown in Figure 12. Regardless,
results from these force-free simulations should be interpreted with caution
and in light of the un-physical dissipation mechanism. In Figure 22, we show
the force-free solution the system attains at late times (with strong
numerical dissipation in the bulk of the cloud). In all cases, the magnetic
Gauss constraint is non-convergent, while the time-averaged outgoing Poynting
flux and energy injection from the superradiance cloud are roughly consistent
across resolutions to within an $\mathcal{O}(1)$-factor (the Poynting flux
estimates is shown in Figure 12). For all the numerical resolutions we
considered, features emerged on the grid scale, suggesting the endstate of the
pair production cascade is a bulk turbulent state. As can be seen in Figure
22, no large scale electric field and charge separation persists, while
$E_{z}/E\sim\mathcal{O}(1)$. The numerical implementation, by construction,
removes any violation of $B^{2}>E^{2}$ at each grid point after each timestep,
such that the ratio $B^{2}/E^{2}$, shown in last panel in Figure 22, is
strictly larger than unity. Similarly, the violations of $E_{i}B^{i}=0$ are at
the level of floating point error.

In the context of resistive magnetohydrodynamics, a macroscopic resistivity is
introduced by means of a suitably chosen Ohm’s law with conductivity $\sigma$.
In order to recover the force-free approximation in the
$\sigma\rightarrow\infty$ limit and to maintain a form for the electrodynamics
that does not require one to also keep track of the fluid dynamics, all
resistive force-free approaches assume the drift velocity $v^{i}_{d}$ of the
charges is altered as Gruzinov (2007) (see also Ref. Komissarov (2006) for a
similar approach)

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

Thus, even for fields with $E^{2}>B^{2}$, the drift velocity $v^{i}_{d}$ is
bounded by the speed of light due to the additional electric field
contribution $E_{0}^{2}$ in the denominator compared with the force-free
prescription (78). Non-vanishing $E_{i}B^{i}$ can only reduce the resulting
drift velocity further. This ensures that around current sheets within a
strongly magnetized plasma, the characteristic speeds remain physical. Three
distinct methods to construct $J^{i}_{\perp}$ have been considered in the
literature. In Ref. Gruzinov (2008), the Ohm’s law was applied in the frame of
vanishing charge density, referred to as (A) in the following. In Ref. Li  _et
al._ (2012) (see also Ref Lyutikov (2003)), the Ohm’s law was applied in the
minimal velocity (with respect to the “lab” or simulation frame) fluid frame,
labeled as (B) in what follows. Lastly, the approach of Ref. Parfrey  _et al._
(2017) introduces resistive effects with a prescription driving $E_{i}B^{i}$
towards $J^{i}B_{i}/\sigma$ over some arbitrary timescale $1/\kappa$, called
approach (C) from here on. Beyond the drift velocity, the three approaches
(A), (B), and (C) differ.

Comparing the three currents, method (A) is manifestly covariant, but lacks a
well-defined vacuum limit, while both (B) and (C) exhibit $J^{i}\rightarrow 0$
as $\sigma\rightarrow 0$. Since the superradiant system is well-understood a
priori only in the vacuum limit, we focus on (B) and (C) in this discussion.
Explicitly, the orthogonal component of the current (B) constructed in Ref.
Lyutikov (2003); Li  _et al._ (2012) reads

$\displaystyle J^{a}_{\perp,(B)}=\sigma
E_{0}\sqrt{\frac{B^{2}+E_{0}^{2}}{B_{0}^{2}+E_{0}^{2}}}\frac{E_{0}E^{a}+B_{0}B^{a}}{B^{2}+E_{0}^{2}}.$
(80)

The prescription, (C) modifies the force-free contribution as Parfrey  _et
al._ (2017)

$\displaystyle
J^{a}_{\perp,(C)}=\frac{\sigma}{(\sigma+\kappa)}\left(J^{a}_{\perp,\text{FF}}+\kappa
E^{i}B_{i}\frac{B^{a}}{B^{2}}\right).$ (81)

Here, the driving timescale $1/\kappa$ can be understood by contracting the
above current by $B^{a}$, using the Maxwell equations to arrive at

$\displaystyle D_{t}(E_{i}B^{i})=-\kappa
N\left(E^{i}-\frac{1}{\sigma}J^{i}\right)B_{i},$ (82)

such that $E_{i}B^{i}$ is driven towards $J_{i}B^{i}/\sigma$. Note, we assume
the resulting system to be hyperbolic (see e.g., Refs. Komissarov (2002);
Pfeiffer and MacFadyen (2013) for a discussion).

Physically, these currents describe the interaction of an effective plasma
with the visible electromagnetic fields, assuming
$E^{2},B^{2}\gg\rho_{p},P_{p}$, where $\rho_{p}$ and $P_{p}$ are the plasma’s
mass density and pressure. The $\sigma\rightarrow 0$ limit corresponds to the
vacuum limit. To understand this, consider the charge conservation
$\nabla_{\mu}I^{\mu}=0$. In the Eulerian frame, together with (31), this leads
to

$\displaystyle
D_{t}\rho_{q}=N\rho_{q}K-D_{i}(NJ^{i})=N\rho_{q}K-D_{i}(N\rho_{q}v_{d}^{i})=0.$
(83)

Hence, if the initial conditions satisfy $\rho_{q}=0$, then the system does
not acquire a non-trivial charge distribution dynamically, i.e., in a medium
with small conductivity charges cannot separate. This implies that the
resistive currents above reduce to $I^{\mu}=0$, assuming the initial data is
neutral. Therefore, the $\sigma=0$ regime is the vacuum limit of the system.
Moving away from this limit to non-zero, but small conductivities,
$\sigma\ll\mu$, the effective fluid coupling to the visible electromagnetic
fields is an efficient insulator. The current $I^{\mu}$ is timelike, and
$J^{i}$ is advection dominated. Due to the residual conductivity, the charges
in the insulating fluid can move along the electric field with mobility
$\sim\sigma/\mu$, i.e., the charge mobility in the fluid frame is conductivity
suppressed. However, since the system is advection dominated, no large charge
gradients can build up, unless the fluid is compressible. In our case, the
fluid velocity in the Eulerian frame is $D_{i}v_{d}^{i}\neq 0$, resulting in
potential charge pile-up in regions of large $v_{d}$-gradients and
compressibility. At intermediate resistivity, $\sigma\sim\mu$, the insulating
fluid transitions to a moderately conducting plasma. Here, the current
$I^{\mu}$ is both locally spacelike and timelike in different places, and the
system has advection and conduction dominated regions. For $\sigma\gg\mu$, the
plasma turns into a highly-conducting plasma with only residual resistivity.
Here, the advection of the fluid is a negligible contribution to the overall
charge distribution. Large scale charge separation is enabled by large
conduction currents along the electromagnetic fields. In this regime, the
conductivity sets the diffusion length scale $\ell=1/\sigma$ that governs
residual resistive features such as current sheets and tearing modes. Finally,
assuming that in the $\sigma\rightarrow\infty$ limit, the system becomes
largely magnetically dominated and $E_{i}B^{i}\rightarrow 0$ while $\sigma
E_{0}$ remains finite, then all three currents reduce to the familiar and
physically well-defined force-free limit. This is discussed further in the
context of the superradiance system in App. E.

Figure 23: We show the rescaled energy emission rates,
$\hat{P}=P(G/\varepsilon^{2})(M/M_{c})$, of the Poynting flux $P_{\text{EM}}$
(extracted at $r=10r_{c}$), the Joule heating $L_{\text{diss}}$, and
$\mathcal{I}_{EB}=\int d^{3}\sqrt{\gamma}|E_{i}B^{i}|$, as functions of
conductivity $\sigma/\mu$ for model (C) with a superradiance cloud of
$\alpha=0.3$ on a Minkowski background with magnetic guide field
$B_{z}=3\times 10^{4}B^{\prime}_{\rm max}$. The corresponding Poynting fluxes
in the vacuum $P_{\rm EM}^{\rm vac.}$ and force-free limits $P_{\rm EM}^{\rm
ff}$ are indicated for reference.

In all cases, the conductivity $\sigma$ is to be understood as a proxy for a
class of dissipative processes and is chosen to be constant in space and time
(primarily due to the lack of a physically motivated prescription for the
spatial dependence of conductivity in this setup), as as typically done, for
instance in Li  _et al._ (2012); Mahlmann  _et al._ (2021) (see Parfrey  _et
al._ (2017) for a notable exception). The advantage of current (C) is its
numerical properties in the high-conductivity limit. There, due to the
prefactor $\sigma/(\sigma+\gamma)$ multiplying the orthogonal component, the
magnitude of the source of the Maxwell equations remains small, ensuring that
the evolution equations do not become stiff. This ultimately allows us to
evolve the system even at relatively large conductivities with moderate
resolution within an explicit forward integration scheme. However, a drawback
of approach (C) is that $J^{i}_{\perp,(C)}$ diverges wherever $B^{2}=0$. In a
magnetically dominated pulsar magnetosphere, this does not lead to problematic
behavior, while in the case of an electrically dominated superradiance cloud,
this causes issues at moderate and high conductivities, since within the
equatorial plane, the magnetic field of the superradiance cloud smoothly
transitions through zero. We tested explicitly, that this magnetic null line
causes the current (81) to diverge in the intermediate and high conductivity
regime, leading to non-convergent features orbiting in the equatorial plane
(particularly on mesh-refinement boundaries). As the resistive methods
outlined above are designed to remove non-converging behavior in, for
instance, current sheets, and we require convergence of our numerical
implementation in order to validate our findings, approach (C) is not well-
suited to tackle the kinetically-mixed superradiance cloud without
modification. Therefore, we resort to approach (B) and current (80) to model
resistive processes and the electromagnetic field geometries throughout the
superradiance cloud. This evolution method has a stiffness problem at large
conductivities, as outlined in App. C, which ultimately limits our ability to
explore the $\sigma/\mu>20$ parameter space.

We briefly illustrate the shortcomings of performing simulations on Minkowski
spacetime, and the extent to which current (C) can be used in the context of a
magnetic guide field removing magnetic null lines. To that end, we consider a
$\alpha=0.3$ superradiance cloud of the form (8) on a fixed Minkowski
background. The constant magnetic guide field is initialized at the beginning
of the simulations as $B^{i}_{0}=(0,0,B_{z})^{i}$ (where $\hat{z}$ is the
spin-direction of the cloud), with magnitude $B_{z}=3\times
10^{4}B^{\prime}_{\rm max}\gg\varepsilon|\textbf{E}^{\prime}|M$. We test that
the following results are independent of the choice of $B_{z}$, as long as the
guide field magnitude is larger than a threshold, $B_{z}>B_{t}$. Below this
threshold, the electric field $\varepsilon\mathcal{E}^{\prime}$ starts
dominating around the origin of the cloud. With this construction, a series of
simulations is performed varying the conductivity from $\sigma/\mu=10^{-3}$ to
$10^{3}$ within the context of the resistive methods (C) introduced above. In
addition, we also study the vacuum limit $J^{a}=0$, as well as the force-free
limit on this flat background.

In Fig. 23, we show the behavior of the total power output of the system as
function of bulk conductivity in the model (C) in (81). Let us compare these
quantities to those obtained on Kerr spacetime without a guide field and using
model (B) [given by Eq. (80)] shown in Figure 12 and Figure 8. The flat
spacetime guide field setup recovers the correct bulk dissipation component
$L_{\rm diss}=L^{\rm bulk}_{\rm diss}$, both in amplitude and in conductivity
dependence, while the turbulent component $L_{\rm diss}^{\rm turb}$ is absent.
The latter is due to the magnetic guide field removing any magnetically
diffusive regions that might form due to turbulence. The behavior of
$\mathcal{I}_{EB}$ in Figure 23 is entirely analogous to the corresponding
quantity in Figure 8. Lastly, the outgoing Poynting flux in Figure 23 is
constant across decades of conductivities, and agrees well with both the
vacuum and the force-free limits. This illustrates the leakage discussed above
due to the lack of gravitational confining potential in flat spacetime,
filtering out the massive propagating states. This demonstrates that, within
the interaction basis, the flat space solution cannot be used to estimate
physical observables associated with the outgoing Poynting flux. Hence, we
cannot use model (C) without a guide field, as discussed above. However, using
a guide field also does not give the correct answer, as this artificially
removes the turbulent dynamics characterizing the high-conductivity limit of
the system.

##  Appendix E Charge distribution and small conductivity regime

Figure 24: The charge distribution $\rho_{q}$ for various small to large
conductivities in the equatorial plane (bottom row) and in a plane spanned by
the BH spin and an arbitrary superradiance cloud phase (top row). The slices
at varying conductivities correspond to the same superradiance cloud phase. We
focus on a $\alpha=0.3$ and $a_{*}=0.86$ BH-cloud system.

In Figure 24, we illustrate the spatial charge distribution $\rho_{q}$ of the
solution at low and intermediate conductivities. At low conductivity, the
largest charge separation occurs along the spin-axis of the BH. This may be
interpreted as follows: Charge separation is suppressed at high resistivity,
$1/\sigma\gg 1$. However, any residual conductivity can separate charges on
scales $\sim\sigma/\mu$ (assuming vanishing charge diffusion). Any separated
charge distribution advects with the drift velocity of the fluid. In the
presence of sufficiently large fluid velocity gradients (with finite fluid
compressibility) large charge densities may build up. In the superradiance
cloud context, regions of high velocity gradients coincide with regions where
the charge density is largest for $\sigma/\mu=10^{-2}$. Once the charge
density is accumulated, and the fluid velocity varies on scales larger than
the charge distribution scale, the latter is frozen into the flow of the
former and is carried away from the BH along the spin axis. At moderate and
large conductivity, $\sigma/\mu\gtrsim 1$, the resistivity is sufficiently
small as to enable large scale charge separation. While for $\sigma/\mu=1$,
the charge density follows roughly the superradiant electric field morphology
(compare with (18)), for $\sigma/\mu=20$, small scale features begin to
appear, likely driven by the turbulent dynamics inside the plasma.

Figure 25: We show the charge separation $\mathcal{I}_{\rho}$ [defined in
(84)] as a function of conductivity in units of $\mathcal{I}_{E^{\prime}}$
[defined in (85)], as well as the global maximum $\max[I_{\mu}I^{\mu}]$ of the
electromagnetic 4-current (37) and global minimum $-\min[I_{\mu}I^{\mu}]$
(recall, we are using the $-+++$ signature). We focus on a $\alpha=0.3$ and
$a_{*}=0.86$ BH-cloud system and consider conductivities
$\sigma/\mu\in\\{0.01,0.1,1,2,5,10,20\\}$.

Let us demonstrate explicitly that at large resistivity, the plasma cannot
charge separate across scales larger than the charge separation scale
$\sim\sigma$ and that the superradiant electric field is screened efficiently
at large conductivity. To that end, we define the quantity

$\displaystyle\mathcal{I}_{\rho}=\int_{D}d^{3}x\sqrt{\gamma}|\rho_{q}|,$ (84)

measuring the charge separation in a coordinate volume $D$ extending out to
$\hat{\rho}=10r_{c}$ within the slice $\Sigma_{t}$ with volume form
$\sqrt{\gamma}d^{3}x$. This is compared with the charge separation required
screen the superradiant electric field $E^{\prime}_{i}$ entirely [see also
(18)]:

$\displaystyle\mathcal{I}_{E^{\prime}}=\int_{D}d^{3}x\sqrt{\gamma}\varepsilon|D_{i}E^{\prime
i}|.$ (85)

We show the behavior of $\mathcal{I}_{\rho}$ as a function of conductivity,
compared with $\mathcal{I}_{E^{\prime}}$, in the right panel of Figure 25
(recall, $\mathcal{I}_{E^{\prime}}$ is conductivity independent). For
$\sigma/\mu\ll 1$, we find that $\mathcal{I}_{\rho}\rightarrow 0$, indicating
that the system tends to the vacuum solution set by the superradiance cloud.
At intermediate conductivity, the charge separation scales roughly as
$\mathcal{I}_{\rho}\sim\sigma/\mu$. At large conductivity,
$\mathcal{I}_{\rho}\sim\mathcal{I}_{E^{\prime}}$, supporting the conclusion
that, for $\sigma/\mu\gg 1$, the solution exhibits large scale charge
separation that screens the field $E^{\prime}_{i}$ efficiently, even in the
turbulent regime. Furthermore, we also show the behavior of the norm
$I_{\mu}I^{\mu}$ of the electromagnetic 4-current (37) in Figure 25. As
outlined in App. D, for $\sigma\ll\mu$, the current is fluid advection
dominated, where a residual charge distribution is flowing with the fluid on
timelike trajectories. Conversely, for large conductivities the conduction
part of the current starts to dominate, the current becomes spacelike, and the
solution begins to asymptote towards a conductivity independent value of
$\max[I_{\mu}I^{\mu}]$.

Figure 26: (top and middle row) We show the energy dissipation density
$\rho_{\rm diss}$, defined in (48), as a function of conductivity in a slice
spanned by the BH spin and an arbitrary superradiance cloud phase (top row),