# Repeating Nuclear Transients as Candidate Electromagnetic Counterparts of
LISA Extreme Mass Ratio Inspirals

Shubham Kejriwal,1 Vojtěch Witzany,2 Michal Zajaček,3 Dheeraj R. Pasham,4 and
Alvin J. K. Chua1,5
1Department of Physics, National University of Singapore, Singapore 117551
2Institute of Theoretical Physics, Faculty of Mathematics and Physics, Charles
University, CZ-180 00 Prague, Czech Republic
3Department of Theoretical Physics and Astrophysics, Faculty of Science,
Masaryk University, Kotlářská 2, 611 37 Brno, Czech Republic
4Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA, USA
5Department of Mathematics, National University of Singapore, Singapore 119076
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(AC)

###### Abstract

Extreme-mass-ratio inspirals (EMRIs) are one of the primary targets for the
recently adopted millihertz gravitational-wave (GW) observatory LISA. Some
previous studies have argued that a fraction of all EMRIs form in matter-rich
environments, and can potentially explain the dozens of soft X-ray band ($\sim
10^{-1}\rm keV$), low-frequency ($\sim 0.1$ mHz) periodic phenomena known as
quasi-periodic eruptions (QPEs) and quasi-periodic oscillations (QPOs). Here,
using a representative EMRI population retrofitted with cutoffs on LISA-band
SNRs and luminosity distances to account for the sensitivity of current
instruments, we estimate the mean frequency band in which QPEs and QPOs
originating from detectable LISA EMRIs may be emitting an X-ray signal “today”
(i.e., in 2024) to be $0.46\pm 0.22$ mHz. We also model the well-known QPO
source, RE J1034+396, which falls in this frequency band, as an EMRI assuming
its primary black hole mass to be $10^{6}-10^{7}M_{\odot}$. Through a prior-
predictive analysis, we estimate the orbiting compact object’s mass to be
$46^{+10}_{-40}M_{\odot}$ and the source’s LISA-band SNR as $\approx 14$,
highlighting it as a candidate multi-messenger EMRI target. We also highlight
the role of current and near-future X-ray and UV observatories in enabling
multi-messenger observations of EMRIs in conjunction with LISA, and conclude
with a discussion of caveats of the current analysis, such as the exclusion of
eccentricity and inclination from the model, and the measurability of sub-
solar mass compact object EMRIs.

###### keywords:

black hole physics – accretion, accretion discs – gravitational waves –
transients: black hole mergers – transients: black hole - neutron star mergers
– galaxies: nuclei

††pubyear: 2024††pagerange: Repeating Nuclear Transients as Candidate
Electromagnetic Counterparts of LISA Extreme Mass Ratio Inspirals–A

##  1 Introduction

Extreme-mass-ratio inspirals (EMRIs) are black hole binary systems in which a
stellar-mass compact object (CO) of mass $\mu\sim 1-10^{2}M_{\odot}$ completes
$\sim 10^{5}$ orbits around a supermassive black hole (MBH) of mass $M\sim
10^{5}-10^{7}M_{\odot}$ in its strong-gravity regime $(r\sim 10~{}GM/c^{2})$
over a period of $\sim 1-10$ years. They emit millihertz gravitational waves
signal (GW) in the process (Amaro-Seoane et al., 2015; Barack & Cutler, 2004),
and are one of the main targets of the recently adopted European Space Agency
(ESA)-led space-based GW observatory, Laser Interferometer Space Antenna
(LISA) (eLISA Consortium et al., 2013; Amaro-Seoane et al., 2017; Colpi et
al., 2024). EMRIs have the potential to be unique scientific laboratories for
conducting tests of general relativity (GR), constraining the astrophysical
MBH mass-function, and inferring cosmological parameters (Amaro-Seoane et al.,
2007; Amaro-Seoane, 2018; Gair et al., 2013; Babak et al., 2017; Berry et al.,
2019; Laghi et al., 2021).

Over the years, various EMRI formation channels have been proposed (Amaro-
Seoane, 2018): (i) Direct capture, where the many-body scattering in the
nuclear star cluster around an MBH pushes a CO into the gravitational-wave
loss cone (Hils & Bender, 1995; Sigurdsson & Rees, 1997); (ii) The tidal
disruption of a CO binary by the MBH (Hills mechanism; Hills, 1988) where one
component is ejected as a hypervelocity star and the other component becomes
more tightly bound to the MBH (Miller et al., 2005); (iii) The capture of COs
and their stellar precursors by the accretion disk in an active galactic
nucleus (AGN) and their subsequent migration to tighter orbits around the MBH
due to disk-CO interactions (Paczynski, 1978; Syer et al., 1991; Rauch, 1995;
Karas & Šubr, 2001; Pan & Yang, 2021); and (iv) The formation of stellar
objects later evolving into COs in situ within the accretion disk and their
continual migration towards the MBH while embedded within the disk (Paczynski,
1978; Collin & Zahn, 1999; Levin & Beloborodov, 2003; Levin, 2007; Sigl et
al., 2007). While the disk-assisted formation channels (iii)-(iv) are only
relevant for the $\sim 10^{-2}-10^{-1}$ fraction of galactic nuclei that are
(or recently have been) AGN, in some population models such as in that of Pan
& Yang (2021) they have significantly higher formation rates than the “dry”
channels (i)-(ii) and, consequently, produce a major fraction of the predicted
LISA EMRIs.

In addition to GW radiation, EMRI sources formed in channels (iii)-(iv) can
trigger high energy electromagnetic (EM)-band emissions owing to CO
interactions with the matter in the MBH accretion disk, making them viable
candidates for multi-messenger detections (Lehto & Valtonen, 1996; Semerák et
al., 1999; Ivanov et al., 1998; Pihajoki, 2016; Suková et al., 2021). Channels
(i)-(ii) can also lead to CO – hot flow interactions in low-luminosity
galactic nuclei. However, because of lower gas densities, the emerging weak EM
signal limits such GW sources to the local Universe and rather fine-tuned
setups with a low likelihood, such as a pulsar on a tight orbit around Sgr A*
(Eckart et al., 2017).

Identification of individual multi-messenger sources would be invaluable for
simplifying data analysis schemes such as the LISA global-fit (Vallisneri,
2009; Cornish & Crowder, 2005; Littenberg & Cornish, 2023). This can be done
by setting tight constraints on the extrinsic parameter set, while lifting
degeneracies arising from the GW-only inference of beyond vacuum-GR EMRI
parameters (Kejriwal et al., 2023). Such observations will also enable
independent constraints on the redshift of the host galaxy and the properties
of the galactic environment, without invoking statistical methods like the
dark siren method (Schutz, 1986a; Del Pozzo, 2012; see also Gray et al., 2023
for a review). Multi-messenger detections of EMRIs would thus help to decouple
the MBH parameters from different galaxy components (Kormendy & Ho, 2013;
Karas et al., 2021) and significantly improve the constraining power of LISA
on the luminosity distance-redshift relation in cosmology (Schutz, 1986b; Holz
& Hughes, 2005). The precise knowledge of the parameters of an accreting MBH
would additionally allow the resolution of long-standing conundrums such as
the radio loud/quiet dichotomy of AGN (Sikora et al., 2007) and their
conjectured relation to MBH spin (see, e.g., Tchekhovskoy et al., 2010).

Independent of the specific model, one can expect the CO-disk interactions to
manifest as regular time variability in the optical to X-ray bands,
corresponding to the temperatures of the accretion disk at various distances
from the MBH (Shakura & Sunyaev, 1973). Another possible signature could come
in the X-ray/radio band in the case of launched outflows (Suková et al., 2021;
Li et al., 2023; Pasham et al., 2024a). Variability in the soft X-ray band
(around $\sim 10^{-1}\ \rm keV$ or $\sim 10^{6}\,\rm K$) corresponds to
regions in the disk within a few gravitational radii near the MBH. At such
distances, GW radiation is also non-negligible. The characteristic period of
this variability should then correspond to the orbital period of the CO
inferred from the GW signal of the source, which will be $\sim$ hours to
minutes for EMRIs relevant to LISA (see Section 2).

In line with this description, the recently discovered quasi-periodic
eruptions (QPEs), which are signals of semi-regular bursts in soft X-ray
energy bands around $\sim 10^{-1}\rm keV$, and which can be traced back to
centers of known galaxies, have been proposed to be EM counterparts of
environment-rich EMRIs (Suková et al., 2021; Chen et al., 2022; Metzger et
al., 2022; Pasham et al., 2024a; Linial & Metzger, 2023; Franchini et al.,
2023; Zhou et al., 2024). Notably, QPEs may also be explained by accretion
disk instabilities around the MBH (Miniutti et al., 2019) or gravitational
lensing during the inspiral of massive black hole binaries (MBHBs) (Ingram et
al., 2021), but such models only explain a subset of the observed signals.
Coincidentally, some QPEs have been found in galaxies hosting MBHs of masses
$\sim 10^{6}M_{\odot}$ (Wevers et al., 2022; King, 2023b), which is consistent
with the mass range of expected LISA-relevant EMRIs ($10^{5}-10^{7}M_{\odot}$)
further motivating a multi-messenger search.

Another type of possibly related phenomenon are the less abrupt quasi-periodic
oscillations (QPOs) in the X-ray fluxes of AGNs, such as the ones reported in
RE J1034+396 (Gierliński et al., 2008), 2XMM J123103.2+110648 (Lin et al.,
2013), Sw J1644+57 (Wang et al., 2014), MS 2254.9-3712 (Alston et al., 2015),
1H 0707-495 (Pan et al., 2016), Mrk 766 (Boller, Th. et al., 2001; Zhang et
al., 2017), ASASSN-14li (Pasham et al., 2019) or MCG-06-30-15 (Gupta et al.,
2018), all with the estimated periods of oscillations falling within few
hours. While many studies have tried to make a connection between AGN QPOs and
their well-known kHz analogues in X-ray binaries (van der Klis, 2000;
Remillard & McClintock, 2006; Ingram & Motta, 2019), their mechanisms may be
completely unrelated. Indeed, it has been suggested in the literature that AGN
QPEs and QPOs may originate from the same source under different
configurations (King, 2023a; Linial & Metzger, 2023).

Hence, in this paper, we proceed with the assumption that AGN QPEs and QPOs
are both caused by EMRIs, jointly referring to them as quasi-periodic
phenomena (QPPs). The relative differences in the intensity of peaks and duty
cycles between QPEs and QPOs will then be assumed to be caused by variations
in the particular parameters of the sources such as the orbital inclination of
the EMRIs or the accretion states of the disks (see Section 4.2.4 for a
detailed discussion).

###  1.1 Executive Summary of Results

An effective multi-messenger study of LISA-relevant EMRIs would require us to
first determine the QPP frequency band in which their EM counterparts will
most likely be observable. We find in Section 2 that EMRIs in the LISA
observation window (assuming an operation time from the year 2037 to 2041)
should emit QPP signals “today”, i.e. in 2024, in a wide frequency range with
a mean of $\approx 0.46$ mHz and a standard-deviation of $\sigma\approx 0.22$
mHz. As a demonstrative case study, in Section 3 we analyze the well-known QPO
source RE J1034+396, which has a suitable oscillatory signature with the
period of $\sim 1$ hr in the relevant EM energy band of $\sim 0.2-10\,\rm keV$
(Gierliński et al., 2008; Jin et al., 2020). Assuming an MBH mass of
$10^{6}-10^{7}M_{\odot}$, we model it as an EMRI to obtain prior-predictive
distributions on the CO’s mass of $46^{+10}_{-40}M_{\odot}$, the source’s
LISA-band SNR with a mean of $\approx 14$, and LISA-band detection probability
of $\approx 0.25$ with an SNR cutoff of $15$, highlighting it as a potential
multi-messenger EMRI candidate. Given the QPP frequency band, in Section 4.1
we then highlight some current and near-future X-ray observatories that are
sufficiently sensitive to signals in this band, and can thus be employed to
obtain a catalog of EM counterparts of LISA-relevant EMRIs. We discuss some
caveats of our analysis in Section 4.2, and conclude in Section 5. We assume
the standard flat $\Lambda$CDM model with $H_{0}=70\,{\rm
km\,s^{-1}\,Mpc^{-1}}$ and $\Omega_{\rm m}=0.3$ throughout the text.

##  2 QPP Frequencies of LISA-relevant EMRIs Today

###  2.1 EMRI catalog

We begin by assuming that the astrophysical EMRI population is described by
the M1 population catalog of Babak et al. (2017)111A more complete analysis
may benefit from analyzing all 12 EMRI catalogs (M1-M12) generated in Babak et
al. (2017) to account for uncertainties pertaining to their underlying
astrophysical models. Our study illustratively examines only the fiducial M1
model.. However, the source parameters in the M1 catalog are given at plunge,
i.e. when the CO is near the MBH’s innermost stable orbit. Pozzoli et al.
(2023) back-evolved these plunge-time parameters to a random time
$t\in[0,T_{\rm back}]$ years before plunge (where $T_{\rm
back}=20(M/10^{4}M_{\odot})$ years) using the leading-order Peters & Mathews
waveform model (Peters & Mathews, 1963), given by:

$\displaystyle\frac{\mathrm{d}f_{\rm
orb}}{\mathrm{d}t}=\frac{96G^{5/3}}{5c^{5}}\left(2\uppi\right)^{8/3}\mathcal{M}^{5/3}f^{11/3}_{\rm
orb}\mathcal{F}(e),$ (1)
$\displaystyle\frac{\mathrm{d}e}{\mathrm{d}t}=-\frac{G^{5/3}}{15c^{5}}\left(2\uppi\right)^{8/3}\mathcal{M}^{5/3}f^{8/3}_{\rm
orb}\mathcal{G}(e),$ (2)

where $f_{\rm orb}$ is the orbital frequency of the EMRI, $e$ is its
eccentricity, and

$\displaystyle\mathcal{M}=\frac{\left(M\mu\right)^{3/5}}{\left(M+\mu\right)^{1/5}},$
(3)
$\displaystyle\mathcal{F}(e)=\frac{1+(73/24)e^{2}+(37/96)e^{4}}{\left(1-e^{2}\right)^{7/2}},$
(4)
$\displaystyle\mathcal{G}(e)=\frac{304e+121e^{3}}{\left(1-e^{2}\right)^{5/2}}.$
(5)

Pozzoli et al. (2023) thus describes a more realistic EMRI catalog such that
only a fraction of EMRIs will plunge in the LISA observation window. We employ
this modified catalog in our study. Additionally, note that the purpose of the
M1 catalog is to synthesize a generous set of EMRIs that is significantly
larger than what will be detected by LISA, which is then used to explore its
detection capabilities. As such, $T_{\rm back}$ for $M\sim 10^{6}M_{\odot}$
can be as large as $\sim 1000$ years. Sources this far away from plunge will
produce negligible $(\ll 1)$ SNRs in the LISA band in general. For simplicity,
we thus only consider sources that plunge within $t\leq 100$ years after the
start of LISA observation window. We additionally require the sources to have
eccentricities, $e\leq 0.8$, mass-ratios, $\mu/M\leq 10^{-4}$, and inclination
angles between the orbital plane and accretion disk, $\iota\geq 0.01$. We
assume the plane of the accretion disk is perpendicular to the spin-axis of
the MBH, but briefly discuss the effect of misaligned planes in Section 4.2.4.
These preliminary cutoffs are based on conservative estimates on the
distribution of environment-rich EMRI parameters and validity of the GW
waveform model (see also discussion in Sections 3 and 4.2.3).

Equations (1)-(5) do not account for contributions from the CO-disk
interactions, due to the unavailability of such models for generic orbit
(eccentric and inclined) Kerr MBH EMRIs. However, their inclusion may
introduce significant dephasing in the long-term inspiral evolution, as hinted
by previous studies (Kocsis et al., 2011; Barausse et al., 2014), and should
be explored in the future. A more realistic analysis should also use higher-
order inspiral evolution models with a larger parameter space (including, for
example, the evolution of inclination, MBH spin, etc.), in contrast to the
leading post-Newtonian (PN) order Peters & Mathews Schwarzschild EMRI model
employed here for its computational feasibility.

###  2.2 Back-evolving to current orbital frequencies

We randomly distribute the M1 catalog parameters in the expected LISA
operation window (Gregorian calendar year 2037 to
2041)222https://sci.esa.int/web/lisa/-/61367-mission-summary, and back-evolve
the EMRI’s orbital frequency and eccentricity parameters to the year 2024,
i.e. “today”, using Eqs. (1) and (2) (see Fig.1). From the EMRI orbital
frequencies “today”, $f_{\rm orb,init}$, we can then calculate the
distribution of QPP frequencies for the population’s EM counterparts, $f_{\rm
QPP}$, as

$\displaystyle f_{\text{QPP}}\equiv n_{\rm col}\times f_{\rm orb,init}$ (6)

where $n_{\rm col}$ is the number of CO-disk collisions per orbit, assumed to
be 2 throughout our analysis and supported by previous estimates (Linial &
Metzger 2023, 2024, see also Suková et al. 2021, Xian et al. 2021, and
Franchini et al. 2023). Note that this approximation is only valid for EMRIs
evolving in circular orbits, which is a suitable assumption for a leading-
order estimate of the distribution of $f_{\rm QPP}$, and since environment-
rich EMRIs are expected to rapidly circularize under the effect of CO-disk
interactions (Tanaka & Ward, 2004; see also Kley & Nelson, 2012), minimizing
the effect of eccentricity and inclination. We also discuss the potential
impact of introducing eccentricity and inclination in the model in Section
4.2.4.

In Fig. 1, we visualize the back-evolved detector-frame orbital frequency and
eccentricity trajectories for the considered EMRI population in the LISA band
till 2001.We also estimate their LISA-band SNRs, $\rho$, using the solar
system barycenter frame point-detector approximation of Bonetti & Sesana
(2020) (see Eq. (17), (18), and (23) in their paper). The cross-section of
these trajectories evaluated for the year 2024 thus gives the distribution of
detector-frame orbital frequencies $f_{\rm orb,init}=f_{\rm orb,source}/(1+z)$
and eccentricities $e_{\rm init}$ “today”, and are summarized in Fig. 2.

Figure 1: Back evolution of the trajectories of detector-frame orbital
frequency, $f_{\rm orb}$, (upper panel) and eccentricity, $e$, (lower panel)
for EMRIs that are uniformly distributed in LISA’s observation window (orange
band), i.e. years 2037-2041 to 2001. Time “today”, 2024, is demarcated by the
black dashed line. The boldness of each trajectory is proportional to the EMRI
SNR in the LISA observation time calculated using the approximate model
described in the text.

###  2.3 Simultaneous EM and GW observability

Name | Type | ref. | $z$ | $L_{x}$ | $f$ | MBH mass range
---|---|---|---|---|---|---
|  |  |  | (erg/s) | ($\times 10^{-4}$ Hz) | ($\times 10^{6}M_{\odot}$)
GSN 069 | QPE | Miniutti et al. (2019) | 0.0181 | $1.1\times 10^{43}$ | 0.3 | 0.4 - 2
eRO-QPE1 | QPE | Arcodia et al. (2021) | 0.0505 | $3.3\times 10^{42}$ | 0.2 | 0.1 - 1
eRO-QPE2 | QPE | Arcodia et al. (2021) | 0.0175 | $1.0\times 10^{42}$ | 1.2 | 0.1 - 1
eRO-QPE3 | QPE | Arcodia et al. (2024) | 0.024 | $4.2\times 10^{41}$ | 0.1 | 0.9 - 5.3
eRO-QPE4 | QPE | Arcodia et al. (2024) | 0.044 | $1.27\times 10^{43}$ | 0.2 | 17 - 68
RX J1301.9+2747 | QPE | Giustini et al. (2020) | 0.0237 | $1.6\times 10^{42}$ | 0.6 | 0.8 - 2.8
RE J1034+396 | QPO | Gierliński et al. (2008) | 0.042 | $4.0\times 10^{43}$ | 2.8 | 1 - 10
2XMMJ123103,2+110648 | QPO | Lin et al. (2013) | 0.13 | $2.5\times 10^{42}$ | 0.7 | 0.09 - 1.1
ASASSN-14li | QPO | Pasham et al. (2019) | 0.0206 | $2.0\times 10^{43}$ | 76.3 | 0.6 - 12.5
ASASSN-20qc | QPOut∗ | Pasham et al. (2024a) | 0.0136 | $(2-4)\times 10^{44}$ | 0.019 | 16-160
Swift J0230+28 | QPE + QPOut(?)∗ | Guolo et al. (2023) | 0.0053 | $(3-6)\times 10^{42}$ | 0.019 | 1.6-10
Table 1: Summarized properties of QPEs and QPOs discovered so far in the
frequency band of $\sim 0.1$ mHz. From left to right, we specify the source’s
name, type (QPO/QPE/QPOut), redshift ($z$), peak soft X-ray (0.2-2 keV or
0.2-10 keV) luminosity, the QPP frequency ($f$), and current lower and upper
bounds on the source’s MBH mass. Asterisk (∗): QPOuts stands for the newly-
discovered phenomenon of quasi-periodic outflows where the X-ray modulation is
driven by a narrow bandpass dominated by an outflow. Recurrent spectral
features reminiscent of outflows were reported from Swift J0230+28, but their
statistical robustness is yet to be verified.

The M1 catalog of Babak et al. (2017) has EMRIs at redshifts up to $z=4.5$ or
luminosity distance $d_{L}\sim 40$ Gpc, which is much larger than the
distances covered by current state-of-the-art X-ray band instruments. The
observable volume is further restricted by the requirement of small
integration times necessary for capturing the $\sim$ hours variability of
relevant QPP sources. Similarly, for LISA-band EMRI observations,
detectability of sources will be constrained by their GW SNRs. A large
fraction of EMRIs in the M1 catalog have SNRs $\leq 1$ even after setting the
cutoffs mentioned in Section 2.1, making them undetectable. For a simultaneous
observation in EM and GW bands, we thus must set realistic cutoffs on the
observable cosmological volume (or equivalently, the maximum luminosity
distance, $d_{L}$) and the lowest detectable LISA-band SNR, $\rho$.

The luminosity distance cutoff can be set by considering both the sensitivity
threshold of current X-ray instruments and the required exposure time to
significantly detect the QPP EM emission. For the mean X-ray luminosity of QPE
sources, $\overline{L}_{\rm X}=2\times 10^{42}\,{\rm erg\,s^{-1}}$ (see Table
1), a luminosity distance of $d_{\rm L}\lesssim 1\,{\rm Gpc}$ results in the
X-ray flux density threshold of $f_{\rm X,sens}\equiv L_{\rm X}/4\pi d_{\rm
L}^{2}\gtrsim 1.7\times 10^{-14}\,{\rm erg\,s^{-1}\,{\rm cm^{-2}}}$. Since
current X-ray detectors (Swift, XMM-Newton, and eROSITA), meet the sensitivity
threshold $f_{\rm X,sens}$ for suitable integration times of a few tens of
minutes, we treat the luminosity distance $d_{L}=1$ Gpc as the cutoff (see
also Section 4.1 for a more detailed discussion of the capabilities of current
and future X-ray detectors).

For the LISA-band SNR cutoff, $\rho$, we choose a value of 15, given the
unavailability of a universally accepted threshold (which should scale with
the information volume of the signal space (Moore et al., 2019)). Our specific
choice is motivated by the estimates produced as part of the mock LISA data
challenges (MLDC), which hinted at an EMRI measurability threshold of
$\rho\approx 20$ for a two-year observation window (Babak et al., 2010), and
accounting for the accumulation of SNR over the full LISA observation window,
assumed here to be four years (Amaro Seoane et al., 2021). The SNR cutoff for
an EMRI’s detectability would also depend on the number of parameters under
simultaneous inference. If, for example, an EM counterpart detection tightly
constrains some extrinsic EMRI parameters, the SNR cutoff can be expected to
further drop, justifying our choice of $\rho_{\rm cutoff}=15$. However, this
cutoff ignores the effect of data gaps and glitches in the LISA data stream
which complicates inference and can thus potentially increase the SNR
requisites. In general, their contributions are expected to be negligible in
the inference of LISA-relevant sources (Dey et al., 2021; Baghi et al., 2019)
and are thus ignored in our analysis.

Fig. 3 describes the distribution of the detector-frame orbital frequencies
and eccentricities “today” (in 2024) after passing the EMRI population through
a cutoff of $d_{\rm L}\leq 1$ Gpc and $\rho\geq 15$.333Note that the large gap
along the $d_{L}$-axis in Fig. 3 between the sources marked with a star is an
artifact of the sampling resolution of Babak et al. (Babak et al., 2017)’s
population, but since we are interested in the span of these sources along the
frequency axis, this artifact can be safely ignored. The corresponding 1-D
sample distribution of $f_{\rm orb,init}$ has mean $\langle f_{\rm
orb,init}\rangle\approx 2.3\times 10^{-4}$ Hz and 1$\sigma$ standard-deviation
of $1.1\times 10^{-4}$ Hz. Using Eq. 6 with $n_{\rm col}=2$, a rough estimate
on the frequency of LISA-relevant EMRI QPPs is obtained as $\langle f_{\rm
QPP}\rangle\approx 4.6\pm 2.2\times 10^{-4}$ Hz. This estimate represents a
generous spread in the $\sim 0.1$ mHz frequency band and numerous candidate
and confirmed QPEs and QPOs in this frequency band have been reported in the
literature (Miniutti et al., 2019; Arcodia et al., 2021, 2024; Giustini et
al., 2020; Chakraborty et al., 2021; Gierliński et al., 2008; Jin et al.,
2020; Alston et al., 2015; Lin et al., 2013; Boller, Th. et al., 2001; Pan et
al., 2016; Pasham et al., 2019; Wang et al., 2014; Gupta et al., 2018).
However, follow-up studies have speculated upon the statistical robustness of
some of these observations (see for example Vaughan (2005) and Zhang et al.
(2023)).

A list of QPPs that fall near the proposed frequency band, and with
statistically significant observations, is presented in Table 1. For
completeness, we additionally list a new class of nuclear transients called
QPOuts (quasi-periodic outflows) discovered by Pasham et al. (2024b) which are
also proposed to be (early) electromagnetic counterparts of large mass-ratio
inspirals. In Pasham et al. (2024b), the authors observed quasiperiodic dips
in the X-ray band caused by the enhanced absorption due to an intermittent
fast outflow. The QPOut was interpreted as possibly caused by an intermediate-
mass black hole orbiting a supermassive black hole as supported by the
simulations of Suková et al. (2021). However, only one such confirmed
transient has been observed so far (but see also Guolo et al. 2023) and is
most likely far from merger. We thus leave its analysis for future studies.

As we discuss in Section 4.1, some current and near-future EM observatories
should be able to detect QPP signals in the estimated frequency band, enabling
potential multimessenger observations of LISA-EMRIs.

Figure 2: The joint distribution (top panel) and 1-D marginals (bottom panel)
of the detector-frame orbital frequency, $f_{\rm orb,init}$ and eccentricity
$e_{\rm init}$ of LISA-band EMRIs “today”, i.e. in 2024. Each marker in the
top panel is darkened in proportion to the source’s SNR. In the bottom panel,
the distributions of sources satisfying the LISA-band detectability cutoff
based on these SNRs are also depicted. Figure 3: Distribution of EMRI sources
with respect to detector-frame $f_{\rm orb,init}$ and their luminosity
distance from the detector, $d_{L}$, “today”, i.e. in 2024. Sources in the
grey band have $d_{L}\geq 1$ Gpc. The stars mark all sources that
simultaneously satisfy $d_{L}\leq 1$ Gpc and the cutoff on the LISA-band SNR,
$\rho\geq 15$. 1D marginal distribution of the detector frame $f_{\rm
orb,init}$ of sources that satisfy the cutoffs is depicted along the
horizontal axis, which thus describes the distribution of orbital frequencies
of candidate multi-messenger EMRIs.

##  3 Modeling RE J1034+396 as an EMRI

RE J1034+396, an AGN located at $z=0.042$, is a confirmed QPO with a frequency
of several times $10^{-4}$ Hz (Alston et al., 2014; Jin et al., 2020). Several
estimates on the mass of its MBH are available, (Bian & Huang, 2010; Chaudhury
et al., 2018; see also Czerny et al., 2016 for the comparison of the inferred
MBH mass using different methods), but it is generally agreed to be in the
range of $10^{6}-10^{7}M_{\odot}$. It has been observed multiple times over
the span of about a decade (2007-2018) — In 2007, the QPO’s time period was
observed to be $T_{2007}\sim 3730\pm 80s$ (Gierliński et al., 2008) which in
2018 reduced to $T_{2018}\sim 3550\pm 60s$ (Jin et al., 2020). If RE J1034+396
is modeled as an EMRI, the reduction in periodicity can be associated with GW-
induced radiation reaction of the system, and/or loss of energy via CO-disk
interactions. However, according to Table 1 in Jin et al. (2020) (see also
Table 2 in Alston et al., 2014), the QPO frequency fluctuated between
$0.25-0.27$ mHz during 2007 - 2011 without any secularly evolving features.
This can still be explained under the environment-rich EMRI model through
eccentric orbit CO-disk interactions, in which the torque forces generated by
the disk’s density wave perturbations both inside and outside of the CO orbit
compete with each other, sometimes canceling each other out to create
“migration traps”, or even exerting a net-outwards force on the CO (Bellovary
et al., 2016; Secunda et al., 2020). A detailed analysis of this effect,
however, is beyond the scope of our study, and we assume a secular drift
between the reported QPO time-periods in 2007 ($3730\pm 80$s) and 2018
($3550\pm 60$s) in our model, which mark the earliest and latest reported
values of the QPO in RE J1034+396.

###  3.1 CO mass estimate for RE J1034+396

Assuming two collisions per orbit, from Eq. (6), the corresponding orbital
frequencies of the EMRI source in 2007 and 2018 are

$\displaystyle f_{{\rm orb},i}=\frac{1}{2T_{i}}=\begin{cases}1.34\times
10^{-4}\,\rm{Hz}&i=2007\\\ 1.41\times 10^{-4}\,\rm{Hz}&i=2018\end{cases}$ (7)

such that the observed rate of change of the orbital frequency, $\dot{f}_{\rm
orb,obs}$, can be approximated to linear-order as

$\displaystyle\dot{f}_{\rm orb,obs}\approx\frac{f_{\rm orb,2018}-f_{\rm
orb,2007}}{2018-2007}\approx 6.18\times 10^{-7}\ \ \text{Hz yr}^{-1}.$ (8)

We can now use the distribution of $\dot{f}_{\rm orb,obs}$ to optimize the
corresponding theoretical prediction,

$\displaystyle\dot{f}_{\rm orb,th}\equiv\dot{f}_{\rm
orb,th}(M,\mu,a,e_{0},\iota,f_{\rm orb,2007})$ (9)

by minimizing the following $\chi^{2}$-error function,

$\displaystyle\chi^{2}\equiv\left(\dot{f}_{\rm orb,obs}-\dot{f}_{\rm
orb,th}\right)^{2}.$ (10)

Here $\mu$ is the unknown mass of the CO, $a$ is the dimensionless spin of the
MBH, $e_{0}$ is the initial eccentricity of the CO’s orbit, and $\iota$ is the
initial inclination of the CO’s orbital plane (which does not usually evolve
significantly in the parameter regime we are considering). The inclusion of
$\\{a,e_{0},\iota\\}$ in the parameter set allows us to gauge their
correlations with the calculated CO mass. A discussion on the choice of prior
on $\\{a,e_{0},\iota\\}$ and their effect on our analysis is presented in
Appendix A. We adopt $\dot{f}_{\rm orb,th}$ from the waveform model of Fujita
& Shibata (2020), as implemented in the 5PNAAK waveform model of the
FastEMRIWaveform (FEW) toolkit (Chua et al., 2021; Katz et al., 2021). Due to
the unavailability of CO-disk interaction models for inclined and eccentric
orbit EMRIs, we exclude their effects from our analysis. However, we stress
that such models will be necessary before drawing strong conclusions.

The first three panels in Fig. 4 describe the distribution on the parameter
set $\\{M,f_{\rm{orb,2007}},\dot{f}_{\rm{orb,obs}}\\}$, from which we infer
the best-fit CO mass of

$\mu\approx 46_{-40}^{+10}M_{\odot},$

as depicted in the figure’s fourth panel. This estimate comfortably falls in
the range of CO masses expected for LISA-relevant EMRIs ($\approx
1-100M_{\odot}$). We also refer to Fig. 11 in the Appendix which explicitly
shows the distribution with respect to all other considered parameters, i.e.
$\\{a,e_{0},\iota\\}$. Note, however, that the correlation coefficient,
$\mathcal{C}_{ij}\equiv\sigma_{ij}/\sqrt{\sigma^{2}_{i}\sigma^{2}_{j}}$,
between $\mu$ and $\\{a,e_{0},\iota\\}$ is relatively small with
$\mathcal{C}_{\mu a}\sim 0.01,\mathcal{C}_{\mu e_{0}}\sim-0.005$, and
$\mathcal{C}_{\mu\iota}\sim 0.02$ which implies that the prior choice of
$\\{a,e_{0},\iota\\}$ does not significantly affect the distribution of $\mu$
in this model. $\mu$’s correlation with $\dot{f}_{\rm{orb,obs}}$ is highest at
$\mathcal{C}_{\mu\dot{f}_{\rm{orb,obs}}}\sim 0.9$, while with the initial
orbital frequency is also high at $\mathcal{C}_{\mu
f_{\rm{orb,2007}}}\sim-0.76$. Comparatively, $\mu$ is moderately correlated
with $M$, $\mathcal{C}_{\mu M}\sim-0.5$.

Figure 4: Marginalized distributions over $\\{M,f_{\rm orb,2007},\dot{f}_{\rm
orb,obs},\mu,\text{SNR}\\}$ of RE J1034+396 modeled as an extreme-mass-ratio
inspiral. $\\{M,f_{\rm orb,2007},\dot{f}_{\rm orb,obs}\\}$ are estimated from
current observational constraints on RE J1034+396. The thus obtained prior-
predictive distribution on $\mu$ is described in the fourth panel, and the
corresponding LISA-band SNR in the final panel. The full 8-dimensional
distribution (including $\\{a,e_{0},\iota\\}$) is available in Appendix A as
Fig. 11. Figure 5: Prior-predictive distribution of LISA-band SNRs of QPO RE
J1034+396 modeled as an EMRI, with an indicative exponential fit. We assume
the MBH mass to be in range $10^{6}-10^{7}M_{\odot}$, mass of the orbiting
body $\mu\approx 46^{+10}_{-40}M_{\odot}$, and conservative prior
distributions on $\\{a,e,\iota\\}$ as described in the text.

###  3.2 LISA band detectability of RE J1034+396

The availability of the full EMRI parameter set,
$(M,\mu,a,e_{0},\iota,f_{\rm{orb,2007}})$, now allows us to evolve RE
J1034+396 to LISA’s operation time, 2037, and calculate the SNR, $\rho$, of
the GW signal, given as (Barack & Cutler, 2004)

$\displaystyle\rho\equiv\sqrt{\left<h|h\right>}$ (11)

where $h\equiv
h(M,\mu,a,e_{0},\iota,f_{\rm{orb,2007}},\boldsymbol{\psi}_{\rm{ext}})$ is the
waveform template corresponding to the signal, and
$\boldsymbol{\psi}_{\rm{ext}}=\\{d_{L}=0.18,\theta_{S}=2.77,\phi_{S}=0.69\\}$
is the extrinsic parameter set defined by the luminosity distance $d_{L}$ and
sky localisation parameters $\theta_{S},\phi_{S}$ of the source. The one-sided
detector noise weighted inner-product $\left<\cdot|\cdot\right>$ between two
discrete time-domain signals $\alpha(t)$ and $\beta(t)$ is given as (Chua &
Cutler, 2022):

$\displaystyle\left<\alpha|\beta\right>\equiv
4\text{Re}\sum_{\chi}\sum_{f=0}^{f_{N}}\delta
f\frac{\tilde{\alpha}^{*}_{\chi}(f)\tilde{\beta}_{\chi}(f)+\tilde{\alpha}_{\chi}(f)\tilde{\beta}^{*}_{\chi}(f)}{S_{{\rm
n},\chi}(f)}.$ (12)

Here $\tilde{\alpha}(f),\tilde{\beta}(f)$ are the frequency domain Fourier
transformations of $\alpha(t),\beta(t)$, $f_{N}$ is the Nyquist sampling
frequency, $\chi\equiv\text{I},\text{II}$ are the two independent data-
channels of LISA in its long-wavelength approximation (LWA) configuration, and
$S_{\rm n}(f)$ is the power spectral density of LISA’s noise channel
$n_{\chi}$ evaluated at detector frequency $f$ (Robson et al., 2019).

The prior-predictive distribution of the source’s LISA-band SNR is
approximately exponential, with a sample mean of $\rho\approx 14$ (see Fig.
5). From this distribution, we can estimate the probability of detecting RE
J1034+396 in the LISA band as ${\rm p}_{\rm detection}(\rho\geq 15)\approx
0.25$ for a LISA-band SNR cutoff of 15. Note that, in our model, the CO is on
average separated from the MBH by about $18.8~{}GM/c^{2}$ when entering the
LISA band and evolves to an average of $18.4~{}GM/c^{2}$ at the end of LISA
observation period. If instead a model accounting for CO-disk interactions in
generic orbits is used, the system will lose energy more rapidly and thus
attain smaller separations when entering the LISA band, boosting the strength
of its GW emission and hence its LISA-band SNR. Additionally, recall that the
SNR cutoff of 15 was set based on conservative estimates for a GW-only
inference. A multiband observation may enable the detection of sources with
smaller SNRs. Illustratively, for smaller LISA-band SNR cutoffs, we have ${\rm
p}_{\rm detection}(\rho\geq 12)=0.33$ and ${\rm p}_{\rm detection}(\rho\geq
10)=0.41$. These estimates suggest that RE J1034+396 may be a potential target
for multimessenger EMRI observations.

##  4 Discussion

###  4.1 Targeted QPP searches using current & planned observatories

####  4.1.1 Detectability threshold in soft X-rays

For the electromagnetic detection of repeating X-ray transients, generally two
aspects should be met: (i) sensitivity close to the eruption peak flux density
and (ii) a sufficient EM-band SNR reached for the integration time that is
shorter than the periodicity of a given source.

Conditions (i) and (ii) can be combined to estimate the minimal count rate
required to detect the flare robustly to constrain its spectral energy
distribution (SED), in particular its temperature. We consider the 5$\sigma$
threshold to be 25 counts. Since the QPE flares last for about an hour and
repeat every few hours, the observation should last for at least a couple of
hours, which results in the count rate of $25/7200\sim 0.003$ counts/s. For
currently operational X-ray telescopes (except for Chandra which has a
deteriorating soft X-ray sensitivity), this count rate translates into the
limiting flux density threshold of $f_{X}\approx 6.3\times 10^{-14}\,{\rm
erg\,cm^{-2}\,s^{-1}}$ for Swift XRT (Burrows et al., 2005), and $f_{X}\approx
3.2\times 10^{-15}\,{\rm erg\,cm^{-2}\,s^{-1}}$ for XMM/PN (Hasinger et al.,
2001; Brunner et al., 2008) 444Count to flux transformation is estimated using
https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl.. In both
cases, we considered the energy range of $0.3-2.0$ keV and the thermal
spectrum with the temperature of $100$ eV. For the NICER telescope, we adopt
the minimal count rate of $0.51$ counts/s in the 0.3-2.0 keV band, which
corresponds to $f_{X}\approx 4.3\times 10^{-13}\,{\rm erg\,cm^{-2}\,s^{-1}}$
(Remillard et al., 2022). Assuming the characteristic X-ray luminosity of the
eruptions to be $L_{X}=5\times 10^{42}\,{\rm erg\,s^{-1}}$, we estimate the
maximum luminosity distance (cosmological volume) which the current detectors
can probe using the simple relation, $d_{\rm L}=[L_{X}/(4\pi f_{X})]^{1/2}$.
Thus, we get $d_{\rm L}\sim 814$ Mpc and $\sim 3613$ Mpc for the Swift and the
XMM-Newton telescopes, respectively. For the NICER X-ray timing telescope, the
limiting luminosity distance is $\sim 312$ Mpc. Overall, this large scatter in
the cosmological volume justifies the distance cut-off of 1 Gpc. The caveat of
the simple estimates above is the assumption that the EM counterparts of gas-
rich EMRIs are all similar to the currently detected soft X-ray QPEs. In
principle, some of them can be more or less luminous and could thus manifest
in different wavebands, in particular harder X-ray or softer bands such as UV,
which directly affects the number of detected sources. Below we discuss some
of those aspects considering the current and future survey telescopes,
including radio and UV surveys in addition to X-ray telescopes.

Figure 6: The flux density or flux sensitivity (in ${\rm
erg\,\,s^{-1}\,cm^{-2}}$) as a function of the cumulative exposure time (in
hours). Left: Sensitivity levels as a function of exposure time (in hours) for
currently used instruments capable to detect soft X-ray emission (NICER,
Einstein Probe, Swift XRT, XMM EPIC, eROSITA). The stars represent already
detected repeating X-ray transients with a peak flux density and periodicity
(QPPs; see Table 1). Right: Same as in the left panel but for some of the
planned X-ray missions (Athena, Lynx, AXIS, STROBE-X). These levels should be
considered as estimates since the actual sensitivity depends on the spectrum
of the source, cosmic and galactic X-ray background, and the real-time
performance of the detector in space.

####  4.1.2 Detectability by specific X-ray missions

Currently detected soft X-ray repeating sources (QPPs) are plotted in Fig. 6
along with the flux sensitivity limits of current (left panel) and future-
planned detectors in the soft X-ray band (right panel). From this, we infer
that the NICER telescope could detect most currently observed sources since
they are above the flux sensitivity limit (blue dotted line in Fig. 6;
Remillard et al., 2022), while the SRG/eROSITA survey in soft X-rays (0.2-2.3
keV) (green solid line; Predehl et al., 2021) and Swift XRT (0.2-10 keV)
(maroon dash-dotted line; Burrows et al., 2005) should be able to detect all
of these sources. XMM-Newton EPIC detector (0.5-2 keV) (cyan dotted line;
Hasinger et al., 2001; Brunner et al., 2008) is the most sensitive in the
X-ray domain from the list of current instruments. The Einstein probe launched
on January 9, 2024 555See e.g.
https://heasarc.gsfc.nasa.gov/docs/objects/heapow/archive/technology/einsteinprobe.html.,
will have a wide field of view ($\sim 3600\,{\rm deg^{2}}$) using the Lobster-
eye technology (the Wide-field X-ray Telescope, WXT) and will be able to
detect about half of the currently detected transients (dashed orange line;
Yuan et al., 2022) at 0.5-4 keV. The future Athena mission (Nandra et al.,
2013), currently under preparation, with a field of view of $0.4\,{\rm
deg^{2}}$ (Wide Field Imager) and angular resolution of 5-10 arcseconds, will
be able to detect even fainter sources in the soft X-ray band with the flux
density at least two orders of magnitude lower (see magenta and orange lines
in the right panel of Fig. 6). The concept of the Lynx mission (Gaskin et al.,
2019) with a smaller field of view of 0.1${\rm deg^{2}}$ and an angular
resolution of 0.5” should reach even lower flux densities (dashed green line).
666The Athena and Lynx sensitivity curves in Fig. 6 were calculated following
the proposed X-ray observatories O1, O2, and O3 according to Lops et al.
(2023). The concept of Advanced X-ray Imaging Satellite (AXIS) mission (0.5-10
keV; Mushotzky et al., 2019; AXIS Time-Domain et al., 2023) with the exquisite
spatial resolution ($<1.5^{\prime\prime}$) across 24’ field of view will reach
the flux sensitivity intermediate between the planned Athena and the current
XMM EPIC detector. The planned STROBE-X Wide Field Monitor as well as its
X-ray Concentrator Array (2-50 keV and 0.2-12 keV, respectively, Ray et al.,
2019) will reach the sensitivity to detect just the brightest QPP sources
(dashed and dotted lime lines in Fig. 6). The detection and classification of
new X-ray QPP sources will benefit from the combination of current and future
survey telescopes, such as eROSITA and Einstein probe, with the dedicated
monitoring of the detected QPPs with more sensitive detectors, e.g. on-board
of XMM-Newton and Athena. We note that Athena and Lynx missions are currently
planned to be launched in the mid-2030s, hence they are not relevant for the
detection of EMRI EM counterparts “today”, though they will be crucial for
such detections and monitoring in the era of future space-borne GW detectors,
possibly including LISA.

####  4.1.3 Dependence on the accretion state and prospects for detectability
in the UV domain

The prospect of detecting EMRIs in the electromagnetic spectrum can also be a
function of the relative accretion rate, as suggested by Linial & Metzger
(2024). Within their model where the QPEs are generated by the ejection of
shocked plasmoids as the CO/star punches through the accretion disk, the
flares can effectively be detected when their radiation is spectrally harder
and hence brighter in the soft X-ray band than that of the underlying
accretion flow. This is not met for accretion rates close to the Eddington
limit, i.e. when the relative accretion rate (or equivalently the Eddington
ratio) is $\dot{m}\equiv\dot{M}/\dot{M}_{\rm Edd}\sim 1$ (where $\dot{M}_{\rm
Edd}$ is the Eddington accretion rate as defined in Linial & Metzger, 2024).
As the relative accretion rate $\dot{m}$ drops to $\dot{m}\sim 0.1$, e.g.
following a tidal disruption event with a decreasing accretion rate, QPE
emission in the soft X-ray band is above the quiescent disc emission due to
the higher temperature associated with QPEs. Linial & Metzger (2024) discuss
that for a further drop in the accretion rate to $\dot{m}\sim 0.01$, the QPEs
become too dim in the X-ray domain due to softer radiation. At this stage, the
eruptions could rather be traceable in the UV domain. An example of such
sources is provided by the AGN QPE source, GSN069, with its generally
decreasing accretion rate following two TDE-like outbursts. Eruptions in the
source were not detected for accretion rates $\dot{m}\geq 0.5$, but
(re)appeared for lower accretion rates (Miniutti et al., 2019, 2023a, 2023b).
Future monitoring of this source by X-ray and UV instruments will reveal how
the eruption SED changes with the drop in the accretion rate. A drop in
eruption amplitudes that can be potentially linked to a decrease in accretion
rate is also traced for eRO-QPE1 (Pasham et al., 2024b; Chakraborty et al.,
2024). Note that a UV-band luminosity of $\sim 10^{41}\rm ergs^{-1}$, typical
for such sources, corresponds to the apparent AB magnitude of $\sim 21.8$ mag
for nearby sources ($z\sim 0.025$). Such sensitivity will be reached by the
future-planned wide field-of-view NUV satellite ULTRASAT (Ben-Ami et al.,
2022; Shvartzvald et al., 2023), which will perform continuous monitoring of
the transient UV sky (230-290 nm) thanks to its wide field-of-view of
$204\,{\rm deg^{2}}$. It will also provide target-of-opportunity triggers for
the planned two-band (FUV, 140-190 nm, and NUV, 260-360 nm) UV photometry
mission QUVIK (Werner et al., 2024; Krtička et al., 2024; Zajaček et al.,
2023), which will be able to perform a high-cadence dedicated monitoring of UV
QPP candidates down to $\sim 22$ and $\sim 20$ AB magnitudes in the NUV and
FUV bands (detector SNR of 5 in 1000 seconds) with the point spread function
of $\lesssim 2.5^{\prime\prime}$. In addition, the UVEX mission (Kulkarni et
al., 2021) with the planned launch in the early 2030s will perform medium-
resolution spectroscopy in the range 115.0–265.0 nm in addition to deep two-
band monitoring of the UV sky with its sensitivity reaching $\sim 24.5$ AB
magnitude (detector SNR of 5 in 900 seconds). Sensitive UV observations may
also reveal long-period UV QPEs tracing EMRIs crossing the disc on wider
orbits, i.e. at lower relative velocities producing softer thermalized shocks.
Alternatively, short-period QPEs around low-mass central black holes of
$\lesssim 10^{5.5}\,M_{\odot}$ are also expected to be detected in the UV
domain (Linial & Metzger, 2024).

####  4.1.4 Prospects for coincident detection in the radio band

The detection of gas-rich EMRI in the radio domain is highly uncertain,
depending mainly on the presence of the magnetized outflow or jet in EMRI
hosts, i.e. whether the source is radio-loud or quiet. The viewing angle also
affects the jet contribution via the Doppler-boosting factor. One of the most
promising sources for hosting a binary MBH – OJ287 – is a blazar that exhibits
a quasiperiodic precession/nutation pattern in the temporal evolution of the
jet position as well as its radio brightness on the timescale of $\sim 20$
years (Britzen et al., 2018, 2023), which is about twice longer than the
typical optical variability timescale in this source (Valtonen et al., 2016).
However, most of the QPP hosts in Fig. 1 are radio quiet and hence the
prospects for tracing the EMRI presence in the radio variability is at best at
the limit of the sensitivity of current radio telescopes. From the repeating
soft X-ray transients (QPE-like) only Swift J0230+28 exhibited one radio flare
coincident with the start of the X-ray flare (Guolo et al., 2023). The radio
flare detected by the Very Large Array (VLA; D configuration) was unresolved
and coincident with the point source in the nucleus. The flux density at 10
GHz was $F_{\nu}=93\pm 7\mu{\rm Jy}$ (13$\sigma$ detection). Given the
distance of the host galaxy ($\sim 165\,{\rm Mpc}$), the radio luminosity at
10 GHz is $L_{\rm 10GHz}\sim 3\times 10^{37}\,{\rm erg\,s^{-1}}$ and the ratio
of the radio luminosity to the corresponding soft X-ray luminosity of the
flare is $L_{\rm 10GHz}/L_{X}=3\times 10^{37}/(10^{42})\sim 3\times 10^{-5}$,
hence the energy going into the radio emission is significantly smaller than
the detected X-ray luminosity. The radio flare could be related to the
synchrotron emission due to star-disc interaction shocks or the magnetic
reconnection event in the magnetized plasma of the hot accretion flow. For the
known QPPs, one can plan observational campaigns of simultaneous X-ray/radio
monitoring with the VLA, whose sensitivity threshold is as low as $10\,{\rm\mu
Jy}$ (RMS noise) at 10 GHz (bandwidth 500 MHz, field of view of 4.2 arcmin)
for the on-source time of 24.9 minutes 777Calculated using the on-source time
calculator https://obs.vla.nrao.edu/ect/.. In the second half of 2020s, the
Square Kilometer Array (SKA-1, mid-frequency) should reach the RMS noise level
of $1.2\,{\rm\mu Jy}$ for 1 hour-long observation at 12.5 GHz (bandwidth 3.75
GHz), which should reveal radio variability of even fainter QPPs. With the
field of view of 6.7 arcmin (SKA1-mid, 12.5 GHz) and a high sensitivity, the
SKA will survey the sky faster than previous radio arrays, making it possible
to visit the sky more frequently (Braun et al., 2019). For detecting radio
flares from previously unidentified QPPs, hence candidates for EMRIs, a large
sky coverage in the radio domain is desirable. The concept of the Canadian
Hydrogen Intensity Mapping (CHIME; CHIME Collaboration et al., 2022)
experiment with the total sky coverage of 31 000 ${\rm deg^{2}}$ for one-year
integration is in principle beneficial for detecting new radio transients such
as fast radio bursts. However, the confusion noise level of $\sim 0.1$ Jy
currently limits the instrument to much brighter transients than those
expected from CO-accretion flow interactions in nominally non-jetted hosts (at
least thousand times brighter than the Swift J0230-like radio eruption).

###  4.2 QPP luminosity and dependence on binary parameters

The dynamics of the putative binaries behind QPPs will depend on a number of
their properties, such as the mass of the perturbing COs, the properties of
the accretion disk, or the orbital parameters. Here we discuss how these will
affect the character and observability of the QPPs.

####  4.2.1 Radius of disk-perturber interaction and QPE luminosity

Using the results on the orbital frequency of EMRI detectable EM candidates
(see Fig. 3) we set the fiducial orbital period to $P_{\rm orb}=1$ hour (QPE-
like eruption occurs every $\sim 0.5$ hours for the CO-disc interaction). The
primary MBH mass is set to $M=10^{6}\,M_{\odot}$ unless stated otherwise. Then
the radius of the circular orbit of the CO is

$\displaystyle r=\left(\frac{Pc^{3}}{GM}\right)^{2/3}\,r_{\rm g}\sim
24\left(\frac{M}{10^{6}M_{\odot}}\right)^{-2/3}\left(\frac{P}{1\,\rm
hour}\right)^{1/3}r_{\rm g}\,,$ (13)

where $r_{\rm g}\equiv GM/c^{2}\sim 1.5\,{\rm km}(M/M_{\odot})$ is the
gravitational radius of the MBH.

The influence radius $R_{\rm inf}$, within which the compact object or the
perturber with the mass $\mu$ affects the surrounding gas gravitationally,
depends on the relative velocity with respect to the accretion disc material
$v_{\rm rel}$ and the disc sound speed $c_{\rm s}$ at the distance $r$, (e.g.
Suková et al., 2021)

$R_{\rm inf}=\frac{2G\mu}{v_{\rm rel}^{2}+c_{\rm s}^{2}}\,.$ (14)

Approximating both the disc velocity and the perturber velocity with Keplerian
radial profiles, and adopting the inclination $\iota$ between the perturber
and the accretion disc, we may estimate the relative velocity as follows,

$v_{\rm rel}=\left[\frac{2GM}{r}(1-\cos{\iota})\right]^{1/2}\,.$ (15)

The sound speed is calculated as $c_{\rm s}=(k_{\rm B}T_{\rm disc}/m_{\rm
p})^{1/2}$, where $k_{\rm B}$ is the Boltzmann constant, $T_{\rm disc}$ is the
temperature of the assumed standard disk, and $m_{\rm p}$ is the proton mass.

Specifically, we can expect the accretion disk sonic speed to be $c_{\rm
s}\lesssim v_{\rm k}\sim\sqrt{GM/r}$, where $v_{\rm k}$ is the Keplerian speed
at radius $r$. For a radiatively efficient accretion disk we will have $c_{\rm
s}/v_{\rm k}\ll 1$ and the height of the disk $H/r\sim c_{\rm s}/v_{\rm k}\ll
1$ (Shakura & Sunyaev, 1973), whereas for hot radiatively inefficient flows we
will have $H/r\sim c_{\rm s}/v_{\rm k}$ of order one (Yuan & Narayan, 2014).
We thus obtain $R_{\rm inf}$ as

$\displaystyle R_{\rm
inf}=r\frac{\mu}{M}\left(1-\cos\iota+\frac{1}{2}\left(\frac{H}{r}\right)^{2}\right)^{-1}$
(16)

This estimate for $R_{\rm inf}$ ignores the gravitational influence of the
primary MBH, which will dominate the gas outside the distance given by the
Hill (tidal) radius $R_{\rm Hill}$ from the CO, where

$R_{\rm Hill}\sim r\left(\frac{\mu}{3M}\right)^{1/3}\,.$ (17)

It is easy to see that for moderate inclinations and small mass ratios
$\mu/M$, the Hill radius will always be larger than $R_{\rm inf}$, so eq. (16)
provides a good estimate for the interaction cross-section. For low
inclinations $\sin\iota\lesssim(\mu/M)^{1/3}$ and
$(H/r)\lesssim(\mu/M)^{1/3}$, however, the Hill radius acts as a cap for the
radius affected by the perturber.

Additionally, we conjecture that a QPE will only be created by an orbiter that
is both not embedded in the disk, and punches through the disk at a relative
speed that is supersonic, both of which lead to the same condition
$\sin\iota\gtrsim H/r$.888Assuming perturbers co-rotating with the gas
direction. There is an edge case of counter-rotating orbits (formally
$\iota\sim 180^{\circ}$) which are embedded in the disk and sweeping through
it at highly supersonic speeds. For $\sin\iota\lesssim H/r$ the orbiter is
embedded in the disk and we surmise that this instead leads to a QPO.

For example, at the distance of the CO corresponding to one hour orbital
period, the sound speed of the accretion-disk gas is $\sim 35\,{\rm
km\,s^{-1}}$ for a fiducial relative accretion rate of $\dot{m}=0.1$ and an
innermost stable circular orbit set to $6\,r_{\rm g}$, while the Keplerian
velocity is $\sim 61\,000\,\rm km\,s^{-1}$. As a result, the orbit will be
non-embedded in the disk and cross it supersonically for $\iota\gtrsim
0.03^{\circ}$. Similarly, the orbit will not reach the cap on the influence
radius set by the Hill radius for $\iota\gtrsim
0.5^{\circ}(\mu/M_{\odot})^{1/3}(M/10^{6}M_{\odot})^{-1/3}$.

To estimate the X-ray luminosity of QPE-like eruptions, we then adopt a toy
model of CO-disc interactions, which lead to shocks and ejection of shocked
spherical, optically thick clouds above the disc plane (see e.g. Franchini et
al., 2023). For the initial cloud radius, we adopt $R_{\rm inf}$ and for its
temperature, we set $T_{\rm cloud}\sim 10^{6}\,{\rm K}$, which was inferred
for QPEs (Miniutti et al., 2019; Arcodia et al., 2021, 2024). Then the soft
X-ray luminosity (in the 0.2-2 keV band) can be estimated as follows,

$L_{\rm X}=4\pi R_{\rm inf}^{2}\int_{0.2\,{\rm keV}}^{2.0\,{\rm
keV}}\frac{2h\nu^{3}}{c^{2}}\frac{\mathrm{d}\nu}{\exp{(h\nu/k_{\rm B}T_{\rm
cloud})-1}}\,,$ (18)

where $\nu$ is the frequency corresponding to the photon energy $E_{\nu}=h\nu$
and $h$ is the Planck constant. We plot examples of the spectral energy
densities corresponding to this model in Figure 7.

Linial & Metzger (2023) present a similar CO-disk interaction model, taking
into account inefficient photon production during the optically thick phase of
an expanding shocked cloud. To reach $L_{X}\sim 10^{42}\,{\rm erg\,s^{-1}}$,
the CO mass should fall into the IMBH range with $\mu\sim 10^{4}\,M_{\odot}$,

$\displaystyle\begin{split}L_{\rm X}&=8.9\times
10^{41}\left(\frac{\dot{m}}{0.1}\right)^{1/3}\left(\frac{M}{10^{6}\,M_{\odot}}\right)^{5/9}\,\\\
&\times\left(\frac{\mu}{10^{4}\,M_{\odot}}\right)^{2/3}\left(\frac{P_{\rm
QPP}}{0.5\,\text{hours}}\right)^{-2/9}\,{\rm erg\,s^{-1}}.\end{split}$ (19)

We indicate the luminosity range according to Eq. (19) in Figure 7 using the
cyan rectangle, which corresponds to the CO masses of
$\mu=10-10^{4}\,M_{\odot}$.

Figure 7: Exemplary spectral energy distributions (SEDs) corresponding to
shocked clouds of material ejected during the passages of a lighter compact
object (CO) around the MBH. We assumed the orbital period of 1 hour and the
MBH mass of $10^{6}\,M_{\odot}$. To reach typical QPE luminosities of
$10^{42}\,{\rm erg\,s^{-1}}$ (red dotted line) in the 0.2-2 keV band (blue
shaded region), the initial cloud size should be $\sim 10^{11}\,{\rm cm}$,
which corresponds to $\sim 0.7$ gravitational radii of the MBH. The compact
objects that are able to eject such clouds fall into the range of $\sim
100-1000\,M_{\odot}$, with the inclination of $5-15$ degrees with respect to
the accretion disc, respectively. Solar and subsolar COs typically generate
X-ray flares of $\lesssim 10^{39}\,{\rm erg\,s^{-1}}$ (grey shaded region),
hence at least three orders of magnitude below the typical QPE luminosities.
The dashed black vertical line represents the peak of the SEDs. The white star
stands for the mean QPE luminosity inferred from seven detected QPE sources
($\log{(L_{\rm x}\,[\rm erg\,s^{-1}])}=42.3\pm 0.6$; see Table 1). The cyan
rectangle represents the X-ray luminosity range for COs of
$\mu=10-10^{4}\,M_{\odot}$ according to the CO-disk interaction model by
Linial & Metzger (2023).

####  4.2.2 Theoretical uncertainties of disk-perturber interaction

The Shakura & Sunyaev (1973) $\alpha$-disk used for various estimates above
corresponds only to the effective structure of an accretion disk averaged over
the orbital period, disk height, and azimuthal angle. As such, it cannot be
used for the dynamical analysis of its interactions with bodies embedded in
the disk (i.e., smaller than the disk height) and/or moving on the orbital
time-scale, since on those scales the structure is dominated by magneto-
rotational-instability-driven turbulence (see Balbus & Hawley (1991) and
especially the discussion in Balbus & Hawley (1998)).

Another difficult aspect of the interaction are the dynamics of the ejected
plasma leading to the observed lightcurve of the QPE. In some models the
ejected material is assumed to expand in the form of a “bubble” for some time
before becoming optically thin and releasing the radiation (Lehto & Valtonen,
1996; Ivanov et al., 1998; Pihajoki, 2016; Linial & Metzger, 2023). However,
the expansion implies significant adiabatic cooling of the plasma, which
requires involving other subtle mechanisms to explain the relatively high
$10^{6}\,\rm K$ effective temperatures observed in QPEs (Linial & Metzger,
2023). Additionally, the CO will eject the material in the form of a trailing
wake as seen in the Hoyle-Lyttleton scenario, which is likely incompatible
with the expanding bubble model. Furthermore, for some configurations, such as
lower-inclination and retrograde orbits, the break-out emission of a bow shock
associated with the supersonic CO can dominate the EM signal (Tagawa & Haiman,
2023). As such, both the action of the perturber on the disk, and the possible
back-reaction of the disk on the perturber is currently uncertain and will
likely be only resolved by first-principle simulations of the interaction of
the plasma of the accretion disk with the orbiting body.

Furthermore, the model of star-disk interactions described above and the
estimated X-ray luminosity rely on the standard geometrically thin and
optically thick accretion disc. However, most of the EMRI hosts in the local
Universe are expected to be low-luminosity sources. Such sources host hot
accretion flows that are advection-dominated and radiatively inefficient (Yuan
& Narayan, 2014). In this case, plasma is more diluted and warmer, hence the
shocks due to interactions with the CO are also weaker because of lower Mach
numbers. The ejected mass and the associated electromagnetic energy is also
lower due to a smaller surface density (Linial & Metzger, 2023). For the hot
flow, the variability could rather be driven by the accretion-rate variations
due to propagating density waves induced by the CO bow shock (see e.g. Suková
et al., 2021).

####  4.2.3 Mass of compact object and detectability

The mass range of orbiting COs in EMRIs/IMRIs is rather broad – from a few
Solar masses (neutron stars, white dwarfs, stellar black holes; Shapiro &
Teukolsky, 1983) to several thousand or even ten thousand in the case of
intermediate-mass black holes (IMBHs; Greene et al., 2020; Pasham et al.,
2024a). The range can be extended to lower, sub-solar masses if we admit the
occurrence of brown dwarfs in galactic nuclei (see, e.g., Amaro-Seoane, 2019;
Gourgoulhon et al., 2019; Zajaček, 2021) or hypothetical light primordial
black holes (see, e.g., Abbott et al., 2022; Bird et al., 2023).

In Fig. 7, we plotted spectral energy distributions (SEDs) of emitting shocked
clouds ejected from the accretion disc due to the passages of lighter compact
objects. To reach typical QPE X-ray luminosities of $L_{X}\sim 10^{42}\,{\rm
erg\,s^{-1}}$ (0.2-2.0 keV; red dotted line), the initial cloud size should be
$R_{\rm in}\sim 10^{11}\,{\rm cm}$. This corresponds to the influence radius
of a CO with mass $\mu\sim 10^{2}-10^{3}\,M_{\odot}$ orbiting around the MBH
at the inclination of $5-15$ degrees with respect to the accretion disc
(prograde sense), respectively. This range corresponds to rather low
inclinations. If the mutual inclination is increased so that the CO and the
disc are approximately perpendicular to each other, then the required CO mass
is in the IMBH range, $\mu\sim 29\,000\,M_{\odot}$, since the relative
velocity increases to $v_{\rm rel}\sim 86\,400\,{\rm km\,s^{-1}}$. If we keep
the same inclination of 15 degrees and decrease the perturber mass by one
order of magnitude, the influence radius decreases by one order and the X-ray
luminosity by two orders of magnitude, see Eqs. (14), (18), and Fig. 7.

For Solar-mass COs, the X-ray luminosity is typically $L_{X}\lesssim
10^{39}\,{\rm erg\,s^{-1}}$ (solid orange line in Fig. 7) for the MBH mass of
$\gtrsim 10^{5}\,M_{\odot}$ and inclination of $\gtrsim 5$ degrees, which
corresponds to the relative velocity of $\gtrsim 2500\,{\rm km\,s^{-1}}$ (the
orbital distance is $111\,r_{\rm g}$ for $M=10^{5}\,M_{\odot}$, one hour
orbital period). For the minimal inclination of $\iota\sim 0.5^{\circ}$ around
$M=10^{6}\,M_{\odot}$, the X-ray luminosity could reach $L_{X}\sim
10^{42}\,{\rm erg\,s^{-1}}$ assuming $T_{\rm cloud}\sim 10^{6}\,{\rm K}$,
which is, however, a special “grazing” setup. Using the CO-disk interaction
model of Linial & Metzger (2023), the luminosity estimate is $L_{X}\sim
2\times 10^{39}\,{\rm erg\,s^{-1}}$ for $M=10^{6}\,M_{\odot}$ and
$\dot{m}=0.1$. For the interactions of a Solar-type star with the disk, the
same model predicts $L_{\star}\sim 1.4\times 10^{42}\,{\rm erg\,s^{-1}}$.
However, Solar-type stars are tidally disrupted at $r\sim 47\,r_{\rm g}$,
hence they are not relevant for QPPs with $f_{\rm QPP}\sim 0.5$ mHz. On the
other hand, brown dwarfs with the mass of $\sim 0.08\,M_{\odot}$ and the
radius of $\sim 0.1\,R_{\odot}$ would disrupt at $\sim r_{\rm t}/r_{\rm g}\sim
11\,(R_{\star}/0.1R_{\odot})(M/10^{6}\,M_{\odot})^{-2/3}(\mu/0.08\,M_{\odot})^{-1/3}$,
hence the brown dwarf-disk interactions with $f_{\rm QPP}\sim 0.5$ mHz are
plausible. The X-ray flares would have luminosities one order of magnitude
below the QPE luminosities (Linial & Metzger, 2023),

$\displaystyle\begin{split}L_{\rm dwarf}&\sim 3\times
10^{41}\left(\frac{R_{\star}}{0.1\,R_{\odot}}\right)^{2/3}\left(\frac{M}{10^{6}\,M_{\odot}}\right)\,\\\
&\times\left(\frac{\dot{m}}{0.1}\right)^{1/3}\left(\frac{P_{\rm
QPP}}{0.5\,\text{hours}}\right)^{-2/3}\,{\rm erg\,s^{-1}}.\end{split}$ (20)

Hence, the prospects for detecting X-ray flares from orbiting and inclined
solar-mass and subsolar-mass COs are currently weak, but they can improve with
the new generation of sensitive X-ray imaging instruments, such as Athena and
Lynx. For instance, for the source at the redshift of $z=0.05$, the eruption
X-ray flux density for the Solar-mass CO is predicted to be $F_{\rm X}\sim
10^{39}\,{\rm erg\,s^{-1}}/[4\pi(222.3\,{\rm Mpc})^{2}]\sim 1.7\times
10^{-16}\,{\rm erg\,s^{-1}\,cm^{-2}}$, which is at the sensitivity limit of
Athena and Lynx for one hour exposure time, see Fig. 6. Solar and subsolar-
mass COs have characteristic influence radii smaller than one gravitational
radius for typical parameters ($M=10^{6}\,M_{\odot}$). Using Eq. (14), we
obtain $R_{\rm inf}\sim 7\times 10^{-4}\,r_{\rm g}$ for $\iota=15^{\circ}$,
and using Eq. (17), we get $R_{\rm Hill}\sim 0.2\,r_{\rm g}$. Hence, this
limits the overall strength of the accretion-disc perturbation and the
resulting emission modulation.

We can also make rough estimates for CO masses causing significant
perturbations of geometrically thick accretion disks based on the simulations
of Suková et al. (2021). In that paper it was found that in order to detect
significant quasiperiodic variability in the inflow (accretion) rate, the size
of the influence radius of the CO should be of the order of $R_{\rm inf}\sim
0.1-1.0\,r_{\rm g}$ for the CO distance of $\sim 10\,r_{\rm g}$. For
$M=10^{6}\,M_{\odot}$, one gravitational radius we obtain the required CO mass
in the range $\mu\sim 160-42\,000\,M_{\odot}$ (inclination between 5 to 90
degrees with respect to the disc midplane) to induce significant quasiperiodic
variability. We can also compare this with the estimate using the Hill radius
in eq. (17); the mass of the CO with the Hill radius equal to $1\,r_{\rm g}$
is

$\displaystyle\begin{split}\mu&=\frac{M^{3}}{P_{\rm
orb}^{2}}\frac{12\pi^{2}G^{2}}{c^{6}}\,,\\\ &\sim
221\,\left(\frac{M}{10^{6}\,M_{\odot}}\right)^{3}\left(\frac{P_{\rm
orb}}{1\,\text{hour}}\right)^{-2}\,M_{\odot}\,,\end{split}$ (21)

which is within the same order of magnitude as the estimate using $R_{\rm
inf}$.

Figure 8: Comparison of different dynamical timescales for COs of different
masses. Left: Timescales expressed in years as a function of the distance from
the MBH of $10^{6}\,M_{\odot}$ (in gravitational radii). Different dynamical
processes (merger, grinding, hydrodynamical drag) are listed in the legend.
Solid lines correspond to the CO mass of $100\,M_{\odot}$, while dashed lines
correspond to the CO mass of $1\,M_{\odot}$. The vertical solid line at the
distance of $\sim 24\,r_{\rm g}$ corresponds to the orbital period of one
hour. The hydrodynamical and griding timescales are estimated for the standard
disk with the relative accretion rate of $\dot{m}=0.1$ and the initial
inclination of the CO set to $15^{\circ}$. The stars represent the positions
of detected QPP sources (see Table 1) assuming that the CO orbital period is
twice as long as the measured QPP period. The distance of the CO is estimated
assuming a circular orbit around a $10^{6}\,M_{\odot}$ MBH. Right: Timescales
(in years) as a function of the CO mass (in Solar masses). The distance from
the MBH of $10^{6}\,M_{\odot}$ is fixed to $\sim 24\,r_{\rm g}$, which
corresponds to the orbital period of one hour. In this panel, we also depict
separately the hydrodynamical-drag and grinding timescales for standard and
hot accretion-flow discs adopting the relative accretion rates of 0.1 and
0.01, respectively.

Another issue with solar-mass and subsolar-mass compact objects is that they
will not merge with the MBH within $\sim 10$ years, even for distances of the
order of 10 gravitational radii. In Fig. 8 (left panel), we compare the merger
timescale using the evolution model of Peters & Mathews (1963) for a CO of
$100\,M_{\odot}$ (solid black line) with the one corresponding to
$1\,M_{\odot}$ CO (dashed black line). For the CO of $100\,M_{\odot}$, the
merger timescale is $\sim 10$ years for the mean distance of $\sim 24r_{\rm
g}$. We compare the merger timescales with the grinding (or alignment)
timescales (Syer et al., 1991, assuming the standard disk with $\dot{m}=0.1$
and the initial inclination of 15 degrees) as well as with the hydrodynamical
drag timescales (Narayan, 2000) when the CO is embedded within the disk. We
can notice that for heavier COs hydrodynamical timescales (alignment and drag)
are shorter by a factor 100. By the same factor, the merger timescale is also
shorter. The interaction of the embedded CO with the realistic accretion disk
is expected to be more complex than the inspiral motion given by the classical
hydrodynamic drag. Density waves can lead to both inward and outward type I
migration and eventually to the formation of “migration” traps (Paardekooper
et al., 2010; Bellovary et al., 2016; Metzger et al., 2022). However, these
dynamical effects are beyond the scope of the comparison of order-of-magnitude
estimates provided here.

The inspiral due to gravitational-wave emission dominates inside $\sim
35\,r_{\rm g}$ (merger timescales is shorter than the hydrodrag timescale),
while the merger timescale is shorter than the characteristic grinding
timescale inside $\sim 350\,r_{\rm g}$. In the right panel of Fig. 8, we fix
the distance of the CO to $\sim 24\,r_{\rm g}$, which corresponds to one-hour
orbital period for $10^{6}\,M_{\odot}$ MBH, and we show the dependency between
the characteristic timescales and the CO mass. We see that the merger
timescale is $\sim 10$ years for $\sim 100\,M_{\odot}$, while heavier COs
would merge too soon (before the LISA launch), while lighter, solar-mass COs
would merge on the timescale of $\sim 1000$ years. Similarly, it should be
noted that the peak gravitational-wave strain of EMRIs scales linearly with
$\mu$, so the GW detection distance of light COs will be smaller as compared
to heavy COs even if they were close to merger within the LISA window.

####  4.2.4 Eccentricity and disk precession

The M1 EMRI catalogue from Babak et al. (2017) used in Figs 1, 2 and 3
corresponds to the MBH mass function of Barausse (2012), with the assumption
that typical MBHs are rapidly spinning, $a_{\rm MBH}\sim 0.98M$, typical mass
of the secondary CO $\sim 10M_{\odot}$, and specific assumptions about the
cusp erosion, regrowth, and structure of the star cluster in the MBH
environments (see Babak et al., 2017, for more details). In particular, the M1
model assumes that the dominant formation channel for EMRIs is the already
mentioned dynamical capture in a “loss-cone” scenario (Hils & Bender, 1995;
Sigurdsson & Rees, 1997). This leads to inspirals that have eccentricities
$e\gtrsim 0.1$ during the entire course of their inspiral as can be seen from
Figs. 1 and 2, and only weak correlation of the orbital plane with the
direction of the central MBH rotation. In other words, the EMRIs in question
have generic inclinations and eccentricities as illustrated in Fig. 9. On the
other hand, the other already mentioned channels including the capture and
separation of a binary by the Hills mechanism (Miller et al., 2005), and disk-
assisted migration (Pan & Yang, 2021) lead to low-eccentricity EMRIs with low
inclinations in the latter cases and generic inclinations in the former.
Hence, it is important to discuss what these distributions of orbital
parameter imply for EM observations. Here we show in particular that generic
orbits will lead to subtle timing features that will make the matching of QPPs
to EMRIs nontrivial.

If the perturbing CO evolves along an eccentric trajectory, this introduces a
more complex two-periodic structure of the time-variability. With the leading
time-scale being the azimuthal orbital period $\Omega_{\phi}\sim 2\pi/P_{\rm
orb}$, the second time-scale is the pericentre precession frequency, which
reads to leading post-Newtonian and mass-ratio order (Robertson, 1938)

$\Omega_{\rm p}=\frac{3(2GM\pi)^{2/3}}{P_{\rm
orb}^{5/3}c^{2}(1-e^{2})}\sim\frac{0.13}{\rm
hour}\frac{1}{1-e^{2}}\left(\frac{M}{10^{6}\,M_{\odot}}\right)^{2/3}\left(\frac{P_{\rm
orb}}{1\,\text{hour}}\right)^{-5/3}\\!,$ (22)

where $e$ is the eccentricity of the orbit. Specifically, if one assumes that
the variability is caused by the CO punching through an accretion disk in the
equatorial plane of the primary, these intersections will occur at different
times and radii for an inclined eccentric trajectory with “beat” patterns in
the time domain. For example, taking a moderate eccentricity $e\sim 0.3$ and
the orbital period and the MBH mass as 1 hour and $10^{6}M_{\odot}$,
respectively, the pericenter precession will complete a cycle in roughly 44
hours, so the beat pattern will span $\sim 88$ eruptions. An example of such
an orbit and its intersection times and radii are plotted in Figs 9 and 10.

Figure 9: The intersections (red) of an eccentric and inclined free test
particle orbit (blue) through an accretion disk (orange) lying in the
equatorial plane of a Kerr black hole with spin $a=0.98$ (black). The orbit is
prograde with the orbital parameters semi-latus rectum $p=24r_{\rm g}$,
eccentricity $e=0.3$ and inclination $i=60^{\circ}$ (defined as in Schmidt,
2002). The interplay of eccentricity, pericenter precession, Lense-Thirring
precession, and other sub-leading relativistic effects lead to a complex
pattern. Figure 10: 30 intersection times and radii of the orbit from Fig. 9
where the times correspond to the fiducial MBH mass $M=10^{6}M_{\odot}$.

Nevertheless, the EM flares are not necessarily produced immediately after the
intersection, so identifying this secondary frequency may be subtle, as
already discussed in Section 4.2.2. In particular, since the eccentric orbits
intersect the disk at different radii and thus interact with a plasma slab
with different temperatures, densities, and geometrical thickness, the delays
between impact and radiation release will vary among impacts. This introduces
additional timing features into the model that depend on the disk properties,
as exemplified in the case of the fit of the parameters of the putative
eccentric MBH binary in OJ 287 (Dey et al., 2018; Zwick & Mayer, 2023).

Another non-trivial timing feature may emerge from the fact that the orbital
plane of the perturber will generally precess due to the spin-orbital
coupling. This is to leading post-Newtonian order and mass ratio given by the
Lense-Thirring frequency (Damour & Schäfer, 1988)

$\displaystyle\begin{split}&\Omega_{\rm LT}=\frac{8\pi^{2}GMa}{P_{\rm
orb}^{2}c^{3}(1-e^{2})^{3/2}}\\\ &\phantom{\Omega_{\rm LT}}\sim\frac{0.11}{\rm
hour}\frac{a}{(1-e^{2})^{3/2}}\left(\frac{M}{10^{6}\,M_{\odot}}\right)\left(\frac{P_{\rm
orb}}{1\,\text{hour}}\right)^{-2}\,,\end{split}$ (23)

where $a\in(-1,1)$ is the dimensionless spin of the MBH. For an accretion disk
aligned with the MBH spin this does not lead to significant effects in the
temporal spacing of the flares. However, if the accretion disk is not aligned
with the equatorial plane, the interplay of the differential precession of the
accretion disk (Bardeen & Petterson, 1975) and the precession of the orbit
will lead to additional intricate variability, a fact that was used by
Franchini et al. (2023) to explain some of the temporal properties of the QPEs
in GSN 069, eRO-QPE1, eRO-QPE2, and RX J1301.9+2747.

##  5 Conclusion

EMRI observations through the LISA GW detector will provide a unique probe for
fundamental physics in our Universe. Environment-rich formation channels of
EMRIs would potentially give rise to their EM counterparts. These may be
observable as soft X-ray ($\sim 10^{-1}$ keV) band, low-frequency phenomena
known as QPEs, and/or the less-abrupt QPOs, which emit semi-regular X-ray
signals from both quiescent and active galactic nuclei. In this paper, we
estimated the frequency band, in which a population of EMRIs in the Universe,
as described in the M1 EMRI catalog of Babak et al. (2017), will most likely
emit a QPE or QPO signal in the soft X-ray band “today” (in 2024). After
establishing this frequency band to be $\approx 0.46\pm 0.22$ mHz, we did a
converse analysis on the well-known QPO source, RE J1034+396, treating it as
an EMRI to estimate its CO’s mass to be $\approx 46M_{\odot}$, and calculated
a prior-predictive distribution on its LISA-band SNR with a best-fit value of
$\rho\approx 14$, highlighting it as a potential multi-messenger EMRI source.
Given the estimate on the frequency band, we explored the possibility of
observing LISA-EMRI relevant QPEs and QPOs in the EM spectrum through a
variety of current and upcoming detectors in the X-ray, UV, and even radio
bands. The impact of several EMRI parameters currently unaccounted for in our
estimates, such as smaller masses of the CO, and the effect of eccentricity
and inclination on the QPE/QPO flares, were also discussed. Our study thus
attempted to strengthen the connection between EMRIs and QPEs/QPOs,
highlighting them as potential multimessenger inference targets. If such
observations become feasible in the future, it would greatly enhance the
science output of the recently adopted LISA observatory.

## Acknowledgements

SK acknowledges the computational resources made accessible by the National
University of Singapore’s IT Research Computing group and the support of the
NUS Research Scholarship (NUSRS). MZ acknowledges the financial support of the
GAČR Junior Star grant no. GM24-10599M. Fig. 9 was generated with the help of
the KerrGeodesics package within the BHPToolkit (bhptoolkit.org). VW and SK
are grateful to the organizers for their invitation to the 1st Trieste meeting
on the physics of gravitational waves where the idea of this paper was
conceived.

## Data Availability

The results of this study are fully reproducible and are made available on a
Github repository called
EMRI_EM_Counterparts999https://github.com/perturber/EMRI_EM_Counterparts.
Babak et al. (2017)’s M1 EMRI catalog as modified by Pozzoli et al. (2023) can
be obtained by contacting the authors101010Federico Pozzoli:
<EMAIL_ADDRESS>but we have also made available the back evolution
trajectories data for our analysis in Section 2 on
Zenodo111111https://tinyurl.com/4s6jfrde.

## References

  * AXIS Time-Domain et al. (2023) AXIS Time-Domain T., et al., 2023, arXiv e-prints, p. arXiv:2311.07658
  * Abbott et al. (2022) Abbott R., et al., 2022, Phys. Rev. Lett., 129, 061104
  * Alston et al. (2014) Alston W. N., Markevičiūtė J., Kara E., Fabian A. C., Middleton M., 2014, Monthly Notices of the Royal Astronomical Society: Letters, 445, L16–L20
  * Alston et al. (2015) Alston W. N., Parker M. L., Markevič iūtė J., Fabian A. C., Middleton M., Lohfink A., Kara E., Pinto C., 2015, Monthly Notices of the Royal Astronomical Society, 449, 467
  * Amaro-Seoane (2018) Amaro-Seoane P., 2018, Living Reviews in Relativity, 21
  * Amaro-Seoane (2019) Amaro-Seoane P., 2019, Phys. Rev. D, 99, 123025
  * Amaro-Seoane et al. (2007) Amaro-Seoane P., Gair J. R., Freitag M., Miller M. C., Mandel I., Cutler C. J., Babak S., 2007, Classical and Quantum Gravity, 24, R113
  * Amaro-Seoane et al. (2015) Amaro-Seoane P., Gair J. R., Pound A., Hughes S. A., Sopuerta C. F., 2015, Journal of Physics: Conference Series, 610, 012002
  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017, Laser Interferometer Space Antenna (arXiv:1702.00786)
  * Amaro Seoane et al. (2021) Amaro Seoane P., et al., 2021, General Relativity and Gravitation, 54
  * Arcodia et al. (2021) Arcodia R., et al., 2021, Nature, 592, 704
  * Arcodia et al. (2024) Arcodia R., et al., 2024, arXiv e-prints, p. arXiv:2401.17275
  * Babak et al. (2010) Babak S., et al., 2010, Classical and Quantum Gravity, 27, 084009
  * Babak et al. (2017) Babak S., et al., 2017, Physical Review D, 95
  * Baghi et al. (2019) Baghi Q., Thorpe J. I., Slutsky J., Baker J., Canton T. D., Korsakova N., Karnesis N., 2019, Phys. Rev. D, 100, 022003
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  * Barack & Cutler (2004) Barack L., Cutler C., 2004, Physical Review D, 69
  * Barausse (2012) Barausse E., 2012, Mon. Not. Roy. Astron. Soc., 423, 2533
  * Barausse et al. (2014) Barausse E., Cardoso V., Pani P., 2014, Phys. Rev. D, 89, 104059
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Bellovary et al. (2016) Bellovary J. M., Mac Low M.-M., McKernan B., Ford K. E. S., 2016, ApJ, 819, L17
  * Ben-Ami et al. (2022) Ben-Ami S., et al., 2022, in den Herder J.-W. A., Nikzad S., Nakazawa K., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 12181, Space Telescopes and Instrumentation 2022: Ultraviolet to Gamma Ray. p. 1218105 (arXiv:2208.00159), doi:10.1117/12.2629850
  * Berry et al. (2019) Berry C. P. L., et al., 2019, The unique potential of extreme mass-ratio inspirals for gravitational-wave astronomy (arXiv:1903.03686)
  * Bian & Huang (2010) Bian W.-H., Huang K., 2010, Monthly Notices of the Royal Astronomical Society, 401, 507
  * Bird et al. (2023) Bird S., et al., 2023, Physics of the Dark Universe, 41, 101231
  * Boller, Th. et al. (2001) Boller, Th. Keil, R. Trümper, J. O’Brien, P. T. Reeves, J. Page, M. 2001, A&A, 365, L146
  * Bonetti & Sesana (2020) Bonetti M., Sesana A., 2020, Physical Review D, 102
  * Braun et al. (2019) Braun R., Bonaldi A., Bourke T., Keane E., Wagg J., 2019, arXiv e-prints, p. arXiv:1912.12699
  * Brenneman (2013) Brenneman L., 2013, Measuring the Angular Momentum of Supermassive Black Holes. Springer New York, doi:10.1007/978-1-4614-7771-6, https://doi.org/10.1007%2F978-1-4614-7771-6
  * Britzen et al. (2018) Britzen S., et al., 2018, MNRAS, 478, 3199
  * Britzen et al. (2023) Britzen S., Zajaček M., Gopal-Krishna Fendt C., Kun E., Jaron F., Sillanpää A., Eckart A., 2023, ApJ, 951, 106
  * Brunner et al. (2008) Brunner H., Cappelluti N., Hasinger G., Barcons X., Fabian A. C., Mainieri V., Szokoly G., 2008, A&A, 479, 283
  * Burrows et al. (2005) Burrows D. N., et al., 2005, Space Sci. Rev., 120, 165
  * CHIME Collaboration et al. (2022) CHIME Collaboration et al., 2022, ApJS, 261, 29
  * Chakraborty et al. (2021) Chakraborty J., Kara E., Masterson M., Giustini M., Miniutti G., Saxton R., 2021, The Astrophysical Journal Letters, 921, L40
  * Chakraborty et al. (2024) Chakraborty J., et al., 2024, arXiv e-prints, p. arXiv:2402.08722
  * Chaudhury et al. (2018) Chaudhury K., et al., 2018, Monthly Notices of the Royal Astronomical Society, 478, 4830
  * Chen et al. (2022) Chen X., Qiu Y., Li S., Liu F. K., 2022, The Astrophysical Journal, 930, 122
  * Chua & Cutler (2022) Chua A. J., Cutler C. J., 2022, Physical Review D, 106
  * Chua et al. (2021) Chua A. J. K., Katz M. L., Warburton N., Hughes S. A., 2021, Phys. Rev. Lett., 126, 051102
  * Collin & Zahn (1999) Collin S., Zahn J.-P., 1999, Ap&SS, 265, 501
  * Colpi et al. (2024) Colpi M., et al., 2024, arXiv e-prints, p. arXiv:2402.07571
  * Cornish & Crowder (2005) Cornish N. J., Crowder J., 2005, Phys. Rev. D, 72, 043005
  * Czerny et al. (2016) Czerny B., et al., 2016, A&A, 594, A102
  * Damour & Schäfer (1988) Damour T., Schäfer G., 1988, Nuovo Cimento B Serie, 101B, 127
  * Del Pozzo (2012) Del Pozzo W., 2012, Phys. Rev. D, 86, 043011
  * Dey et al. (2018) Dey L., et al., 2018, The Astrophysical Journal, 866, 11
  * Dey et al. (2021) Dey K., Karnesis N., Toubiana A., Barausse E., Korsakova N., Baghi Q., Basak S., 2021, Physical Review D, 104
  * Eckart et al. (2017) Eckart A., et al., 2017, Foundations of Physics, 47, 553
  * Franchini et al. (2023) Franchini A., et al., 2023, A&A, 675, A100
  * Fujita & Shibata (2020) Fujita R., Shibata M., 2020, Physical Review D, 102
  * Gair et al. (2013) Gair J. R., Vallisneri M., Larson S. L., Baker J. G., 2013, Living Reviews in Relativity, 16
  * Gaskin et al. (2019) Gaskin J. A., et al., 2019, Journal of Astronomical Telescopes, Instruments, and Systems, 5, 021001
  * Gierliński et al. (2008) Gierliński M., Middleton M., Ward M., Done C., 2008, Nature, 455, 369
  * Giustini et al. (2020) Giustini M., Miniutti G., Saxton R. D., 2020, A&A, 636, L2
  * Gourgoulhon et al. (2019) Gourgoulhon E., Le Tiec A., Vincent F. H., Warburton N., 2019, A&A, 627, A92
  * Gray et al. (2023) Gray R., et al., 2023, Joint cosmological and gravitational-wave population inference using dark sirens and galaxy catalogues (arXiv:2308.02281)
  * Greene et al. (2020) Greene J. E., Strader J., Ho L. C., 2020, ARA&A, 58, 257
  * Guolo et al. (2023) Guolo M., et al., 2023, arXiv e-prints, p. arXiv:2309.03011
  * Gupta et al. (2018) Gupta A. C., Tripathi A., Wiita P. J., Gu M., Bambi C., Ho L. C., 2018, A&A, 616, L6
  * Hasinger et al. (2001) Hasinger G., et al., 2001, A&A, 365, L45
  * Hills (1988) Hills J. G., 1988, Nature, 331, 687
  * Hils & Bender (1995) Hils D., Bender P. L., 1995, ApJ, 445, L7
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, ApJ, 629, 15
  * Ingram & Motta (2019) Ingram A. R., Motta S. E., 2019, New Astronomy Reviews, 85, 101524
  * Ingram et al. (2021) Ingram A., Motta S. E., Aigrain S., Karastergiou A., 2021, MNRAS, 503, 1703
  * Ivanov et al. (1998) Ivanov P. B., Igumenshchev I. V., Novikov I. D., 1998, ApJ, 507, 131
  * Jin et al. (2020) Jin C., Done C., Ward M., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3538
  * Karas & Šubr (2001) Karas V., Šubr L., 2001, A&A, 376, 686
  * Karas et al. (2021) Karas V., Svoboda J., Zajaček M., 2021, in RAGtime: Workshops on black holes and netron stars. p. E1 (arXiv:1901.06507), doi:10.48550/arXiv.1901.06507
  * Katz et al. (2021) Katz M. L., Chua A. J. K., Speri L., Warburton N., Hughes S. A., 2021, Phys. Rev. D, 104, 064047
  * Kejriwal et al. (2023) Kejriwal S., Speri L., Chua A. J. K., 2023, Impact of Correlations on the Modeling and Inference of Beyond Vacuum-GR Effects in Extreme-Mass-Ratio Inspirals (arXiv:2312.13028)
  * King (2023a) King A., 2023a, Monthly Notices of the Royal Astronomical Society: Letters, 523, L26
  * King (2023b) King A., 2023b, Monthly Notices of the Royal Astronomical Society: Letters, 526, L31
  * Kley & Nelson (2012) Kley W., Nelson R., 2012, Annual Review of Astronomy and Astrophysics, 50, 211–249
  * Kocsis et al. (2011) Kocsis B., Yunes N., Loeb A., 2011, Phys. Rev. D, 84, 024032
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Krtička et al. (2024) Krtička J., et al., 2024, Space Sci. Rev., 220, 24
  * Kulkarni et al. (2021) Kulkarni S. R., et al., 2021, arXiv e-prints, p. arXiv:2111.15608
  * Laghi et al. (2021) Laghi D., Tamanini N., Pozzo W. D., Sesana A., Gair J., Babak S., Izquierdo-Villalba D., 2021, Monthly Notices of the Royal Astronomical Society, 508, 4512
  * Lehto & Valtonen (1996) Lehto H. J., Valtonen M. J., 1996, ApJ, 460, 207
  * Levin (2007) Levin Y., 2007, MNRAS, 374, 515
  * Levin & Beloborodov (2003) Levin Y., Beloborodov A. M., 2003, ApJ, 590, L33
  * Li et al. (2023) Li R.-N., Zhao Z.-Y., Gao Z., Wang F.-Y., 2023, The Astrophysical Journal Letters, 956, L2
  * Lin et al. (2013) Lin D., Irwin J. A., Godet O., Webb N. A., Barret D., 2013, The Astrophysical Journal, 776, L10
  * Linial & Metzger (2023) Linial I., Metzger B. D., 2023, ApJ, 957, 34
  * Linial & Metzger (2024) Linial I., Metzger B. D., 2024, ApJ, 963, L1
  * Littenberg & Cornish (2023) Littenberg T. B., Cornish N. J., 2023, Phys. Rev. D, 107, 063004
  * Lops et al. (2023) Lops G., Izquierdo-Villalba D., Colpi M., Bonoli S., Sesana A., Mangiagli A., 2023, MNRAS, 519, 5962
  * Metzger et al. (2022) Metzger B. D., Stone N. C., Gilbaum S., 2022, The Astrophysical Journal, 926, 101
  * Miller et al. (2005) Miller M. C., Freitag M., Hamilton D. P., Lauburg V. M., 2005, ApJ, 631, L117
  * Miniutti et al. (2019) Miniutti G., et al., 2019, Nature, 573, 381
  * Miniutti et al. (2023a) Miniutti G., Giustini M., Arcodia R., Saxton R. D., Read A. M., Bianchi S., Alexander K. D., 2023a, A&A, 670, A93
  * Miniutti et al. (2023b) Miniutti G., Giustini M., Arcodia R., Saxton R. D., Chakraborty J., Read A. M., Kara E., 2023b, A&A, 674, L1
  * Moore et al. (2019) Moore C. J., Gerosa D., Klein A., 2019, Monthly Notices of the Royal Astronomical Society: Letters, 488, L94–L98
  * Mushotzky et al. (2019) Mushotzky R., et al., 2019, in Bulletin of the American Astronomical Society. p. 107 (arXiv:1903.04083), doi:10.48550/arXiv.1903.04083
  * Nandra et al. (2013) Nandra K., et al., 2013, arXiv e-prints, p. arXiv:1306.2307
  * Narayan (2000) Narayan R., 2000, ApJ, 536, 663
  * Paardekooper et al. (2010) Paardekooper S. J., Baruteau C., Crida A., Kley W., 2010, MNRAS, 401, 1950
  * Paczynski (1978) Paczynski B., 1978, Acta Astron., 28, 91
  * Pan & Yang (2021) Pan Z., Yang H., 2021, Physical Review D, 103
  * Pan et al. (2016) Pan H.-W., Yuan W., Yao S., Zhou X.-L., Liu B., Zhou H., Zhang S.-N., 2016, ApJ, 819, L19
  * Pasham et al. (2019) Pasham D. R., et al., 2019, Science, 363, 531
  * Pasham et al. (2024a) Pasham D. R., et al., 2024a, arXiv e-prints, p. arXiv:2402.10140
  * Pasham et al. (2024b) Pasham D. R., et al., 2024b, The Astrophysical Journal Letters, 963, L47
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Pihajoki (2016) Pihajoki P., 2016, MNRAS, 457, 1145
  * Pozzoli et al. (2023) Pozzoli F., Babak S., Sesana A., Bonetti M., Karnesis N., 2023, Computation of stochastic background from extreme mass ratio inspiral populations for LISA (arXiv:2302.07043)
  * Predehl et al. (2021) Predehl P., et al., 2021, A&A, 647, A1
  * Rauch (1995) Rauch K. P., 1995, MNRAS, 275, 628
  * Ray et al. (2019) Ray P. S., et al., 2019, arXiv e-prints, p. arXiv:1903.03035
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Remillard et al. (2022) Remillard R. A., et al., 2022, AJ, 163, 130
  * Reynolds (2013) Reynolds C. S., 2013, Classical and Quantum Gravity, 30, 244004
  * Robertson (1938) Robertson H. P., 1938, Annals of Mathematics, pp 101–104
  * Robson et al. (2019) Robson T., Cornish N. J., Liu C., 2019, Classical and Quantum Gravity, 36, 105011
  * Schmidt (2002) Schmidt W., 2002, Classical and Quantum Gravity, 19, 2743
  * Schutz (1986a) Schutz B. F., 1986a, Nature, 323, 310
  * Schutz (1986b) Schutz B. F., 1986b, Nature, 323, 310
  * Secunda et al. (2020) Secunda A., et al., 2020, The Astrophysical Journal, 903, 133
  * Semerák et al. (1999) Semerák O., Karas V., de Felice F., 1999, PASJ, 51, 571
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Shapiro & Teukolsky (1983) Shapiro S. L., Teukolsky S. A., 1983, Black holes, white dwarfs and neutron stars. The physics of compact objects (A Wiley-Interscience Publication, New York: Wiley), doi:10.1002/9783527617661.
  * Shvartzvald et al. (2023) Shvartzvald Y., et al., 2023, arXiv e-prints, p. arXiv:2304.14482
  * Sigl et al. (2007) Sigl G., Schnittman J., Buonanno A., 2007, Phys. Rev. D, 75, 024034
  * Sigurdsson & Rees (1997) Sigurdsson S., Rees M. J., 1997, MNRAS, 284, 318
  * Sikora et al. (2007) Sikora M., Stawarz Ł., Lasota J.-P., 2007, ApJ, 658, 815
  * Suková et al. (2021) Suková P., Zajaček M., Witzany V., Karas V., 2021, The Astrophysical Journal, 917, 43
  * Syer et al. (1991) Syer D., Clarke C. J., Rees M. J., 1991, MNRAS, 250, 505
  * Tagawa & Haiman (2023) Tagawa H., Haiman Z., 2023, MNRAS, 526, 69
  * Tanaka & Ward (2004) Tanaka H., Ward W. R., 2004, ApJ, 602, 388
  * Tchekhovskoy et al. (2010) Tchekhovskoy A., Narayan R., McKinney J. C., 2010, ApJ, 711, 50
  * Vallisneri (2009) Vallisneri M., 2009, Classical and Quantum Gravity, 26, 094024
  * Valtonen et al. (2016) Valtonen M. J., et al., 2016, ApJ, 819, L37
  * Vaughan (2005) Vaughan S., 2005, A&A, 431, 391
  * Volonteri et al. (2005) Volonteri M., Madau P., Quataert E., Rees M. J., 2005, The Astrophysical Journal, 620, 69
  * Wang et al. (2014) Wang J.-Z., Lei W.-H., Wang D.-X., Zou Y.-C., Zhang B., Gao H., Huang C.-Y., 2014, The Astrophysical Journal, 788, 32
  * Werner et al. (2024) Werner N., et al., 2024, Space Sci. Rev., 220, 11
  * Wevers et al. (2022) Wevers T., Pasham D. R., Jalan P., Rakshit S., Arcodia R., 2022, A&A, 659, L2
  * Xian et al. (2021) Xian J., Zhang F., Dou L., He J., Shu X., 2021, ApJ, 921, L32
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Yuan et al. (2022) Yuan W., Zhang C., Chen Y., Ling Z., 2022, in , Handbook of X-ray and Gamma-ray Astrophysics. p. 86, doi:10.1007/978-981-16-4544-0_151-1
  * Zajaček (2021) Zajaček M., 2021, ApJ, 915, 111
  * Zajaček et al. (2023) Zajaček M., et al., 2023, arXiv e-prints, p. arXiv:2306.15082
  * Zhang et al. (2017) Zhang P., Zhang P.-f., Yan J.-z., Fan Y.-z., Liu Q.-z., 2017, ApJ, 849, 9
  * Zhang et al. (2023) Zhang H., Yang S., Dai B., 2023, ApJ, 946, 52
  * Zhou et al. (2024) Zhou C., Huang L., Guo K., Li Y.-P., Pan Z., 2024, Probing orbits of stellar mass objects deep in galactic nuclei with quasi-periodic eruptions (arXiv:2401.11190)
  * Zwick & Mayer (2023) Zwick L., Mayer L., 2023, MNRAS, 526, 2754
  * eLISA Consortium et al. (2013) eLISA Consortium T., et al., 2013, The Gravitational Universe (arXiv:1305.5720)
  * van der Klis (2000) van der Klis M., 2000, ARA&A, 38, 717

##  Appendix A Additional priors on RE J1034+396 parameters

Here, for the QPO source RE J1034+396, we determine the priors on the MBH
dimensionless spin parameter, $a$, and the eccentricity and inclination of the
CO’s orbit, i.e. $e_{0}$ and $\iota$ respectively. We use a right skewed-
normal prior on $a$ with a mean of 0.9. This is based on evidence suggesting
that MBHs in AGNs are generally expected to be endowed with a large spin (see,
for example, Volonteri et al., 2005; Reynolds, 2013; Brenneman, 2013). Our
model assumes prograde orbits, i.e. orbits in the same direction as the MBH
spin. This is expected for EMRI sources, in which CO is either captured by the
accretion disk surrounding the MBH or forms in-situ within the disk and
migrates towards tighter orbits (see references for formation channels
(iii)-(iv) as described in Section 1, and Bellovary et al., 2016). For the
eccentricity and inclination of the CO’s orbit, we choose exponentially
decaying distributions with means 0.05 and 0.01 for eccentricity and
inclination respectively, allowing for only slightly inclined and eccentric
orbits. This builds on the evidence that small objects orbiting a massive body
in the presence of a disk lose their energy and angular momentum through
interactions with the disk, circularizing their orbits and decreasing their
inclination (Tanaka & Ward, 2004; see also Kley & Nelson, 2012 for a review).
Energy and angular momentum loss through GW radiation in EMRIs further assists
this circularization, such that environment-rich EMRIs can be expected to have
small eccentricities and inclinations (Pan & Yang, 2021). As noted in Section
3, the correlations of the estimated parameters ($\mu,\rho$) with the set
$\\{a,e_{0},\iota\\}$ is small $(|C_{ij}|\lesssim 0.02)$, and thus the choice
of prior does not change our results significantly. The full 8-dimensional
prior-predictive estimates on the RE J1034+396 parameter set
$\\{M,a,e_{0},\iota,f_{\rm orb,2007},\dot{f},\mu,\rm{SNR}\\}$ are provided in
Fig. 11.

Figure 11: Distributions of $\\{M,a,e_{0},\iota,f_{\rm orb,2007},\dot{f}_{\rm
orb,obs},\mu,\text{SNR}\\}$ of RE J1034+396 modeled as an extreme-mass-ratio
inspiral. $\\{M,f_{\rm orb,2007},\dot{f}_{\rm orb,obs}\\}$ are estimated from
current observational constraints on RE J1034+396 and distributions of
$\\{a,e_{0},\iota\\}$ are based on expected properties of EMRIs as described
in the text. The thus obtained prior-predictive distribution on $\mu$ is
described in the seventh panel, and the corresponding LISA-band SNR in the
final panel.