Black holes and gravitational waves in models of minicharged dark matter

In viable models of minicharged dark matter, astrophysical black holes might be charged under a hidden U(1) symmetry and are formally described by the same Kerr-Newman solution of Einstein-Maxwell theory. These objects are unique probes of minicharged dark matter and dark photons. We show that the recent gravitational-wave detection of a binary black-hole coalescence by aLIGO provides various observational bounds on the black hole's charge, regardless of its nature. The pre-merger inspiral phase can be used to constrain the dipolar emission of (ordinary and dark) photons, whereas the detection of the quasinormal modes set an upper limit on the final black hole's charge. By using a toy model of a point charge plunging into a Reissner-Nordstrom black hole, we also show that in dynamical processes the (hidden) electromagnetic quasinormal modes of the final object are excited to considerable amplitude in the gravitational-wave spectrum only when the black hole is nearly extremal. The coalescence produces a burst of low-frequency dark photons which might provide a possible electromagnetic counterpart to black-hole mergers in these scenarios.


Introduction
Astrophysical black holes (BHs) are considered to be electrically neutral due to quantum discharge effects [1], electron-positron pair production [2][3][4], and charge neutralization by astrophysical plasmas. These arguments rely -one way or the other -on the huge chargeto-mass ratio of the electron, 1 e/m e ≈ 10 21 . Together with the celebrated BH no-hair theorems (cf. ref. [5] for a review), these arguments imply that -within Einstein-Maxwell theory -vacuum astrophysical BHs are described by a special case of the Kerr-Newman metric [6], namely the Kerr solution [7]. The latter is characterized only by its mass M and angular momentum J := χM 2 , since the electric BH charge Q em is assumed to be negligible in astrophysical scenarios.
On the other hand, models of minicharged dark matter (DM) predict the existence of new fermions which possess a fractional electric charge or are charged under a hidden U(1) symmetry [8][9][10][11][12][13]. Their corresponding charge is naturally much smaller than the electron charge and their coupling to the Maxwell sector is suppressed. These minicharged particles are a viable candidate for cold DM and their properties have been constrained by several cosmological observations and direct-detection experiments [12,[14][15][16][17][18][19][20][21][22]. In some other models JCAP05(2016)054  Figure 1. The parameter space of a minicharged DM fermion with mass m and charge q = e 2 h + 2 (see main text for details). The left and right panels respectively show the planes h = 0 and = 0 of the three dimensional parameter space (m, h , ). As a reference, an electron-like particle (m ∼ 0.5 MeV, = 1) is denoted by a black marker. In each panel, the red and blue areas below the two threshold lines denote the regions where charged BHs with a charge-to-mass ratio Q/M > 10 −3 and Q/M = 1 can exist [cf. eq. (2. 23)]. The region above the black dashed line is excluded because in this region extremal BHs would discharge by plasma accretion within less than the Hubble time [cf. eq. (2. 20)]. Left panel: the hatched region is excluded by the effects of the magnetic fields of galaxy clusters [21] and it is the most stringent observational constraint on the model (we also show the region excluded by the direct-detection experiment LUX [24], cf. ref. [21] for details and other constraints). Right panel: when = 0 our model reduces to that of DM with dark radiation [23] and the region above the solid black line is excluded by soft-scattering effects on the galaxy dynamics [23]. In the region above the dark red dot-dashed line hidden photons emitted during the ringdown of a M ∼ 60M would be absorbed by hidden plasma of density ρ DM ∼ 0.4 GeV/cm 3 [cf. eq. (3. 42)]. dark fermions do not possess (fractional) electric charge but interact among each other only through the exchange of dark photons, the latter being the mediators of a long-range gauge interaction with no coupling to Standard-Model particles [23].
Although rather different, both these DM models introduce fermions with a small (electric and/or dark) charge. It is therefore natural to expect that minicharged DM can circumvent the stringent constraints that are in place for charged BHs in Einstein-Maxwell theory. In this work we will show that this is indeed the case and that even extremal Kerr-Newman BHs are astrophysically allowed in the presence of minicharged DM.
Charged BHs are remarkably sensitive to the presence of even tiny hidden charges but are otherwise insensitive to the details of their interaction. This is a consequence of the equivalence principle of general relativity. We shall take advantage of this universality and discuss BHs charged under a fractional electric charge or under a hidden dark interaction on the same footing. Figure 1 summarizes the main results of section 2, showing the parameter space of a minicharged fermion with mass m and charge q = e 2 h + 2 in which astrophysical charged BHs can exist (here and in the following h and are the fractional hidden and electric charges of the dark fermion, respectively). Interestingly, such region does not overlap with the region excluded by direct-detection experiments and by cosmological observations.
Having established that charged astrophysical BHs can exist in theories of minicharged -2 -

JCAP05(2016)054
DM, we proceed to study their gravitational-wave (GW) signatures in section 3. We consider the coalescence of a binary BH system similar to GW150914, the event recently detected by the LIGO/Virgo Collaboration using the GW interferometers aLIGO [25], and show that different phases of the coalescence can be used to constrain the U(1) charge of the BHs in the binary and of the final BH produced in the post-merger phase. Finally, we explore the excitation of the ringdown modes and the total radiated energy during the collision of unequal-mass BHs, in a perturbative approach, and show that they are in good agreement with previous, restricted results in Numerical Relativity, as well as with simplified flat-space calculations.
2 Charged BHs in minicharged DM models

Setup
We consider the following classical Lagrangian [10] where F µν := ∂ µ A ν − ∂ ν A µ and B µν := ∂ µ B ν − ∂ ν B µ are the field strengths of the ordinary photon and of the dark photon, respectively, j µ em and j µ h are the electromagnetic (EM) and the hidden number currents, e is the electron charge, and e h is the gauge coupling of the hidden sector. The model (2.1) describes a theory in which a charged fermion is coupled to ordinary photons with coupling 2 e 2 and to dark photons with coupling e 2 h := 2 h e 2 . The parameters and h are free. 2 The field equations arising from the Lagrangian (2.1) read

3)
G µν = 8π(T em µν + T h µν ) , (2.4) where we defined the effective stress-energy tensors for the standard Maxwell field and for the dark photon, respectively. The continuity equations ∇ ν j ν em = 0 = ∇ ν j ν h for the standard and hidden currents follow directly from eqs. (2.2) and (2.3). Note that a hidden electron carries both electric charge e and hidden charge e h .
In the absence of currents (j µ em = j µ h = 0), the most general stationary solution [5] of the above field equations is the Kerr-Newman [6] metric with total charge Q = Q 2 em + Q 2 h , 2 In the model studied in refs. [10,18], / h = tan θ defines the kinetic mixing angle θ of the photon fields before the diagonalization leading to the Lagrangian (2.1) in which Fµν and Bµν are decoupled (cf. refs. [10,18] for details). In this case, the effective electron charge is e cos θ. This kinetic coupling is constrainted by arguments related to Big Bang nucleosynthesis, cosmic microwave background, and large-scale structure formation (cf., e.g., ref. [26]). For simplicity we neglect such coupling here and consider the Lagrangian (2.1) as fundamental. This corresponds to the model studied in refs. [10,18] when θ → 0. where the (EM and hidden) BH charges are defined from the solution of the (standard and hidden) Maxwell equations in static equilibrium, When the BH is spinning the presence of a charge induces a magnetic field along the angular directions [6]. Note that the hidden current j ν h provides both electric charge and hidden charge to the BH. Thus, a standard electron and a dark electron feel respectively the force If the BH acquires charge only through the accretion of a hidden current (i.e., if j ν em = 0, a condition that is enforced by several considerations as discussed in section 2.2), then Q = 0 and the EM and hidden charge are proportional to each other, Q em = h Q h . In this case the BH charge reads Q = Q h 1 + 2 / 2 h , and the force felt by a hidden electron reduces to where we have defined the effective charge of a dark fermion in this model, q := e 2 h + 2 . Observational bounds on minicharged DM typically constrain the coupling to Standard-Model particles, especially the coupling to ordinary photons. As such, they are insensitive to the coupling e h which is indeed typically neglected. 3 On the other hand, the dark coupling e h plays a crucial role in models of dark radiation [23]. Gravitational tests do not require ordinary photons as mediators and are indeed sensitive to the entire parameter space ( h , ). In particular, the effects we are going to discuss are present also when = 0, i.e. when DM does not couple to Standard-Model particles, as in dark radiation models [23].

Theoretical bounds on the charge-to-mass ratio of astrophysical BHs
There are several mechanisms that conspire to limit the electric charge of astrophysical BHs. One is purely kinematical. Take a BH with mass M and electric charge 4 Q em and throw in a low-energy electron of charge e and mass m e . For the electron to be absorbed by the BH, then the (classical) electric force must satisfy which can be written in terms of the dimensionless charge-to-mass ratio of the BH as (2.13) 3 When h → 0 the model (2.1) simply describes a dark fermion coupled to ordinary photons with coupling 2 e 2 e 2 . 4 We are now considering the standard scenario in which the electric charge is produced by an ordinary current j µ em . In other words, Q h = 0 and Qem = Q as defined in the previous section.

JCAP05(2016)054
These numbers change if the particle is thrown at large velocities, but show that the maximum charge-to-mass ratio Q em /M is typically very small. In addition, BHs can be neutralized by surrounding plasma. If the eletrical force overwhelms the gravitational force, charge separation can occur and the BH charge can be neutralized by particles of opposite charge.
For an extremal BH with Q em = M , the total number N of elementary charges (each with charge q and mass m) that it needs to accrete from the surrounding plasma to be neutralized is [27] N ∼ 10 39 e q M M . (2.14) This corresponds to a plasma mass of with m p being the proton mass. This plasma mass is easily available under the form of interstellar matter within a small region surrounding the BH [27]. To estimate the time needed to accrete an amount of plasma such that an extremal BH is discharged, let us assume that plasma accretion occurs at the same rate of gas accretion from an ordinary accretion disk. In the most conservative scenario mass accretion occurs at the Eddington rate,Ṁ Edd = 2.2 × 10 −8 (M/M )M yr −1 , corresponding to a discharge time scale For q = e and m = m e , τ discharge ∼ 8 × 10 −7 s. The above results show that -in Einstein-Maxwell theory -it is difficult to charge a BH past Q em /M ∼ 10 −21 (in geometric units), and that -if such BH ever acquires a charge -it discharges very quickly. Another possibility to form charged BHs is if the latter were born through the gravitational collapse of (charged) stars. Self-gravitating stars are globally charged due to pressure effects, as shown by Eddington, Rosseland and others [28][29][30]: in a nutshell, lighter charges (electrons) are easily kicked out of the star by pressure effects whereas heavy ions (protons) are stuck in the interior. The calculation, which proceeds by assuming thermal equilibrium for positive and negative charges, yields the following result for the charge of a star [28][29][30], where the stellar material is assumed to be composed mainly of two species of ions with charge and mass (m 1 , q 1 ) and (m 2 , q 2 ), respectively. For standard Maxwell theory, the chargeto-mass ratio in the star is of the order of eq. (2.13). Thus, stars are typically charged with ≈ 100 C. It is reasonable to expect that, if they collapse to a BH, this small charge will remain hidden behind the horizon, giving rise to a (very weakly) charged BH with Q em /M ≈ 10 −19 (M /M ). The above discussion imposes an extremely stringent upper bound on the EM charge of astrophysical BHs. To fulfill this bound, we assume that astrophysical BHs do not accrete ordinary charges and, therefore, we set the current j ν em to zero. Thus, in the rest of this work we can safely assume that the force felt by a putative hidden electron is given by eq. (2.11).
The above discussion also shows that the bounds on the charge-to-mass ratio of astrophysical BHs become much less stringent in minicharged DM models. For a hidden electron which can be less stringent than unity 5 depending on the parameters m and q. Likewise, a BH with hidden charge surrounded by a plasma of hidden electrons can be neutralized by accreting a mass M plasma given by eq. (2.15). If q e or m m p , M plasma might be a considerable fraction of the BH mass so that a charged BH would be difficult to discharge in this scenario. Furthermore, DM being almost collisionless it does not form accretion disks around compact objects so that dark plasma is accreted at the Bondi rate for collisionless fluids,Ṁ where ρ is the local DM density and v is the relative velocity between DM and the BH. Using the above equation and eq. (2.15), we estimate the discharge time which is much longer than eq. (2.16) since the DM density is low. In the above expression we have normalized ρ and v by their typical local values. Thus, the discharge process is much slower than for ordinary plasma. The black dashed line in figure 1 shows the threshold τ discharge = 1.4 × 10 10 yr obtained from eq. (2.20) in the (m, q) parameter space. In the region below this line, the time scale to discharge an extremal BH through accretion would be much longer than the Hubble time, which we assume as a very conservative limit. Furthermore, BHs with a large charge also have a large electric field close to the horizon. Such electric field is prone to produce spontaneous pair production via the Schwinger mechanism [31]. This effect becomes relevant when the work done by the electric field on a Compton wavelength is of the order of the rest mass of the lightest particle [1], In the rest of this work we shall focus on this regime, i.e. we do not consider ultralight DM whose Compton wavelength is larger than the gravitational radius of the BH.
To summarize, in minicharged DM models eq. (2.18) provides an intrinsic constraint on the maximum charge-to-mass ratio Q/M of astrophysical BHs. This constraint can be also turned around to determine the region in the (m, q) plane in which BHs with a certain Q/M can exist, namely q e < 2 × 10 −18 m GeV M Q . (2.23)

JCAP05(2016)054
This region is shown in figure 1 for two different values of Q/M . Interestingly, existing constraints on minicharged DM models do not rule out charged BHs. On the contrary, even extremal BHs with charge-to-mass ratio Q/M ∼ 1 are allowed, even when e h = 0. We stress that EM-based constraints on minicharged DM models are insensitive to the coupling e h and that charged BHs can exist also when = 0, a regime which cannot be ruled out by EM observations. Finally, the colored regions shown in figure 1 lay well below the black dashed threshold line for plasma discharge [cf. eq. (2.20)], so that charged BHs should not easily discharge in this scenario.

Gravitational-wave tests of charged BHs in minicharged DM models
The recent GW detection of a binary BH coalescence by aLIGO [25] has given us access to the strong-field/highly-dynamical regime of the gravitational interaction. In this regime, precision GW measurements can be used to develop BH-based tests of fundamental physics [32,33]. It is therefore natural to investigate whether present and upcoming GW observations can be used to constrain a putative hidden charge of astrophysical BHs in minicharged DM scenarios. The GW-driven coalescence of a compact binary can be characterized by three phases [34][35][36]: the inspiral, the merger and the ringdown. The inspiral corresponds to large orbital separations and is well approximated by post-Newtonian theory; the merger phase corresponds to the highly nonlinear stage right before and after merger, and can only be described accurately through numerical simulations of the full Einstein's equations; finally, the ringdown phase corresponds to the relaxation of the highly-deformed end-product to a stationary, equilibrium solution of the field equations, and can be described by BH perturbation theory [36][37][38]. In the rest of this paper we will focus on each of these phases, with various degrees of approximation, in order to estimate the effect of a hidden BH charge in the GW signal.

The inspiral of two charged BHs in minicharged DM models
We model the initial inspiral of the BH binary by considering two point charges, q 1 = e 2 h,1 + 2 1 and q 2 = e 2 h,2 + 2 2 , with masses m 1 ≤ q 1 / √ G and m 2 ≤ q 1 / √ G, respectively. 6 To leading order, by using eq. (2.11), the motion is governed by the equation where the upper (lower) sign refers to i = 1 (i = 2), r i is the position vector of the mass m i , r := r 2 − r 1 is the relative position of the bodies. If we define q i := λ i m i √ G, from the equation above it is clear that the problem can be mapped into a standard Keplerian motion of two uncharged particles where G → G eff := G(1 − λ 1 λ 2 ).
Let us assume circular orbits for simplicity. 7 An (electric or hidden) charge in circular motion emits dipolar radiation governed by Larmor's formula. The dipolar energy flux dominates over the quadrupolar GW flux at large distances, so it might play an important role in the early inspiral. Because both particles are charged under the standard Maxwell

JCAP05(2016)054
field and under the hidden field, the dipolar flux consists of two copies of Larmor's energy dissipation [40] with different couplings, in addition to the standard GW energy flux [41]. Namely, where R is the radius of the orbit, M = m 1 +m 2 and η = m 1 m 2 /M 2 . Note that if the chargeto-mass ratio of the two objects is the same, the corresponding dipole term is zero and (only in this case) the leading term would be quadrupolar like the GW flux. It is convenient to write the total dipolar flux, From the above equation it is clear that the dipolar emission in the inspiral phase depends only on the combination ζ. Note that the latter is nonvanishing even when i = 0, provided e h,1 or e h,2 are nonzero. On the other hand, ζ ∼ 0 if the two BHs have similar charge-to-mass ratios.
The inspiral phase of GW150914 was compatible with the prediction of post-Newtonian theory [25,42]. This implies that any putative dipolar contribution to the energy flux must be small. In the limit ζ 1, the condition can be written as where Ω is the orbital frequency and we assumed m 1 ∼ m 2 . This simple estimates suggest that the dipolar correction might be small when the binary enters the aLIGO band.
To quantify the effect of the dipolar energy loss more precisely, we can compute the GW phase associated with such effect through a simple quasi-Newtonian evolution [43]. The binding energy of the two-body system is E = −G eff m 1 m 2 /R. By assuming an adiabatic approximation,Ė = −dE GW /dt − dE dip /dt, we can obtain a differential equation for the orbital radius R = R(t). This equation can be solved analytically in the limit |dE dip /dt| |dE GW /dt|, i.e. when the dipolar loss is a small correction compared to the quadrupolar GW flux. In this case, a standard procedure [43] allows us to compute the amplitude and the phase of the quadrupolar GWs emitted by the system through an adiabatic evolution. The phase of the "+" GW polarization reads  [25]. Note that when λ 1 λ 2 < 0 the total mass of the system can be significantly smaller than in the uncharged case.
where t c and Φ c are the time and phase at coalescence, f is the GW frequency and M := M η −3/5 is the chirp mass. When i = e h,i = 0, then G eff = G, ζ = 0, and eq. (3.8) yields the standard leading-order result, to which we added the first next-to-leading order post-Newtonian term in the second line of the above equation.
However, when at least one of the parameters i , e h,i are nonzero, we obtain two types of corrections. The first one is a rescaling of Newton's constant, which affects also the Newtonian result. This correction is present even if ζ = 0, provided λ i = 0. Because to leading order the GW phase depends on G eff M, a rescaling of Newton's constant is degenerate with the measurement of the chirp mass. Extracting the latter from the Newtonian GW phase obtained by neglecting charge effects would yield a result that is rescaled by a factor G eff /G relative to the real chirp mass of the system, namely Figure 2 shows the total mass M as a function of m 1 for fixed M measured ≈ 30M [25] and for different values of λ i . When λ i = 0 we recover the standard result, namely a minimum total mass M min ≈ 69M . However, significant changes occur if 0.1 |λ i | < 1. In particular, when λ 1 λ 2 < 0 the effective Newton's constant is larger and the real total mass of the system can be significantly smaller than in the uncharged case. This property is intriguing since it shows that neglecting charge effects might systematically lead to overestimate the measured BH masses. The other correction appearing in eq. (3.8) is the second term inside the square brackets in eq. (3.8). The latter is larger at small frequencies, as expected, and in fact reassembles the leading-order correction for neutron-star binaries in scalar-tensor theories. 8 The ratio between the second term and the third term inside the square brackets in eq. (3.8), namely is shown in figure 3 as a function of the GW frequency for a typical inspiral. For f > 30 Hz, the charge-induced corrections are at least ∼ 0.03 times smaller than the first post-Newtonian term when ζ 0.1, but they can be as large as 25% when ζ ≈ 0.3. Another relevant quantity is the number of cycles spent in the detector bandwidth, N := fmax f min df fḟ . In our case we obtain and we have also expanded for f max f min to simplify the final expression. In the smallcharge limit, for f max ∼ 100 Hz and f min ∼ 30 Hz, we obtain Therefore, dipolar effects change the number of cycles relative to the Newtonian case by a few percent when ζ ≈ 0.1 and by less than 0.01% when ζ < 0.01. On the other hand, these corrections become important at smaller frequencies and might produce detectable effects for space-based interferometers such as eLISA [44].

JCAP05(2016)054
Very recently, ref. [45] performed a detailed analysis to derive GW-based constraints on generic dipolar emissions in compact-binary inspirals (see also ref. [46]). It is straightforward to map eq. (3.5) into this generic parametrization. In our case the parameter B defined in refs. [45,46] reads B = 5 24 ζ 2 . The analysis of ref. [45] shows that GW150914 sets the upper bound |B| 2 × 10 −2 , whereas a putative eLISA detection of a GW150914-like event with an optimal detector configuration or a combined eLISA-aLIGO detection set a projected bound as stringent as |B| 3 × 10 −9 . In our case these bounds translate into for the combination ζ defined in eq. (3.6). Note that the bound derived from the aLIGO detection is roughly consistent with our simplified analysis.
Finally, we stress that if the two BHs have similar charge-to-mass ratios (which might be the case if their formation mechanisms are similar), then ζ ≈ 0 and the dipolar emission is suppressed. In this case the zeroth-order corrections due to G eff in eq. (3.8) would still be present [cf. figure 2]. In addition, the first nonvanishing radiative effect would be a quadrupolar term that will also modify the standard Newtonian quadrupole formula.

Ringdown phase and bounds on the BH charge
Once the two BHs merge, they form a single deformed charged and spinning BH which will relax to its final stationary (Kerr-Newman [6]) state by emission of GWs, EM and dark radiation. The final, "ringdown" stage of this process is well-described by the superposition of exponentially damped sinusoids, called quasinormal modes (QNMs), which are the characteristic oscillation modes 9 of the final BH [37,[48][49][50]. Here, l and m are angular indices describing how radiation is distributed on the final BH's sky (|m| ≤ l), and n is an overtone index. Usually, the modes excited to larger amplitudes 10 are the (2, 2, 0) and (3, 3, 0) gravitational modes [34,35,51]. Because the final state only depends on three parameters, the knowledge of the (ω 220 , ω 330 , τ 220 ) triplet, for example, allows us to invert the problem and to determine the mass M , spin J := χM 2 and charge Q of the final BH.
To complete this program we must first know the QNMs of Kerr-Newman BHs, which has been an open problem for more than 50 years. Fortunately, this problem was recently solved in a series of papers [52][53][54][55][56] and the QNMs of a Kerr-Newman BH can be now computed numerically.
For simplicity, here we focus on the small-charge case. The l = m = 2 ringdown frequencies of Kerr-Newman BHs in the small-charge limit are well approximated by the 9 Very recently, ref. [47] showed that compact objects without an event horizon but with a light ring might display a ringdown signal similar to BHs even when their QNM spectrum is completely different. Such effect is due to the different boundary conditions that occur for a horizonless ultracompact object and does not play any role for charged BHs. For this reason in the rest of the paper we will refer to QNMs or to ringdown modes without distinction. 10 As discussed in section 3.3, when the BHs are highly charged the final ringdown might also depend on the l = m = 2 Maxwell modes which might be excited to considerable amplitude. We neglect this possibility here for simplicity.
As discussed in appendix A, there is a tight relation between the BH QNMs and the some geodesic properties associated to the spherical photon orbits, as established in the eikonal limit [59][60][61]. Indeed, it turns out that the geodesics correspondence provides an estimate for the relations (3.16)-(3.17), as well as for the correction for the (3, 3, 0) mode. Because numerical data for the l = m = 3 modes are not available, in the rest we will estimate δω 330 from the geodesic correspondence presented in appendix A.
Let us then assume that the two dominant modes of the gravitational waveform were extracted from a GW detection, so that the two frequencies ω KN 220 , ω KN 330 and the damping time τ KN 220 were measured. In principle, using the formulas above, the mass, spin and charge of the BH could be determined precisely. Unfortunately, detection is always done in the presence of noise, which introduces some uncertainty in the determination of the ringdown frequencies. The proper way to handle noise is by either using Monte-Carlo simulations and a Bayesian analysis or by approximating the process through a Fisher-matrix analysis [62]. A Fisher-matrix study of multi-mode ringdown is done in ref. [57], which we follow. For a single mode, the relevant entries are shown in the appendix B.
Consider now two modes, mode 1 (l = m = 2) with amplitude A 1 , frequency ω 1 = ω KN 220 and damping time τ 1 = τ KN 220 and mode 2 (l = m = 3) with amplitude A 2 , frequency ω 2 = ω KN 330 and damping time τ 2 = τ KN 330 . Define also the quality factor for i = 1, 2. For detection of multi-modes having "orthogonal" angular structure (i.e., the two modes are characterized by different l, m indices), the Fisher matrix is simply an addition of matrices of different modes. In this case the errors σ ω i and σ τ i associated with frequency and damping time measurements read [63] ρσ ω 1 = 1 2 √ 2 where ρ is the signal-to-noise ratio (SNR) of the ringdown phase. 11 We thank Aaron Zimmerman for useful correspondence on this issue and for sharing some data of ref. [55]. The small-charge results are available in ref. [55] and agree very well with the full numerical results of ref. [56] and with the small-spin expansion of refs. [52,53]. Our fit is accurate to within 0.5% in the region j ≡ a/M ∈ [0, 0.99].  Figure 4. Projected bound on the BH charge-to-mass ratio Q/M as a function of the BH spin and for two different SNR of the ringdown phase, ρ = 100 and ρ = 7. We assume that the dominant mode has an amplitude approximately three times larger than the second, sub-dominant mode as is the case for many BH coalescences [63]. As a reference, the vertical band denotes the final BH spin (J/M 2 = 0.67 +0.05 −0.07 ) of GW150914 [25] for which ρ ≈ 7 in the ringdown part [46]. We stress that our results neglect terms of the order Q 4 /M 4 or higher and, therefore, are only qualitative in the nearly-extremal limit.
We can convert the errors on the frequency and damping time to errors on physical quantities by using a simple propagation of errors (this procedure yields the correct analytic result in the single-mode case and we expect it to be also accurate generically, since correlations between different physical quantities are small). Specifically, we impose where X = (ω 1 , ω 2 , τ 1 ). It is straightforward to solve the system of three equations above for σ M , σ χ and σ Q ; this yields where f i are cumbersome analytical functions. Finally, we can now view σ Q as an upper bound on Q and use it to estimate the minimum charge that can be measured by a ringdown detection with a certain SNR ρ. Because ρσ Q ∼ 1/Q, the condition σ Q = Q gives a minimum detectable charge that scales as 1/ √ ρ. This value is shown in figure 4, for A 2 /A 1 = 1/3, which is appropriate for a wide range of BH binaries [51,63]. Interestingly, the upper bound on Q/M becomes more stringent as the final BH spin increases and it can improve by some orders of magnitude between χ = 0 and χ ∼ 0.9. The bound scales as ∼ ρ −1/2 so that it becomes more stringent for higher SNR. Thus, for a final BH with χ ∼ 0.9, our simplified analysis suggests that ringdown tests can set an upper constraint of the order |Q| M 0.1 100 ρ . Our analysis is valid up to O(Q 2 /M 2 ) and should be extended to include the nearly extremal case. Nonetheless, it provides an indication for the SNR necessary to constrain the BH charge with a given SNR ratio. For example, the SNR of GW150914 is roughly ρ ≈ 7 in the ringdown part [46]. From the spin measurement J/M 2 = 0.67 +0.05 −0.07 of GW150914 [25], our figure 4 suggests that the ringdown phase of GW150914 does not exclude that the final BH was nearly extremal. (We note that, for a relatively large Q/M ∼ 0.7, our analysis neglects terms of the order Q 4 /M 4 ∼ 0.24, which should modify the final result by roughly a factor of 25%.)

Excitation of gravitational and (hidden) EM modes in dynamical processes
The discussion of the previous section relies on the fact that the gravitational QNMs of a Kerr-Newman BH are affected by the charge Q, cf. eq. (3.16). However, in addition to the shift of these modes, another feature of the QNM spectrum in the presence of charge is the appearance of a new family of modes, which reduce to the standard Maxwell modes of a Reissner-Nordstrom BH when the spin vanishes. As a reference, the fundamental gravitational and Maxwell mode of a (neutral) Kerr BH and of a static Reissner-Nordstrom BH are shown in figure 5. In the Kerr-Newman case, the modes of the Kerr BH acquire charge corrections proportional to Q 2 when Q M [55], whereas the modes of the Reissner-Nordstrom BH acquire corrections proportional to χ in the small-spin case [52,53]. The general case of the gravito-EM modes of a Kerr-Newman BH for arbitrary spin and charge was recently discussed in ref. [56]. Figure 5 shows the well-known fact that the l = 2 Maxwell modes of a BH are well separated from the GW modes. As we now discuss, in dynamical processes like a BH merger these modes are coupled to the gravitational sector and -if excited -they can in principle contribute to the ringdown phase of the GW. In the ringdown analysis of the previous section we have neglected such possibility and assumed that the two dominant modes were the l = m = 2 and the l = m = 3 gravitational modes. In this section we discuss under which conditions the extra Maxwell modes can be neglected.
For this purpose, we consider a simplified model in which a point particle (modelling a small BH) with charge q and mass µ falls radially into a static BH with charge Q and mass M . This model does not capture the effects of the angular momentum of the small -14 -

JCAP05(2016)054
inspiralling BH. However, in the last stages of coalescence, once the orbiting particle reaches the innermost stable circular orbit, it will plunge into the BH. Since the ringdown is excited as the particle crosses the light ring, we suspect that a radial infall will yield a good estimate of the effect that we're trying to study, which is the relative excitation of different QNMs. Our model also does not capture individual spin effects. We can only hope that these effects are subdominant.
For simplicity, we consider that the charge q is either a fractional electric charge e or a hidden charge e h and that the BH is charged accordingly (namely, either Q = Q em or Q = Q h ). In this case, we can set either B µν or F µν to zero and the problem is effectively mapped into an electric charge q plunging onto a Reissner-Nordstrom BH. The general case in which the particle and the BH have both electric and hidden charges is a simple extension of our computation. We consider the charge and the mass of the particle to be small (µ M , q ≤ µ) so that the effect of the particle can be treated perturbatively (cf. appendix C for details). The metric perturbations can be analyzed through a harmonic decomposition, by separating the angular dependence of the perturbations in spherical harmonics. In the frequency domain, the EM and gravitational radiation due to a charged particle falling radially into a charged BH are described by the following coupled radial equations where ω is the frequency, and ψ g and ψ e denote the gravitational and EM master functions, respectively, and r * is the tortoise coordinate of a Reissner-Nordstrom BH, dr/dr * = 1 − 2M/r + Q 2 /r 2 . The potentials V g,e , the coupling functions I i , the source terms S g,e , and a detailed derivation of the above equations are given in appendix C. The source terms depend on the charge q, on the mass µ and on the initial Lorentz factor γ of the particle at infinity. The coupling functions I i are proportional to the BH charge Q, so only when Q = 0 the gravitational and the EM perturbations are decoupled. In the general case we expect a contamination of the EM modes into the gravitational sector and viceversa. We employed two different methods to solve the above equations: the first is a standard Green's function technique which makes use of the solutions of the associated homogeneous systems, whereas the other is a direct integration of the full inhomogeneous system through a shooting method. We explain both procedures in appendix C. The two methods agree with each other within numerical accuracy, and as we explain below they reproduce earlier results in the literature. With the solutions ψ g,e (ω, r) at hand, the GW and EM energy spectra at infinity for each multipole read [65] respectively, where the wavefunctions are evaluated at spatial infinity. The energy spectra defined above are shown in figure 6 for some representative cases. The most salient (and more relevant for this discussion) feature of the GW spectrum (top left panel) is that it has a cutoff at roughly the lowest QNM frequency of the central BH. A  cut is also present for the EM energy spectrum with l = 2 (top right panel) and with l = 1 (bottom panels). The vertical lines in the panels of figure 6 denote the fundamental QNMs, which are presented in table 1 for completeness. The largest contribution to the EM energy spectrum comes from the dipole (bottom panels) which, for sufficiently large values of |Q| ∼ |q|, is even larger than the GW energy flux.
We are now in a position to consider our initial question on the excitation of the EM modes in dynamical processes. Figure 6 shows that the flux for Q = −q = 0.9M displays two peaks: the first one corresponds to the excitation of the l = 2 gravitational mode, whereas the second peak corresponds to the l = 2 EM mode. The latter is excited due to the coupling between gravitational and EM perturbations [cf. eqs. (3.27)-(3.28)], which is a feature of charged BH spacetimes [66,67]. In other words, the GWs emitted in the process contain information about the EM modes of central BH, not only about its gravitational modes.
The ratio between the flux at the two peaks depends on the charge and on the initial Lorentz factor γ of the particle, and is well fitted by where ω e,g are the frequencies at the peaks, and the coefficients (a, b, c) depend on γ. In ap-  Table 1. Fundamental QNMs of a Reissner-Nordstrom BH computed with continued fractions [37] for different values of the BH charge. Gravitational-led and EM-led modes are denoted by s = 2 and s = 1, respectively. Total energy radiated to infinity, summing the gravitational multipoles up to l = 3 and the EM multipoles up to l = 2. We consider γ = 1.1 and the charge-to-mass ratio of the particle q/µ = |Q|/M . pendix C we show some supplemental results about the energy fluxes emitted in the process. In particular, we show the ratio (3.31) for some values of γ [cf. figure 9 in appendix C]. Thus, the relative amplitude of the EM peak is larger when γ ∼ 1, at least for γ 1.2. Furthermore, R is nonnegligible only when the BH is near extremality and only when Qq < 0.

Emission of (hidden) EM radiation in a binary BH merger
Let us now discuss the total energy radiated in the collision. This is shown in figure 7 for the representative case of an initial Lorentz factor γ = 1.1 of the point particle. An interesting aspect of this figure is that in the extremal limit, Q = M , both the EM and gravitational radiation are suppressed. The reason is that in this limit the gravitational attraction cancels the EM repulsion, leading to constant dipole and quadrupole moments (there is in fact a static solution of the field equations describing two maximally charged BHs in equilibrium, the Majumdar-Papapetrou solution [68,69]).

JCAP05(2016)054
For collisions from rest (γ ≈ 1), assuming Qq > 0, we obtain This result can be compared with flat-space estimates using simple quadrupole and dipole Larmor formulae for the GW and EM emission, respectively. For the GW flux we get Using energy conservationṙ = 2M/r for infalls from rest, and we obtain The same procedure for the flat-space Larmor formula gives Thus, a flat-space Newtonian calculation agrees with our numerical results to within a factor two. An interesting quantity is the ratio E e /E g which is shown as a function of Q in the inset of figure 7. From the above relations, in the small-Q limit we obtain which is of the order of Q 2 /M 2 as expected.
Up to now, our results are formally valid only in a perturbative scheme; however, decades of work shows that even fully nonlinear results for equal-mass BHs can be recovered with point-particle calculations if one replaces µ with the reduced mass of the system [70,71]. This substitution would immediately recover (within a factor of two) the results in ref. [72] for the gravitational radiation produced during the head-on collisions of two equal-mass, equalcharge BHs. However, our point particle results overestimate the amount of EM radiation during that same process. The reason is that, as we remarked earlier, for two equal mass-tocharge ratio objects, dipole emission is suppressed and one only gets quadrupole emission. In fact, our l = 2 results for the EM channel agree with the numbers reported in ref. [72] once the extrapolation to equal-mass is done.
In a binary BH coalescence, the total energy loss in GWs can be enormous. For example, the GW luminosity of GW150914 was dE g /dt ≈ 3.6 × 10 56 erg/s [25]. Equation (3.37) shows that, even for BHs with small Q, the EM luminosity can still be very large. In this model the EM luminosity of a GW150914-like event is roughly Even for weakly charged BHs with Q ∼ 10 −4 M , a GW150914-like event would produce an EM luminosity dE e /dt ≈ 10 48 erg/s, comparable to the luminosity of the weakest gamma-ray bursts [73].

JCAP05(2016)054
Nonetheless, the spectra presented in figure 6 show that there exists a cutoff frequency associated to the fundamental EM QNM of the final BH. This mode has a typical frequency of the order of [cf. (the precise value depends on the spin and on the charge [cf. figure 5]). If the BH is electrically charged, the energy is released in GWs and ordinary photons with frequency f f QNM . Lowfrequency photons are absorbed by the interstellar medium if their frequency is smaller than the plasma frequency f em plasma = (2π) −1 4πn e e 2 m e ∼ 10 4 n e 1 cm −3 Hz , (3.40) where n e is the electron number density. Photons with frequency f < f em plasma do not propagate in the plasma. In the interstellar medium n e ≈ 1 cm −3 , so most of the EM energy released in the process is absorbed by the plasma.
On the other hand, if the BH is charged under a hidden U(1) charge, a sizeable fraction of the luminosity is released in dark photons. The latter do not interact with ordinary electrons, but would nevertheless interact with hidden plasma, whose typical frequency reads

Discussion and final remarks
We have shown that, in models of minicharged DM and dark radiation, astrophysical BHs can have large (electric and/or hidden) charge and are uniquely described by the Kerr-Newman metric. In these models the standard arguments that prevent astrophysical BHs to have some electric charge within Einstein-Maxwell theory can be circumvented and, in particular, nearly-extremal BHs with Q ∼ M are also phenomenologically available. Charged BHs in these scenarios are interesting GW sources. The GW signal from the coalescence of two charged BHs contains a wealth of information about the properties of the system. In particular, we have shown that the inspiral and the post-merger ringdown stages provide complementary information. The initial inspiral can constrain a combination of the -19 -

JCAP05(2016)054
initial BH charges [cf. eq. (3.6)] and a rescaling of Newton's constant, whereas ringdown tests are only sensitive to the final BH total charge.
Combined inspiral and ringdown tests might also be performed provided the energy lost in GWs and EM waves during the coalescence is known. Therefore, it would be very interesting to perform fully numerical simulations of charged-BH binary systems close to coalescence (extending the work of ref. [72]) and to estimate the mass, spin and charge of the final Kerr-Newman BH formed after the merger. The combined information from the inspiral and ringdown phases (together with estimates of the mass and charge loss during the merger phase) can be used to disentangle part of the degeneracy appearing in the dipole formula. Such information would be crucial for cross checks similar to those performed for the masses and spins of GW150914 [25,42]. A more detailed analysis in this direction is left for future work.
Likewise, our analysis of the upper bounds derived from ringdown detections can be improved in several ways, for example by considering multiple detections, multiple modes, and a more sophisticated statistical analysis. It is also reasonable to expect that GW150914 was not a statistical fluctuation and that even louder GW events (with higher SNR in the ringdown phase) might be detected in the near future, when second-generation detectors will reach their design sensitivity. In this case our results suggest that an analysis of the entire GW signals can provide stringent constraints on the charge of BHs in minicharged DM models; we hope that this exciting prospect will motivate further studies on this topic.
BH-based tests provide a unique opportunity to constrain the hidden coupling e h , which is otherwise challenging to probe with EM-based tests. An interesting prospect in this direction is the burst of low-frequency dark photons emitted during the merger [cf. eq. (3.38)]. As we have shown, these dark photons are not absorbed by the DM plasma and their luminosity can be extremely high. In some models of minicharged DM, dark photons are coupled to ordinary photons through a kinetic mixing term [10,18,26] proportional to ∼ sin θ (where tan θ := / h ), so that conversion of dark photons to ordinary photons might occur when = 0. The frequency of dark photons emitted in BH mergers are typically smaller than the kilohertz, and therefore next-to-impossible to detect with ordinary telescopes. However, there might be mechanisms in which dark photons can convert to higher-frequency ordinary photons and to ordinary fermions. This conversion might provide an exotic EM counterpart of BH mergers and might leave a detectable signal in current experiments [74,75]. Futhermore, it is in principle possible that the (gigantic) burst of dark photons affects nearby (hidden-) charged stars, via the same mechanisms we described. In such a case, a passing burst of hidden radiation would cause nearby stars to oscillate, with r.m.s fluctuations that could be measurable and extracted from asteroseismology studies, in a phenomena similar to that described recently for GWs [76,77].
We have focused on models of massless dark photons. Massive dark photons are an appealing candidate to explain the muon g −2 discrepancy [78]. If ultralight, these bosons are known to turn spinning BHs unstable [79][80][81] due to a superradiant instability (cf. ref. [82] for an overreview). A more rigorous analysis of this instability is a further interesting application of BH physics in the context of dark-radiation models.
An important open question concerns the formation of charged BHs in these scenarios. Accretion of charged DM particles is a natural charging mechanism. Charged BHs might also form in the gravitational collapse of charged compact stars, the latter might acquire (electric or hidden) charge by DM capture in their interior. To the best of our knowledge there are no studies on the DM accretion by BHs in models of minicharged DM and -in light of -20 -

JCAP05(2016)054
our results -it would be very interesting to fill this gap (see refs. [83,84] for some related work). GW-based bounds on the BH charge might be combined to realistic accretion models to constrain the parameter space of minicharged DM models.

A Quasinormal modes of Kerr-Newman BHs from the geodesic correspondence
In this appendix we discuss the correspondence between the BH QNMs and the properties of spherical photon orbits [59][60][61]. In the static case, the real part of the QNM frequency is related to the azimuthal orbital frequency, whereas the imaginary part of the frequency corresponds to the Lyapunov exponent of the orbit [60]. In the rotating case the relation between modes with generic (l, m) and some geodesic properties is more involved [61]. For simplicity, here we focus on the l = m case in which the analysis is remarkably straightforward, since these modes are associated only with equatorial motion [61]. Let us start with the stationary and axisymmetric line element ds 2 = g tt dt 2 + g rr dr 2 + g θθ dθ 2 + 2g tφ dtdφ + g φφ dφ 2 , where all metric coefficients are functions of r and θ only. The radial motion of null particles on the equatorial plane is governed bẏ where E and L are the (conserved) specific energy and angular momentum of the geodesic, the metric coefficients are evaluated on the equatorial plane, and a dot denotes derivative with respect to the affine parameter of the null geodesic. The light ring and the corresponding ratio L/E are defined by V = 0 = V , where a prime denotes a radial derivative. The orbital frequency of the light ring is simply the azimuthal frequency,

JCAP05(2016)054
evaluated at the light-ring location. Other relevant orbital frequencies are where U(r, θ) = g tt − 2(L/E)g tφ + (L/E) 2 g φφ and the expressions above are evaluated at the light-ring location on the equatorial plane. Ω θ represents the frequencies of small oscillations around quasicircular equatorial orbits along the angular direction, respectively, whereas Ω pre is the precession frequency of the orbital plane.
Reference [61] shows that the QNM frequency in the eikonal limit reads where in the last step we used l = m 1. Interestingly, in the l = m limit the above expression coincides with that derived for static spacetimes in ref. [60], i.e. ω R ∼ lΩ.
Although not relevant for our analysis, for completeness we discuss the geodesic correspondence of the damping time of the modes. In the eikonal limit the latter is related to the Lyapunov coefficient of the orbit [60] again evaluated at the light-ring location on the equatorial plane. Thus, for l = m 1 the complex QNM can be written as [59][60][61] ω R + iω I ∼ Ωl − i(n + 1/2)|λ| , where n is the overtone number. Note that the above expression formally coincides with that obtained in refs. [60] for static spacetimes and it extends the results of ref. [59] which are valid only for slowly-rotating BHs. A more involved result for QNMs with generic (l, m) is derived in ref. [61]. Strictly speaking, the geodesic prediction (A.8) should only be valid in the eikonal limit, i.e. when l 1. However, in figure 8 we show that the analytical result (A.8) agrees remarkably well with the exact numerical results for the QNMs of a Kerr-Newman BH even when l = m = 2. Relative errors are always smaller then ≈ 4% for any spin, both in the neutral case (top panels of figure 8) and in the Kerr-Newman case with Q/M = 0.2 (bottom panels of figure 8). In the latter case the exact results are only available for l = 2 [55], but their deviation from the geodesic predition is always smaller than 3%.
In the main text, we have used this striking agreement to estimate the l = m = 3 modes of a weakly-charged Kerr-Newman BH. Note that the deviations from the geodesic predictions are likely smaller than the observational errors on these modes, therefore this approximation should not affect our analysis significantly.

B Fisher matrix analysis
We follow ref. [57] for the analysis of uncertainties associated with measurements of ringdown waveforms in noise. We assume that the GW signal during the ringdown phase can be In this expansion the spheroidal functions S lmn = S lm (aω lmn ) are evaluated at the (complex) QNM frequencies, so they are complex numbers (henceforth we drop the angular dependence on the S lmn ). The waveform measured at a detector is given by where F +,× are pattern functions that depend on the orientation of the detector and the direction of the source, namely F + (θ S , φ S , ψ S ) = 1 2 (1 + cos 2 θ S ) cos 2φ S cos 2ψ S − cos θ S sin 2φ S sin 2ψ S , (B.3a) F × (θ S , φ S , ψ S ) = 1 2 (1 + cos 2 θ S ) cos 2φ S sin 2ψ S + cos θ S sin 2φ S cos 2ψ S . (B.3b) -23 -

JCAP05(2016)054
We will follow the prescription outlined in ref. [57] to compute the SNR ρ. We assume large quality factor Q lmn and average the source over sky position and over detector and BH orientations, making use of the angle averages: F 2 + = F 2 × = 1/5, F + F × = 0, and |S lmn | 2 = 1/4π.
With a given noise spectral density for the detector, S h (f ), one defines the inner product between two signals h 1 (t) and h 2 (t) by whereh 1 (f ) andh 2 (f ) are the Fourier transforms of the respective gravitational waveforms h(t). The components of the Fisher matrix Γ ab are then given by where θ a are the source parameters. In the limit of large SNR, if the noise is stationary and Gaussian, the probability that the GW signal s(t) is characterized by a given set of values of the source parameters θ a is where δθ a = θ a −θ a , and p (0) (θ) represents the distribution of prior information. An estimate of the rms error, ∆θ a = ( (δθ a ) 2 ) 1/2 , in measuring the parameter θ a can then be calculated, in the limit of large SNR, by taking the square root of the diagonal elements of the inverse of the Fisher matrix, The Fisher matrix components in the parameter basis of (A + , A × , φ + lmn , φ × lmn , f lmn , Q lmn ), were computed and presented in ref. [57] [here f lmn := ω lmn /(2π)]. In the large Q lmn limit, they read Here, .

C Technical details on the radial infall of a point charge into a charged BH
In this appendix we give some details about the derivation and integration of the coupled system (3.27)-(3.28). We follow the procedure outlined in ref. [65], but correcting typos and possible errors. Our integration technique is fully consistent, and has been validated by two independent codes, as well as with previous results in the literature for uncharged BHs. One of the codes, written in Mathematica, is freely available online [85].

C.1 EM and metric perturbations
The spacetime metric due to a point charge falling into a charged BH can be written as g ab is a small perturbation to the background Reissner-Nordstrom geometry, g (0) ab = −e v dt 2 + e −v dr 2 + r 2 dθ 2 + r 2 sin 2 θdφ 2 , (C.1) where e v = 1 − 2M/r + Q 2 /r 2 . Here we consider the Regge-Wheeler decomposition of h ab . In the case of a radially falling particle, the metric perturbations have even (polar) parity and h ab can be written as functions of (t, r) only. The stress-energy tensor of the particle can also be decomposed in terms of spherical harmonics, where, once again, the A's are functions of (t, r) only. In a stationary background, we can explicitly eliminate the time dependence by a Fourier transformation of the perturbation functions: we define the Fourier transform of a function ψ(t, r) as Below we shall consider already the transformed Fourier quantities, using the same notation for each perturbation function. By substituting the expressions for the metric and the vector potential into the Einstein equations and expanding up to first order, we obtain the following set of equations

JCAP05(2016)054
In the above equations we already used the fact that H 2 = H 0 , required by the (θ, φ) component of the perturbed Einstein equations. Due to the Bianchi identities, not all of the above equations are independent. Let us first look into the Maxwell equations (C.11)-(C.13). Equation (C.11) is automatically satisfied, and is a consequence of (C.11)-(C.12) and the continuity equation for the currents. Indeed, it is easy to see that, by defining f 12 = e v ψ e and manipulating eqs. (C.12) and (C.13), we obtain eq. (3.28).
Simplification of the gravitational sector is more involved [65]. By manipulating eqs. (C.6)-(C.9) we can obtain a system of two differential equations for H 1 and K, and an algebraic relation between H 0 , H 1 , and K. By defining the vector ψ = (K, −iω −1 H 1 ), the differential equations can be written as where A = α β γ δ , S = (S 1 , S 2 ), and the coefficients (α, β, γ, δ) are linear in ω 2 , e.g. α = α 0 + ω 2 α 2 . To simplify the system further, we wish to find a transformation ψ = Fψ +S , (C. 15) with F = f g h k andψ = (ψ g , ψ g ), such that the function ψ g obeys the following differential equation Note that S z will also involve EM perturbations, since these appear as a source terms in the first-order equations. By comparing the coefficients of different order in ω in the differential equations it is possible to obain a relation involving (f, g, h, k, n, V g , S z ) and their derivatives. We obtain g = 1, n = e v , k = −re −v , and The remaining functions are given below.

C.3 Numerical procedure
To solve the perturbation equations (3.27) and (3.28) we employed two different methods. The latter agree with each other within numerical accuracy. The first method relies on the Green's function approach, also called method of variation of parameters [87]. Let Ψ = (ψ g , ψ e , dψ g /dr * , dψ e /dr * ). The system (3.27)-(3.28) can be written as d dr * Ψ + V Ψ = S. (C.30) We start by constructing the fundamental 4 × 4 matrix X, whose columns are independent solutions of the associated homogeneous differential equations. The independent solutions -28 -

JCAP05(2016)054
can be obtained in the following way. We notice that at the horizon, the required solutions have the following form: ψ g,e ∼ A r + g,e e −iωr * , (C. 31) and at infinity we have ψ g,e ∼ A ∞ g,e e iωr * . (C.32) In the above equations A r + ,∞ g,e are constants. The first two columns of X can be obtained by integrating the homogeneous equations from the horizon outwards with boundary conditions (C.31), by choosing two independent solutions, say (A Let us now expose the second method that we used to solve the system (3.27)-(3.28). The shooting method relies on an integration of the full system of inhomogeneous equations. First, we impose that the solutions near the horizon are of the form (C.31). We then integrate the full equations up to infinity where the general solution will be a superposition of the ingoing and outgoing modes, namely ψ g,e ∼ A in g,e e −iωr * + A out g,e e iωr * . (C.36) The physical solutions corresponding to a particle falling into the BH require A in g,e = 0, and, therefore, this becomes a two-parameter shooting problem for the amplitudes A r + g,e . With the proper values of A r + g,e , we can compute the amplitude of the GW and EM waves at infinity, which enable us to compute the GW and EM energy spectra through eqs. (3.29)-(3.30).

C.4 Supplemental results
In this section we present some supplemental results on the GW and EM emission in the radial infall of a point charge into a Reissner-Nordstrom BH. Figure 9 shows the ratio between the two peaks of the quadrupolar GW flux [cf. eq. (3.31)] as a function of the BH charge. As discussed in the main text, the relative amplitude of the EM peak is nonnegligible only when Qq < 0 and when the BH is highly charged, cf. figure 6. Note that the ratio depends on the boost factor γ.
Finally, for completeness in figure 10 we present some representative cases for the GW and EM energy spectra for the radial infall of a high-energy point charge, i.e. γ → ∞. In this regime our results for the EM flux are in perfect agreement with those presented in ref. [88].  Figure 10. Quadrupolar GW (left panel) and EM (right panel) energy spectra for a high-energy charged particle plunging radially on a RN BH. The particle and the BH have the same charge-to-mass ratio. In the legend of both panelsQ := Q/M andq := q/µ.