Constrains on the electric charges of the binary black holes with GWTC-1 events

Testing black hole's charged property is a fascinating topic in modified gravity and black hole astrophysics. In the first Gravitational-Wave Transient Catalog (GWTC-1), ten binary black hole merger events have been formally reported, and these gravitational wave signals have significantly enhanced our understanding of the black hole. In this paper, we try to constrain the amount of electric charge with the parameterized post-Einsteinian framework by treating the electric charge as a small perturbation in a Bayesian way. We find that the current limits in our work are consistent with the result of Fisher information matrix method in previous works. We also develop a waveform model considering a leading order charge effect for binary black hole inspiral.

Testing black hole's charged property is a fascinating topic in modified gravity and black hole astrophysics. In the first Gravitational-Wave Transient Catalog (GWTC-1), ten binary black hole merger events have been formally reported, and these gravitational wave signals have significantly enhanced our understanding of the black hole. In this paper, we try to constrain the amount of electric charge with the parameterized post-Einsteinian framework by treating the electric charge as a small perturbation in a Bayesian way. We find that the current limits in our work are consistent with the result of Fisher information matrix method in previous works. We also develop a waveform model considering a leading order charge effect for binary black hole inspiral.

I. INTRODUCTION
The prominent black hole (BH) no-hair theorems [1,2] imply vacuum astrophysical black hole can be described by a special case of Kerr-Newman metric, which could be characterized by its mass, spin and charge only [3][4][5]. Although astrophysical BHs are usually considered as electrically neutral due to charge neutralization by astrophysical plasma, quantum discharge effects, and electron-positron pair production [6][7][8], an accurate upper limit on the amount of the charges of BHs is still absent. Therefore, it is important and meaningful to give a quantitative measurement on the amount of charges of BHs. Besides, many novel charging mechanisms has been discussed theoretically, such as primordial black hole [9], and BHs formed by millicharge dark matter [10]. In the later case, the discharge process is much slower than for ordinary plasma so the charged BHs are viable. So it is meaningful to study the effect on gravitational wave if electric charges are presented in these special theories.
At present, using the shadow of supermassive BH to estimate its charge has been developed and widely practiced [11], for example Sgr A* and M87* measured by VLBI [12][13][14][15]. These proposals are either far from being accurate or model-dependent. Therefore, developing new technique to model-dependently test the charged property is needed, and gravitational wave has always been expected, which encodes the BH's information [16]. Fortunately, after decades of hard work, LIGO-Virgo announced the first detection of the gravitational wave (GW) signal GW150914 in 2015, generated by the merger of a binary BH [17]. To characterize GW, we often use inspiral, merger and ringdown to describe the whole coalescence of a binary BH, where the inspiral stage is generally described by post-Newtonian theory [18], merger stage is generally approximated by the numerical simulation [19], and ringdown stage is described by the BH perturbation theory [20,21]. Some waveform models, such as IMRPhenomPv2 (IMR) waveform model, have incorporated these three stages [22].
The GW waveform will be affected if the two BHs are electrically charged. For simplicity, consider the electric dipole radiation during their orbital motions, which will be reflected in the phase of inspiral stage, and thus provides us a new trick to detect charged binary BH. Depending on the amount of the charge, there are two different treatments. When the charge of the black hole is small, the dipole radiation can be taken as a perturbation or correction to the phase of the inspiral stage of the GW waveform, which can be incorporated into the parameterized post-Einsteinian (ppE) framework [16]. In this case, the correction due to the electric dipole radiation is completely described by the coefficient of the −1 post-Newtonian (PN) order in the waveform, the rest part of the waveform is the same as that of the accurate waveform model describing the coalescence of two neutral black holes, such as the well-known phenomenological waveform model [22]. This is similar to the dipole radiation caused by the scalar charges carried by black holes in the scalar-tensor theory [23]. However, if the charge of the black hole is large, the ppE framework may not be applicable anymore, as the latter requires the expansion in PN is always linear, which may not be assured as a prior.
We also consider only the leading order gravitational quadrupole radiation and the electric dipole radiation, based on which we call such waveform the leading order charged (LOC) one. It is apparent that the LOC waveform has a bad accuracy when applied to the analysis of GW data unless the inspiral duration is long enough and the orbit of the two black holes decays very slow. Despite this shortcoming, the LOC waveform provides us a toy model such that the effect of the charge can be demonstrated explicitly without the assumption that the charge of the black hole is small. The detailed calculation of the LOC waveform is shown in Sec. V.
In this work, based on IMR waveform model in LIGO-Virgo Algorithm Library [24], we use the Bayesian method to test the dipole radiation of GW signals [25,26] with the ppE framework. As the main conclusion, we find no visible charge taking by astrophysical BH.
The rest of this paper is organized as follows. In Sec. II, we introduce the ppE framework. In Sec. III, we briefly introduce the Bayesian method for the GW data processing and the results are presented in Sec.IV. In order to verify the consistency of the results obtained from ppE method, in Sec. V, we introduce the toy model, the LOC waveform, to study the effect of the amount of the charge on the waveform. In Sec. VI, we summarize the calculation results and present the conclusion, and discuss the deficiencies of this work as well as what can be done further. We assume c = 1 throughout the paper unless otherwise specified.

II. EFFECTS ON GRAVITATIONAL WAVE SIGNALS
In this section, we introduce the waveform models adopted in this work. We treat the charge effect as a small perturbation, which can be well described by the parameterized post-Einsteinian framework. In this case, the phase of the inspiral part of the charged binary BHs can be incorporated into the ppE formalism, which describes the gravitational waveforms of theories alternative to general relativity (GR) in a model-independent way [16,25]. Formally, the effect of the electric dipole radiation can be captured by the ppE parameter entering at −1 post-Newtonian order, similar to the dipole radiation from the scalar-tensor theory [23].
Assuming the early-inspiral stage of the IMR waveform model is reduced to the form of h e−ins (f ) = A e−ins (f )e iΨe−ins , then the waveform model due to modifications from different modified gravity effects can be written as In ppE framework [25,26], where β is the amplitude coefficient and b is the exponent coefficient. Note that the modification enters at (b + 5)/2 PN order for a related modified gravity. M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 is the chirp mass with component masses m 1 and m 2 , and f is the frequency of the GW. The study in Tahura et al. [27] shows that at least in theories where the leading corrections enter at negative PN orders, the phase-only analyzes can produce sufficiently accurate constraints. In this analysis, we only consider the correction on the phase and neglect the correction to the amplitude. The frequency at the end of this stage is given by f c = 0.018/[G(m 1 + m 2 )] [26], above which this model is calibrated with numerical-relativity data and can not be applied to ppE formalism.
The effect of the charge on the waveform can be taken as a perturbation term in the phase [10,28], where η = m 1 m 2 /(m 1 + m 2 ) 2 is symmetric mass ratio, ζ represents the difference between the charges of the binary BH and is defined by is the charge-to-mass ratio and q i is the electric charge of a BH. Here κ = 1 − λ 1 λ 2 ≈ 1 since we consider the effects of λ 1 and λ 2 as small perturbations. In this case, β = − 5 3584 η 2/5 ζ 2 and b = −7, corresponding to Eq. (2). In TABLE I. of Yunes et al. [25], Yunes et al. obtained constraint |β| < ∼ 1.6 × 10 −4 at −1 PN order for GW150914. The relationship between β and ζ is: β = −5/3584η 2/5 ζ 2 . Thus, one can derive the corresponding constraint ζ < ∼ 0.45. There are two ways to enhance the constraint on the charge. Firstly, one may expect the increase of the sensitivities of GW detectors, and a more stringent constraint on the charge should appear. Secondly, a full gravitational waveform model of charged binary BH is excepted to give a more stringent constraint on the charge. This full GW waveform should apply to the coalescence of two black holes carried with arbitrary amount of charge. It is a challenging job and in this work we move forward by developing the leading order charged (LOC) waveform model, which is introduced in Sec. V.

III. BAYESIAN INFERENCE METHODS
To infer the uncertainty of the source parameters θ, which are quantified by the posterior probability distribution p( θ|d, M ), we perform Bayesian analysis with the prior p( θ|M ) and the likelihood with Gaussian noise assumption for the GW data, where d i is the data of the i-th instrument, M is the model assumption, N is the number of detectors in the network of Advanced LIGO and Advanced Virgo, and h i is the waveform model calculated with θ for the i-th detector. The noise weighted inner product a(f )|b(f ) is defined by where f low and f f igh are the high and low pass cut-off frequencies respectively, S n (f ) is the power spectral density of the detector noise. The Bayesian theorem is described by where p(d|M ) = p(d| θ, M )p( θ|M )d θ is the evidence of a specific model assumption.
There are fifteen parameters in the IMR waveform model, including the redshifted chirp mass (M z ), mass ratio (q), luminosity distance (d L ), inclination angle (θ jn ), the reference orbital phase (φ c ), the geocentric time (t c ), the polarization angle (ψ), the two dimensional sky location, and six spin parameters. And the early-inspiral stage of the IMR waveform model for ppE framework (hereafter e-insp-ppE waveform model) has an additional ppE parameter ζ. We marginalize over the reference phase φ c and the geocentric time t c , thus we have 14 free parameters for e-insp-ppE waveform model. The GW data and power spectral density for each event are downloaded from LIGO-Virgo GW Open Science Center [29]. To estimate parameters with data from the first two observation runs (O1 and O2), we carry out Bayesian inference with Bilby [30], using Pymultinest [31] as our sampler.
The prior on the redshifted chirp mass is chosen to be uniform in the range of 5 M ≤ M z ≤ 20 M for GW151226 and GW170608 while 5 M ≤ M z ≤ 50 M for the other events, the prior on the mass ratio is chosen to be uniform in the range of 0.25 ≤ q ≤ 1. We apply comoving uniform prior on the luminosity distance 50Mpc ≤ d L ≤ 4000Mpc, while the prior on the inclination angle is chosen to be uniform in the range of −1 ≤ cos ι ≤ 1. The prior on the polarization angle is chosen to be uniform in the range of 0 ≤ ψ ≤ π, and the prior of the sky location is chosen to be uniform in spherical coordinates. The prior on |ζ| is chosen to be uniform in the range of 0 ≤ |ζ| < √ 2. For other parameters in the e-insp-ppE waveform model, we use the same priors presented in Abbott et al. [32].

IV. RESULTS OF BAYESIAN INFERENCES
In this section, on the promise that the BH has a small amount of charge, we set up constraints on the dipole radiation with the ppE framework. We analyze with e-insp-ppE waveform model to constrain the dipole radiation, with a high frequency cut-off at f c = 0.018/[G(m 1 + m 2 )]. Specifically, this analysis is applied to GW150914, GW151226, GW170104, GW170608, and GW170814. We neglect other events since the early-inspiral signal-to-noise ratios (SNRs) of them are all less than 6. As shown in Table I, for GW170608, the GW data can constrain |ζ| to be less than about 0.21 at 90% credible level, while the loosest case is given by GW170104, |ζ| ≤ 0.70 (at 90% credible level) due to its lowest SNR of the inspiral signal among all these five events. In Yunes et al. [25], the authors got constraints on scalar dipole radiation through a Fisher parameter estimation study, using a fitted spectral noise sensitivity curve. As a result, the constraint on GW150914 is |ζ| < ∼ 0.45, and the constraint on GW151226 is |ζ| < ∼ 0.24. Their results agree well with ours shown in Table I. One can also get constraints on |ζ| by reweighting the posteriors of parameters δφ −2 from results in Abbott et al. [26], where δφ −2 = −5/84ζ 2 is the parameterized violation of GR at −1 PN. Recently, a similar analysis was performed in Wang et al. [33] with both reweighting method and Bayesian inference method, and they find that the Bayesian analysis is more reliable than the reweighting analysis.
To further discuss the reliability of our results, we apply Bayesian inference on GW data of GW150914 with a set of low-pass cut-off frequencies f high . And we show the chirp mass distribution in Figure 1, in which the distribution becomes wider as f high less than 400 Hz. Nevertheless, the result is not biased even the cut-off frequency is as low as 50 Hz, although the uncertainty is a bit larger.
Actually, the charge will significantly affect the space-time metric of the binary system if the charge is as large as the maximum allowable value. This effect will be fully reflected in waveform model when higher PN order terms are taken into account. However, an analytic solution will not be available in the forthcoming future and one can only get it numerically, i.e. Bozzola and Paschalidis [34]. In Bozzola and Paschalidis [34], the binary charged black holes were studied via numerical general relativity method, where the higher PN order effects of the charge are included and the amount of charge is set free. One of the main conclusions from this paper is that the greatest difference between charged and uncharged black holes arises in the earlier inspiral. This provides a convincing evidence to support the reasonability of our work.
For future space-based GW detectors such as Laser Interferometer Space Antenna (LISA) [35] and TianQin [36], the duration of GW strain from BH binaries are much longer [37,38]. Since the dipole gravitational radiation dominates at −1 PN order, it is more suitable for detecting by LISA and TianQin. Barausse et al. [39] shows that joint observations of GW150914-like systems by LIGO-Virgo and LISA will improve bounds on dipole emission from BH binaries by several orders of magnitude relative to current constraints. We expect that with multi-band GW the constraint on the electric charge of BH can be improved as well.

V. THE LEADING ORDER CHARGED WAVEFORM
In this section, we try to learn more about the effect of the charge by developing a LOC waveform model. For this waveform model, we do not consider the spin of binary black holes since the spin evolution of the binary BH only occurs at high orders of the PN expansion [40,41]. Therefore, we only need to analyze the orbit (circular) evolution in the inspiral phase due to the energy loss. For two point particles with mass m 1 and m 2 and charge q 1 and q 2 respectively, we define λ i = q i /( √ Gm i ). They orbit each other with orbital radius R, and the orbit decays with time t. The dissipation of total energy can be divided into two parts, one is the emission of GW, the other is electromagnetic dipole radiation, which can be written as [10] where M = m 1 + m 2 is total mass of the binary, G eff = G(1 − λ 1 λ 2 ) is the effective Newton constant, ζ = |λ 1 − λ 2 |/ √ 1 − λ 1 λ 2 is used to characterize the difference between the charges of two BHs, and η = m 1 m 2 /(m 1 + m 2 ) 2 is symmetric mass ratio. The evolution equation of the orbital radius arising from Eq. (7) is where A = 64GG 2 eff M m 1 m 2 /5 and B = 4ζ 2 G 2 eff m 1 m 2 /3. When B is not equal to zero, the relation between R and t can be parameterized as where τ = t c − t and t c is the coalescence time. It is clear that an analytical solution of R(t) is not possible in Eq. (9) , so we have to solve it numerically. Due to the last part of Eq. (9), the numerical solution is not accurate enough when BR/A ∼ 0, thus we expand Eq. (9) as a series of BR/A . According to the Kepler's law, the orbital frequency is G eff M/R 3 and the gravitational frequency is twice as much as the orbital frequency, so we have ω gw = 2πf gw = 2 G eff M/R 3 . The waveform in time domain contains two parts ( [41,42]) where ι is inclination angle, d L is luminosity distance and the phase of waveform is where t ISCO is the time when the BH reaches the innermost stable circular orbit (ISCO). We cut the phase before t ISCO for three reasons: first, the phase before t ISCO can be included in φ c and the effect is equal to a time shift, which means this does not affect the results of parameter estimation; second, ω gw is infinity when τ = 0, and this can not be integrated; last but not least, the LOC waveform model cannot truly describe the motion of the binary BH when the orbital distance is too small. According to Eqs. (8) and (9), t ISCO is fixed for a given R ISCO . In reality, it is difficult to know the exact values of M and λ for the final BH. Because with a given λ 1 λ 2 and |ζ|, we still cannot uniquely determine the respective charge of the two BHs. For example if λ 1 = a, λ 2 = b is the solution then λ 1 = b, λ 2 = a could still be the solution. Approximately we take M = m 1 + m 2 , λ = min{| m1λ1+m2λ2 m1+m2 |, | m1λ2+m2λ1 m1+m2 |}. then in the similar way to Pugliese et al. [43] we obtain where C = −(9 − 8λ 2 − 4 √ 4λ 4 − 9λ 2 + 5) 1/3 . The ISCO of a charged BH decreases with the charge, and we have R ISCO = 4GM for |λ| = 1 and R ISCO = 6GM for the non-rotating uncharged BH, respectively. Based on the discussion above, we obtain the LOC waveform model. Obviously, it reduces to the 0PN when both BHs are uncharged. As we have already pointed out, under the ppE framework, LIGO-Virgo has very weak restrictions on the charges of GW150914, with |ζ| < ∼ 0.45. Here we show that such large amount of charge could have non-negligible effect on the waveform. We choose parameters similar to GW150914, where m 1 = m 2 = 30 M , d L = 540Mpc, ι = 0, φ c = 0, and we set low cut-off frequency to be f low = 20 Hz. For the charge of GW150914, if |ζ| = 0.45, then possible allowed range of λ 1 λ 2 is (−0.053, 0.80). For proper comparison, we have chosen four special values with λ 1 λ 2 = {−0.05, 0.05, 0.15, 0.25} to study the effect of charge. As shown in Figure 2, both the amplitude and phase of the GW signal will be affected significantly if GW150914 is endowed with a large amount of charges. The effect of charge on the early inspiral of the gravitational waveform is more significant, which is consistent with previous works [10,16,25]. It is also worth noting that G eff and M have similar effects on shaping the GW signals, the smaller G eff , the lower GW amplitude, as described in Eq. (10). We also find that the larger value of λ 1 λ 2 , the slower frequency evolution of GW signals.

VI. SUMMARY
In this work, we study the charge effect on the inspiral stage of the GW waveform. By considering the charge effect as an perturbation in the inspiral stage, we perform Bayesian inference on five events of the first Gravitational-Wave Transient Catalog (GWTC-1), GW150914, GW151226, GW170104, GW170608, and GW170814, whose SNRs are larger than 6 after applying a low-frequency cut off f c . Based on the ppE formalism, the most stringent limit from FIM up to now is |ζ| ∼ 0.24 [10,25], while our Bayesian based result is |ζ| ∼ 0.21 at the 90% credible level, which is given by GW170608. The analysis on GW151226 and GW170814 give similar constraints, i.e., ζ < 0.26 and ζ < 0.27 at 90% credible level respectively. For GW150914, the constraint is ζ < 0.42 at 90% credible level, which is similar to the analysis before [25].
To explicitly show the effect of the charge on the waveform, we developed the LOC waveform model in Sec. V. Although due to the lack of accuracy, this waveform cannot be used to analyze the realistic GW data, this toy model has the advantage of being adjustable and intuitive. The LOC waveform is obtained by taking into account the dissipation effect caused by the electric dipole radiation and the quadrupole gravitational radiation on the circular orbit of the charged binary BHs. Like the 0PN waveform model, the spin of BH is not considered for the LOC waveform model as which emerges in the waveform only at higher orders. Besides, the analytical LOC waveform is achievable only if the charge is treated as a correction, so the general charge case is gotten numerically. We find that both the amplitude and phase of the GW signal will be strongly affected if the BHs are endowed with a large amount of charge.
The work in this paper can be improved in several aspects. For example, in the employment of the ppE framework, we only consider the electric dipole radiation, so when the two black holes have the same charge-to-mass ratio, the charge effect disappears in the phase of the ppE waveform. To overcome this problem, one can consider the electric quadrupole radiation. Moreover, as we have shown above, the constraints on the electric charges of the black holes are still not very stringent, one of the possible ways to improve this is to add higher order corrections to the waveform. However, this may not be helpful, because the parameterized constraints on the PN coefficients obtained by LIGO-Virgo show that the -1 PN gets the most stringent constraint than other higher PN orders [26], which we expect also applies to the study of the charged black holes. The contributions from the higher orders corrections of some other specific theories has been studied in Yunes et al. [25].
Above all, we give the constraints on the dipole radiations of GW events from GWTC-1. These constraints can be also converted to the constraints on radiation effects in other forms, such as the radiation caused by the magnetic or other U(1) dark charges carried by black holes. In the future, similar analysis can be applied to GW data that published in GWTC-2 [44].

VII. ACKNOWLEDGMENTS
We sincerely thank Vitor Cardoso and Laura Bernard for their kindhearted help. Hai-Tian Wang also appreciates Yi-Fu Cai, Jian-dong Zhang, Bing-bing Zhang, Si-ming Liu, Can-ming Deng, Ming-Zhe Han, and Shao-Peng Tang for insightful comments and discussions. This work has been supported by NSFC under grants of No. 11525313 (i.e., Funds for Distinguished Young Scholars), No. 11921003, No. 11335012, No. 11325522, No. 11735001, No. 11703098, No. 11847241 and No. 11947210. Peng-Cheng Li is also funded by China Postdoctoral Science Foundation Grant No. 2020M670010. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is