The exploration of dark sector interactions via gravitational waves GWs from binary inspirals has been a subject of recent interest. We study dark forces using extreme mass ratio inspirals EMRIs pointing out two issues of interest. Firstly the innermost stable circular orbit ISCO of the EMRI which sets the characteristic length scale of the system and hence the dark force range to which it exhibits enhanced sensitivity probes force mediator masses that complement those studied with supermassive black hole SMBH or neutron star binaries. The LISA mission the proposed Ares detector will probe mediators with masses m SUBVSUB  10SUP16SUP eV m SUBVSUB  10SUP18SUP eV corresponding to ISCOs of 10SUP6SUP M SUBSUB 10SUP8SUP M SUBSUB central SMBHs. Secondly while the sensitivity to dark couplings is typically limited by the uncertainty in the binary component masses independent mass measurements of the central SMBH through reverberation mapping campaigns or the motion of dynamical tracers enable one to break this degeneracy. Our results therefore highlight the necessity for coordinated studies loosely referred to as multimessenger between future Hz mHz GW observatories and ongoing and forthcoming SMBH mass measurement campaigns including OzDESRM SDSSRM and SDSSV Black Hole Mapper.
['astrophysical black holes', 'gravitational waves / sources', 'massive black holes', 'Astrophysics - High Energy Astrophysical Phenomena', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Phenomenology']
Exploring dark forces with multimessenger studies of extreme mass ratio inspirals
{'Badal Bhalla, a Kuver Sinha, a and Tao Xu a': 'a Homer L. Dodge Department of Physics and Astronomy, University of Oklahoma, Norman, OK 73019, USA  \nE-mail: badalbhalla@ou.edu, kuver.sinha@ou.edu, tao.xu@ou.edu  \nAbstract. The exploration of dark sector interactions via gravitational waves (GWs) from binary inspirals has been a subject of recent interest. We study dark forces using extreme mass ratio inspirals (EMRIs), pointing out two issues of interest. Firstly, the innermost stable circular orbit (ISCO) of the EMRI, which sets the characteristic length scale of the system and hence the dark force range to which it exhibits enhanced sensitivity, probes force mediator masses that complement those studied with supermassive black hole (SMBH) or neutron star binaries. The LISA mission (the proposed µ Ares detector) will probe mediators with masses m V ∼ 10 -16 eV ( m V ∼ 10 -18 eV ), corresponding to ISCOs of 10 6 M ⊙ ( 10 8 M ⊙ ) central SMBHs. Secondly, while the sensitivity to dark couplings is typically limited by the uncertainty in the binary component masses, independent mass measurements of the central SMBH through reverberation mapping campaigns or the motion of dynamical tracers enable one to break this degeneracy. Our results therefore highlight the necessity for coordinated studies, loosely referred to as \'multimessenger", between future µ Hz -mHz GWobservatories and ongoing and forthcoming SMBH mass measurement campaigns, including OzDES-RM, SDSS-RM, and SDSS-V Black Hole Mapper.', '1 Introduction': '1  \nA wide range of physics scenarios beyond the Standard Model (SM) and General Relativity (GR) have as their template the following potential generated by an object of mass M :  \nU = -GM r [ 1 + ˜ α \' e -r λ ] , (1.1)  \nwhere the extra contribution to the Newtonian potential, sometimes dubbed a \'fifth force" or \'dark force" depending on the context, is controlled by an effective coupling strength parameter ˜ α \' and a characteristic length scale λ . The Yukawa contribution can also be interpreted as the static limit of an interaction mediated by bosons of mass  \nm V ∼ 1 λ . (1.2)  \nThis widely used framework has historically served to parametrize deviations from GR and the inverse square law [1-3], the effects of light string moduli [4-7] and, in more recent times, as a template to study light mediators in dark sector models [8-13]. Interpreted purely in terms of deviations from GR (i.e., without reference to dark sector physics), this formalism describes, for example, models of large extra dimensions (for small values of λ ≪ 10 -3 m) [14, 15] as well as braneworld [16-18], and the linearized regime of massive gravity theories (for large values of λ ) [1]. In the study of dark sectors, on the other hand, ˜ α \' takes on the role of the dark sector coupling and m V the role of the dark mediator mass, while the Yukawa interaction is assumed to be sourced by the accumulation of dark charge on the central and secondary objects. Indeed, a host of interaction potentials (not just the Yukawa potential) are possible, depending on the nature of the mediator being exchanged (we refer to [19, 20] for a unified treatment of such exchange interactions, and [21, 22] for recent work along these lines). We will mainly restrict our attention to long-ranged Yukawa potentials in this work.  \nThe phenomenological framework of Eq. (1.1) encapsulates a wide variety of theories. For the purposes of this work, it is useful to distinguish between ( i ) theories of modified gravity or fifth  \nforces, which would apply to all bodies; and ( ii ) theories of dark sector physics, where Eq. (1.1) only applies to bodies with dark charge. While our main interest is in the second category, we briefly summarize the various systems that have been used to constrain fifth forces in general. The parametrization in Eq. (1.1) dictates the use of diverse experimental settings, where the characteristic scale of the probe should be λ probe ∼ λ . The pursuit of probing λ through measurements at laboratory scales has a long history [23]. In the large λ (light mediator) regime, competitive bounds arise from planetary motion and the motion of stars around Sgr A ∗ (we refer to [23] and [24] for comprehensive treatments of these two frontiers). The most stringent Solar System constraints on the (˜ α \' , λ ) plane come from modifications of Kepler\'s Third Law. On the other hand, broadly three avenues have been pursued with respect to constraints from Sgr A ⋆ [25]: ( i ) Orbital precession data by near infrared monitoring of stars that are sufficiently close to Sgr A ⋆ ; ( ii ) high-precision timing of pulsars by existing and future radio telescopes; and ( iii ) space- and time-resolved probes of the accretion flow by the Event Horizon Telescope. The advent of gravitational wave (GW) astronomy has opened up new opportunities in this frontier (we refer to [26] for a recent review) 1 .  \nCompact object mergers in particular present unique opportunities: the population distributions and the GW waveforms emanating from these systems are sensitive to modifications of the type parametrized in Eq. (1.1). Different classes of compact binary systems enable different scales λ to be probed, depending on the characteristic frequency of the GW emission. The dark sector interpretation of Eq. (1.1), which is our main interest, works under the assumption that both the interacting objects accumulate dark charge. In the case of inspiraling neutron stars, such accumulation could result from dark matter capture; in the case of solar or super massive black holes, such accumulation could result from spike formation, or from the fact that such black holes are charged under millicharged dark fermions [28]. Here, too, the recent literature is substantial, and we refer to [8-12] for some representative studies.  \nThe purpose of this paper is to explore the possibility of probing dark forces with extreme mass ratio inspiral (EMRI) systems consisting of a central object - the super-massive black hole (SMBH) - and a stellar mass black hole in a decaying orbit around the central object [29, 30]. These systems produce millions of cycles in the GW sensitivity band of LISA, enabling probes of the background geometry, as well as the detection of small changes to the radiation-reaction force. The inspiral process is driven mainly by the emission of gravitational radiation, thus providing an excellent laboratory in which to probe any altered emission due to dark forces. EMRI systems have been extensively used as laboratories of GR and modified gravity: we refer to [31, 32] for papers relevant to our study. Their use as probes of dark forces and dark sectors is somewhat more limited and recent: we refer to [31, 32] for scalar-mediated dark force studies, and [9] for somewhat related dark sector studies of SMBH binaries.  \nThere are two main reasons why EMRIs should be studied in this context. The first is that the innermost stable circular orbit (ISCO) of the EMRI, which determines the characteristic length scale of the system and thus the characteristic wavelength of the mediator, probes a dark force range that is complementary to that probed by SMBH or neutron star binaries. Moreover, since the ISCO scales linearly with the mass of the central SMBH, different EMRI systems are sensitive to different dark force ranges corresponding to λ ∼ R ISCO ∝ M 2 , where M 2 is the mass of the SMBH. In particular,  \nwe show our results for LISA [33, 34] (SMBH masses ∼ 10 6 M ⊙ with corresponding λ ∼ O (10 10 ) m and m V ∼ O (10 -16 ) eV) and for the proposed µ Ares detector [35] (SMBH masses ∼ 10 8 M ⊙ with corresponding λ ∼ O (10 12 ) mand m V ∼ O (10 -18 ) eV).  \nThe second reason EMRIs are interesting in this context is that mass measurements of the central SMBH, performed by methods that are independent of the primary EMRI system (a combination that we somewhat loosely call \'multimessenger") can enable us to probe small values of the dark sector coupling. The main challenge in probing small couplings ˜ α \' is that the correction in Eq. (1.1) is parametrically limited by the uncertainty in the binary component masses 2 . Independent mass measurements can break the degeneracy between gravity and the dark force, if the dark force is assumed to be subdominant or \'switched off" between the SMBH and ordinary matter, while being operational between the SMBH and the black hole comprising the EMRI system. Stated in another way, the central SMBH is envisioned to have ordinary matter in its vicinity, which allows for the determination of the mass of the SMBH, without the confounding effect of the dark force. This knowledge of the mass of the SMBH can then be used to constrain the dark force using the EMRI between the SMBH and stellar black hole. Although we will be somewhat agnostic to the underlying particle physics model that achieves this, we will indicate the most natural setting where our scenario may be applicable: a dark sector fermion χ with very feeble millicharge (near the weak gravity conjecture limit) interacting via an Abelian gauge boson V µ .  \nMeasuring the mass of SMBHs is a topic of intense research in astronomy. In quiescent galaxies, the methods typically rely on the motion of dynamical tracers - stars or gas - combined with sophisticated dynamical models. On the other hand, for active galaxies with prominent broad emission lines, the most promising method is reverberation mapping [37-39]. While measuring the motion of dynamical tracers requires high resolution and is applicable for nearer galaxies, reverberation mapping can be applied to galaxies at greater distances. We will also consider a third method that is applicable to nearby galaxies: GW signals from extremely large mass-ratio inspirals (XMRIs) formed by the central SMBH and a brown dwarf. We note that there are other emerging proposals for measuring the mass of SMBHs that we do not consider, such as X-ray flux variability of active galactic nuclei [40]. Our work motivates cross-pollination between the exciting field of SMBH mass determination and the future field of EMRI studies with future GW detectors, all within the context of dark sector physics.  \nThe rest of the paper is organized as follows. In Section 2, we discuss multimessenger measurements of the SMBH mass in different astrophysical systems. In Section 3, we show the dark force effect in the GW emission from EMRIs and the resulting sensitivity to the dark force with future LISA and µ Ares observations. We conclude in Section 4.', '2 EMRIs, Multimessenger Studies, and Ultralight Mediators': 'For ultralight mediators, the prospects of distinguishing the dark force from pure gravity is challenging. Compared to equal mass binaries, however, EMRIs offer certain advantages in addressing this important problem. The key point is the determination of the mass of the SMBH by a method for which the dark force can be switched off, followed by the determination of the mass of the mediator by an EMRI between the SMBH and a stellar black hole.', '2.1 Multimessenger Studies with XMRI-EMRI Combinations': 'The combination of an extremely large mass-ratio inspiral (XMRI) and an EMRI with the same central SMBH but two different smaller objects (a brown dwarf for the XMRI [41-43] and a black hole for the EMRI) constitutes our first example of a multimessenger signal. XMRIs have been recently explored as interesting sources for LISA or TianQin. With mass ratios as small as ∼ 10 -8 , a brown dwarf would complete ∼ 10 8 cycles before the merger, spending time on band for millions of years. These are expected to be abundant and loud systems: 10 4 ( 10 3 ) years before merger, the SNR is expected to be O (10 3 ) ( O (10 4 ) ) [43]. It should be noted that the high SNR of XMRI systems implies that they can be detected in nearby galaxies as well ([43] calculates that XMRI systems with SNR = 10 can be detected out to 50 Mpc 200 yrs before plunge). It is instructive to consider the case of Sgr A ⋆ . XMRIs between Sgr A ⋆ and brown dwarfs can constrain the mass of Sgr A ⋆ at the level of O (10 -2 ) M ⊙ [42] 3 . For the event rate of eccentric XMRIs, [42] obtain an average of N = 8 +9 -3 for the case of low spin of Sgr A ⋆ ( a = 0 . 1 ) and N = 12 +6 -4 for the case of high spin of Sgr A ⋆ ( a = 0 . 9 ). The numbers quoted above signify sources that emit GWs in the LISA and TianQin detection band with SNR ∼ 30 . The brown dwarfs are assumed to have mass 0 . 05 M ⊙ in prograde orbits and have an orbital inclination of 0 . 1 rad. The initial eccentricity of the XMRI system is taken to be between 0 . 99 and ≥ 0 . 999 with semimajor axis 4 × 10 -4 pc ≤ a ≤ 8 × 10 -3 pc. Similar numbers are reported by [43].  \nWe now turn to a discussion of EMRI formation. Generally, the rate of EMRI formation in a galaxy is determined by two timescales. The first (the two-body relaxation timescale T rlx ) is the time it takes to produce high eccentricity EMRI candidates from the dense cluster of stars and stellar remnants in the nuclear cluster surrounding the SMBH. The timescale T rlx is determined by the number density of compact objects in the nuclear cluster, their velocity dispersion, steepness profile, and half mass radius [44]. The second timescale pertains to the time it takes for the EMRI candidate to decouple from the rest of the nuclear cluster, and is determined by the rate of GW emission. Typical estimates for the relaxation time are in the range of O (0 . 1 -1) Gyr. Combining the estimates of these timescales, the benchmark figure for EMRI event rates from studies in the literature is ∼ 10 -6 per galaxy per year. This translates to tens to hundreds of EMRI events per year at LISA, which probes the Universe out to redshift z ∼ 1 [45]. It should be noted that the event rate for EMRIs is a subject of ongoing research, and the benchmark rate quoted above can be boosted by other mechanisms. While the ∼ 10 -6 per galaxy per year is based on the \'loss cone" scenario where a candidate is scattered by interactions with other bodies in the nucleus into a highly eccentric orbit around the SMBH, the presence of an accretion disk can significantly enhance this process. For example, [46] find that accretion disk-assisted EMRIs form ∼ O (10 1 -10 3 ) times faster than the \'loss cone\' mechanism.  \nThe possibility of a concurrent observation of an XMRI-EMRI combination at Sgr A ⋆ is pessimistic: while the expected event rate of an XMRI is relatively high, the event rate of an EMRI is more suppressed, and the joint probability is further suppressed. Nevertheless, given that XMRI signals may be detected up to 50 Mpc, and EMRIs even further, the possibility of a coincident XMRIEMRI signal may be significantly boosted if one considers other galaxies. Since such a concurrent multimessenger observation represents one of the best case scenarios for constraining ultralight mediators, we depict results corresponding to this benchmark scenario in our plots in Section 3.', '2.2 Multimessenger Studies with Stellar Kinematics - EMRI Combinations': 'We now turn to another possibility: multimessenger combinations from stellar kinematics around SMBHs on the one hand, and EMRIs of stellar black holes on the other. While the precision of mass measurements using stellar orbits is far lower than that using the XMRI-EMRI combination, the event rate of such a possibility is more optimistic.  \nFor the case of Sgr A ⋆ , the mass measurement by the GRAVITY Collaboration [47-49] relies on stellar orbits at several thousands of gravitational radii and yields a value of M 2 = (4 . 297 ± 0 . 012) × 10 6 M ⊙ . On the other hand, [50] provide six SMBH mass measurements at 3 σ confidence level from SINFONI observations of nearby fast rotating early-type galaxies: NGC 584 ( M 2 = (1 . 3 ± 0 . 5) × 10 8 M ⊙ ), NGC 2784 ( M 2 = (1 . 0 ± 0 . 6) × 10 8 M ⊙ ), NGC 3640 ( M 2 = (7 . 7 ± 5) × 10 7 M ⊙ ), NGC 4570 ( M 2 = (6 . 8 ± 2 . 0) × 10 7 M ⊙ ), NGC 4281 ( M 2 = (5 . 4 ± 0 . 8) × 10 8 M ⊙ ) and NGC 7049 ( M 2 = (3 . 2 ± 0 . 8) × 10 8 M ⊙ ). Other SMBH mass measurements are also available: the MASSIVE Survey [51] is a volume-limited, multi-wavelength, integral-field spectroscopic and photometric survey of O (100) of the most massive early-type galaxies within 108 Mpc. A sample of 72 SMBHs and their host galaxies based on kinematic data and modeling efforts is provided in [52]. We note that the general level of the uncertainty in the mass measurement using stellar kinematics is ∼ O (10 -30)% . Thus, the observation of an EMRI between the SMBH and a stellar black hole in any of these galaxies would result in multimessenger combinations of stellar kinematics-EMRI.', '2.3 Multimessenger Studies of Spectroscopic Reverberation Mapping and EMRIs': 'The third multimessenger combination scenario that we will consider is the measurement of the SMBH mass with spectroscopic reverberation mapping. In this method, the time delay in changes in the continuum emission arising from the accretion disk near the SMBH and the response of these changes on the emission from the (farther away) broad line region are measured. This time delay effectively measures the distance between the accretion disk and the broad line region (more precisely, the responsivity-weighted average radius of the broad line region). Finally, from this information, the mass of the SMBH can be extracted following [53]:  \nM 2 = fcτ cent σ 2 line G . (2.1)  \nHere, f ∼ O (1) is a parameter that is specific to the broad line region and depends on its structure, orientation, and kinematics. The parameter τ cent is the measured time delay (as measured from the centroid of the cross-correlation function around the emission peak), while σ line denotes the rootmean square width of emission in the variable part of the spectrum. A major source of uncertainty is the value of f . Mass measurements of SMBHs using this method can be obtained from the AGN Black Hole Mass Database [54]. Most of the measured SMBH masses lie at redshifts z < 0 . 2 , although there are a few at z ∼ O (2 -3) . The general level of uncertainty in the mass measurement is ∼ O (10%) .  \nIn this context, we take particular note of the OzDES and SDSS reverberation Mapping Programs [39, 55, 56], which have pioneered large-scale reverberation mapping campaigns [55]. These programs have observed hundreds of active galactic nuclei, probing a wide range luminosities and redshifts. The OzDES-RM program recently published data from successful campaigns for eight active galactic nuclei between 0 . 12 < z < 0 . 71 [55]. In previous work, they measured the mass of two SMBHs at some of the highest redshifts, z = 1 . 905 and z = 2 . 593 . The measured masses are 4 . 4 +2 . 0 -1 . 9 × 10 9 M ⊙ for DES J0228-04 and 3 . 3 +1 . 1 -1 . 2 × 10 9 M ⊙ for DES J0033-42, calculated with f = 4 . 47 ± 1 . 25 . The SDSS-RM Program published results for 849 broad line quasars between  \n0 . 1 < z < 4 . 5 . The observation of an EMRI between the SMBH and a stellar black hole in any candidate observed by these campaigns would result in multimessenger combinations of reverberation mapping - EMRI, which constitutes our third benchmark scenario.', '3 Results': "In this section, we obtain constraints on the mediator mass versus coupling plane ( m V , ˜ α ' ) , assuming that the central SMBH and a second stellar black hole constitute an EMRI system. Our strategy will be to compute the temporal evolution of the angular frequency ˙ ω , including the correction due to the dark force in Eq. (1.1), and subsequently compute the GWs from perturbations of the Kerr metric following [57] (we refer to [58, 59] for classic papers and [60, 61] for earlier work). We will then compare our results to the case of pure gravity, where the dark force is absent. Our calculation most closely resembles that of [11], who performed an analogous calculation for regular equal-mass binary inspirals. We will assume an uncertainty of O (10)% in the mass of the SMBH, corresponding to the case where the SMBH mass measurements comes from either stellar kinematics or reverberation mapping. On the other hand, for the case of an XMRI-EMRI multimessenger combination, we will assume that the mass of the SMBH is determined up to negligible fractional uncertainty ∼ O (10 -3 )% .  \nWenote that several prior studies of dark sectors with binary inspirals have concentrated on constraining either the strength of the dipole radiation term, or the charge of the remnant. For example, [11, 28, 62-64] consider black hole binary inspirals with massless mediators and constrain the dipole radiation term from the early pre-merger inspiral stage. The charge of the merger remnant can be constrained from the detection of the quasinormal modes. Dipole radiation in the context of EMRIs has been considered by [65-67]. Our focus, as stated before, will instead be on constraining the dark force in the limit of vanishing dipole radiation, relying instead on independent mass measurements of the SMBH.  \nThe use of independent mass measurements to constrain dark forces works under the assumption that the dark force is only operational between the SMBH and the stellar black hole, while being switched off between the SMBH and surrounding baryonic matter, including stars, other objects such as brown dwarfs, and matter in the accretion disk. This is possible in models where the black holes are charged under a dark force acting only between particles in the dark sector, and not between the dark and visible sectors. We discuss a simple model that achieves this as a proof of principle, in Appendix A. The model consists of millicharged dark fermions with mass ∼ 1 GeV and charge approximately equal to their mass, near the weak gravity conjecture limit. The dark force carrier is a vector V . The discharge time of a black hole initially charged with such dark fermions can be shown to be large under the most conservative assumptions. We remain somewhat agnostic about the question of how the black holes acquire their initial distribution of dark fermions, indicating some possibilities in Appendix A. For the case of the SMBH, the question reverts to the open question of how seeds of such black holes are formed.", '3.1 Effect of Dark Forces': "To set our notation, we consider an EMRI system consisting of two masses M 1 and M 2 , with the mass ratio, η ≡ M 1 /M 2 ≪ 1 . The heavy object ( M 2 ) is a SMBH. In addition to the gravitational attraction between the two objects, we introduce a dark force operating between them. For simplicity, we consider a Yukawa force mediated by a mediator particle V  \n| ⃗ F Yukawa | = α ' Q 1 Q 2 r 2 e -m V r (1 + m V r ) . (3.1)  \nHere α ' is the dark fine structure constant, m V is the mass of the dark force mediator, r is the separation between the two objects constituting the EMRI system, and Q i is the dark charge of the object i , with i = 1 , 2 . We will consider both the case when α ' > 0 (attractive force between charges of opposite signs) as well as α ' < 0 (repulsive force between charges of the same sign) 4 . Although the Yukawa force is typically used to describe interactions mediator by scalar particles, we do not distinguish between scalar and vector mediators in this study when introducing dark forces in the form of Eq. (3.1). The dark force, together with the gravitational attraction, gives the total force acting on each object  \n| ⃗ F total | = GM 1 M 2 r 2 [1 + ˜ α ' e -m V r (1 + m V r )] , (3.2)  \nwhere it is convenient to normalize the strength of the dark force to the gravitational force by introducing  \n˜ α ' ≡ α ' Q 1 Q 2 GM 1 M 2 . (3.3)  \nThe introduction of the dark force modifies the EMRI system in two ways. Firstly, the angular frequency ω ( r ) is determined by the net force on the lighter object, which is modified by the dark force, resulting in a departure from the inverse-square law as given in Eq. (3.2). Secondly, dark (dipole) radiation in the form of emission of the light force mediator will make the inspiral faster, affecting the temporal dependence of the separation r ( t ) . Both ω ( t ) and r ( t ) , and therefore ultimately the GW frequency f GW ( t ) and strain h c ( t ) are modified from the pure gravity case in the presence of a dark force. In this work, we will assume that the dipole radiation is negligible.  \nGiven a centripetal force in Eq. (3.2), the relation between the orbital frequency and the separation r in the presence of a dark force can be obtained as  \nω 2 = G ( M 1 + M 2 ) r 3 (1 + ˜ α ' e -m V r (1 + m V r )) . (3.4)  \nThe GW frequency of a harmonic m is determined with the orbital angular frequency  \nf GW ,m = m 2 π ω. (3.5)  \nNear the ISCO, GW radiation drives orbits to be circular; therefore the m = 2 mode dominates and we focus on f GW ≡ f GW , 2 in the rest of the paper. We determine r ( ω ) numerically from Eq. (3.4). The EMRI losses energy mainly through GW radiation when dipole radiation of the mediator is turned off. The GW emission rate is  \ndE GW dt = 32 5 Gµ 2 ω 6 r 4 . (3.6)  \nHere, µ = M 1 M 2 / ( M 1 + M 2 ) is the reduced mass.  \nThe time evolution of the orbital frequency due to GW radiation is [10]  \ndω dt = -32 5 Gµω 5 r 2 g N -1 . (3.7)  \nThe g factor is from the dark force correction to the gravitational attraction,  \ng = -3 + ˜ α ' e -m V r (3 + m V r (3 + m V r )) 1 + ˜ α ' e -m V r (1 + m V r (1 -m V r )) . (3.8)  \nWhile a precise determination of the gravitational wave emission from the EMRI systems remains an ongoing effort, several sets of waveforms have been provided for data analysis. We adopt the result provided by Ref. [57], which was obtained by numerically solving the Teukolsky equation [58, 68] for quasi-circular orbit cases. The result for each physical quantity is presented as the usual equation but with an additional factor added which captures the relativistic correction. In the case of GW frequency evolution in Eq. (3.7), we include the correction N -1 obtained from Table VIII of [57].  \nThe strain amplitude of the GW signal can be calculated from the gravitational energy radiation rate. For the characteristic strain in the harmonic m , we follow the definition of [57],  \nh c,m ≡ h o,m √ 2 f 2 GW ,m / ˙ f GW ,m , (3.9)  \nwith the dominant contribution coming from the harmonic m = 2 . The f 2 GW / ˙ f GW factor accounts for the time the GWs are produced at a fixed frequency for detection. The rms amplitude is defined as  \nh o,m = 2 √ ˙ E ∞ m mωd L , (3.10)  \nwhere ˙ E ∞ m is the total energy emission rate into the m -th harmonic to infinity. d L is the luminosity distance from the EMRI system to Earth. The full expression for h o,m appears in Eq. (3.11b) of [57]. We reproduce the expression for the m = 2 mode here for completeness  \nh o, 2 = √ 32 5 η G M 2 d L ( GM 2 ω ) 2 / 3 H o, 2 . (3.11)  \nThe relativistic correction to the strain H o, 2 can be obtained from the relativistic correction to the energy loss rate H o, 2 = √ ˙ E ∞ 2 , and we use numerical values of ˙ E ∞ 2 reported in Table IV of [57]. We will focus on the m = 2 harmonic of the characteristic strain amplitude h c ≡ h c, 2 , which is calculated through Eq. (3.9) with numerical results of ˙ f GW from Eq. (3.7) and h o, 2 from Eq. (3.11).", '3.2 Temporal Evolution of Gravitational Waves': "We use formulas introduced in the last section to calculate GWs emitted by the EMRI system. Fig. 1 and Fig. 2 illustrate the evolution of the GW frequency f GW and characteristic strain h c for different choices of the dark force strength parameter ˜ α ' . The mass of the SMBH is chosen as M 2 = 10 6 M ⊙ , and the mass of the inspiralling black hole is M 1 = 10 M ⊙ . We assume for simplicity both black holes are Schwarzschild black holes, though dark forces can also impact EMRIs involving rotating black holes. The mass of the mediator is m V = (10 7 km) -1 ≃ 1 . 98 × 10 -17 eV , such that the force range is comparable to the ISCO radius of 8 . 9 × 10 6 km for a 10 6 M ⊙ SMBH. The associated GW frequency emitted at the ISCO is about 4 . 4 mHz . We evolve the EMRI system using Eq. (3.7) from the initial condition f GW = 1 . 34 mHz at t = 0 , and calculate the GW emission power with Eq. (3.9). The corresponding initial EMRI components separation is about 2 × 10 7 km with small difference depending on the type of interactions.  \nThe colored curves in Fig. 1 show f GW as a function of time in the cases with the absolute value of the dark force strength | ˜ α ' | = 0 . 5 (purple), 0 . 25 (red), and 0 . 1 (blue). The attractive force case  \n<!-- image -->  \nFigure 1 : Evolution of the EMRI GW frequency f GW of a 10 M ⊙ object inspiraling a 10 6 M ⊙ SMBH for a range of dark force strength values ˜ α ' in the attractive case (solid) and repulsive case (dashed). The mediator mass is (10 7 km) -1 . ( left ): the case of multimessenger EMRI signals with either reverberation mapping or stellar kinematics data, resulting in an uncertainty (displayed by the gray band) of 10% in determining both of the masses M 1 and M 2 . ( right ): the case of multimessenger EMRI-XMRI signals, resulting in a fractional uncertainty (displayed by the gray band) of 10 -3 % in the SMBH mass M 2 and 10% in M 1 .  \n<!-- image -->  \n<!-- image -->  \nFigure 2 : The temporal dependence of the characteristic strain h c of the EMRI system for a range of dark force strength values in the attractive (solid) and repulsive (dashed) cases. The black hole and dark force mediator masses are chosen to be consistent with Fig. 1. ( left ): multimessenger measurement with reverberation mapping. The uncertainty of masses M 1 and M 2 are both 10% . ( right ): multimessenger measurement with XMRI. The uncertainty of masses M 1 and M 2 are 10% and 10 -3 % respectively.  \n<!-- image -->  \n( ˜ α ' > 0 ) is shown with solid curves while the repulsive force case ( ˜ α ' < 0 ) is shown with dashed curves. As a comparison, we show the case in the absence of dark force with the black curve, together with a gray band representing the typical uncertainty from the determination of the black hole masses. In the left panel, we show the uncertainty of ∆ M 2 /M 2 = 10% corresponding to reverberation mapping campaigns. In the right panel, the uncertainty is ∆ M 2 /M 2 = 10 -3 % corresponding to  \ng  \nFigure 3 : The coupling g defined in Eq. (3.8) as a function of the dimensionless force range m V r with different choices of ˜ α ' for attractive dark force (solid) and repulsive dark force (dashed). The largest deviations from the gravitation only case, g = -3 , occur when m V r ∼ 1 .  \n<!-- image -->  \nthe mass determination with an XMRI measurement. In both panels, a 10% uncertainty for the mass of the inspiraling black hole ∆ M 1 /M 1 = 10% is incorporated within the gray bands. Evidently, the attractive dark force accelerates the system towards ISCO, while the repulsive dark force slows down the system's evolution towards ISCO. Larger values of the effective dark force strength induce more significant deviations from the gravitation-only scenario. With a more precise multimessenger determination of M 2 in the right panel, a smaller | ˜ α ' | = 0 . 1 can be distinguished from the uncertainty band compared to the sensitivity of about | ˜ α ' | = 0 . 25 in the left panel.  \nIn Fig. 2, we show the strain amplitude h c as a function of time for the same EMRI event consisting M 1 = 10 M ⊙ and M 2 = 10 6 M ⊙ located at d L = 1 Gpc , in the attractive dark force case (solid) and the repulsive dark force case (dashed). The inverse mediator mass is chosen to be 1 × 10 7 km . In all figures, the pure gravity case is shown with the black curve, and the benchmarks with non-zero dark force couplings are shown for | ˜ α ' | = 0 . 5 (purple), | ˜ α ' | = 0 . 25 (red), and | ˜ α ' | = 0 . 1 (blue), corresponding to the same benchmark values in Fig. 1. The SMBH mass uncertainty is shown with the gray band for 10% in the left panel, and for 10 -3 % in the right panel, while the mass uncertainty of the lighter component is kept to be 10% in both panels. We note that for an attractive (repulsive) dark force, the characteristic strain h c is weaker (stronger) than the pure gravity case at a given time t > 0 . This contrasts with the evolution of rms amplitude h o, 2 of GWs emitted toward infinity. The dependence is reversed as a result of the effect of the dark force on the GW frequency: when the EMRI system evolves faster to ISCO, it spends less time f GW / ˙ f GW in the fixed observation frequency window, such that the h c in Eq (3.9) is suppressed. Similar to Fig. 1, we see that the effect of the dark force can be discerned best with a stronger dark force strength | ˜ α ' | .  \nThe strain amplitude of a resolved EMRI signal is inversely proportional to the luminosity distance. Therefore, EMRI events that are closer can be measured with greater precision; indeed, as mentioned before, multimessenger studies involving EMRIs rely on stellar kinematics or reverberation mapping that are only effective at much smaller distances. Nevertheless, d L = 1 Gpc represents a conservative choice.  \nIn addition to the dark force strength, EMRIs are also sensitive to the mass of the mediator m V , which determines the range of the dark force. The dependence on the mediator mass is from the g factor in Eq. (3.8). In the scenario where there is no dark force, as well as in the scenario where  \n]  \nz  \nH  \nm  \n[  \nW  \nG  \nf  \n3.5  \n3.0  \n2.5  \n2.0  \n1.5  \n0  \n2  \n4  \n6  \n8  \nM  \n1  \n=  \n10  \nM  \nα  \n'  \n=  \nmV  \nmV  \nmV  \n0  \n=(  \n=(  \n=(  \n|  \nα  \n'  \n| =  \n☉  \n,  \nM  \n×  \n5  \n×  \n1  \n×  \n1  \n0.5  \n2  \n6  \n10  \n7  \n10  \n9  \n10  \n=  \n10  \nkm  \nkm  \nkm  \n)  \n)  \n)  \n6  \nM  \n-  \n-  \n-  \n1  \n1  \n1  \n☉  \nt  \n[  \nyears  \n]  \nFigure 4 : GW frequency versus time plot of a 10 M ⊙ object inspiraling a 10 6 M ⊙ SMBH for a range of the dark force mediator mass m V . ( left ): the case of multimessenger EMRI signals with reverberation mapping data, resulting in an uncertainty of 10% in determining each of the masses M 1 and M 2 . ( right ): the case of multimessenger EMRI-XMRI signals, resulting in an uncertainty of 10 -3 % in the SMBH mass M 2 and 10% in M 1 .  \n<!-- image -->  \nthe range of the dark force is much longer or shorter than the EMRI separation, Eq. (3.8) reduces to g = -3 . Conversely, the correction factor deviates from -3 when m V is comparable to the inverse mutual separation of the EMRI components. We show the dependence of the g factor on the relative values of m V and r in Fig. 3, for the case of attractive (denoted by solid lines) as well as repulsive (denoted by dashed lines) dark force interactions. It is clear that significant deviations from the gravity only case occur when m V r ∼ 1 . The deviation of the g from -3 introduces an effect on the EMRI evolution that agrees with the attractive or repulsive nature of the dark force. Specifically, the system evolves more quickly for an attractive dark force interaction when g is more negative around m V r ∼ 1 , while it evolves more slowly for a repulsive dark force interaction when g is less negative around m V r ∼ 1 .  \nThe effect of the mediator mass on the EMRI is displayed in Fig. 4 with f GW as a function of time for a fixed dark force strength but different mediator masses. We choose ˜ α ' = 0 . 5 in the left panel, and ˜ α ' = 0 . 2 for the right panel. The other EMRI parameters are kept to be consistent with Fig. 1. We also choose the initial frequency to be f GW = 1 . 34 mHz for all EMRIs. The corresponding initial black hole separation is about 2 × 10 7 km with small difference depending on the type of interactions. As a comparison, the ISCO radius is about 8 . 9 × 10 6 km for M 2 = 10 6 M ⊙ . The mediator mass values are chosen to represent three different force range scenarios: short force range m -1 V = 5 × 10 6 km (red); intermediate force range comparable to the ISCO radius m -1 V = 10 7 km ∼ R ISCO (yellow); and the very long-range force limit m -1 V = 10 9 km (blue). As is expected from the dependence of the g factor on m V r , the dark force effect of yellow curves exhibit more significant distinction from the gravity-only case (black), compared to the short-range and long-range force limits.  \nIn the regime where m V ≪ R -1 ISCO , EMRI GWs can still discern between the gravitational attraction and an additional dark force if one advocates, as we have, for an independent mass measurement for which the dark force is switched off. This is different from the case of GW measurements from equal-mass binary systems, where the effect of the dark force becomes entirely degenerate with gravity: in that case, a rescaling of the binary component mass by a factor of (1 + ˜ α ' ) 2 / 5 mimics the presence of the dark force. The primary challenge in distinguishing the effect of the dark force stems  \ng  \ni  \nn  \np  \na  \np  \nM  \nn  \ni  \no  \nt  \nr  \na  \ne  \nr  \nb  \ne  \ne  \nv  \nR  \n/  \nΔ  \n,  \n0  \n%  \n1  \n=  \n1  \nM  \n1  \nΔ  \nM  \n/  \n0  \n=  \n%  \n1  \n2  \nM  \n2  \nM  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 5 : Strain h c as a function of the GW frequency f GW near ISCO. Color curves show the EMRI evolution with a range of dark force strengths. Black curve shows the evolution in a gravitation-only scenario. Brown dots indicate the evolution time is 5 years from ISCO while cyan dots represents 1 year left to ISCO. ( upper ): Dark force effect compared to the gravitation-only case, with mass uncertainties ∆ M 1 /M 1 = 10% and ∆ M 2 /M 2 = 10% , depicted on the left, and zoom-in of specific regions on the right. ( lower ): Dark force effect compared to the gravitation-only case with ∆ M 1 /M 1 = 10% and ∆ M 2 /M 2 = 10 -3 % on the left, and zoom-in plot on the right.  \n<!-- image -->  \nfrom uncertainties in the pure gravitational contribution, which ultimately arises due to uncertainties in M 1 and M 2 (illustrated by the gray band). For instance, for a benchmark value of ˜ α ' = 0 . 5 , the rescaling factor on the mass, (1 + ˜ α ' ) 2 / 5 ≃ 1 . 18 , exceeds the uncertainty in the reverberation mapping and the XMRI measurements. Consequently, the degeneracy is broken for EMRI GWs in the ultra-long-range force limit.  \nIn the end, we show the EMRI GWs on the plane of strain versus frequency in Fig. 5, with the black and cyan dots indicating the remaining time from ISCO is 5 years and 1 year respectively. The mediator mass is set to m -1 V = 10 7 km to demonstrate the maximum effect. The upper panels illustrate the effects of the dark force compared to uncertainties determined by reverberation mapping, with a detailed zoom-in on the region where the dark force leads to a significant deviation from the gravity-only scenario. The lower panels show the case when uncertainties are determined by the XMRI measurement. The GW data collected over the last few years from the EMRI events in future  \nobservations will be utilized to precisely track the evolution of the strain versus the frequency, crucial for distinguishing the existence of the dark force. We will use the h c -f GW relation to calculate the sensitivity to the dark force in the next section.", '3.3 Sensitivity with Multimessenger Measurements': "We perform a likelihood analysis to test the future sensitivity to ˜ α ' with the EMRI GWs and the additional independent measurement of the black hole masses. The likelihood analysis is based on the strain-frequency data of the GW signal. The likelihood is defined as [69]  \nΛ( s | ˜ α ' ) = K exp[( h ˜ α ' | s ) -1 2 ( h ˜ α ' | h ˜ α ' ) -1 2 ( s | s )] (3.12)  \nHere K is a normalization constant, s ( t ) is the observed GW waveform in the time domain, and h ˜ α ' ( t ) is the waveform corresponding to a selected value of ˜ α ' . The scalar product between the two waveform is given by  \n( h ˜ α ' | s ) = 4 Re ∫ f max f min df GW ˜ h ˜ α ' ( f GW ) ˜ s ( f GW ) S n ( f GW ) (3.13)  \nwhere S n ( f GW ) represents the noise spectral density of the GW observatory and has a dimension of Hz -1 . Quantities denoted with a tilde represent the Fourier transform of the time domain data. Following the definitions in [70], we can write the scalar product as  \n( h ˜ α ' | s ) = Re ∫ f max f min df GW h c, ˜ α ' ( f GW ) s c ( f GW ) f 2 GW S n ( f GW ) (3.14)  \nThe characteristic strain h c, ˜ α ' ( f GW ) is calculated with the chosen dark force strength. The s c ( f GW ) is obtained from the observed GW waveform. In the EMRI measurement, the upper limit f max and the lower limit f min of the integral in Eq. (3.14) is determined by the frequency range over which the waveform is observed, with the upper limit set to f ISCO GW and the lower limit set to the EMRI GW frequency at 10 years before ISCO 5 . As long as the EMRI signal is stronger than the detector noise, the likelihood is mostly determined by the deviation of the signal spectrum in the presence of the dark force from the gravity-only spectrum (with mass uncertainties). We calculate the value of ˜ α ' for a given m V that maximizes the likelihood. This is equivalent to maximizing the log-likelihood, which we can read from Eq. (3.12) as  \nlog Λ( s | ˜ α ' ) = log K +( h ˜ α ' | s ) -1 2 ( h ˜ α ' | h ˜ α ' ) -1 2 ( s | s ) = ( h ˜ α ' | s ) -1 2 ( h ˜ α ' | h ˜ α ' ) (3.15)  \nFor an observed waveform s ( t ) , the term log K 1 2 ( s | s ) is a constant independent of the dark force effect and can be omitted from the definition of log-likelihood in the rest of the analysis.  \nTo estimate the sensitivity to ˜ α ' , we assume that the observed signal corresponds to the dark force scenario. The likelihood of a gravity-only scenario is calculated with the assumed uncertainties in the masses of the EMRI components and ˜ α ' = 0 . A scan across the parameter space of ˜ α ' is performed to assess the likelihood with non-vanishing ˜ α ' . We further assume that the test statistic,  \ndenoted as τ , is χ 2 -distributed under the null hypothesis. The two likelihoods are then compared at a 95% confidence level. The likelihood test statistics is  \nτ = -2 ( ln Λ( s | 0) -max ˜ α ≥ 0 [ln Λ( s | ˜ α )] ) (3.16)  \nWe will use the likelihood test statistics to calculate the sensitivity to the dark force strength and the force range. While the precision of GW measurements improves the sensitivity to ˜ α ' , the sensitivity to the dark force range depends on the mass of the SMBH involved in the EMRI event. We discuss the relation between the dark force range and the SMBH mass in the following. EMRI GWs are observed mostly when the orbits of the inspiraling object have decayed close to the ISCO radius, which can be written as R ISCO = 6 GM 2 for Schwarzschild SMBHs. The frequency of GWs emitted when the separation is equal to the ISCO radius, r = R ISCO , is  \nf ISCO GW = 1 6 √ 6 π GM 2 ≃ 4 . 35 ( 10 6 M ⊙ M 2 ) mHz . (3.17)  \nThis determines the observed GW frequency for EMRIs with different center SMBH masses. The frequency window of LISA is ideal to be used for observing GWs from ∼ 10 6 M ⊙ SMBHs. On the other hand, EMRI GWs from heavier SMBHs appear at lower frequencies that are outside the LISA range. For example when M 2 ∼ 10 9 M ⊙ , the GW frequency lies in the unexplored µ Hz range. Future observations of µ Hz GWs have been proposed with space-based interferometers µ Ares [35], coldatom interferometers [71], and using Solar System asteroids as test masses [72]. With a larger SMBH mass and consequently a wider ISCO radius, we can probe the dark force of a lighter mediator mass m V , which would be indistinguishable from a rescaling of gravitational force with measurements at smaller separations. Therefore, the mediator mass of optimal sensitivity, the SMBH mass and the GWfrequency can be correlated as follows,  \nm -1 V ≃ 8 . 96 × 10 8 ( M 2 10 8 M ⊙ ) km ≃ 3 . 90 × 10 8 ( 10 -4 Hz f GW ) km (3.18)  \nWe present the sensitivity to ˜ α ' across a range of the mediator mass m V in Fig. 6. In particular, we choose two sets of representative GW event parameters: (1) M 2 = 10 6 M ⊙ , d L = 10 Mpc for lighter SMBHs; and (2) M 2 = 10 8 M ⊙ , d L = 0 . 1 Gpc for heavier SMBHs, to demonstrate their sensitivity to different mediator masses. In both benchmarks, the results are calculated by requiring τ = 4 in Eq. (3.16) with different multimessenger measurements. In both cases, we will show that the optimal sensitivity is achieved when the mediator mass approximates the inverse of the EMRI separation, as can be seen from the minimum of the sensitivity curve at m V r ∼ 1 .  \nIn the first benchmark, the SMBH mass uncertainty can be independently constrained to within 10% by spectroscopic reverberation mapping or stellar kinematics, with which the sensitivity of EMRI GW measurement is shown with the red dotted curve in Fig. 6. Additionally, the SMBH mass can be measured within 10 -3 % uncertainty with XMRIs. The corresponding sensitivity is shown with the cyan dotted curve in Fig. 6. The EMRI GWs from ∼ 10 6 M ⊙ SMBHs will be measured precisely with the future LISA observatory. We use noise power spectral densities reported in [29] for the LISA instrumental noise and the galactic confusion noise to obtain the sensitivity. We assume M 1 = 10 M ⊙ with 10% uncertainty. The results show that LISA has optimized sensitivity to m V ∼ 10 -17 -10 -16 eV , while it can still detect effects from longer mediator wavelengths. To illustrate sensitivities achieved through further reduction of uncertainties associated to the black hole masses, we present the sensitivity corresponding to 10 -3 %precision for both M 1 and M 2 in the black dashed curve.  \nFigure 6 : Projected 95% C.L. sensitivity of LISA and µ Ares to an attractive dark force on its strength ˜ α ' and force mediator mass m V . The cyan, red, and black curves assume the EMRI system consists of a stellar black hole of mass M 1 = 10 M ⊙ and a SMBH of mass M 2 = 10 6 M ⊙ . The cyan dotted curve corresponds to the case of EMRI-XMRI signals, resulting in an uncertainty of 10% in determining the value of M 1 , and 10 -3 % in determining the value of M 2 . The red dotted curve displays the sensitivity of EMRI signals complemented with spectroscopic reverberation mapping measurements, leading to a 10% uncertainty in both M 1 and M 2 . The black dashed curve depicts the sensitivity of an EMRI waveform observed with an uncertainty of 10 -3 %in determining both masses. We also show the case of a heavier SMBH mass with the blue dashed curve for an IMRI-XMRI multimessenger measurement. The black hole masses are M 1 = 10 4 M ⊙ and M 2 = 10 8 M ⊙ , with 10% uncertainty in M 1 and 10 -3 % uncertainty in M 2 .  \n<!-- image -->  \nIn the second benchmark, we assume the sensitivity is set through GWs measured by the future µ Ares observation, combined with an independent SMBH mass measurement ∆ M 2 /M 2 = 10 -3 % from the XMRI. µ Ares is capable of detecting GWs with a sensitivity of h c ≈ 10 -18 at f GW = 10 -6 Hz . With such sensitivities, it will be able to observe Intermediate Mass Ratio Inspirals (IMRIs) involving 10 3 -10 4 M ⊙ intermediate-mass black holes inspiralling onto 10 7 -10 8 M ⊙ SMBHs out to z ≥ 6 [35]. An IMRI system will respond to the dark force similarly to an EMRI, and the GW emission will be altered in a comparable manner. According to Eq. (3.18), an IMRI signal from a 10 8 M ⊙ SMBH is able to probe m -1 V ≳ 10 9 km . We assume M 1 = 10 4 M ⊙ for the intermediatemass black hole, and use the µ Ares noise level reported in [35] to calculate its sensitivity. The final sensitivity is depicted by the blue dashed curve in Fig. 6. In the regions where the blue curve overlaps with the red and cyan curves, the dark force effect can be measured for both SMBHs of different masses, allowing the parametric dependence on ˜ α ' and m V to be disentangled.  \nOur study concentrates on the role of independent mass measurements; however, previous studies employing Fisher information matrix analysis demonstrated that uncertainties in EMRI parameters could potentially be reduced to an order of magnitude of O (10 -4 ) [73], benefiting from the information in the GW waveform (see also applications to dark force effects in neutron star binaries [11]). In order to estimate the prospective sensitivity of an EMRI observation to the strength of the dark force utilizing the maximum likelihood criterion, we use the mass determination uncertainty of 10 -3 % to demonstrate the improvement of the sensitivity with both the multimessenger measure-  \nments and the dedicated waveform analysis. While a comprehensive waveform analysis that takes into account the additional dark force coupling within the model parameters is outside the scope of this study, our findings indicate that an independent measurement of the SMBH mass-one not influenced by the dark force- could resolve the degeneracy between the EMRI model parameters.", '4 Conclusions': 'Our focus in this paper has been to probe dark forces using EMRI systems. In contrast to other studies using compact object mergers, we have stressed on the importance of independent mass measurements of the central SMBH that do not depend on the interaction between the SMBH and stellar black hole partner of the EMRI system. Such independent mass measurements break the degeneracy of the effect of the dark force with a simple rescaling of the binary component masses. We have discussed two classes of such independent mass measurements, based on multimessenger studies involving other gravitationally bound systems and on photons. In the first class are XMRI systems involving the SMBH and a smaller object, such as a brown dwarf. In the second class are various existing methods of determining the SMBH mass, such as stellar kinematics and reverberation mapping campaigns. The SMBH mass determination in the first case can be very precise, with an uncertainty of ∼ 10 -3 % , while a standard benchmark of the mass uncertainty using reverberation mapping is ∼ 10% . We have depicted both benchmarks in our results.  \nThe question of localization of EMRIs using future GW detectors is important for our work. The LISA-TianQin or TaiJi-TianQin networks have the possibility of localizing GW events; such localization is frequency-dependent and has been estimated to be 200 , 3 , and 0 . 005 deg 2 at 1 , 10 , and 100 mHz, respectively [74]. A typical benchmark for the localization of EMRIs is ≲ 10 deg 2 [30]. Multimessenger observations involving electromagnetic follow-ups, with the accretion disk or the broad line region, would further localize EMRIs (their potential as standard candles has been discussed by [75]).  \nIn this study, we have assumed that the corrections from environmental effects are negligible. SMBHsareknown to exist in dense astrophysical environments, which may influence the dynamics of EMRIs [76-78]. Key processes such as gravitational pull from matter, accretion onto EMRI components, dynamical friction, and planetary migration can impact the gravitational waveform, potentially leading to biased parameter estimation and reduced detectability of EMRIs. We direct readers to [79] for information on the maximum level of corrections due to various environmental effects over a typical LISA mission. The environmental corrections are expected to be minimal for SMBHs surrounded by thick disks, but more substantial for those in thin disks, which are predominantly found in active galactic nuclei at high redshifts. Although SMBHs in thin disks contribute up to only a few percent of the EMRIs detectable by LISA with z ≲ 1 , potential uncertainties from environmental effects should also be considered when applying our method to these events.  \nOur results motivate multimessenger studies with coordination between GW detectors and reverberation mapping campaigns. In particular, SMBH mass determination by the OzDES-RM and SDSS-RM programs, and the AGN Black Hole Mass Database serve as rich targets for future EMRI studies. The upcoming SDSS-V Black Hole Mapper Reverberation Mapping Project is also important in this regard [80]. Similarly, SMBH mass determination by the SINFONI and MASSIVE Surveys using dynamical tracers also serve as possible targets for future EMRI sources. The best case scenario is the concurrent observation of an EMRI and an XMRI for a given SMBH.', 'Acknowledgments': 'The authors would like to thank Huai-Ke Guo and Chen Sun for collaboration in the early stages of the project. KS would like to thank the Aspen Center for Physics for hospitality during the course of this work. TX would like to thank the Center for Theoretical Underground Physics and Related Areas (CETUP*), The Institute for Underground Science at Sanford Underground Research Facility (SURF), and the South Dakota Science and Technology Authority for hospitality during a period part of this work was completed. The work of KS and TX is supported by the U.S. Department of Energy Grant DE-SC0009956.', 'A Dark Charged Black Holes': "In this appendix, we discuss the existence of black holes with dark charge. Several mechanisms of endowing black holes with long-lived charge content are mentioned by the authors of [9, 81] one interesting possibility is to introduce millicharged fermions in the dark sector [28]. Denoting the mass of the SMBH by M 2 and its charge by Q 2 , one can obtain Q 2 /M 2 ∼ 1 by choosing the charge-to-mass ratio of dark fermions to be ∼ 1 (close to the limit imposed by the weak gravity conjecture). In this case, one could consider a SMBH-brown dwarf XMRI system (where the brown dwarf is the uncharged Object A alluded to above) which determines the mass of the SMBH; concurrently, a SMBH-black hole EMRI system (where the stellar black hole is the dark-charged Object B) constrains the light mediator mass. We outline the model below.  \nOur model will consist of a millicharged dark fermion χ with charge q χ and mass m χ , interacting by the exchange of a vector mediator V µ . The Lagrangian is  \nL = -1 4 V µν V µν + 1 2 m V V µ V µ + ¯ χ ( iγ µ D µ -m χ ) χ . (A.1)  \nHere, the covariant derivative is given by D µ = ∇ µ + iq χ V µ , and we normalize the dark fine structure constant α ' = 1 6 . The total charge of the SMBH can then be written as  \nQ 2 = N 2 q χ , (A.2)  \nwhere N 2 gives the number of fermions of opposite charge it needs to accrete to become neutral. We now discuss in turn estimates of the discharge time for a dark charged black hole, and then how it could acquire dark charge in the first place. These issues all hinge on the possibility that a charge q χ with mass m χ is able to be absorbed by a black hole with charge Q and mass M . A simple (classical) comparison of the associated gravitational attraction and possible dark repulsion yields the condition  \nq χ Q ∼ m χ M ⇒ Q M ∼ m χ q χ . (A.3)  \nIn the Standard Model, this condition is badly violated unless Q ∼ 0 . However, for a dark sector with  \nm χ ∼ m p , and q χ ∼ 10 -18 e (A.4)  \nwhere e is the charge of the electron and m p is the mass of a proton, one can obtain Q ∼ M .  \nAssuming that the SMBH is in an environment with density ρ of oppositely charged dark fermions moving with velocity v , and that the accretion follows the Bondi rate for collisionless fluid, [28] estimates the discharge time to be  \nτ discharge ∼ ( v 220 km / s ) 3 ( 0 . 4 GeV / cm 3 ρ ) ( 10 M ⊙ M 2 )( m χ m p )( e q χ ) yr . (A.5)  \nWith the choice of dark sector parameters given in Eq. (A.4), this results in a large discharge time. Another choice of the mass and charge could be m χ ∼ 10 -3 m p and q χ ∼ 10 -21 e , which leaves the discharge time unchanged. If the black hole is located in an environment where ρ ∼ 0 , the discharge time can be even longer. We note that the above estimate for the discharge time assumes an accretion rate given by  \ndM 2 dt ∼ M 2 2 Q ρ , (A.6)  \nwhere Q ρ ∼ ρ/v 3 is the dark fermion phase space density [82]. Moreover, it assumes that the entire density ρ is constituted by oppositely charged dark fermions, which is the most conservative assumption.  \nThe question of how a SMBH can acquire dark charge in the first place is a difficult one and is tied to the formation mechanism of SMBHs, which is unknown. It may be easier to answer this question for the case of stellar black holes first. Classic calculations [83] indicate that the charge of a typical star at birth may be given by  \n( Q M ) birth ∼ m χ q χ (A.7)  \nunder the assumption that the black hole is formed by the collapse of dark charged particles. For the Standard Model, the charge obtained from Eq. (A.7) is miniscule, and is promptly discharged as well. For dark sectors with m χ q χ ∼ 1 , however, one can obtain ( Q M ) birth ∼ 1 , and the discharge time can also be large.  \nFor a SMBH, the situation is less clear. The proliferation of SMBH observations with z > 5 raises the question about their origin, which is a mystery (we refer to [84] for a review). One possibility is the existence of intermediate mass black hole seeds with M 2 ∼ 10 3 -5 M ⊙ at z ∼ 20 , along with a population of stellar black holes. If these SMBH seeds are formed from the collapse of a dark matter halo, calculations similar to the one sketched above may apply to their initial charge. The subsequent accretion history of the seed is expected to follow mean Eddington rates with episodes of super-Eddington accretion [85-88]:  \nM 2 ( t ) ∼ M seed × e t/τ Edd , (A.8)  \nwhere τ Edd = 45 Myr is the Eddington timescale. Accretion under these circumstances can grow ∼ O (10 -100) M ⊙ seeds into black holes with mass ∼ 10 3 M ⊙ by z = 18 and ∼ 10 8 M ⊙ by z = 7 .  \nAssuming that the growth of the black hole mass between z ∼ 30 to z ∼ 7 proceeds by the accretion of equal numbers of positively and negatively charged dark particles, the final charge of the SMBH would be equal to the initial charge, while the mass is enhanced by orders of magnitude. This would yield, for example, ( Q/M ) final ∼ 10 -5 , assuming that ( Q/M ) birth ∼ 1 . On the other hand, preferential accretion of one charge over the other could produce ( Q/M ) final ∼ 1 .", 'References': "- [1] C. de Rham, J.T. Deskins, A.J. Tolley and S.-Y. Zhou, Graviton Mass Bounds , Rev. Mod. Phys. 89 (2017) 025004 [ 1606.08462 ].\n- [2] E.G. Adelberger, B.R. Heckel and A.E. Nelson, Tests of the gravitational inverse square law , Ann. Rev. Nucl. Part. Sci. 53 (2003) 77 [ hep-ph/0307284 ].\n- [3] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, Modified Gravity and Cosmology , Phys. Rept. 513 (2012) 1 [ 1106.2476 ].\n- [4] D. Brzeminski, Z. Chacko, A. Dev, I. Flood and A. Hook, Searching for a fifth force with atomic and nuclear clocks , Phys. Rev. D 106 (2022) 095031 [ 2207.14310 ].\n- [5] T. Damour and J.F. Donoghue, Phenomenology of the Equivalence Principle with Light Scalars , Class. Quant. Grav. 27 (2010) 202001 [ 1007.2790 ].\n- [6] B.S. Acharya, A. Maharana and F. Muia, Hidden Sectors in String Theory: Kinetic Mixings, Fifth Forces and Quintessence , JHEP 03 (2019) 048 [ 1811.10633 ].\n- [7] G. Kane, K. Sinha and S. Watson, Cosmological Moduli and the Post-Inflationary Universe: A Critical Review , Int. J. Mod. Phys. D 24 (2015) 1530022 [ 1502.07746 ].\n- [8] D. Croon, A.E. Nelson, C. Sun, D.G.E. Walker and Z.-Z. Xianyu, Hidden-Sector Spectroscopy with Gravitational Waves from Binary Neutron Stars , Astrophys. J. Lett. 858 (2018) L2 [ 1711.02096 ].\n- [9] J.A. Dror, B.V. Lehmann, H.H. Patel and S. Profumo, Discovering new forces with gravitational waves from supermassive black holes , Phys. Rev. D 104 (2021) 083021 [ 2105.04559 ].\n- [10] J. Kopp, R. Laha, T. Opferkuch and W. Shepherd, Cuckoo's eggs in neutron stars: can LIGO hear chirps from the dark sector? , JHEP 11 (2018) 096 [ 1807.02527 ].\n- [11] S. Alexander, E. McDonough, R. Sims and N. Yunes, Hidden-Sector Modifications to Gravitational Waves From Binary Inspirals , Class. Quant. Grav. 35 (2018) 235012 [ 1808.05286 ].\n- [12] T.K. Poddar, A. Ghoshal and G. Lambiase, Listening to dark sirens from gravitational waves:\\itCombined effects of fifth force, ultralight particle radiation, and eccentricity , 2302.14513 .\n- [13] M.I. Gresham, V.S.H. Lee and K.M. Zurek, Astrophysical observations of a dark matter-Baryon fifth force , JCAP 02 (2023) 048 [ 2209.03963 ].\n- [14] A. Kehagias and K. Sfetsos, Deviations from the 1/r**2 Newton law due to extra dimensions , Phys. Lett. B 472 (2000) 39 [ hep-ph/9905417 ].\n- [15] D.B. Kaplan and M.B. Wise, Couplings of a light dilaton and violations of the equivalence principle , JHEP 08 (2000) 037 [ hep-ph/0008116 ].\n- [16] P. Callin and F. Ravndal, Higher order corrections to the Newtonian potential in the Randall-Sundrum model , Phys. Rev. D 70 (2004) 104009 [ hep-ph/0403302 ].\n- [17] S. Fichet, E. Megias and M. Quiros, Continuum effective field theories, gravity, and holography , Phys. Rev. D 107 (2023) 096016 [ 2208.12273 ].\n- [18] S. Fichet, E. Megias and M. Quiros, Cosmological dark matter from a bulk black hole , Phys. Rev. D 107 (2023) 115014 [ 2212.13268 ].\n- [19] P. Brax, S. Fichet and G. Pignol, Bounding Quantum Dark Forces , Phys. Rev. D 97 (2018) 115034 [ 1710.00850 ].\n- [20] H. Banks and M. Mccullough, Charting the Fifth Force Landscape , Phys. Rev. D 103 (2021) 075018 [ 2009.12399 ].\n- [21] A. Costantino, S. Fichet and P. Tanedo, Exotic Spin-Dependent Forces from a Hidden Sector , JHEP 03 (2020) 148 [ 1910.02972 ].  \n- [22] X.-j. Xu and B. Yu, On the short-range behavior of neutrino forces beyond the Standard Model: from 1/r 5 to 1/r 4 , 1/r 2 , and 1/r , JHEP 02 (2022) 008 [ 2112.03060 ].\n- [23] C. Talmadge, J.P. Berthias, R.W. Hellings and E.M. Standish, Model Independent Constraints on Possible Modifications of Newtonian Gravity , Phys. Rev. Lett. 61 (1988) 1159.\n- [24] T. Johannsen, Sgr A* and General Relativity , Class. Quant. Grav. 33 (2016) 113001 [ 1512.03818 ].\n- [25] A. Hees et al., Testing General Relativity with stellar orbits around the supermassive black hole in our Galactic center , Phys. Rev. Lett. 118 (2017) 211101 [ 1705.07902 ].\n- [26] N. Yunes and X. Siemens, Gravitational-Wave Tests of General Relativity with Ground-Based Detectors and Pulsar Timing-Arrays , Living Rev. Rel. 16 (2013) 9 [ 1304.3473 ].\n- [27] C.M. Will, Bounding the mass of the graviton using gravitational wave observations of inspiralling compact binaries , Phys. Rev. D 57 (1998) 2061 [ gr-qc/9709011 ].\n- [28] V. Cardoso, C.F.B. Macedo, P. Pani and V. Ferrari, Black holes and gravitational waves in models of minicharged dark matter , JCAP 05 (2016) 054 [ 1604.07845 ].\n- [29] S. Babak, J. Gair, A. Sesana, E. Barausse, C.F. Sopuerta, C.P.L. Berry et al., Science with the space-based interferometer lisa. v. extreme mass-ratio inspirals , Phys. Rev. D 95 (2017) 103012.\n- [30] C.P.L. Berry, S.A. Hughes, C.F. Sopuerta, A.J.K. Chua, A. Heffernan, K. Holley-Bockelmann et al., The unique potential of extreme mass-ratio inspirals for gravitational-wave astronomy , 1903.03686 .\n- [31] S. Barsanti, A. Maselli, T.P. Sotiriou and L. Gualtieri, Detecting Massive Scalar Fields with Extreme Mass-Ratio Inspirals , Phys. Rev. Lett. 131 (2023) 051401 [ 2212.03888 ].\n- [32] A. Maselli, N. Franchini, L. Gualtieri, T.P. Sotiriou, S. Barsanti and P. Pani, Detecting fundamental fields with LISA observations of gravitational waves from extreme mass-ratio inspirals , Nature Astron. 6 (2022) 464 [ 2106.11325 ].\n- [33] LISA collaboration, Laser Interferometer Space Antenna , 1702.00786 .\n- [34] M. Colpi et al., LISA Definition Study Report , 2402.07571 .\n- [35] A. Sesana et al., Unveiling the gravitational universe at µ -Hz frequencies , Exper. Astron. 51 (2021) 1333 [ 1908.11391 ].\n- [36] B.C. Seymour and K. Yagi, Probing massive scalar and vector fields with binary pulsars , Phys. Rev. D 102 (2020) 104003 [ 2007.14881 ].\n- [37] R.D. Blandford and C.F. McKee, Reverberation mapping of the emission line regions of Seyfert galaxies and quasars. , The Astrophysical Journal 255 (1982) 419.\n- [38] B.M. Peterson, Reverberation Mapping of Active Galactic Nuclei , Publications of the Astronomical Society of the Pacific 105 (1993) 247.\n- [39] J.K. Hoormann, P. Martini, T.M. Davis, A. King, C. Lidman, D. Mudd et al., C IV black hole mass measurements with the Australian Dark Energy Survey (OzDES) , Monthly Notices of the Royal Astronomical Society 487 (2019) 3650 [ 1902.04206 ].\n- [40] A. Akylas, I. Papadakis and A. Georgakakis, Black hole mass estimation using X-ray variability measurements in Seyfert galaxies , Astron. Astrophys. 666 (2022) A127 [ 2208.12490 ].\n- [41] P. Amaro-Seoane, The gravitational capture of compact objects by massive black holes , 2011.03059 .\n- [42] V. Vazquez-Aceves, Y. Lin and A. Torres-Orjuela, Sgr A* Spin and Mass Estimates through the Detection of an Extremely Large Mass-ratio Inspiral , Astrophys. J. 952 (2023) 139 [ 2206.14399 ].\n- [43] P. Amaro-Seoane, Extremely large mass-ratio inspirals , Phys. Rev. D 99 (2019) 123025 [ 1903.10871 ].\n- [44] J. Binney and S. Tremaine, Galactic dynamics (1987).  \n- [45] M.C. Miller, T. Alexander, P. Amaro-Seoane, A. Barth, C. Cutler, J. Gair et al., Probing Stellar Dynamics in Galactic Nuclei , in astro2010: The Astronomy and Astrophysics Decadal Survey , vol. 2010, p. 209, Jan., 2009, DOI [ 0903.0285 ].\n- [46] Z. Pan and H. Yang, Formation Rate of Extreme Mass Ratio Inspirals in Active Galactic Nuclei , Phys. Rev. D 103 (2021) 103018 [ 2101.09146 ].\n- [47] GRAVITY collaboration, Polarimetry and astrometry of NIR flares as event horizon scale, dynamical probes for the mass of Sgr A* , Astron. Astrophys. 677 (2023) L10 [ 2307.11821 ].\n- [48] GRAVITY Collaboration, Detection of orbital motions near the last stable circular orbit of the massive black hole SgrA* , Astron. Astrophys. 618 (2018) L10 [ 1810.12641 ].\n- [49] GRAVITY collaboration, Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole , Astron. Astrophys. 615 (2018) L15 [ 1807.09409 ].\n- [50] S. Thater, D. Krajnovi'c, M. Cappellari, T.A. Davis, P.T. de Zeeuw, R.M. McDermid et al., Six new supermassive black hole mass determinations from adaptive-optics assisted SINFONI observations , Astron. Astrophys. 625 (2019) A62 [ 1902.10175 ].\n- [51] C.-P. Ma, J.E. Greene, N. McConnell, R. Janish, J.P. Blakeslee, J. Thomas et al., The massive survey. i. a volume-limited integral-field spectroscopic study of the most massive early-type galaxies within 108 mpc , The Astrophysical Journal 795 (2014) 158.\n- [52] N.J. McConnell and C.-P. Ma, Revisiting the scaling relations of black hole masses and host galaxy properties , The Astrophysical Journal 764 (2013) 184.\n- [53] B.M. Peterson, L. Ferrarese, K.M. Gilbert, S. Kaspi, M.A. Malkan, D. Maoz et al., Central masses and broad-line region sizes of active galactic nuclei. ii. a homogeneous analysis of a large reverberation-mapping database , The Astrophysical Journal 613 (2004) 682.\n- [54] M.C. Bentz and S. Katz, The AGN Black Hole Mass Database , Publications of the Astronomical Society of the Pacific 127 (2015) 67 [ 1411.2596 ].\n- [55] U. Malik, R. Sharp, A. Penton, Z. Yu, P. Martini, C. Lidman et al., OzDES Reverberation Mapping Program: H β lags from the 6-yr survey , Mon. Not. Roy. Astron. Soc. 520 (2023) 2009 [ 2210.03977 ].\n- [56] Y. Shen et al., The Sloan Digital Sky Survey Reverberation Mapping Project: Key Results , 2305.01014 .\n- [57] L.S. Finn and K.S. Thorne, Gravitational waves from a compact star in a circular, inspiral orbit, in the equatorial plane of a massive, spinning black hole, as observed by LISA , Phys. Rev. D 62 (2000) 124021 [ gr-qc/0007074 ].\n- [58] S.A. Teukolsky, Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations , Astrophys. J. 185 (1973) 635.\n- [59] S.L. Detweiler, Black holes and gravitational waves. I. Circular orbits about a rotating hole. , The Astrophysical Journal 225 (1978) 687.\n- [60] M. Shibata, Gravitational waves by compact star orbiting around rotating supermassive black holes , Phys. Rev. D 50 (1994) 6297.\n- [61] K. Glampedakis and D. Kennefick, Zoom and whirl: Eccentric equatorial orbits around spinning black holes and their evolution under gravitational radiation reaction , Phys. Rev. D 66 (2002) 044002.\n- [62] O. Christiansen, J. Beltrán Jiménez and D.F. Mota, Charged Black Hole Mergers: Orbit Circularisation and Chirp Mass Bias , Class. Quant. Grav. 38 (2021) 075017 [ 2003.11452 ].\n- [63] P.K. Gupta, T.F.M. Spieksma, P.T.H. Pang, G. Koekoek and C.V.D. Broeck, Bounding dark charges on binary black holes using gravitational waves , Phys. Rev. D 104 (2021) 063041 [ 2107.12111 ].\n- [64] L. Liu, Z.-K. Guo, R.-G. Cai and S.P. Kim, Merger rate distribution of primordial black hole binaries with electric charges , Phys. Rev. D 102 (2020) 043508 [ 2001.02984 ].  \n- [65] L. Liu, O. Christiansen, Z.-K. Guo, R.-G. Cai and S.P. Kim, Gravitational and electromagnetic radiation from binary black holes with electric and magnetic charges: Circular orbits on a cone , Phys. Rev. D 102 (2020) 103520 [ 2008.02326 ].\n- [66] L. Liu and S.P. Kim, Merger rate of charged black holes from the two-body dynamical capture , JCAP 03 (2022) 059 [ 2201.02581 ].\n- [67] C. Zhang and Y. Gong, Detecting electric charge with extreme mass ratio inspirals , Phys. Rev. D 105 (2022) 124046 [ 2204.08881 ].\n- [68] M. Sasaki and T. Nakamura, Gravitational Radiation From a Kerr Black Hole. 1. Formulation and a Method for Numerical Analysis , Prog. Theor. Phys. 67 (1982) 1788.\n- [69] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments , Oxford University Press (2007), 10.1093/acprof:oso/9780198570745.001.0001.\n- [70] C.J. Moore, R.H. Cole and C.P.L. Berry, Gravitational-wave sensitivity curves , Class. Quant. Grav. 32 (2015) 015014 [ 1408.0740 ].\n- [71] AEDGE collaboration, AEDGE: Atomic Experiment for Dark Matter and Gravity Exploration in Space , EPJ Quant. Technol. 7 (2020) 6 [ 1908.00802 ].\n- [72] M.A. Fedderke, P.W. Graham and S. Rajendran, Asteroids for µ Hz gravitational-wave detection , Phys. Rev. D 105 (2022) 103018 [ 2112.11431 ].\n- [73] L. Barack and C. Cutler, LISA capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy , Phys. Rev. D 69 (2004) 082005 [ gr-qc/0310125 ].\n- [74] C. Zhang, Y. Gong, H. Liu, B. Wang and C. Zhang, Sky localization of space-based gravitational wave detectors , Physical Review D 103 (2021) .\n- [75] LIGO SCIENTIFIC, VIRGO, 1M2H, DARK ENERGY CAMERA GW-E, DES, DLT40, LAS CUMBRES OBSERVATORY, VINROUGE, MASTER collaboration, A gravitational-wave standard siren measurement of the Hubble constant , Nature 551 (2017) 85 [ 1710.05835 ].\n- [76] E. Barausse, V. Cardoso and P. Pani, Can environmental effects spoil precision gravitational-wave astrophysics? , Phys. Rev. D 89 (2014) 104059 [ 1404.7149 ].\n- [77] E. Berti et al., Testing General Relativity with Present and Future Astrophysical Observations , Class. Quant. Grav. 32 (2015) 243001 [ 1501.07274 ].\n- [78] S. Kejriwal, L. Speri and A.J.K. Chua, Impact of Correlations on the Modeling and Inference of Beyond Vacuum-GR Effects in Extreme-Mass-Ratio Inspirals , 2312.13028 .\n- [79] E. Barausse, V. Cardoso and P. Pani, Environmental Effects for Gravitational-wave Astrophysics , J. Phys. Conf. Ser. 610 (2015) 012044 [ 1404.7140 ].\n- [80] L.B. Fries, J.R. Trump, M.C. Davis, C.J. Grier, Y. Shen, S.F. Anderson et al., The sdss-v black hole mapper reverberation mapping project: Unusual broad-line variability in a luminous quasar , The Astrophysical Journal 948 (2023) 5.\n- [81] G. Bozzola and V. Paschalidis, General Relativistic Simulations of the Quasicircular Inspiral and Merger of Charged Black Holes: GW150914 and Fundamental Physics Implications , Phys. Rev. Lett. 126 (2021) 041103 [ 2006.15764 ].\n- [82] S. Peirani and J.A. de Freitas Pacheco, Dark Matter Accretion into Supermassive Black Holes , Phys. Rev. D 77 (2008) 064023 [ 0802.2041 ].\n- [83] A.S. Eddington, The Internal Constitution of the Stars (1926).\n- [84] K. Inayoshi, E. Visbal and Z. Haiman, The Assembly of the First Massive Black Holes , Ann. Rev. Astron. Astrophys. 58 (2020) 27 [ 1911.05791 ].\n- [85] M.A. Abramowicz, B. Czerny, J.P. Lasota and E. Szuszkiewicz, Slim Accretion Disks , The Astrophysical Journal 332 (1988) 646.  \n- [86] M. Volonteri, J. Silk and G. Dubus, The Case for Supercritical Accretion Onto Massive Black Holes at High Redshift , Astrophys. J. 804 (2015) 148 [ 1401.3513 ].\n- [87] K. Inayoshi, Z. Haiman and J.P. Ostriker, Hyper-Eddington accretion flows on to massive black holes , Mon. Not. Roy. Astron. Soc. 459 (2016) 3738 [ 1511.02116 ].\n- [88] M.J. Hardcastle and J.H. Croston, Radio galaxies and feedback from AGN jets , New Astron. Rev. 88 (2020) 101539 [ 2003.06137 ]."}