Constraints on Relic Magnetic Black Holes

We present current direct and astrophysical limits on the cosmological abundance of black holes with extremal magnetic charge. Because they don't Hawking radiate, much lighter primordial black holes could exist today if they are extremal. The dominant constraints come from white dwarf destruction for intermediate masses, and intergalactic gas heating for heavier black holes. Extremal magnetic black holes may catalyze proton decay, and thus we derive robust limits -- independent of the catalysis cross section -- from the above as well as from white dwarf heating. We discuss other bounds such as those from neutron star heating, solar neutrino production, binary formation and annihilation into gamma rays, and magnetic field destruction. We note that stable magnetically charged black holes can assist in the formation of neutron star mass black holes.


Introduction
A simple, well-motivated extension of the Standard Model of particle physics (and Maxwell's equations) is the inclusion of magnetic monopoles. Magnetic monopoles are a natural consequence of grand-unified theories (GUTs) [1,2], have been searched for by a number of different experiments [3][4][5], and were an original motivation for cosmic inflation [6,7]. One possible manifestation of magnetic charges in our universe is magnetically charged black holes. Primordial black holes could have formed in the early Universe from, for example, large-amplitude, short-wavelength density perturbations or early phase transitions [8,9]. If monopoles existed (or perhaps dominated) at that time, one might expect √ N fluctuations to leave primordial black holes with magnetic charges. Light enough black holes could eventually Hawking radiate until they become extremal, leaving a magnetically charged remnant today, a fate different from that of small uncharged black holes. This is touched upon in Ref. [10]. Additionally near Planck mass primordial black holes could acquire magnetic charge by spontaneously emitting magnetic monopoles, leading their charge to fluctuate and perhaps get caught at a stable extremal value. This formation process is explored for extremal electrically charged black holes in Ref. [11]. More exotic formation scenarios may exist, but we will not address those here. EMBHs are interesting because they are not expected to Hawking radiate, making them stable dark matter candidates at masses ruled out for uncharged black holes. Their large mass compared to that predicted for magnetic monopole also allows many EMBHs to avoid the constraints on magnetic monopoles, effectively providing a way to "hide" magnetic charge. Understanding the behavior of EMBHs gives insight into other exotic heavy magnetic objects, such as some higher dimensional black holes [12].
In this paper, we put constraints on the abundance of primordial black holes with magnetic charge. We focus on extremal magnetic black holes (EMBHs), but point out the applicability and strength of the bounds for magnetic black holes (MBHs) heavier than 10 15 g with non-extremal charges. The dominant bounds come mainly from two sources -Galactic gas cloud heating at large masses (> 4 × 10 12 g) and destruction of white dwarfs (WDs) at lower masses. We present novel constraints from gamma-ray emission for EMBHs at low masses, though they are never dominant. In addition, we touch upon bounds from the destruction of magnetic fields and direct searches at MACRO.
Magnetic monopoles derived from some GUTs are predicted to catalyze proton decay with significant cross sections [13,14]. The bounds on EMBHs change significantly when one assumes catalyzed proton decay in stellar objects -neutron stars, white dwarfs, and the Sun -due to heating and neutrino emission. We scan over cross sections for this catalysis to generate robust bounds on relic abundances of EMBHs, independent of the size of this effect. The dynamics of nucleon decay near EMBHs are highly non-trivial, and we leave a detailed analysis for future work.
This paper is organized as follows. In Section 2, we briefly discuss the relevant physical properties of EMBHs. In Section 3, we describe bounds due to the annihilation of binary EMBHs, the heating of the intergalactic medium, the destruction of WDs, destruction of large scale magnetic fields, and non-detection of magnetic charges neglecting effects from catalysis of proton decay. In Section 4, we include the catalysis of proton decay for all possible cross sections. This weakens the bounds in some mass regimes by, for example, preventing the destruction of WDs. For illustration, we also show, for a fixed catalysis cross section, bounds which arise from neutron star (NS) heating and proton decay in the Sun. In Section 5, we explore bounds on MBHs heavy enough to be stable over cosmological timescales even when sub-extremal. Section 6 includes comments on how to generalize the bounds on EMBH abundance when EMBHs have a non-monochromatic mass spectrum, how EMBHs may seed large scale magnetic fields. Section 7 concludes. We compare our constraints with those from Refs. [15,16] in Figs.1 and 3. Throughout this paper, unless units are explicitly specified otherwise, we take c = = k = 1 and 0 = 1 4π .

Extremal Magnetic Black Holes (EMBHs)
In general relativity, MBHs are described by the Reissner-Nordstrom metric, one that is spherically symmetric and carries an outer event horizon and an inner Cauchy horizon. The extreme limit of the magnetic charge is one in which these horizons merge, producing the metric: where M is the ADM (asymptotic) mass of the black hole and G is Newton's constant. The EMBH has, in natural units, a magnetic charge where m P = G −1/2 = 1.22 × 10 19 GeV is the Planck mass. In these units, Q = 10 5 corresponds to M ∼ 1 g. The minimum value for Q is set by the Dirac quantization Q min = 1 2e ∼ 6, where e is the electron charge [17]. Note that this is a different definition of Q than is used in other works [15,18] on EMBHs.
EMBHs have some interesting and unresolved properties. Being extremal, they will not Hawking radiate [19]. Uncharged black holes with M 5 × 10 14 g are predicted to decay within the lifetime of the Universe. EMBHs are expected to be stable and could, in principle, be cosmological relics. 1 Ref. [18] indicates that there are large regions around EMBHs where the magnetic fields are strong enough to condense electroweak bosons and restore the Higgs field minimum to its origin, making Standard Model fermions massless. Interestingly, one would also expect the sphaleron barrier for baryon number violation to disappear [20] in this region. In addition, the interior of non-extremal magnetically charged black holes will have an inner Cauchy horizon with known instabilities [21]. In general relativity, the interior metric can only be extended in physically nonsensical ways (requiring the existence of an infinite number of other universes). It is thus expected that new physics lives at this horizon, and in the extremal limit, this would be at r = GM . The physics at this surface, for example, could involve the Planckian density [22] and violate baryon number in a significant way or source other fields (black hole 'hair'). Thus significant model dependence dictates the physics at close ranges.
A possibly important feature of EMBHs is their potential for catalyzing proton decay, akin to that in GUT-induced magnetic monopoles [13,14]. It is a non-trivial question here -not only is there model dependence with respect to the properties of the core of the EMBH, but also the path to the core for a charged particle, which will involve gauge-boson interactions at medium and long distances and gravitational and (potentially) black-hole hair interactions at short distances 2 . A full analysis for specific models of EMBHs would be interesting. However, to put a robust bound on these objects, we will allow the proton catalysis cross section to be a free parameter (within reason), and will see it is possible to put more model-independent limits from various astrophysical phenomena.
A number of other interesting properties of EMBHs were explored in Ref. [18]. One is that charged Fermions near an EMBH will experience radically different dynamics. In regions where the magnetic field is larger than the fermion's squared mass, the low-lying energy quantum states are Landau levels and should correspond to two-dimensional fermion states which move radially and have degeneracy of order Q. This is predicted to enhance the Hawking decay rate of non-extremal black holes by a factor of Q. This will have some effect on the discussion of bounds on merger rates.

Bounds on EMBHs without Catalysis
Relic EMBHs are heavy magnetic charges expected to be non-relativistic at all relevant times and follow the dark matter distribution once they decouple from the background plasma. Their physical effects include gas heating, WD destruction, magnetic field destruction, and high-energy particle emission from annihilation. We describe the bounds from these effects in the subsections below and summarize them as constraints on the average relic density of EMBHs as a fraction of the average dark matter energy density, f ≡ ρ bh /ρ dm . We take the dark matter density in the Milky Way to be ρ M W dm = 0.3 GeV cm −3 . In this section, we are ignoring the possibility of proton-decay catalysis. We plot the bounds (without catalysis) for a range of masses in Fig.1.

Gas Cloud Heating
When an EMBH passes through an ionized plasma, it loses energy by exciting long range Eddy currents. This heats the plasma via the kinetic energy of the accelerated electrons. We constrain f by requiring that EMBHs not change the observed temperature of clouds of warm ionized medium (WIM) in the Milky Way. The heating rate from friction must be less than the average cooling rate determined in a survey of the WIM [23].  [15] in black. The dashed bound comes from the Parker effect in the Andromeda Galaxy (see Section 3.4), and the dotted one comes from EMBH annihilation in the Sun producing excess neutrinos (see Section 4.4). We cut off the latter bound at M > 300 g because convection in the Sun can suppress the annihilation rate for these EMBHs.
The friction between a magnetic monopole and a plasma is described in Ref. [24]. We adapt it here for an EMBH with charge Q.
where n e , m e are the electron number density and mass, respectively, T is the gas temperature, v is the velocity of the EMBH, λ D = T 4πnee 2 is the plasma Debye length, and l = 2T πme 1/4 1 v 1/2 ωp is the attenuation length of the plasma with plasma frequency ω p = 4πe 2 ne me . Ionized gas clouds observed on this cooling survey [23] have an average n e ∼ 0.08 cm −3 and T ∼ 6000 K. Using Eq.(3.1) we find that a single EMBH moving through at virial velocity v ∼ 10 −3 will deposit energy into the plasma at a rate of dE dt = 9 × 10 6 M 10 9 g 2 erg s −1 . (3.2) Unregulated by cooling, this will significantly change the temperature and behavior of the gas within δt ∼ T dT /dt ∼ 1 Myr 10 9 g M 1 f . A wide range of EMBHs can, therefore, change the properties of the WIM in the Milky Way in a time comparable to that of other less exotic heating processes [25].
The temperature of WIM clouds is regulated by a variety of heating and cooling processes. The primary one is electron and hydrogen atom collisions with singly ionized carbon atoms. These excite fine structure transitions in the ground state of carbon atoms from the 2 P 1/2 to 2 P 3/2 state. Observations of infrared emission from de-excitations of these carbon atoms are used to estimate the cooling rate in ionized gas clouds [23], which is sensitive to the temperature. The heating rate of the gas is estimated from and should not exceed the observed cooling rate [23]. Doing so would disturb the observed dynamics of the gas, changing both its temperature and cooling rate.
Cooling rates for a variety of WIM clouds are reported in Ref. [23] and grouped by the velocity of the cloud. We compare the EMBH gas heating rate to the cooling rate for low velocity clouds because this gives the strongest bounds. 3 The average cooling rate is (dE/dt) ∼ n H × 10 −25.65± 0.11 0.15 erg s −1 , where n H ∼ 0.25 cm −3 is the number density of hydrogen atoms in a typical cloud [23,25]. The rate is given with 3σ uncertainties. Requiring that EMBH WIM heating not exceed the 3σ upper limit in the cooling rate gives the following constraint: This bound is only sensible if there are enough EMBHs, their density in the Galaxy being ∼ 10 31 f (M/g) −1 pc −3 , in a warm gas cloud to heat it. Warm ionized gas clouds have a typical radius of ∼2 pc and would thus be well populated by EMBHs with masses at least up to 10 15 g [25]. Gas heating effects from EMBHs with M > 10 15 g is discussed in Section 5.1. The above constraint is plotted in Fig.1.

White Dwarf Destruction
"Small" energy injections in the WD core of a Carbon-Oxygen (CO) WD can initiate runaway fusion, which leads to Type Ia supernovae [26]. For example, depositing ∼ 6 × 10 21 GeV within a region of radius of ∼ 6 × 10 −3 cm in less than ∼ 3 × 10 −12 s leads to runaway fusion and supernovae in WDs with masses above 1M [26][27][28]. EMBHs collect inside WDs, and can eventually merge with each other and annihilate. The merger event releases a large amount of electromagnetic (EM) energy, typically more than enough to destroy a WD. We use this phenomenon and the observation of old (>1 Gyr) WDs in the Galaxy to constrain the EMBH abundance. The energy needed to destroy a WD depends on the its mass, with lighter WDs requiring larger energy deposits. Not all EMBHs are heavy enough to supply the energy needed to destroy all WDs. We account for this by considering bounds from two different sets of WDs: a collection of known non-magnetic 1M WDs, and WD J0551+4135, a 1.14M WD thought to have formed from the collision of two smaller CO WDs [29]. We use 1M WDs for multiple reasons. Hundreds such WDs have been observed [30]. Their average age is 1.9 Gyr, though some ages over ∼ 10 Gyr [30]. Even allowing for some uncertainty in their cooling ages, more than half would not exist if EMBHs could destroy them within 1 Gyr. WDs with starting masses above 1.05M are expected to be primarily composed of Oxygen and Neon, instead of Carbon and Oxygen [31], and require much larger energy depositions to trigger runaway fusion [26]. These do not give useful bounds. 1M is about as heavy (having the smallest trigger energy) as CO WDs can initially form.
We also consider bounds from WD J0551+4135, a rare CO WD heavier than 1.05M thought to have formed from the collision of two smaller CO WDs just over 1 Gyr ago [29]. Its high mass of 1.14M makes its trigger energy low enough that EMBHs too light to trigger runaway fusion in 1M WDs can destroy it. Key physical properties of these WDs used for calculations are shown in Table 1. Both sets of WDs accumulate EMBHs, which 4 × 10 12 4 × 10 12 Table 1: Above are the relevant properties of the WDs used to derive bounds in this work. The radius is R wd , ρ c is the density at the core, E T is the threshold energy needed to initiate runaway fusion, τ dif f is the diffusion time, and B c is the estimated magnetic field at the core. The EMBH masses M min , M T h , and M max are the minimum mass needed to trigger runaway fusion, the threshold mass required to overcome the WD core magnetic field, and the maximum mass that can sink to the core in less than a Gyr, respectively. The first column is for generic non-magnetic 1 M CO WDs and the second is for WD J0551+4135. The radii, densities, and diffusion times are derived using WD profiles from Ref. [28].
sink to the core, merge and initiate runaway fusion. WDs are composed of a degenerate Fermi gas of electrons and an ion plasma of nuclei. Friction with the Fermi gas and ion plasma traps traversing EMBHs. The stopping power of a degenerate Fermi gas at zero temperature for a monopole is given in Ref. [32]. Adapting this for an EMBH, we find: where v f is the Fermi velocity, E f is the Fermi energy, and ν ei σ ∼ 10 17 s −1 is the electron ion scattering frequency [33] at relevant WD densities and temperatures [34] (we use 10 6 g/cm 3 and 10 6.8 K respectively). For these parameters the energy loss rate is where ρ wd is the average white dwarf density. The velocity is scaled to the escape velocity at the surface of a 1M WD, v esc ∼ 2.3 × 10 −2 , for illustrative purposes. The above description of the friction is appropriate when the energy gained by electrons passing the EMBH, ∼ pv, is large compared to the temperature of the Fermi gas. Here, p is the momentum of passing electrons. When the energy exchanged is small compared to the temperature of the Fermi gas, finite temperature effects must be considered. The stopping power for a magnetic charge in a Fermi gas with a finite temperature is shown in Appendix A to be is the fraction of states occupied for a given electron energy E , and is the energy gained by a passing electron. E f is the Fermi energy of the gas, and x and y are respectively the angles between the electron's incoming and outgoing momentum and v. The total friction in a WD can be estimated by combining the contributions from (3.6) and (3.1) for a plasma of carbon nuclei. The lightest EMBHs we consider, M ∼ 2 × 10 −4 g ∼ 10 20 GeV, have kinetic energies within the Galaxy of ∼ 6 × 10 13 GeV. Using Eq.(3.5) to make a simplified estimate of the energy losses, we find that the lightest EMBHs get trapped after traversing about ∼ 1000 km inside a WD. As heavier EMBHs lose kinetic energy over even shorter distances, it is clear that every one which passes through a WD becomes gravitationally bound to it. WDs accumulate EMBHs at a Sommerfeld enhanced rate: where R wd and v esc are the WD's radius and escape velocity, respectively.
Once inside, EMBHs must sink to the WD core in order to pair and annihilate. The EMBH mass and WD the core density and magnetic field determine an EMBH's behavior. The core magnetic field separates lighter EMBHs of opposite charge, preventing them from annihilating initially. Heavier EMBHs can overcome the magnetic field, and friction with the WD directs their behavior instead.
To annihilate, light EMBHs must first sink to where the magnetic and gravitational forces from the WD balance the attractive force toward other EMBHs. The strength, prevalence, and structure of WD magnetic fields remain an area of active research. The surface magnetic field of WD J0551+4135 is ∼ 50 kG [29]. Current measurement techniques only detect surface magnetic fields and are only sensitive to those 1 kG [35]. When considering 1M WDs, we only use those with magnetic fields too small to measure ( 1 kG). We argue that the core magnetic field is equal to the observed surface field. This is justified because magnetic fields in WDs are thought to either be fossil fields left over from the progenitor star or the result of differential rotation between proto-WDs in a binary system and the common envelope they formed in [36]. In the former case, the WD magnetic field arises from conserving the magnetic flux of the larger progenitor. Modeling indicates that the poloidal fossil field can only persist over long time scales if it is supported by a twisted toroidal field outside of the core [37]. In the latter case, the magnetic field comes from a magnetized accretion disk, which settles onto the surface of the WD [38]. In both cases, the magnetic field is generated well outside of the core, leaving no reason to expect the core magnetic field to be much larger than that at the surface.
A sufficiently strong poloidal magnetic field will separate oppositely charged EMBHs and prevent them from interacting. Interactions between EMBHs of the same sign charge are comparatively weak and neglected in this estimate. One can estimate where EMBHs will settle by balancing the attractive forces of gravity and the opposite-sign EMBHs with the separating force of the magnetic field: (3.10) Here, B int is the magnetic field at the center, ρ c is the local density, r is the radial distance from the core and N is total number of EMBHs in the core, which we assume are split evenly by charge. Note that the attraction between EMBHs and the opposite charge cluster is doubled due to the combined magnetic and gravitational attraction. If only one EMBH is present in the star, it will settle at r = 3 4π B int m P ρc , which for a 1 M WD with a core density of ∼ 3 * 10 7 g/cm 3 and a core magnetic field strength of 1 kG is r ∼ 0.03 cm. We determine the time needed to reach this region using Eq.(3.4) and a density profile derived in Ref. [28] and find all EMBHs light enough to be separated by the core magnetic field sink to their equilibrium point in well under 1 Gyr.
As N increases, the stable solution to Eq.(3.10) moves closer to the center and then disappears. At this point, EMBHs fall to the center of the WD instead of separating into polarized regions. This happens when which is met when N satisfies Here 'Min' refers to the minimum as a function of r. EMBHs in separated clusters moving at thermal velocity do not have enough kinetic energy to overcome the magnetic field and reach each other until Eq.(3.12) is satisfied. When it is satisfied, mutual attraction between the EMBH clusters dominate over other forces, and annihilation proceeds quickly. Only a single annihilation is needed to destroy the WD. EMBHs heavy enough to satisfy Eq.(3.12) when N = 1 are not affected by the core magnetic field. They fall toward the core and interact with opposite sign EMBHs when their mutual attraction becomes stronger than their gravitational attraction to the WD core. This happens once they fall within a radius r * . For both the 1 M and 1.14 M WD, the time EMBHs need to sink to r * scales approximately linearly with their mass, with 4 × 10 12 g being the heaviest EMBHs that can reach the core within a Gyr. The large friction with the WD should dissipate the angular momentum and kinetic energy of infalling EMBHs, leading them to sink directly to r * . They then only need to traverse r * to meet and annihilate. We can conservatively estimate the time to reach each other as r * /v * , where v * is the EMBH velocity at the core, which scales roughly as v * ∼ 10 −21 (M/10 12 g) −1 in the 1.14 M WD case and faster in the 1 M WD case due to their lower core density. Even the heaviest, slowest EMBHs that make it to the core in under a Gyr annihilate in well under a Gyr.
When two opposite charge EMBHs collide, they release energy through a magnetic "chirp" during the collision and through Hawking radiation of the non-extremal black hole left behind afterward. The magnetic dipole they represent vanishes in about t s ∼ 2r s = 4GM ∼ 10 −38 M g s, the time it would take an EM signal to traverse the Schwarzschild radii, r s = 2GM , of the two EMBHs. One can expect an EM "chirp" from the sudden changing and rearranging of the magnetic field, akin to the gravitational chirp that happens during the final moments of an uncharged binary merger and the subsequent ringdown. A magnetic dipole, m, will radiate energy at a rate dE dt = 2 3 d 2 m dt 2

2
. Just before colliding, two EMBHs separated by a distance 2r s will form a dipole with a dipole moment 2r s * Q ∼ 4 M 2 If this were to suddenly disappear in t s , a burst of 2 3 M ∼ 3.7 × 10 23 M g GeV of EM energy would be emitted. The energy needed to destroy each of the WDs is shown in Table 1. The "chirp" alone allows any EMBH with M 0.02 g to furnish the energy needed to destroy a 1M WD and with M 8 × 10 −5 g to supply the energy needed to destroy WD J0551+4135.
EMBHs can also supply enough energy to destroy a WD through the Hawking radiation emitted by the non-extremal remnant left behind by the merger. Ref. [28] suggests a black hole of mass M inside of a WD will Hawking radiate at a rate of ∼ 1.4 × 10 49 GeV If we assume the EMBHs involved in the merger have approximately equal masses, then the remnant left behind would then be non-extremal and need to lose ∼ M of mass to become extremal again. These radiate at nearly the same rate as uncharged black holes, and may radiate even faster if we consider Q enhancements to the decay rate described in Refs. [15,18]. Assuming the merger remnant radiates at least ∼ 1.4 × 10 49 GeV EMBHs with masses as low as 0.012 g for 1 M WDs and 2 × 10 −5 g for WD J0551+4135 will be able to supply the energy needed to destroy their host WD. Works such as Ref. [39] suggest that the semi-classical description of Hawking emission is not correct and that black holes may radiate different amounts of energy, with different characteristics and over different time scales than predicted above. If true, this changes the mass range of EMBHs able to destroy a white dwarf. If the emission rate slows down, then the heaviest EMBHs may not radiate quickly enough to initiate runway fusion. If the amount of radiation is reduced, then the lightest EMBHs may not emit enough energy to trigger fusion. Precisely how this mass range would evolve depends on the details of how Hawking emission changes outside of the semi-classical approximation, which remains uncertain [39]. This does not particularly affect the bounds here, as most EMBHs able destroy a WD through Hawking radiation can also do so through the magnetic chirp produced when they annihilate.
We limit the abundance of EMBHs so that WDs do not accumulate enough to overcome the magnetic field barrier and annihilate within 1 Gyr (see Eq.(3.12)). The bound for EMBHs too heavy to be affected by the WD core magnetic field is set by requiring that WDs not accumulate 6 EMBHs in 1 Gyr. In a group of 6 EMBHs there is a > 95% chance that at least one has the opposite charge as the others and is able to annihilate. The existence of hundreds of non-magnetic WDs with masses around 1 M and cooling ages longer than 1 Gyr, without any apparent suppression based on age or mass [30], rules out EMBH parameter space where annihilation within 1 Gyr would be expected. For the 1 M WDs this bound is and for WD J0551+4135 this is (3.14)

Binary Mergers and Gamma-Ray Emission
Oppositely charged EMBHs can form binary systems in the early Universe. These binaries spin down and eventually merge. Assuming the EMBHs involved have approximately equal and opposite charges, the merger leaves behind a new non-extremal black hole. The new black hole will have lost a "chirp mass" fraction of its mass to EM and gravitational radiation and then, unprotected by an extremal charge, would rapidly Hawking radiate down to its new extremal mass. Conservatively, we assume it will radiate ∼ M , the mass of the EMBHs. The hot radiation from these merger remnants gets reprocessed by the thermal bath, and reaches Earth as a gamma-ray signal. We can place constraints on the EMBH abundance by requiring this to be less than the 2σ upper limit on the diffuse extra-galactic gammaray flux reported in the most recent isotropic background analysis of data collected by the Fermi-LAT Telescope [40].
A binary decouples from the Hubble flow and begins spiraling inward when two opposite charge EMBHs' free-fall velocity toward each other exceeds the Hubble velocity pulling them apart. The unique redshift, z d , when this happens is determined by M and by the co-moving separation between the EMBHs, x. The time for a binary to merge, t m , is determined by x, z d , and the distance y to the next nearest EMBH. EMBHs do not interact strongly with others of the same sign charge, so only one EMBH in the binary will interact with the next nearest charge. The next nearest EMBH will tug the in-falling opposite sign one out of its free-fall path, causing the newly decoupled system to form an eccentric binary instead of immediately annihilating in a head-on collision. The closer the next nearest EMBH, the smaller the eccentricity of the binary system. Less eccentric binaries take longer to spin down. This all affects the merger rate per black hole for a given redshift, Γ m (z), which is derived in full in Appendix B.2.
When a binary merges, it produces a sub-extremal MBH that begins to decay. An uncharged black hole decays in [19,28] which is a reasonable estimate for the time a non-extremal MBH needs to decay down to it's near-extremal mass (this is, in fact, an over estimate, as we only included degrees of freedom up to MeV). This is fast compared to the age of the Universe for all black holes constrained by radiation from this decay (roughly 10 7 g), so it is reasonable to treat the entire decay as simultaneous with the merger. Recent works [15,18] suggest that the decay rate for non-extremal magnetic black holes is enhanced by a factor of the black hole charge, Q, This would reduce the decay time by 1 Q , further justifying our decision to treat all decays as instantaneous.
Sub-extremal black holes in this mass range radiate both electromagnetically and hadronically. The hadronic radiation quickly fragments into more stable particles, including neutral pions. These then decay into 2 photons. A gamma-ray signal comes from both the photons radiated directly by the black hole and those resulting from pion decays. We estimate the combined photon spectrum d 2 Nγ dEγ dt using the derived analytic expressions for a decaying black hole's photon spectra found in equations (30)-(37) of Ref. [41]. While the merger remnant's mass still exceeds its charge by an O(1) factor it will radiate with approximately the same surface temperature as a non-magnetic black hole. We estimate the total photon spectrum produced during this decay as dNγ dEγ t bh d 2 Nγ dEγ dt , as the black hole will radiate most of it's energy in this time and around its initial mass. Ref. [39] suggests that the semi-classical description of Hawking radiation may be inaccurate. If true, one would estimate the EMBH bounds using a modified photon emission spectrum derived from the new description of black hole emission. The bounds on EMBH abundance primarily depend on the total energy of radiated photons able to scatter with the thermal bath, those with E th = 36000 1+z GeV, where z is the emission redshift. These bounds will scale inversely with the fraction of hawking radiation emitted above E th .
Photon-photon interactions with CMB photons and Inverse Compton scattering off electrons reprocess the initial photon energy distribution from the binary annihilation. The reprocessed gamma-ray spectral energy distribution has the form [42,43] where E th = 36000 1+z GeV is the threshold energy a photon must have to be reprocessed, E γ is the energy of the photon, and z is the redshift at the time of emission. Photons with initial energies below E th free stream to Earth [43], while those at higher energies are reprocessed. A total reprocessed energy spectrum can be found by multiplying L(E γ , z) by the amount of energy radiated above E th . The shape of the spectrum is independent of this energy. The parts of the initial spectrum, N Einit (E γ , z), that were and were not reprocessed combine to give a total visible spectrum (3.17) Here, θ(x) represents the Heavyside function. The flux per unit energy of gamma-rays that reach Earth can be found by convoluting the reprocessed spectrum with the merger rate as a function of redshift: where H(z) is the Hubble rate at z, and E γ0 is the photon energy today, and Γ m (z) is the merger rate per black hole at z. z m = E th (z=0) E γ0 − 1 is the highest redshift at which a photon with observed energy E γ0 could be emitted and still reach Earth. Any photon with E γ0 that reaches Earth today must have been emitted at a redshift z < z m because those emitted at earlier times will be reprocessed and redshifted down to lower energies. The total flux per Fermi energy bin is found by integrating F γ (E γ 0 ) over the energy range for that bin. We require the integrated gamma-ray fluxes to be less than the 2σ upper bound observed by Fermi. The strongest constraints come from the highest energy bin, 580 GeV-820 GeV.
We have, so far, neglected friction with the thermal bath. As explained in Section 3.1, EMBHs interact strongly with plasmas such as the thermal bath. Friction slows infalling EMBHs in binaries, causing them to decouple at later times and get stretched further apart by the Hubble flow. The stretching reduces the number of binaries that form and causes those that do to merge at later times. The result is that EMBHs that are still strongly interacting with the thermal bath at z d , their friction-free decoupling redshift, do not significantly contribute to either the merger gamma-ray signal or to the EMBH abundance bound. We consider an EMBH to be strongly interacting with the thermal bath when friction more strongly influences its velocity than Hubble acceleration, dE dx f /(M V ) > H(z). The frictional force grows as M 2 , making heavier EMBHs more sensitive to it. Friction begins to weaken the bound for EMBHs with M ∼ 10 6 g and removes it entirely for those with M > 10 7 g. The green region shows the bounds on the EMBH abundance, f , based on gamma-ray emission from mergers of EMBH binaries not strongly influenced by friction with the thermal bath. The dotted orange line gives a rough estimate of the bounds that arise from gamma-ray emission from binaries whose initial infall dynamics were strongly influenced by friction with the thermal bath. The friction-free bounds weaken for high f for EMBHs with 10 6 g < M < 10 7 g because less eccentric binaries are more sensitive to frictional effects and raising f lowers the typical eccentricity.
The bounds are shown in Fig.2. The most notable feature is that they weaken at higher f for EMBHs with masses 10 6 − 10 7 g. Increasing f decreases the typical spacing between EMBHs, which reduces the distance between a binary and the next nearest EMBH, lowering the eccentricity of the most eccentric binaries. The merger rate is dominated by these binaries, so decreasing the maximum eccentricity slightly lowers the merger rate. This is normally insignificant, but becomes important when friction suppresses binary formation.
Weaker bounds may be found for heavy EMBHs (M > 10 7 g) by accounting for friction. These EMBHs form bound pairs if they can traverse the distance separating them within one Hubble time while infalling at a terminal velocity set by friction, their charge, and their initial separation. Friction reduces the number of bound pairs, so the constraint from gamma-ray emission for EMBHs with M > 10 7 g is weak compared to the other new constraints presented in this paper. For comparison, we show a rough constraint for EMBHs strongly affected by friction in Fig.2, but do not include this in our final set of constraints.
Ref. [18] suggests that fermionic radiation from near extremal MBHs is be enhanced by a factor of Q, reducing the black hole lifetime by a factor of 1 Q , without changing the surface temperature. This should not significantly change the the photon spectrum that reaches Earth from mergers of EMBHs with M < 10 7 g. The low energy portion of the spectrum comes from hadronic decays. If the surface temperature remains unchanged, the hadronic spectrum should have the same form, so the low energy part of the spectrum is unchanged.
Reprocessed photons contribute to the high energy part of spectrum. The Q enhancement may suppress the number of direct photons produced due to the reduction in the black hole lifetime. However, electrons and positrons would be produced at a Q enhanced rate, and produce the same spectrum as photons after reprocessing by the thermal bath [43].

Destruction of Large Scale Magnetic Fields
EMBHs interact with and can destroy large scale coherent magnetic fields. The bounds from this are subdominant to those derived in this paper, so we only touch on them briefly.
Parker Bound There are multiple versions of the Parker bound [44][45][46] on the abundance of magnetic monopoles based on their interaction with the Milky Way's magnetic field. This field accelerates magnetic charges, and drains its energy in the process. The magnetic charges must not destroy the magnetic field or prevent it from forming. The protogalactic Parker bound offers the most stringent monopole abundance constraints [44]. It limits the magnetic charge density present during the collapse of the protogalaxy, when the Galactic magnetic field was only a small seed field unsupported by a Galactic dynamo. When applied to EMBHs in the Milky Way this gives Stronger constraints, f < 4 × 10 −4 [15] and f < 1.7 × 10 −3 [16], were found by applying a Parker type bound to the Andromeda Galaxy whose magnetic field is coherent over larger length scales than those of the Milky Way. Uncertainty in this bound comes from uncertainty in the properties of the magnetic field in Andromeda [15]. We show the bound from Ref. [15] in Fig.3 for comparison.
Cluster Magnetic Fields The largest scale magnetic fields observed in galaxy clusters are in radio relics. These are coherent over ∼ 2 Mpc [47,48]. EMBHs in these clusters would form a diffuse magnetic plasma, suppressing magnetic fields that extend over scales larger than its Debye length, λ d . Requiring λ d > 2 Mpc for an EMBH plasma in a typical galaxy cluster gives the bound f < 0.7. (3.20)

MACRO
MACRO was a dedicated monopole detector that ran from 1989 through 2000 and relied on ionization effects to detect passing magnetic charges [3]. While the detectors and analyses focused on Dirac monopoles with Q = g D [3], they were at least as sensitive to higher charge magnetic particles [49], and have been used to constrain a variety of similar objects, including extremal electrically charged black holes [11]. MACRO flux bounds should therefore be applicable to EMBHs. Assuming EMBHs do not cause nucleon decays, to be discussed in Section 4, MACRO places an upper limit on the EMBH flux of [3] This corresponds to a bound

The Effects of Proton Decay Catalysis on EMBH Bounds
As discussed in the introduction, EMBHs might catalyze proton decay, perhaps analogous to the effect in monopoles from grand unified theories [13,14]. However, uncertainty around what lies at the 'core' of the EMBH (the strength of B-violation, the existence of shortdistance hair etc.) makes the cross section for this process model dependent.
We circumvent many of these complications by looking for robust bounds independent of the cross section for catalysis. Inside of WDs, heat produced by catalysis can generate convection currents that prevent the EMBHs from annihilating, weakening the bounds from WD destruction. On the other hand, catalysis can cause anomalous heating of NSs causing x-ray emission. In addition, excess heating of WDs can set a (weaker) bound, while catalysis in the Sun can be constrained by neutrino detectors, such as Super Kamiokande. EMBH abundance constraints from the MACRO detector are weakened but still present when accounting for catalysis effects. Whatever the catalysis cross section, we can place a robust bound on the EMBH abundance at any mass.
We estimate the most robust bound by finding the threshold catalysis cross section needed to activate convection in WDs, σ th . This minimal cross section is then used to find bounds from WD heating. Increasing the cross section above σ th increases the heating in WDs, strengthening the heating bound. Reducing the cross section below σ th will allow EMBH annihilations to resume in WDs. For a set cross section, NS heating generates stronger bounds than WD heating. However, many estimates of the catalysis cross section depend on the velocity of passing nuclei [50][51][52], which are quite different in WD (∼ 10 −4 ) and NS (∼ 0.3) cores, making it unreasonable to assume EMBHs will have the same effective cross sections in both environments. We only compare bounds that rely on catalysis in WDs, where we can give a more consistent description of the cross section. Catalysis bounds in the Sun and in NSs are described and shown for a specific cross section for comparison and future application to specific models of EMBHs.
Whether σ th is physically reasonable depends on two different length scales set by the EMBH and by the properties of the WD core: a geometric scale and a scattering length. One might expect a catalysis cross section at the geometric scale where the energy scale of the magnetic field is larger than the proton mass and the dynamics are better described by two-dimensional degenerate radial fermions [18]. This cross section will be referred to as the QCD cross section.
where Λ QCD ≡ 218 MeV is the approximate QCD scale. All nuclei entering this region will have dynamics dominated by the background magnetic field and could plausibly have non-negligible overlap with a proton-decaying core. We note that this cross section can be altered by an angular momentum barrier between the strong magnetic field of an EMBH and the electric charge and magnetic moment of an incoming nuclei [50]. While σ QCD defines the region where nucleon decay may occur, the strong magnetic field around EMBHs allows them to have even larger effective catalysis cross sections. For example, charged nuclei with non-zero magnetic moments accelerate and radiate photons when passing through the magnetic field near an EMBH, causing them to fall into a bound state which overlaps σ QCD [51,52]. Any enhancement to the effective catalysis cross section is limited by scattering in the WD medium. Nuclei falling toward the EMBH can scatter and regain enough kinetic energy to escape before reaching the inner decay region. This introduces the second important length scale: the scattering length of the nuclei in the local medium, Here, m n is the mass of the nucleus involved, and λ d , T , and ρ c are the Debye length, temperature, and density in the WD core, respectively. The length scales described above set effective upper limits on σ th . Given that σ QCD describes a region where nucleon decay may occur, any value of σ th ≤ σ QCD is physically plausible. Mechanisms exist to allow σ th to grow larger than σ QCD [51,52], but this is limited by nuclear scattering. σ th is only physically plausible if where we define σ M as the maximal allowed catalysis cross section within the WD core. How we determine σ th , the EMBH annihilation rate, and the WD heating bounds is elaborated on in the sections below. Fig.3 shows the minimal bounds once catalysis effects are considered, along with the bounds not impacted by catalysis for a full picture of the constraints if catalysis is non-negligible. We have made the conservative simplifying assumption that once convection starts in a WD, the EMBHs will be spread throughout the star and unable to annihilate. In truth, they may still interact and annihilate under these conditions; however, calculating how the annihilation rate and WD convective currents evolve with the catalysis cross section requires detailed modeling of convection in WD cores, which is beyond the scope of this paper and will only yield more stringent bounds.
For illustration only, we display all of the bounds that arise if we assume the catalysis cross section is fixed at σ QCD in all environments in Fig.4. Here, additional bounds arise from physical processes, such as NS heating and neutrino production from proton decay in the Sun.

Convection in White Dwarfs
As explained in Section 3.2, EMBH abundances are strongly constrained by their annihilation destroying WDs. Catalysis heating bounds are significantly weaker but still relevant. When an EMBH falls into a WD, or any medium, it can cause nearby nuclei to decay into lighter particles, releasing the nucleon mass energy and making the EMBH a source of radiation. The energy radiated can create a severe enough temperature gradient to initiate convection in a WD core [53]. Friction between the EMBHs and convecting material drags them away from center, potentially preventing them from annihilating within 1 Gyr and undermining the WD annihlation bound described in Section 3.2. The catalysis cross section determines the amount of heat generated in the WD core, which sets the temperature  Figure 3: The colored regions denote the minimum bounds on the EMBH abundance when catalysis effects are considered. The WD overheating and annihilation bounds are related by catalysis effects and marked as one in light blue. The bounds from overheating or destroying 1 M WDs dominate for EMBHs with M > 0.012g. Lighter EMBHs cannot destroy a 1 M WD when they annihilate so overheating WD J0551+4135 is used to constrain those with M < 0.012 g. No physically reasonable catalysis cross section could prevent EMBHs with 7 g M 2 × 10 5 g from annihilating in 1 M WDs at all times, so the original annihilation bound remains in a somewhat weakened form. The bounds appear jagged because the catalysis effects were considered at discrete cooling temperatures and would vary smoothly if these effects were considered at arbitrary cooling temperatures and times. There are plausible catalysis cross sections that prevent EMBH annihilation outside of this mass range, and so bounds come from WD heating instead. Bounds from gas heating and binary mergers are unaffected by catalysis and are described in Sections 3.1 and 3.3, respectively. Constraints derived in Ref. [15] are shown in black for comparison and labeled Bai et al 2020. The dashed line comes from a Parker bound applied to the Andromeda Galaxy (see Section 3.4). The dotted line comes from energetic neutrinos produced from EMBH annihilations in the Sun (see Section 4.4). We cut this bound off for EMBHs with M > 300 g because convection in the Sun suppresses their annihilation rate. gradient, which itself determines whether convection occurs. If convection does not happen, then the bound determined in Section 3.2 stands with some modifications. Once convection starts, we assume the EMBHs, now spread throughout the WD, do not annihilate. Estimating the annihilation rate once convection begins requires detailed modeling of how EMBHs move within WD convection currents and would only yield stronger bounds.
We determined whether EMBH annihilations could happen by comparing the minimum luminosity from EMBH induced nuclei decays needed to initiate convection, L th , to L m , the maximum luminosity that can be radiated by EMBHs in the WD core assuming realistic catalysis cross sections and that enough EMBHs are present to annihilate under no catalysis conditions. We determined L th for a variety of WD ages because younger WDs are hotter, and can accommodate more EMBHs before convection turns on. This was done separately for both sets of white dwarfs considered. Throughout the WD's conductive region, the temperature gradient is where L tot is the combined luminosity from all of the EMBHs in the WD, κ is its conductivity, and r is the distance from the center of the WD. We use the expression for κ presented in Ref. [33] for a WD made of equal parts carbon and oxygen with the density profiles for 1.14 M and 1 M WDs derived in Ref. [28]. Two conditions must be met for convection to occur. First, there must be a region in the star which satisfies where T is the temperature, is the adiabatic temperature gradient, is the pressure, and dP dr = GM int (r)ρ(r) r 2 (4.8) is the pressure gradient. M int (r) is the mass contained within the radius r, and ρ(r) is the mass density at r. Starting with a surface temperature, T s , we calculate T and dT dr as a function of r in km steps moving inward from the surface and in m steps for the inner-most kilometer of the WD. At each step we check if the convective condition has been met. WD interiors cool over time, causing older ones to develop more extreme temperature gradients and leading convection to start at lower decay luminosities than in younger stars. Some EMBHs will start convection in cool Gyr-old WDs but not hotter 100 Myr-old ones. To account for this, we considered the temperature gradient for WDs at cooling ages {10 9 , 10 8 , 10 7 , 10 6 , 10 5 } years and their corresponding internal temperatures T s = {10 6.8 , 10 7.3 , 10 7.7 , 10 7.85 , 10 7.95 } K, as estimated in Ref. [34]. If Eq.(4.5) is satisfied near the core, we check for the second condition needed to initiate convection. Buoyant forces acting on the WD material must be greater than the viscous forces resisting deformation. The Rayleigh number, Ra, parameterizes the relationship between these forces. For spherical geometries where is the adiabatic excess at a given radius r [53], χ is the thermometric conductivity, and ν is the kinematic viscosity [54], which we take from Ref. [55]. Ra must be 10 3 for convection to occur in a spherical system where heat is generated in the convection medium [56]. For each cooling temperature and set of WDs, we determined L th , the minimum luminosity where Eq.(4.5) and Ra 10 3 , are satisfied. To determine if L th is physically plausible, we also calculated L m , the maximal luminosity that could be radiated from all of the accumulated EMBHs if they each had the largest possible cross section: Here, N wd is the number of EMBHs present in the WD core, which we take to be the minimum number needed to destroy a WD if catalysis did not happen, set by Eq.(3.12), and v n is the thermal velocity. Convection only starts if L th ≤ L m .
The original annihilation bounds remain unchanged if, at every temperature checked, L th > L m . If L th > L m at some, but not all, temperatures, then the annihilation bound is modified to prevent WDs from accumulating enough EMBHs for an annihilation to occur within the maximum time period for which L th < L m . Finally, if L th ≤ L m at every temperature scale checked, then there exists a physically reasonable catalysis cross section that prevents annihilation at all times. To find this, we first determine the catalysis cross section, σ that corresponds to L th for each mass and cooling age: For a given mass, σ th is the maximum σ found among all of the different cooling temperatures considered. σ th is determined separately for 1 M WDs and WD J0551+4135. The maximum of the two is used to set the catalysis heating bounds in WDs and is referred to as σ H . When M < 0.012 g, the annihilation bound only applies to WD J0551+4135 and so only its value of σ th is used to set σ H . The calculated values of σ th and σ H can be found in Fig.5.

White Dwarf Heating
The heat produced from EMBH catalyzed nucleon decays can make WDs older than a few billion years more luminous than expected. Observations confirm the existence of WDs with luminosities below 10 −4 L [57, 58] and ages above 8 Gyr [30,59,60]. The EMBH abundance must be low enough to allow these old WDs to cool to their observed luminosities and WD J0551+4135 to cool to its observed luminosity, 3×10 −4 L [29], within its estimated lifetime of 1.3 Gyr [29]. The luminosity of EMBHs in a WD core is with v n = 3T wd mn and ρ c as the thermal velocity of nuclei and density in the core, respectively. We take T wd , the temperature of the WD core, to be 10 6.5 K in a WD with L = 10 −4 L and 10 6.8 K in WD J0551+4135, based on modeling in Ref. [34]. The cross sections used for different EMBH masses, σ H , are presented in Fig.5. Some works [50][51][52] suggest that the catalysis cross section can grow as the velocity of passing particles decreases. If such a velocity dependence is present, then σ H for a young hot WD may correspond to a larger cross section in an older cooler one with lower thermal velocities. We use σ H in our estimate because the velocity dependencies are model dependent, the thermal velocities of nuclei in the hottest and coolest WDs considered vary by less than an order of magnitude, and accounting for velocity dependencies will only strengthen the heating constraint.
A WD that has accumulated EMBHs for t acc years without annihilation will radiate with a combined luminosity L cat = LΓ wd t acc , (4.14) where Γ wd can be found in Eq.(3.9). One must keep L cat below 10 −4 L for 1M WDs older than t acc ∼ 8 Gyr and 3 × 10 −4 L for WD J0551+4135, assuming t acc ∼ 1.3 Gyr. For σ th for EMBHs in White Dwarfs Figure 5: The black and gray lines above show the minimum catalysis cross sections, σ th , needed to start convection at all cooling times considered in 1 M WDs and WD J0551+4135, respectively. The red region excludes implausible cross sections in the WDs, with the red line marking σ M . No physically viable value for σ th for 1 M WDs exists for EMBH masses where the black line crosses into the red region. Here, some form of the original annihilation bound stands. Heating bounds are set wherever the annihilation bounds break down (M < 8 g and 2 × 10 5 g< M ) using σ H , which is marked by the green dotted line and follows whichever value of σ th for the two sets of WDs is larger. It follows σ th1.14M at low masses where no annihilation bound exists for 1 M WDs. σ M and σ th are temperature sensitive and are set for WDs at T = 10 7.95 K which corresponds to a cooling time of 10 5 years. This is when the star is hottest and least sensitive to energy injections from catalysis. WD J0551+4135 this age is a conservatively low estimate, as it is thought to result from the merger of two older WDs (1.3 Gyr ago) that would themselves have had more time to accumulate EMBHs [29]. The bounds on EMBH abundances from WD heating can be found in Fig.3.

Neutron Star Heating and Background X-Rays
The effects of magnetic monopoles falling into NSs was explored thoroughly in Ref. [61]. Magnetic monopoles cause nearby neutrons to decay and release their mass energy as lighter particles. The power radiated by all the monopole-catalyzed decays heat the NS, increasing its x-ray radiation. Constraints come from limiting the diffuse x-ray signal to be below the observed soft x-ray background. This same analysis can be applied to EMBHs. This picture is complicated by annihilation effects and uncertainty around what lies at the NS core. EMBH self annihilation reduces the number able to catalyze nucleon decay, and, therefore, the resulting x-ray signal. EMBHs cannot annihilate in NSs with a type II superconducting core, whereas a core that transitions to a type I superconductor or to an exotic state of matter may offer little resistance to annihilation. This uncertainty will be discussed below.
The NS must first capture passing EMBHs. They lose energy due to friction with the NS's electron Fermi gas [62]. We estimate the stopping power using Eq. (3.6), where E f ∼ 100 MeV is the electron Fermi energy, n e ∼ 10 36 cm −3 is the electron density, and v is the velocity of the infalling EMBH which we take to be the escape velocity ∼ 0.3 [63]. Numerically we find that a NS with a 10 km radius can capture the entire relevant mass range of EMBHs.
Once inside, an EMBH radiates with a luminosity where v n ∼ 0.3 is the neutron velocity, ρ ns ∼ 5 × 10 14 g is the core density, and σ cat is the catalysis cross section. For illustrative purposes, we show this bound assuming σ cat = σ QCD . When annihilation effects are not relevant, the EMBHs have a combined luminosity of L tot = N acc L ns , where N acc is the total number of EMBHs accumulated. When annihilation is fast compared to the accumulation rate, the EMBHs that actually remain in the NS all carry the same sign charge. Their number is set by a random walk with N acc steps. In this case L tot = √ N acc L ns . We take N acc = t ns Γ ns , where t ns is the age of the NS, and the EMBH accumulation rate is for an NS with r ns = 10 km and an escape velocity v esc = 0.3. The energy output from all of the EMBHs give the NS a surface temperature where σ b is the Stefan-Boltzmann constant.

Bounds from Overheating Neutron Stars
The bounds on EMBH abundance come from limiting the x-ray radiation from catalysis heated NSs. We follow the analysis in Ref. [61]. A NS heated to a temperature T ns radiates photons with energy E like a black body with a differential luminosity of [61] dL dE = 2 × 10 36 r ns 10km The differential flux that reaches Earth is where l is the distance to the NS and τ (l, E) is the absorption length for photons with energy E in the Galaxy and can be found in Ref. [61].
The catalysis decay luminosity of an individual NS depends on the number of EMBHs it has accumulated, which depends on its age. To get the total x-ray flux from all of the NSs, one must integrate their luminosities as a function of time over the range of possible ages. For simplicity, we assume NSs in the Galaxy are produced at a constant rate over ∼ 10 10 years, and have a local number density of n ns ∼ 10 −4 pc −3 [61]. As we will explain in Section 4.3.2, EMBHs may be free to annihilate in up to 90% of NSs. Therefore when calculating the total radiation from nearby ones, we assume 10% accumulate all EMBHs that pass through them, while the other 90% only retain EMBHs that do not annihilate.
Accounting for absorption in the ISM, the total differential flux from all NSs out to a distance of d ns ∼ few kpc, about as far as the x-ray signal is expected to travel without scattering [61], is where t ns is the age of the NS, t 0 is 10 Gyr, and B(E) is the integrated fraction of photons with energy E that reach Earth: The total flux that reaches Earth will be compared against sounding rocket soft x-ray observations from Ref. [64]. The counting rate per frequency band for a given flux is where A band (E) is the detector response function for each of the four different x-ray detection bands and can be found in Ref. [61].
In large regions of the sky the counting rate measured by Ref. [64] was 1 3 the all-sky average [61]. As in Ref. [61], we set limits by requiring the count rate of soft x-rays produced from radiating NSs to be less than this. We use the weakest of the limits derived from the four bands. The resulting bounds are shown in Fig.4.

Annihilation Effects and the Superconducting Core
Whether EMBHs can annihilate in a NS depends on the properties of its inner core. The inner ∼5 km is thought to be occupied by a type II superconductor [65]. When an EMBH falls into a NS, it sinks into the superconducting surface, which produces 2 flux tubes around the EMBH for each unit of dirac magnetic charge, q d = 1 2e , it carries. Each flux tube carries a magnetic flux of π e and an energy per unit length of ∼ 2 × 10 11 GeV cm −1 [66]. This acts as a tension pulling the EMBH back toward the surface of the superconducting core.
The EMBHs form a suspended surface where gravitational and tension forces balance, ∼ 0.6 km above the center of the NS independent of the EMBH mass. In agreement with Ref. [15], we find that annihilation between suspended EMBHs is too slow to alter the radiation bounds. However, NSs are poorly understood at these depths and may transition to a new state which does not support the EMBH surface and prevent annihilation (such as a type I superconductor [67,68] or some exotic state of matter [69]). Any such transition would reduce the number of EMBHs available to catalyze neutron decay and severely weaken the bounds.
NS cores require densities above ∼ 7×10 14 g cm −3 to transition [67]. While the equation of state and mass range that avoid exotic transitions remains an area of active research [69], lighter NSs are expected to have less dense cores that do not transition [70]. Based on modeling in Ref. [71], we estimate that NSs with masses 1.1M can reasonably be treated as having a core density below the transition threshold. Based on the available mass measurements for NSs [72], we estimate that ∼ 10% have masses below ∼ 1.1M . These NSs do not permit EMBH annihilation and so keep those they accumulate. The remaining 90% amass EMBHs at a random walk rate proportional to the square root of the number that pass through them. We make the conservative assumption that annihilation is fast whenever it is possible, as this gives the lowest estimate for the number of EMBHs accumulated, the smallest NS x-ray luminosity, and the weakest bounds.
At high masses (M 10 9 g), a single EMBH can cause a NS to overheat, so annihilation suppression no longer matters. Only few EMBHs with 10 7 g < M < 10 9 g are needed to overheat a NSs, so N acc is not much smaller than √ N acc , and all NSs contribute to the bound. The bound on EMBHs with M < 10 7 g is dominated by radiation from nonannihilating NSs.

Proton Decay in the Sun, Neutrinos and SuperK
Proton decay in the Sun has been constrained by neutrino observations at Super Kamiokande [73]. Most of the EMBHs we focus on would become trapped while passing through the Sun and cause proton decay once inside. Heavier EMBHs get caught in the convection zone, while those lighter than about 300 g make it to the core and eventually annihilate. Bounds come from the steady-state number in the Sun. These bounds are weaker than those derived for NSs and WDs when the catalysis cross section is held constant. However, some enhancements to the cross section from EMBH-proton bound states [51,52] or angular momentum effects [50] can make the solar neutrino bound stronger than the WD heating bound for some masses. Details of this constraint are included for reference when considering specific EMBH models with specific σ cat . As a toy example, we will assume σ cat = σ QCD for the remainder of this section. We also note that Ref. [15] derives stronger solar bounds by considering neutrino emission from EMBHs that annihilate in the solar core. These bounds should apply to EMBHs which are light enough; however, Ref. [15] does not consider how those heavier than 300 g get trapped in the convection zone which suppresses their annihilation rate and weakens these constraints.
Super Kamiokande is a water Cherenkov detector sensitive to neutrinos with energies above 5.5 MeV [74]. About half of proton decays produce charged pions, which themselves dominantly decay via π + → µ + + ν µ → e + + ν e +ν µ + ν µ , producing three neutrinos with O(10 MeV) energies [73]. The analysis in Ref. [73] limits the flux of neutrinos resulting from proton catalysis in the Sun with 90% certainty to I 90 < 166.6 cm −2 s −1 . This limits the proton decay rate in the Sun to be less than where b π + ∼ 0.5 is the branching fraction of protons into π + , a π + ∼ 0.2 is the probability that a pion gets absorbed by the solar core [73], and d is the Earth-Sun distance. An EMBH passing through the Sun will lose energy as described in Eq.(3.1). Since the average velocity is v 10 −3 and the escape velocity from the solar surface is v esc 2×10 −3 , the EMBH needs to lose much of its kinetic energy to become gravitationally bound to the Sun. Taking an average internal temperature for the Sun of T = 2 × 10 6 K in Eq.(3.1), one finds this requires M 0.5 g. The Sun accumulates such EMBHs at the (slightly) Sommerfield enhanced rate of where r 7 × 10 5 km is the solar radius. Heavier EMBHs will get trapped in the convection zone, which makes up the outer third of the Sun. The temperature varies from 2 × 10 6 K at the base of the convective zone to 5700 K at the surface. The convective motions have velocities around 10 −6 [75]. Using the solar mass profile M (x) [76], we add the term GM M (x)/(r − x) 2 to Eq.(3.1) to account for the change in gravitational potential and numerically integrate the energy loss as an EMBH traverses the Sun. The velocity drops below the convection value if M 300 g. Thus, this range of masses will be trapped in the convection region.
One can check that annihilation within the convection region is unimportant. EMBHs trapped in the convection region will move with stellar convection currents and annihilate at a rate Γ ∼ N n σv th , where N and n are the number and number density of EMBHs accumulated in the convection region, v th is the thermal velocity, and σ is a cross section defined by the region where two EMBHs approaching each other at free-fall velocity can meet within the time a convection flow would cross the convection region ∼ 7×10 5 s (a rough 'coherence' time for the flow). Γ is never large enough to significantly alter the number of EMBHs accumulated in the Sun, so we neglect its contribution to these constraints.
The proton catalysis rate is where t = 4.6 Gyr is the age of the Sun, n p is the proton number density set by the local density, 0.2 g cm −3 , and v p ∼ 10 −4 is the local proton thermal velocity. For illustrative purposes, we set σ cat = σ QCD . Comparing this rate with the current bound (4.23) gives our constraint for M > 300 g. EMBHs with M < 300 g pass through the convection region and sink to the core. The vast majority that reach the core will ultimately annihilate. Only those that have not yet done so can contribute to proton catalysis. The total number in the Sun at any time, N tot , is made up of those that have not yet annihilated because they are either in the process of sinking to the core (N sink ) or part of a charge imbalance due to the random walk of charge accumulation (N ran ). There could additionally be a delay in annihilation if there is a strong enough magnetic field in the core, but we will ignore this possibility to set a conservative bound. Thus, N tot = N sink + N rand .
We estimate the number of EMBHs in transit as: where the sinking time t sink can be computed numerically. Using estimated temperature and density profiles in the Sun [76], we find for M > 0.5 g, implying N sink ∼ 10 9 f . On the other hand, the number accumulated from random asymmetries in the signs of the EMBHs passing through the Sun over its lifetime t is . (4.28) The EMBHs actively sinking toward the core and those in the innermost region of the star are in different environments and have different catalysis rates. We estimate the outer core where EMBHs are sinking to have T out ∼ 7 × 10 6 K and ρ out ∼ 50 g cm −3 , and the inner core where the random walk EMBHs reside to have T in ∼ 16 × 10 6 K and ρ in ∼ 150 g cm −3 . Using this information and Eq.(4.25), one can find the catalysis rates for both sinking and random walk EMBHs. Combining these gives the total proton decay rate, which must be less than F max p given by Super K. The resulting bounds are shown in Fig.4. One might worry that the accumulated random walk magnetic charge could be large enough to alter the dynamics in the stellar core. The maximum possible random walk charge for EMBHs that accumulate here would produce a magnetic field ∼ 10 −6 G in the lower convection region. This is negligible compared to the magnetic field already present, which could be as strong as 100 G [77], and is thus unlikely to influence stellar dynamics.

Direct Constraints
The magnetic monopole detector MACRO, described in Section 3.5, is sensitive to nucleon decay effects. Nuclear decay catalyzed by a traversing monopole produces extra radiation which changes the detector's sensitivity. MACRO produced separate analyses and flux bounds for monopoles that do [78] and do not catalyze nucleon decay [3]. In the catalysis case, the flux limit from MACRO data depends on the catalysis cross section. The flux limit for a magnetic charge moving at v ∼ 10 −3 , F (σ cat ), is presented for magnetic charges with catalysis cross sections up to σ cat = 10 −24 cm 2 in Ref. [78]. Increasing the catalysis cross section weakens MACRO's ability to verify EMBH events [78]. We take this cross section and the corresponding flux limit F 10 −15 cm −2 s −1 sr −1 , as it is the most conservative limit considered in Ref. [78]. This gives the below constraint, which is presented in Fig.4.

Sub-Extremal Magnetic Black Holes
This work so far has focused on EMBHs in a mass range, (M < 10 15 g), in which MBHs must be extremal to stable over the lifetime of the Universe. If the Q enhancement to the radiation rate of leptons and fermions described in Ref. [18] is present, then sub-extremal MBHs only survive to the present if their surface temperature is too low to effectively radiate electrons and positrons (M 10 17 g) [18]. The sections below explore the behavior of and abundance constraints on MBHs heavy enough to be stable over the lifetime of the Universe even when sub-extremal (M > 10 17 g). Some of the previously discussed constraints on EMBHs continue to apply at these high masses and also apply to sub-extremal MBHs. Heavy MBHs are constrained by their interactions with warm ionized gas clouds in the Milky Way and with WDs. Some can cause NSs to collapse into black holes. For convenience, we describe the charge of an MBH using q, a fractional value of the charge compared to the extremal charge, Q ext , of a black hole with mass M .

Gas Heating
Generalizing the constraints on EMBH abundance due to overheating warm gas clouds (see Section 3.1) for MBHs gives: This bound cuts off when MBHs either pass through WIM clouds too rarely to alter their behavior or lose a significant fraction of their kinetic energy to friction with the interstellar medium. Based on Eq.(3.1), an MBH with M 3 × 10 21 g may deposit the total energy needed to overheat the entire cloud in a single passage. The energy is initially deposited in a small region defined by the plasma attenuation length l ∼ 16 km. Depositing the energy needed to heat a ∼ 2 pc cloud into a cylindrical cross section with a radius of ∼16 km will heat the region to extreme temperatures, causing it to expand supersonically. This would significantly alter the cloud's structure and behavior within its sound crossing time, ∼ 5 × 10 5 yr for typical warm clouds [25]. Natural processes such as cooling and supernova shocks evolve the WIM over timescales of 10 6 yr [25]. MBHs must not disrupt WIM clouds more quickly than naturally observed processes. We limit the abundance of MBHs with M > 3 × 10 21 g, so that typical WIM clouds encounter them less than once per 10 5 yr.
This bound also cuts off when MBHs lose an O(1) fraction of their kinetic energy to friction with the WIM within ∼ 10 Gyr. Losing this much kinetic energy causes them to sink toward the Galactic center. There may be interesting phenomenology around this, but we have not found new constraints that arise from it. Modifying Eq. for MBHs passing through warm ionized regions. They lose energy most efficiently in the WIM, which fills 1% of the Galactic disk [79]. Most dark matter in the Galaxy, ∼ 90%, sits outside of the Galactic disk in the dark matter halo. We assume then that MBHs moving along virial orbits around the Galactic center spend ∼ 10% of their time in the Galactic disk. Using all of this, we estimate that MBHs lose an O(1) fraction of their kinetic energy within 10 Gyr if M 10 19 g q 2 1.6.
The gas heating bound cuts off where this is saturated. The bounds for MBHs with different parameters are shown in Fig.6.

White Dwarf Consumption
MBHs with M > 10 17 g can get trapped inside of a 1.2M WD and consume it within 1 Gyr. Their abundance is constrained by the existence of hundreds of WDs older than a Gyr and with masses ∼ 1.2M [30]. These WDs yield the strongest bounds because they are both abundant and dense enough to capture a wide range of MBHs. Uncharged black holes simply pass through a WD, while MBHs are trapped by friction with the Fermi gas. MBHs experience the same frictional forces in a WD as EMBHs (described in Eqs.(3.6) and (3.1)) but reduced by q 2 . Using this and the 1.2 M WD density profiles derived in [80], we find that all black holes that satisfy will get trapped in any 1.2M WD they pass through.
Estimates in Ref. [28] indicate that an uncharged black hole should be able to consume an entire WD within where c s ∼ 3 × 10 −2 is the speed of sound in a 1.2 M WD, λ is an O(1) constant presented in Ref. [28], ρ c ∼ 10 8 g cm −3 is the density at the core of a 1.2 M WD, and m W D is the mass of said WD [28]. One might worry that interactions with the large magnetic field around MBHs slows accretion or that MBHs cannot absorb electrically charged particles without becoming super-extremal. While the magnetic field around an MBH can repel charged particles before they reach the event horizon [81], the outward flow of repelled material would be resisted by pressure from the WD. The pressure from ρ c ∼ 10 8 g cm −3 of material moving away from the MBH at c s is approximately four times greater than the ambient pressure in the WD core. Balancing forces, ρ wd is reduced by at most a factor of ∼4 in the area around the MBH, increasing the accretion time by up to a factor of ∼4. All MBHs with M > 10 17 g will still be able to consume an entire 1.2 M WD within a Gyr. EMBHs in this mass range do not become super-extremal upon absorbing a single electrically charged particle because the magnetic and electric charges add in quadrature.

Sub-Extremal Black Holes Unstable
Uncharged Black Hole Bounds for all MBHs that satisfy Eq.(5.5). This constraint is shown in Fig.6

Neutron Star Collapse
The recent LIGO/Virgo observation of a 2.6 M compact object occupying the mass-gap between the theoretical maximum mass for NSs and minimum mass for stellar black holes [83] has raised questions about how such an object could form. Ref. [84] suggests that a low mass black hole, perhaps the object observed by LIGO/Virgo, could result from a NS absorbing non-annihilating dark matter and collapsing into a black hole. Some MBHs can get trapped in NSs due to friction with the Fermi gas and cause this collapse.
The minimum accretion rate for a black hole in a degenerate neutron gas with a stiff equation of state (the conditions of a NS core) was found [85] to be: An MBH can therefore consume a ∼ 2M NS within (5.9) Those with masses above ∼ 10 22 g can easily consume a NS within 1 Gyr, causing it to collapse into a mass-gap black hole.
Only MBHs that get caught in NSs can consume them. Using Eq.

Generalized Mass Functions
Throughout this paper we have assumed that EMBHs have a monochromatic mass function. It is worth exploring how the bounds described in the above sections transform for EMBHs with non-monochromatic mass functions, as these commonly arise from models of primordial black hole formation. We outline how one would apply the constraints described in this paper to EMBHs with a general mass function below. We describe the mass function using dn dM (M ), which gives the differential number density of EMBHs as a function of mass.
The constraint from WIM cloud heating in Section 3.1 is estimated by comparing the total heat generated by many EMBHs passing through WIM clouds in the Milky Way to the observed heating and cooling rates of said clouds. The heating rate for EMBHs with a general mass function can be found by multiplying the heating rate per EMBH, Eq.(3.2), and dn dM , and then integrating over all relevant masses. This can be compared against the observed cooling rate to obtain bounds.
The constraints from EMBH annihilation in WDs leading to runaway fusion and supernovae is described in Section 3.2. The total mass of EMBHs accumulated in a WD depends on the local density of dark matter and not the mass of the EMBHs themselves. Additionally, a set accumulated mass of EMBHs, M T h , is needed to overcome the core magnetic field and start annihilation (see Table 1). The annihilation bound is sensitive to the whether EMBHs have enough mass to provide the energy needed to initiate runaway fusion, and whether they are light enough to sink to the core within 1 Gyr. The bounds from WD J0551+4135 can be straightforwardly adapted for general EMBH mass functions because all EMBHs lighter than 4 × 10 12 g can sink to the core and trigger runaway fusion within 1 Gyr. WD J0551+4135 must not accumulate a group of at least six EMBHs, each lighter than 4 × 10 12 g, whose combined mass exceeds M T h within 1 Gyr. Not all EMBHs supply enough energy when they merge to trigger runaway fusion in 1 M WDs. This complicates estimating the bounds because the non-extremal black hole produced in an EMBH merger will only radiate off the mass of the smaller EMBH. One can simplify this estimate by considering two different sets of EMBHs. The first is those with 4 × 10 12 > M M T h . These will trigger runaway fusion once the WD accumulates 6 EMBHs whose combined mass exceeds M T h . The second set is EMBHs with M M T h . If we make the simplifying assumptions that all EMBHs in this group pair and annihilate exactly once, that annihilations are fast compared to the time needed to accumulate M T h of EMBHs, and that the distribution of EMBHs in the WD matches the true distribution, then one could estimate the number of times the WD would need to accumulate M T h of EMBHs before it would be expected that one of the annihilating pairs contains two EMBHs heavy enough to trigger runaway fusion. The faster of the two processes can be used to set the bound.
Generalizing the constraints on EMBH abundance due to gamma-rays produced when EMBH binaries annihilate (Section 3.3) requires modifying the derivation for the EMBH merger rate presented in Appendix B.2. The details of the original and modified derivation are described there. Broadly, the merger rate comes considering a collection of cold EMBHs whose separations in the early Universe are described by a Poisson distribution. These initial separations dictate which EMBHs form binaries, the semi-major axes and eccentricities of these binaries, and the time they take to merge. The merger rate at any time t comes from integrating over all binaries whose semi-major axis and eccentricity lead them to merge after t has passed. When the mass spectrum is not monochromatic, one must include terms in the distribution of binary parameters that account for the distribution of EMBH masses. The merger rate then comes from integrating this over the space of EMBH masses, binary semi-major axes, and eccentricities which produce binaries that merge after t has passed. Each merger remnant will radiate off the mass of the lighter EMBH in the binary before becoming extremal again. This amount of radiation per merger, along with the modified merger rate described in Appendix B.2 can be used to estimate the total diffuse gamma-ray signal produced by binary mergers, which can be compared to the observed diffuse extra galactic gamma-ray signal observed by FERMI-LAT.
The existence of large magnetic fields in both the Milky Way and in galaxy clusters was used to limit EMBH abundance in Section 3.4. Too many EMBHs would damage the Galactic magnetic field by removing energy from it and would screen the large coherent magnetic fields observed in galaxy clusters. Both effects are mass independent, and limit the total fraction of dark matter that can be made of EMBHs of any mass.
The MACRO monopole detector set upper limits on the flux of magnetic charges near the Earth (see Section 3.5). One can use this to constrain EMBHs with arbitrary mass distributions by comparing the expected flux of EMBHs, given by 1 4π v dn dM dM , where v ∼ 10 −3 is the local EMBH velocity, to the MACRO flux. The way EMBH catalysis of nucleon decay affects WD annihilation bounds is explored in Section 4.1. This is somewhat more complicated for EMBHs with non-monochromatic distributions, as the catalysis cross section typically varies with EMBH mass. The heat produced by catalysis can cause convection in WD cores which prevents EMBH annihilation. If the annihilation bounds break down, the EMBH abundance is constrained by comparing the predicted luminosity from EMBHs catalysing nucleon decay in WDs to that observed in old dim WDs. One can determine whether the annihilation bounds stand by comparing L th , the luminosity from catalysis decays needed to start convection in WDs, to L M , the maximum plausible luminosity from a collection of EMBHs in the WD. When the typical EMBH masses are small compared to M T h , L M can be determined by integrat- Here n is the number density of EMBHs, σ M is the maximum plausible catalysis cross section, defined in Eq.(4.3), ρ c is the density at the WD core, and v th is the thermal velocity of nuclei there. As in the monochromatic case, the annihilation bounds remain in tact unless L th < L M at all WD cooling temperatures. When annihilation bounds break down, one can estimate the average cross section per unit of EMBH mass using L th v th M T h ρc . This can be used to estimate the luminosity EMBHs would produce in an old dim WD and compared to the observed luminosities to get constraints. If the typical EMBH mass is large compared to M T h , then only 6 are needed to cause the WD to explode. In this case, one can estimate L M as the expectation value for the luminosity of 6 randomly selected EMBHs radiating at σ M . If the annihilation bound breaks down, one can find the heating bound by using the same general method described for light EMBHs, though using L th /6v th ρ c as the average catalysis cross section. Section 4.3 describes the EMBH abundance constraints from nucleon catalysis in many nearby NSs producing too much x-ray radiation. The catalysis cross section we used to demonstrate this constraint, σ QCD , is proportional to the EMBH mass. With this cross section, the nucleon decay rate and the resulting bound relies on the total mass of EMBHs accumulated in NSs and not generally on the mass of the specific EMBHs involved. The bound becomes mass dependent once EMBHs are heavy enough (M > 10 7 g) that one or a few of them cause NSs to produce too much radiation. At higher masses (M > 10 9 g) individual EMBHs captured in a fraction of all NSs can saturate the limits on xray radiation. One can assume the EMBH mass distribution across all NS reflects the true distribution. In this case, one can estimate the fraction of NSs which would acquire a given luminosity based on the probability that they accumulate a certain number of EMBHs of a given mass. Integrating this over all masses would give the expected combined luminosity from all NSs, which could then be compared to observed x-ray measurements to get constraints. Section 4.4 describes the constraints that emerge from EMBH induced proton decay in the Sun producing too many energetic neutrinos. As in Section 4.3, we demonstrate this bound using σ QCD as the catalysis cross section, which makes it dependent on the total mass of EMBHs trapped in the Sun, and not the mass of said EMBHs. Those lighter than 300 g annihilate too quickly to produce observable decays, while those heavier than 300 g get trapped in the convection zone where annihilation is strongly suppressed. Generalized bounds come from limiting the fraction of EMBHs heavier than 300 g so that the proton decay rate stays below observed limits.

Seeding Magnetic Fields
One possibly interesting area we did not address is the connection between EMBHs and primordial magnetic fields. It has been suggested that primordial magnetic fields present just before recombination could reduce the Hubble tension by altering the way charged particles cluster [86]. Additionally, many models of galactic magnetic field formation rely on the presence of a small seed field around the time of galaxy collapse [87]. The seed fields needed to generate the µG magnetic fields observed in galaxies today vary with the particular model, but reach as low as 10 −30 G fields coherent over distances of 100 pc [88]. EMBHs and MBHs carry high magnetic charges and can generate reasonably strong magnetic fields coherent over the average separation between them. Taking the typical EMBH separation, one can estimate the scale of the typical magnetic fields they produce: (1 + z), (6.1) where ∆ is the local EMBH overdensity. Such a field would be coherent over length scales There are some unconstrained regions of parameter space where MBHs could play a role in galactic magnetic field formation. This idea is explored further in Ref. [10].

Conclusion
In this work, we have extensively explored the phenomenology of EMBHs and stable MBHs.
We have placed stringent bounds on their abundances based on their ability to destroy WDs either by initiating supernovae (EMBHs with M < 4×10 12 g) or by consuming them (MBHs with M> 10 17 g). We also derived bounds which apply to both extremal and non-extremal MBHs based on overheating WIM clouds in the Milky Way. We described how EMBHs merge and the observable emission that results from these mergers. We provided minimal model-independent constraints for EMBHs that catalyze nucleon decay, while leaving a precise calculation of the catalysis cross section for future work. As in Refs. [15,16], we do not consider a detailed exploration of how EMBHs form. There is still much to understand about EMBHs. How exactly low mass EMBHs accrete, whether nucleon catalysis happens, and the exact details of such a process are not yet understood. While it is quite clear that EMBHs cannot constitute a meaningful fraction of the dark matter, they remain long-lived, well-motivated, very simple extensions of the Standard Model. They offer new mechanisms for heating Galactic plasma, producing gamma-rays, initiating supernovae, and generating NS mass black holes. Their behavior also gives us insight into more exotic objects, such as higher dimensional black holes, which carry magnetic charge [12]. EMBHs remain interesting tools to consider for many current and future astrophysical questions.
Acknowledgements We would like to thank Emanuele Berti, Rob Farmer, Michael Fedderke, Ely Kovetz, Michael Shara, and Raman Sundrum for their insightful conversations and guidance. Special thanks to Leandro Althaus for sharing their detailed analysis of WD interior compositions. This work was supported in part by the National Science Foundation under grant PHY-1818899.

A Magnetic Friction with a Fermi Gas at Finite Temperatures
The friction between a magnetic charge and an electron Fermi gas at finite temperature can be found by modifying the derivation presented in Ref. [32]. The energy of electrons that scatter off a passing magnetic charge will change by ∆E = pv(cosθ − cosθ ), where θ is the angle between the velocity of the incoming electron and the magnetic charge velocity v, θ is the angle between the scattered electron's velocity and v, and p is the momentum of the electron. The energy given to a magnetic charge from all of the scattering electrons can be found by integrating: Here is the number density of electrons with energy E in a Fermi gas at temperature T and with a Fermi energy E f .
represents the fraction of occupied states at energy E. The scattering cross section dσ dΩ is analogous to a Rutherford scattering cross section [32] dσ dΩ = Q 2 e 2 4p 2 sin 4 (Ψ/2) , where cosΨ = sinθsinθ cos∆φ + cosθcosθ . Ψ is the angle between the velocities of the incoming and scattered electron. The integral over Ψ diverges as a log as Ψ approaches zero. We cut it off at Ψ min = ν ie σ 2E f , which limits the scattering angle by the scattering length in the medium [32]. Integrating Eq.(A.1) over Ψ and φ gives 16π log 1 Ψ . Substituting x = cos θ and y = cos θ gives: where G is the gravitational constant, e is the electron charge, m e is the elctron mass, and ν ei σ ∼ 10 17 s −1 is the electron ion scattering frequency [33].

B.1 Derivation of the Merger Time for an EMBH Binary
The time it takes an EMBH binary to merge can be found by considering the energy lost to EM radiation emitted during orbit. EMBHs moving through an orbit radiate as two perpendicular oscillating magnetic dipoles. Energy losses via gravitational radiation come from the binary's quadropole moment, and are subdominant. The binary emission rate for black holes of arbitrary electric charge and mass in circular orbits is presented in [89,90]. The emission rate for black holes with arbitrary electric and magnetic charge can be found in [91]. The merger rate for charged black holes in a circular orbit is also considered in [16]. The electric and magnetic fields of two perpendicular magnetic dipoles which vary arbitrarily with time can be described in spherical coordinates by orienting them in the θ = π/2 plane.
where m 1 (τ ) and m 2 (τ ) are the dipole moments as a function of the retarded time τ = t−r, and r is the radius from the center of the dipole. The energy loss is found by taking dE dt = 1 4π (E × B) · dA at r = ∞. Changes to the shape of the orbit are described by changes to the angular momentum along the θ = 0 axis, which can be described by dL dt = 1 4π cos θ (r × (E × B)) dA integrated over a sphere at r = ∞. Using the terms from Eq.(B.1), this gives: To use these expressions to get a merger time, one must describe the dipole moments in terms of orbital parameters of the binary system. Each EMBH moves about the center of mass of the system with an angular velocity: where a is the semi-major axis of the orbit, e is the eccentricity, and M 1 and M 2 are the masses of each EMBH. EMBHs feel an equal mutual magnetic and gravitational attractive force, which is accounted for by the extra factor of 2 in the above expression. The distance of an EMBH from the center of mass of the system, d, is set by: In terms of these parameters the dipole moments are: where m p is the Plank mass and m −2 P = G. Using these expressions, all of the time derivatives of the dipole moments can be found in terms of a, e, and φ. Putting these into Eq.(B.2) gives: Averaging the energy and angular momentum losses over a single orbit gives: For comparison to the energy losses from gravitational radiation, we rewrite these assuming  [92], very small compared to the power from EM radiation. The extra factor of 16 comes from enhancements to the energy and angular momentum of the binary due to the magnetic attraction between EMBHs. Power loss from gravitational radiation will only match EM radiation once a drops down to 12 5 r s , where r s = 2GM is the Schwarzschild radius, and is thus unimportant for estimating the time needed for EMBH binaries to merge.
The energy of the system is related to the semi-major axis by E = − GM 2 a . This is twice the energy of a neutral black hole system of the same mass due to the extra attractive force between the EMBHs. The average rate of energy loss can be related to the average change in the semi-major axis by dE dt = GM 2 a 2 da dt , which gives; For e = 0 one can solve for the merger time analytically where a 0 is the initial semi-major axis. For comparison, the merger time if only the gravitational radiation were considered would be (B.11) When e = 0 the merger time also depends on how the eccentricity evolves. Using L 2 = GM 3 a(1 − e 2 ) gives the relation which can be rearranged to find The EMBH merger rate can be found by modifying derivations for the merger rate of uncharged black hole binaries [93,94] to account for the extremal magnetic fields. The initial orbital parameters of binaries, their semi-major axes, a 0 , and eccentricities, e 0 , are characterized by their co-moving separation, x, the co-moving distance to the next nearest EMBH, y, and the redshift, z d , when the binary decouples from the Hubble flow. Co-moving distances will correspond to distances in the present, when the scale factor equals 1, if the separated objects never decouple from the Hubble flow. When friction with the thermal bath is not important, EMBH binaries decouple from the Hubble flow and form a stable binary when the time needed for them to free-fall into each other becomes less than the Hubble time [93,94]. During radiation domination this corresponds to a decoupling redshift 4πρ dm , ρ dm is the dark matter density, and z eq is the redshift of matter radiation equality.
The infall dynamics change when friction with the thermal bath becomes large compared to Hubble acceleration. In this scenario, the binary decouples when the EMBH terminal infall velocity exceeds the Hubble velocity v H = H(z)x (1+z) , where H(z) is the Hubble expansion rate at a redshift z. The terminal infall velocity is found by balancing the attractive force between the EMBHs, 2GM 2 (1+z) 2 x 2 , with the frictional force exerted by the plasma on them. The frictional force between a magnetic monopole and a plasma is described in Ref. [24]. We adapt it here for an EMBH with mass M : where n e , is the electron number density in the thermal bath which we take to be ∼ (1 + z) 3 ρ c0 Ω b mpr , ρ c0 is the critical density today, Ω b is the baryon fraction, m pr is the proton mass, e and m e are the electron charge and mass, respectively, T is the gas temperature which we take to be ∼ 2.7(1 + z)K, V is the velocity of the EMBH, λ D = T 4πnee 2 is the plasma Debye length, l = 2T πme 1/4 1 V 1/2 ωp is the attenuation length of the plasma and ω p is the plasma frequency. To simplify comparing forces we will represent ( dE dx ) f as M 2 V (1 + z) 5/2 F. F represents the mostly z, M and V independent component of ( dE dx ) f If Fz 5/2 M > H(z) at the regular decoupling redshift, z d , then the EMBHs will continue moving apart with the Hubble flow until 2G x 2 (1+z) 1/2 F = H(z)x (1+z) . During radiation domination this gives a decoupling redshift , Ω r is the fraction of the critical energy density composed of radiation at z = 0, and H 0 is the Hubble constant. Binaries must decouple by matter radiation equality because the ratio between the free-fall time and the Hubble time becomes constant during matter domination.
The semi-major axis of the binary is a 0 = x 2(1+z d ) . The eccentricity, e 0 , of the binary is set by interactions with the next nearest EMBH. It pulls the opposite charge EMBH in the binary out of it's free-fall path, without affecting the same charge one. The displacement from its free-fall path sets the semi-minor axis for the binary, b. This can be estimated by The probability that two oppositely charged EMBHs are separated by a co-moving distance between x and x + dx and that the next nearest one is between y and y + dy can be described by a Poisson distribution. dp = 1 2 (4πn) 2 x 2 y 2 dxdy, (B. 24) where n is the co-moving number density of EMBHs. The factor of 1 2 arises because each EMBH can only pair with an oppositely charged partner. x is by definition less than y. y can be at most the average separation between EMBHs. This defines y max ≡ 3M 4πρ dm f All binaries that merge at t m have parameters that lie along the curve in e − a space a 0 (t m , e 0 ). With a bit of algebra, one finds when e ∼ 1 the probability that an EMBH binary merges between t m and t m + dt is dp = (4πn) 2  The merger rate per EMBH, Γ m , at a specific t m can be found by integrating dp dt over all possible starting eccentricities for a binary merging at t m . f is the fraction of dark matter composed of EMBHs. Multiplying Γ m by n bh , the proper number density of EMBHs at t m , would give the rate at which EMBH mergers happen at t m . The bounds used in these integrals e max and e min are the largest and smallest eccentricities a binary merging at t m could have. We define e upper as the eccentricity a binary with semi-major axis a 0 and y = y max , the highest eccentricity possible for that binary given the average separation between EMBHs. e max is set by where the a 0 (e 0 , t m ) curve intersects e upper or a eq , the semi-major axis of a binary that decouples at the latest allowed decoupling redshift z eq . Finally, e min is the lowest eccentricity a binary that merges in t m and decouples after z f can have. z f denotes the redshift at which frictional and Hubble acceleration are equal: (1 + z f ) = Here e min is determined in the limit that e min ∼ 1 because the merger rate is dominated by binaries with e ∼ 1.
If one wanted to generalize the merger rate to EMBHs with a non-monochromatic mass specturm then one could follow the same general procedure outlined above with some modified values of dp, z d , and e. Instead of simply integrating over the range of eccentricities that will lead a binary to merge in the present, one must integrate over all possible eccentricities and EMBH masses for a given mass distribution. One can approach this by rewriting Eq.(B.24) as , and the range of possible masses. This will likely need to be done numerically for arbitray mass distributions.