Prospects for EMRI/MBH parameter estimation using Quasi-Periodic Eruption timings: short-timescale analysis

TL;DR

This paper demonstrates that QPE timing analysis can effectively constrain key parameters of EMRIs and MBHs, providing a new method for probing galactic nuclei dynamics through Bayesian inference and analytic Kerr geodesics.

Contribution

It introduces QPE-FIT, a GPU-accelerated Bayesian inference framework using analytic Kerr geodesics to extract EMRI and MBH parameters from QPE timings.

Findings

Mild-eccentric EMRIs can constrain MBH mass to 10% within tens of orbits.

MBH spin remains unconstrained at larger semimajor axes and observation baselines.

Disk precession properties can be constrained within 10-50% despite parameter degeneracies.

Abstract

Quasi-Periodic Eruptions (QPEs) are luminous, recurring X-ray outbursts from galactic nuclei, with timescales of hours to days. While their origin remains uncertain, leading models invoke accretion disk instabilities or the interaction of a massive black hole (MBH) with a lower-mass secondary in an extreme mass ratio inspiral (EMRI). EMRI scenarios offer a robust framework for interpreting QPEs by characterizing observational signatures associated with the secondary's orbital dynamics. This, in turn, enables extraction of the MBH/EMRI physical properties and provides a means to test the EMRI scenario, distinguishing models and addressing the question: what can QPE timings teach us about massive black holes and EMRIs? In this study, we employ analytic expressions for Kerr geodesics to efficiently resolve the trajectory of the secondary object and perform GPU-accelerated Bayesian inference to assess the information content of QPE timings. Using our inference framework, referred to as QPE-FIT (Fast Inference with Timing), we explore QPE timing constraints on astrophysical parameters, such as EMRI orbital parameters and MBH mass/spin. We find that mild-eccentricity EMRIs ($e\sim0.1-0.3$) can constrain MBH mass and EMRI semimajor axis/eccentricity to the 10% level within tens of orbital periods, while MBH spin is unconstrained for the explored semimajor axes $\geq 100R_g$ and monitoring baselines $\mathcal{O}(10-100\rm)$ orbits. Introducing a misaligned precessing disk generally degrades inference of EMRI orbital parameters, but can constrain disk precession properties within 10-50%. This work both highlights the prospect of QPE observations as dynamical probes of galactic nuclei and outlines the challenge of doing so in the multimodal parameter space of EMRI-disk collisions.

Tables

Table 1. Table 1: Typical speed-up from GPU-accelerating one evaluation of the likelihood function. We used an NVIDIA A100 GPU, compared with an AMD EPYC 7763 256-core CPU operating with a 2652.814 MHz clock speed.

	100 walkers				1000 walkers
	CPU	GPU	Gain	Gain	CPU	GPU	Gain	Gain
Computing step	(sec)	(sec)	(wall time)	(per core)	(sec)	(sec)	(wall time)	(per core)
EMRI trajectory	0.232	0.0119	$19.5\times$	$4992\times$	2.37	0.0119	$199\times$	$50944\times$
Disk crossings	0.0257	0.00199	$12.9\times$	$3302\times$	0.352	0.0136	$25.9\times$	$6630\times$
Timing residuals	0.00395	0.00142	$2.78\times$	$712\times$	0.0439	0.00452	$9.71\times$	$2486\times$
Total log-likelihood	0.2622	0.0154	$17.0\times$	$\mathrm{𝟒𝟑𝟓𝟐}\times$	2.77	0.0263	$105\times$	$\mathrm{𝟐𝟔𝟖𝟖𝟎}\times$

Table 2. Table C.1: Recovered maximum a posteriori values and deviations from the true values for the semi-major axis, eccentricity, and black hole mass, comparing different observing inclination angles θ obs \theta_{\rm obs} .

Parameter	${\theta }_{\mathrm{obs}}=\pi /4$	${\theta }_{\mathrm{obs}}=\pi /2$	${\theta }_{\mathrm{obs}}=3\pi /4$
Semi-major axis $a$	103 (0.15$\sigma$)	103 (0.14$\sigma$)	95 (-0.23$\sigma$)
Eccentricity $e$	0.298 (-0.04$\sigma$)	0.315 (0.26$\sigma$)	0.300 (0.02$\sigma$)
BH Mass	4.98 (-0.15$\sigma$)	4.98 (-0.14$\sigma$)	5.03 (0.26$\sigma$)

Table 3. Table E.1: Average 16/50/84th percentile of posteriors, marginalized over other parameters, for our large grid of simulations.

		Static aligned disk		Misaligned precessing disk
Parameter	Observing strategy	True val.	Average est.	True val.	Average est.
EMRI semimajor axis ($R_g$)	$1\times 100$ ks	100	${114}_{-18}^{+16}$	100	${178}_{-73}^{+91}$
		200	${211}_{-29}^{+23}$	200	${312}_{-84}^{+100}$
		400	${405}_{-35}^{+30}$	400	${354}_{-122}^{+98}$
	$4\times 100$ ks (200 ks gaps)	100	${108}_{-16}^{+17}$	100	${187}_{-80}^{+102}$
		200	${209}_{-23}^{+21}$	200	${268}_{-86}^{+99}$
		400	${402}_{-25}^{+25}$	400	${342}_{-107}^{+94}$
	$4\times 100$ ks (1 Ms gaps)	100	${110}_{-20}^{+16}$	100	${227}_{-92}^{+123}$
		200	${206}_{-23}^{+20}$	200	${261}_{-96}^{+109}$
		400	${400}_{-24}^{+29}$	400	${295}_{-108}^{+115}$
EMRI eccentricity	$1\times 100$ ks	0.01	${0.0151}_{-0.0045}^{+0.0026}$	0.01	${0.031}_{-0.01}^{+0.03}$
		0.1	${0.105}_{-0.021}^{+0.011}$	0.1	${0.097}_{-0.03}^{+0.03}$
		0.3	${0.304}_{-0.027}^{+0.024}$	0.3	${0.278}_{-0.04}^{+0.04}$
		0.5	${0.526}_{-0.032}^{+0.020}$	0.5	${0.474}_{-0.06}^{+0.07}$
	$4\times 100$ ks (200 ks gaps)	0.01	${0.0104}_{-0.0020}^{+0.0044}$	0.01	${0.033}_{-0.01}^{+0.03}$
		0.1	${0.128}_{-0.017}^{+0.014}$	0.1	${0.167}_{-0.06}^{+0.05}$
		0.3	${0.319}_{-0.049}^{+0.092}$	0.3	${0.302}_{-0.05}^{+0.06}$
		0.5	${0.521}_{-0.043}^{+0.040}$	0.5	${0.512}_{-0.08}^{+0.10}$
	$4\times 100$ ks (1 Ms gaps)	0.01	${0.0323}_{-0.011}^{+0.016}$	0.01	${0.080}_{-0.04}^{+0.08}$
		0.1	${0.0989}_{-0.017}^{+0.091}$	0.1	${0.162}_{-0.07}^{+0.09}$
		0.3	${0.292}_{-0.028}^{+0.11}$	0.3	${0.294}_{-0.08}^{+0.11}$
		0.5	${0.44}_{-0.078}^{+0.084}$	0.5	${0.461}_{-0.10}^{+0.14}$
EMRI inclination (degrees)	$1\times 100$ ks	60	${58.9}_{-18.5}^{+20.2}$	60	${57.480}_{-24.90}^{+19.81}$
	$4\times 100$ ks (200 ks gaps)	60	${59.6}_{-19.3}^{+20.1}$	60	${56.280}_{-25.05}^{+20.24}$
	$4\times 100$ ks (1 Ms gaps)	60	${60.3}_{-19.3}^{+19.0}$	60	${52.845}_{-23.84}^{+20.73}$
MBH spin	$1\times 100$ ks	0.5	${0.478}_{-0.32}^{+0.33}$	0.5	${0.505}_{-0.31}^{+0.31}$
		0.9	${0.512}_{-0.34}^{+0.31}$	0.9	${0.507}_{-0.31}^{+0.32}$
	$4\times 100$ ks (200 ks gaps)	0.5	${0.471}_{-0.31}^{+0.31}$	0.5	${0.496}_{-0.31}^{+0.31}$
		0.9	${0.497}_{-0.32}^{+0.30}$	0.9	${0.492}_{-0.31}^{+0.32}$
	$4\times 100$ ks (1 Ms gaps)	0.5	${0.462}_{0.34}^{+0.31}$	0.5	${0.501}_{-0.30}^{+0.30}$
		0.9	${0.43}_{-0.28}^{+0.33}$	0.9	${0.486}_{-0.28}^{+0.30}$
log(MBH mass/$M_{\odot }$)	$1\times 100$ ks	5	${4.96}_{-0.058}^{+0.075}$	5	${4.870}_{-0.20}^{+0.29}$
		6	${5.93}_{-0.086}^{+0.14}$	6	${5.680}_{-0.23}^{+0.26}$
		7	—	7	—
	$4\times 100$ ks (200 ks gaps)	5	${4.97}_{-0.055}^{+0.065}$	5	${4.885}_{-0.23}^{+0.30}$
		6	${5.97}_{-0.055}^{+0.064}$	6	${5.880}_{-0.23}^{+0.29}$
		7	${6.86}_{-0.14}^{+0.22}$	7	${6.060}_{-0.22}^{+0.31}$
	$4\times 100$ ks (1 Ms gaps)	5	${4.97}_{-0.059}^{+0.07}$	5	${4.881}_{-0.24}^{+0.34}$
		6	${5.92}_{-0.10}^{+0.12}$	6	${5.818}_{-0.30}^{+0.36}$
		7	${6.92}_{-0.13}^{+0.17}$	7	${6.332}_{-0.29}^{+0.33}$
Disk inclination (degrees)	$1\times 100$ ks	5	—	5	${8.698}_{-3.61}^{+5.52}$
		20	—	20	${21.716}_{-7.27}^{+6.33}$
	$4\times 100$ ks (200 ks gaps)	5	—	5	${8.232}_{-3.99}^{+3.49}$
		20	—	20	${20.720}_{-7.93}^{+6.68}$
	$4\times 100$ ks (1 Ms gaps)	5	—	5	${13.108}_{-5.95}^{+7.81}$
		20	—	20	${22.822}_{-8.86}^{+8.57}$
Disk precession period ($T_{\mathrm{disk}}/P_{\mathrm{orb}}$)	$1\times 100$ ks	5	—	5	${46.659}_{-19.31}^{+23.79}$
		20	—	20	${46.053}_{-18.70}^{+26.64}$
		50	—	50	${63.156}_{-23.38}^{+23.17}$
	$4\times 100$ ks (200 ks gaps)	5	—	5	${46.190}_{-20.93}^{+27.24}$
		20	—	20	${45.028}_{-20.32}^{+27.17}$
		50	—	50	${56.790}_{-19.85}^{+23.39}$
	$4\times 100$ ks (1 Ms gaps)	5	—	5	${47.691}_{-21.52}^{+26.43}$
		20	—	20	${45.416}_{-20.14}^{+25.79}$
		50	—	50	${57.260}_{-23.39}^{+23.83}$

Equations

$$d{s}^2$$

$$+\frac{1-z^2}{\Sigma }\left(2{\chi }_{\bullet }^{2}r(1-z^2)+({\chi }_{\bullet }^2+r^2)\Sigma \right)d{\varphi }^2$$

$$-\frac{4{\chi }_{\bullet }r(1-z^2)}{\Sigma }dtd\varphi ,$$

$$r(\lambda )$$

$$z(\lambda )$$

$$k_r$$

$$k_z$$

$$q_r$$

$${\left(\frac{dr}{d\lambda }\right)}^2$$

$$-\Delta [r^2+({\mathcal{L}}_z-{\chi }_{\bullet }\mathcal{E}{)}^2+\mathcal{Q}],$$

$${\left(\frac{dz}{d\lambda }\right)}^2$$

$$-\left[(1-z^2)(1-{\mathcal{E}}^2){\chi }_{\bullet }^2+{\mathcal{L}}_z^2\right]z^2.$$

$$\frac{d\varphi }{d\lambda }$$

$$+\frac{{\mathcal{L}}_z}{1-z^2}-{\chi }_{\bullet }\mathcal{E},$$

$$\frac{dt}{d\lambda }$$

$$-a^2\mathcal{E}(1-z^2)+{\chi }_{\bullet }{\mathcal{L}}_z.$$

$$r(\lambda )$$

$$\cos \theta (\lambda )$$

$$t(\lambda )$$

$$\varphi (\lambda )$$

$$t^{(A)}(\lambda )$$

$${\varphi }^{(A)}(\lambda )$$

$$\displaystyle\hat{d}(t)=\begin{bmatrix}\sin\theta_{\rm disk}\cos\Big{(}2\pi t/T_{\rm disk}+\phi_{{\rm disk},0}\Big{)}\\ \sin\theta_{\rm disk}\sin\Big{(}2\pi t/T_{\rm disk}+\phi_{{\rm disk},0}\Big{)}\\ \cos\theta_{\rm disk}\end{bmatrix},$$

$$D(t)$$

$$\displaystyle+\sin\theta_{\rm disk}\sin\bigg{(}2\pi t/T_{\rm disk}+\phi_{{\rm disk},0}\bigg{)}y(t)$$

$$+\cos {\theta }_{\mathrm{disk}}z(t).$$

$${\Delta }_{\mathrm{Shap}}$$

$${\hat{n}}_{\mathrm{obs}}=\sin {\theta }_{\mathrm{obs}}\hat{x}+\cos {\theta }_{\mathrm{obs}}\hat{z}.$$

$${\Delta }_{\mathrm{geom}}=-{\vec{r}}_{\mathrm{cross}}\cdot {\hat{n}}_{\mathrm{obs}}.$$

$$\vec{\xi }=\{$$

$$q_{r,0},q_{z,0},q_{\varphi ,0},{\varphi }_{\mathrm{disk},0}\},$$

$$r_i(\vec{\xi },t_{\mathrm{obs},i})=\left|t_{\text{obs},i}-t_{\text{model},i}(\vec{\xi })\right|,$$

$${\chi }^2(\vec{\xi })=\sum _j^{N_{\text{obs}}}\sum _{i=1}^{n_{\text{qpe}},j}{\left(\frac{r_i(\vec{\xi })}{{\sigma }_i}\right)}^2.$$

$${\mathcal{L}}_z(\vec{\xi })\propto \exp \left(-\frac{1}{2}{\chi }^2(\vec{\xi })\right),$$

$$T_{\rm aps}\approx\frac{1}{3}\bigg{(}\frac{a}{R_{g}}\bigg{)}(1-e^{2})P_{\rm orb}.$$

$$\Phi (M_{\bullet },z)=\frac{dN}{d\log M_{\bullet }dV},$$

Full text

††thanks: Both authors contributed equally to this work††thanks: Both authors contributed equally to this work††thanks: NHFP Einstein Fellow

1 Introduction

Quasi-periodic eruptions (QPEs) are a recently identified class of high-amplitude soft X-ray flares recurring on timescales of hours to days, typically observed from galactic nuclei hosting low-mass massive black holes (MBHs). The known QPE sample has grown to around eleven sources (Miniutti et al., 2019; Giustini et al., 2020; Arcodia et al., 2021, 2024a; Chakraborty et al., 2021, 2025a; Quintin et al., 2023; Nicholl et al., 2024; Hernández-Garc´ıa et al., 2025), enabled primarily by high-cadence X-ray monitoring from XMM-Newton and NICER. These systems exhibit repeating outbursts with recurrence times that are generally stable over many cycles. The flare energetics and recurrence behavior point to a compact, coherent engine that is both localized (operating in the immediate environment of the MBH) and temporally stable across multiple eruption episodes. However, the physical origin of QPEs remains under active investigation, with three classes of theoretical frameworks emerging in recent years.

One class of models attributes QPEs to intrinsic variability in the accretion flow. Proposed mechanisms include radiation-pressure instabilities and the relativistic precession of a warped or geometrically thick inner disk around a spinning MBH (Raj & Nixon, 2021; Pan et al., 2021a, 2022, 2023; Kaur et al., 2023; Śniegowska et al., 2020, 2023; Middleton et al., 2025). These scenarios can naturally produce variable accretion flows, but face challenges explaining the QPE recurrence times, spectral evolution, and long viscous timescales compared to the observed eruption durations. Continued investigation is nevertheless warranted to further explore QPE observables using disk instability models.

A second class of interpretations focuses on mass transfer in extreme mass-ratio inspirals (EMRIs). EMRIs are gravitationally bound binaries which consist of a stellar-mass object (e.g. a star or stellar-remnant black hole) orbiting a supermassive black hole (MBH), with the inspiral driven predominantly by gravitational-wave emission. We refer to the stellar-mass object as the “secondary”. In this class of models, the eruptions are powered by episodic accretion from the secondary onto the MBH (King, 2020, 2022, 2023; Metzger et al., 2022; Zhao et al., 2022; Chen et al., 2022; Krolik & Linial, 2022; Linial & Sari, 2023; Yang et al., 2025; Lau & Yu, 2025). These models naturally yield periodicity tied to the orbital motion, but still require rapid accretion compared to the viscous timescale, raising questions about the feasibility of this scenario. Additionally, these models often produce clocks that are too stable, in contrast to the slight quasi-periodicity and frequency wandering observed in some QPE sources.

A third possibility, which is gaining increasing attention, proposes that QPEs are powered by repeated collisions between a compact secondary and the accretion disk of the MBH. In this picture, the secondary plunges through the disk typically twice per orbit, generating shocks, outflows, or thermal transients that give rise to the observed X-ray flares (Dai et al., 2010; Xian et al., 2021; Suková et al., 2021; Lu & Quataert, 2023; Franchini et al., 2023; Tagawa & Haiman, 2023; Linial & Metzger, 2023, 2024; Zhou et al., 2024a, 2025b, b, 2025a; Tsz-Lok Lam et al., 2025; Yao et al., 2025; Yao & Quataert, 2025; Huang et al., 2025). This model combines the advantages of a stable orbital clock from EMRIs with a physically motivated and localized energy release mechanism. Crucially, the eruption timings correspond to the approximately geodesic motion of the secondary and are largely independent of viscous timescales, allowing for a robust link between the flare cadence and underlying orbital dynamics. We caution that the EMRI collisional model is far from confirmed as the true physical origin of QPEs—see e.g. Mummery (2025), Guo & Shen (2025), and Guolo et al. (2025b) for further discussion on the shortcomings of these models in adequately producing the energetics, temperatures, and scaling relations predicted by the simple-case pictures. Additional tests across a variety of observables, e.g. timing residuals (Chakraborty et al., 2024; Arcodia et al., 2024b; Pasham et al., 2024; Giustini et al., 2024; Miniutti et al., 2025; Zhou et al., 2024b) and spectral/SED analysis (Wevers et al., 2025; Guolo et al., 2025b; Chakraborty et al., 2025b; Guo & Shen, 2025) are well-motivated, and are indeed an increasing focus in the literature.

In the EMRI model scenario, some studies also propose that QPEs are formed when the orbiter begins interacting with the compact accretion disks formed in the aftermath of tidal disruption events (TDEs) (Linial & Metzger, 2023; Franchini et al., 2023). In this EMRI+TDE framework, the secondary predominantly inspirals via GW emission and two-body scattering, until it reaches $\mathcal{O}(100R_{g})$, where it intersects the TDE disk. This scenario is especially compelling for QPEs where a TDE was explicitly observed several years prior to the onset of the eruptions (Quintin et al., 2023; Nicholl et al., 2024; Bykov et al., 2025; Chakraborty et al., 2025a), and is bolstered by mounting evidence for long-lived precursor flares in other systems without a definitive TDE association (Chakraborty et al., 2021; Miniutti et al., 2023; Arcodia et al., 2024a; Guolo et al., 2025a).

QPEs may represent the first electromagnetic detection of EMRIs, one of the key anticipated gravitational-wave (GW) sources targeted by the upcoming space-based millihertz observatories such as LISA (Amaro-Seoane et al., 2023) and TianQin (Luo et al., 2016). EMRIs offer extremely precise tests of strong-field gravity and detailed probes of massive black hole environments, motivating their importance in the LISA science case (Barack & Cutler, 2004). It has been proposed that “wet” EMRIs (those formed in the presence of an accretion disk) may constitute a common formation channel for LISA sources (Pan et al., 2021b; Wang et al., 2024), as their inward migration rate may be significantly enhanced via dissipative interactions with the accretion flow. EMRIs formed in these matter-rich environments can carry imprints of the accretion disk in their gravitational wave signals (Derdzinski et al., 2021; Speri et al., 2023; Copparoni et al., 2025; Duque et al., 2025), and may also produce quasi-periodic electromagnetic signatures during their inspiral.

Wet EMRIs have therefore been proposed as a possible source of periodic soft X-ray variability at millihertz frequencies, including quasi-periodic oscillations (QPOs) and QPEs, and may enable multi-messenger observations of individual systems (Kejriwal et al., 2024; Suzuguchi et al., 2025). While the population of thus-far known QPEs is at periods too long to be detectable by mHz observatories (Kara & Garc´ıa, 2025), a subset of LISA-detectable EMRIs may emit X-rays in the same frequency band as observed QPOs (Gierliński et al., 2008; Masterson et al., 2025), making them possible multi-messenger candidates. There is therefore growing potential for current and future X-ray observatories to identify EMRIs in advance of LISA’s launch.

In principle, the EMRI collisional model encodes information about the orbiter-disk system via the QPE timings (assuming a constant delay between collision and peak light). In this work, we address the question: how much information about the central MBH and the orbit of the secondary can be extracted from QPE timing data alone? We assess how different observing strategies, as well as different astrophysical parameters, impact inference. Using an efficient Post-Newtonian expansion for describing Kerr geodesics and including the impact of disk precession and geometric/Shapiro delays, we develop a computationally optimized, GPU-accelerated framework for simulating the timing of eruptions and assessing the prospects for astrophysical inference from existing and future observations. We make our QPE timing calculation code, dubbed QPE-FIT (Fast Inference with Timing), publicly available on GitHub111https://github.com/joheenc/QPE-FIT, to support future observational data fitting and serve as a resource for the community. QPE-FIT will undergo active development, and the most up-to-date features will be described at the provided link.

In Section 2 we describe the features of our model, including the EMRI orbital trajectory calculation, disk prescription, QPE timing method, and likelihood function for Bayesian inference. In Section 3 we show the results of applying our inference algorithm across a large grid of simulated QPE timings to assess how EMRI/disk properties affect parameter inference efficacy, and fit the timing data of eRO-QPE1 to show that the precessing disk model is able to describe the erratic timing behaviors observed in some QPEs. We discuss the possibility of further astrophysical inference from QPE timings in Section 4, and make concluding remarks in Section 5.

2 Methods

In this paper, we design our timing model to achieve the following objectives:

Develop an efficient framework for QPE source analysis that can be applied broadly and evaluated quickly. This is accomplished by (i) using a Post-Newtonian (PN) model for Kerr geodesics, which is computationally efficient and described in detail in Sec. 2.1; and (ii) using GPU-acceleration. Evaluation speed is crucial to allow Bayesian parameter inference, which typically entails $>$millions of model evaluations. We note that the PN formalism provides limited accuracy for low orbital separations ($\lesssim 50R_{g}$) and large eccentricities ($e\gtrsim 0.5$); we explore its limitations in Appendix A. 2. 2.

Create a flexible model that accounts for diverse QPE timing behaviors, including sources with irregular timing patterns. To achieve this, we model flare timings as interactions between the secondary object and a rigidly precessing accretion disk. The disk precession model is discussed in 2.2. 3. 3.

Coherently resolve flare timings while accounting for observational gaps. This is achieved by constructing the flare timing likelihood to match observed flares to model timings, enabling accurate parameter estimation even when large sections of the dataset are omitted. We use this framework to evaluate different observing strategies and assess their relative abilities to recover EMRI and disk parameters. We note that in this paper, we only model short observational gaps/baselines, of order tens to hundreds of consecutive QPEs, and defer long-term monitoring including secular evolution effects to future work.

This approach ensures our timing model is both adaptable and practical for analyzing QPE sources with varied timing behaviors and observational challenges. We build heavily on the existing work of Franchini et al. (2023), with the key addition that our computational accelerations allow for full Bayesian inference. To illustrate the efficiency gains, we benchmark the speed of one evaluation of the log-likelihood function with and without a GPU. We used an NVIDIA A100 GPU and an AMD EPYC 7763 256-core processor operating with a 2652.814 MHz clock speed; the resulting wall-clock times of the computation steps are shown in Table 1. We performed two separate benchmarks using 100 and 1000 parallel walkers. As is typical for GPU-accelerated programs, the best improvements are seen for higher parallelism, with a speed-up of $26880\times$ comparing one GPU to one CPU for the 1000-walker case. These performance gains enable rapid posterior exploration using ensemble samplers, underscoring the practical utility of QPE-FIT for future population-level studies.

2.1 Trajectory of the orbiting body

The trajectory of the secondary object in the QPE timing model is governed by bound geodesics in Kerr spacetime. The Kerr metric, describing a rotating black hole with spin parameter $\chi_{\bullet}$, is expressed in Boyer-Lindquist coordinates $\{t,r,\theta,\phi\}$ as follows:

[TABLE]

where $z=\cos\theta$, $\Delta=r(r-2)+\chi_{\bullet}^{2}$ and $\Sigma=r^{2}+\chi_{\bullet}^{2}z^{2}$. In the above equations, the black hole’s mass is normalized to unity ($M_{\bullet}=1$).222Throughout this work, $\chi_{\bullet}$ denotes the dimensionless spin of the massive black hole, with the mass normalized to $M_{\bullet}=1$. We restore physical units, including the mass scale, when computing observables such as flare timings. For reference, $GM_{\bullet}/c^{3}=49.2\,{\rm sec}(M_{\bullet}/10^{7}M_{\odot})$.

The Kerr metric admits four constants of motion that facilitate analytic solutions for geodesic motion. The first conserved quantity, $u_{\mu}u^{\mu}$, arises from the definition of proper time, $\tau$. The specific energy $\mathcal{E}$ and axial angular momentum $\mathcal{L}_{z}$ are conserved due to the spacetime’s Killing vector fields, while the existence of a Killing tensor leads to the Carter constant $\mathcal{Q}$. In the Schwarzschild limit, $\mathcal{Q}$ is simply related to the total angular momentum of an orbit: $\mathcal{Q}(\chi_{\bullet}=0)=|\boldsymbol{\mathcal{L}}|^{2}-\mathcal{L}_{z}^{2}$. For $a_{\bullet}\neq 0$, this correspondence does not hold precisely, but provides useful intuition about the meaning of $\mathcal{Q}$, especially for orbits with wide separation.

The radial and polar motion of the geodesic can be described analytically using elliptic functions Schmidt (2002); Fujita & Hikida (2009); van de Meent (2020):

[TABLE]

Here, $r_{1}$, $r_{2}$, and $z_{1}$ are the radial and polar turning points, while $k_{r}$ and $k_{z}$ are given by

[TABLE]

The roots $r_{3}$, $r_{4}$ and $z_{2}$ are determined by the radial and polar potentials, as described by Eqs. (13)–(14) in van de Meent (2020). The motion is parameterized in terms of the Mino time parameter $\lambda=\int d\tau/\Sigma$ (Mino, 2003) which decouples the radial and polar motion, thus simplifying the geodesic equations.

The functions $\mathsf{K}(k)$ and $\text{sn}(u|k)$ represent the complete elliptic integral of the first kind and the Jacobi elliptic sine function, respectively. The phases $q_{r}$ and $q_{z}$ evolve linearly with Mino time as:

[TABLE]

where $\Upsilon_{r}$ and $\Upsilon_{z}$ are the Mino-time frequencies, and $q_{r,0}$ and $q_{z,0}$ are initial phases. Each bound Kerr geodesic can be uniquely mapped to a set of quasi-Keplerian orbital parameters that offer a more intuitive description of the motion. These include the black hole spin $\chi_{\bullet}$, semi-latus rectum $p$, eccentricity $e$, inclination angle $i$, and three phase constants ($q_{r,0},q_{z,0},q_{\phi,0}$) that specify the initial conditions. This seven-dimensional parameterization fully characterizes the geodesic and forms the basis of our trajectory model.

While the orbital inclination $i$ provides an intuitive measure of tilt relative to the black hole spin, it is important to note that generic bound geodesics in Kerr spacetime do not lie in a fixed plane. Instead, the motion fills a two-dimensional torus due to independent radial and polar oscillations. The inclination angle $i$ used here refers to a conserved parameter derived from the orbital constants of motion via $\cos i=\mathcal{L}_{z}/\sqrt{\mathcal{L}_{z}^{2}+Q}$ and should not be interpreted as the geometric angle between a unique orbital plane and the black hole spin axis.

The first derivatives of $r$ and $z$ with respect to Mino time are determined by the potentials $R(r)$ and $Z(z)$:

[TABLE]

The geodesic motion also includes azimuthal and temporal components, which are expressed as:

[TABLE]

The temporal component relates Mino time $\lambda$ to the coordinate time $t$.

While the exact Kerr geodesic equations presented above enable high-precision modeling of the secondary’s trajectory, their numerical implementation can be computationally intensive, particularly when incorporated into a Bayesian inference pipeline which requires repeated trajectory evaluations. The bottleneck arises from the need to evaluate special functions such as elliptic integrals which can be slow to compute.

To mitigate this cost, we employ the post-Newtonian (PN) expansion of Kerr geodesics. In particular, we use the 3PN Fourier-series eccentricity expansions for bound orbits derived by Sago & Fujita (2015), which express the coordinate motion $r(\lambda)$, $\theta(\lambda)$, $\phi(\lambda)$, and $t(\lambda)$ as analytic functions of Mino time, without requiring computation of special functions. These expansions exploit the integrable structure of Kerr geodesics and are valid in the weak-field limit for moderate eccentricities and inclinations. This post-Newtonian framework offers orders-of-magnitude speedup over the exact method, making it well-suited for use in Bayesian inference algorithms. Each coordinate is given by a Fourier series expansion:

[TABLE]

where $Y=\cos{i}$ and the coefficients $t^{(A)},\phi^{(A)}$ are given by:

[TABLE]

The full expressions for the Fourier coefficients $\alpha_{n_{r}}$, $\beta_{n_{\theta}}$, $\tilde{t}_{n_{A}}^{(\{r,\theta\})}$, $\phi^{(\{r,\theta\})}(\lambda)$, and Mino-time frequencies $\Upsilon_{\{t,r,\theta,\phi\}}$, in terms of the orbital parameters $\{p,e,\chi_{\bullet},i\}$, are given in Appendices A and B of Sago & Fujita (2015) up to sixth order in $e$.

The increased computational efficiency comes at the cost of accuracy in certain regions of parameter space, particularly for highly eccentric orbits or those deep in the strong-field regime (i.e., small semi-latus rectum $p$). The PN formalism remains sufficiently accurate for orbits with $p\gtrsim 50$ and moderate eccentricity $\lesssim 0.5$, which encompass the parameter range relevant for most QPEs. We explore the accuracy of our 3PN computations, and the limits in which their validity breaks down, in Appendix A.

2.2 Rigidly precessing disk

Equipped with the orbital trajectory of the secondary, the natural next step is to compute its crossing times for a static accretion disk, under the simplifying assumption that it is aligned with the MBH spin axis. In fact, this was explored in several early QPE timing studies (Xian et al., 2021; Zhou et al., 2024a, b). However, after the first few observational studies of QPE timings, it quickly became evident that the simple-case static disk model is an incomplete description in at least some QPEs. Some systems display significant, unexpected irregularity in their burst timings (Arcodia et al., 2022; Giustini et al., 2024) not immediately reconciled with model predictions. Closer looks in Chakraborty et al. (2024) and Miniutti et al. (2025) observed super-periodic structure in the QPE arrival times of eRO-QPE1 and GSN 069, respectively, suggesting the presence of a further modulating timescale not present in the static disk model.

One promising avenue to explain these super-orbital modulations is rapid disk precession. As QPEs are sometimes associated with the aftermath of TDEs (Chakraborty et al., 2021, 2025a; Quintin et al., 2023; Miniutti et al., 2023; Nicholl et al., 2024), and TDE accretion disks are generally formed with arbitrary initial inclination, a significant disk inclination with respect to the MBH spin axis is not surprising. Moreover, the chaotic fluid dynamics of tidal stream circularization, stream-stream collisions, and accretion flow formation in Kerr spacetime further complicate the relationship between disk, stellar orbit, and MBH inclinations, further motivating the presence of misaligned disks. Such disks will be subject to rapid nodal precession via the frame dragging effect (Stone & Loeb, 2012; Franchini et al., 2016), which in turn modulates the orbiter-disk collision timings. The precessing disk is bound to align on a timescale dependent on black hole mass and effective viscosity, which can be as long as several years (Franchini et al., 2016). Previous modeling efforts have already laid the groundwork for modeling the effect of an inclined, rigidly precessing accretion disk on QPE timings (Franchini et al., 2023; Zhou et al., 2025b); we build upon it here with our efficient Bayesian inference code.

We adopt a coordinate system with orthonormal basis $(\hat{x},\,\hat{y},\,\hat{z})$ where the spin axis of the MBH, $\hat{a}_{\bullet}$, is aligned with the $\hat{z}$-axis ($\hat{z}=\hat{a}_{\bullet}$). The origin of the coordinate system is also chosen to coincide with the MBH. In this configuration, the normal vector of the accretion disk, $\hat{d}$, is inclined at an angle $\theta_{\rm disk}$ relative to $\hat{a}_{\bullet}$. Similarly, the normal vector to the secondary’s orbital plane, $\hat{E}$, is inclined relative to $\hat{a}_{\bullet}$ at an angle $i$. This coordinate system is particularly convenient because both the secondary’s orbit and the accretion disk undergo nodal precession about the $\hat{z}$-axis. As a result, only two components of $\hat{d}$ evolve with time.

The disk normal vector as a function of time is given by:

[TABLE]

where $T_{\rm disk}$ is the Lense-Thirring (nodal) precession period of the accretion disk, and $\phi_{{\rm disk},0}$ is an arbitrary initial phase.

The EMRI orbit is parameterized as $\vec{r}(t)=(x(t),y(t),z(t))$, which represents the full geodesic trajectory. This trajectory incorporates the complexities of apsidal precession, nodal precession, and orbital evolution. The orbit of the secondary in relation to the primary black hole and the accretion disk is shown in Fig. 1. The distance $D$ between the disk and the EMRI as a function of time is given by the dot product $\hat{d}\cdot\vec{r}(t)$:

[TABLE]

The zeros of $D(t)$ correspond to the QPE timings. $D(t)$ can take positive or negative values, indicating which side of the disk the EMRI is located on at a given time.

We note that a limitation of our timing model is the assumption of a constant offset between the orbiter-disk collision and the QPE peak timing. It is certainly possible that variations in the accretion disk properties, mutual inclination between orbiter and disk, and/or instantaneous orbital velocity may result in differing rise timescales between eruptions. Moreover, recent work has suggested QPE may be powered by stellar debris-disk collisions (Yao et al., 2025; Linial et al., 2025) which would further complicate the timings. However, given the considerable uncertainty regarding which physical effects dominate and to what extent, the lack of a complete theoretical understanding of how flare profiles relate to collision dynamics, and the generally constant rise times observed in any given QPE (within a factor of a few), we assume a constant offset for simplicity.

2.3 Likelihood for QPE timings

The likelihood function quantifies how well the predicted disk-crossing times match the observed QPE timings. These predicted crossings are determined by the zeros of the disk-crossing distance function $D(t,\vec{\xi})$, previously defined in Eq. 2.2, where $\vec{\xi}$ denotes the full set of model parameters. The disk precession parameters $(\theta_{\rm disk},T_{\rm disk},\phi_{{\rm disk},0})$ describe the time-dependent orientation of the accretion disk: $\theta_{\rm disk}$ is the inclination of the disk relative to the MBH spin axis, $T_{\rm disk}$ is the disk precession period, and $\phi_{{\rm disk},0}$ is the initial phase of the disk orientation.

The components of the secondary’s trajectory, $x(t)$, $y(t)$, and $z(t)$, are determined by bound Kerr geodesics and encapsulate relativistic effects such as apsidal and nodal precession. These trajectories depend on the orbital parameters $(\chi_{\bullet},p,e,i,q_{r,0},q_{z,0},q_{\phi,0})$, which uniquely characterize the geodesic. Given a realization of model parameters $\vec{\xi}$, the function $D(t,\vec{\xi})$ can be evaluated over time, and its zeros correspond to the predicted flare times $t_{\text{model},j}(\vec{\xi})$. For our code, we evaluate the zeros by identifying values of $t$ for which $D(t)$ changes sign, meaning the sampling resolution at which we evaluate $D(t)$ is important for inference accuracy. These computed disk crossings are then compared to the observed QPE timings $t_{\text{obs},i}$ to compute the likelihood. Figure 2 illustrates how the structure of $D(t,\vec{\xi})$ determines the timing of observable eruptions.

The observed QPE timings, compared to the orbital disk crossing times, are modified by the Shapiro and geometric (Roemer) delays. The Shapiro delay is given by:

[TABLE]

where $\hat{n}_{\rm cross}$ is the unit vector pointing from the black hole to the orbiter–disk collision location, as defined in Eq. 2.2, and $\hat{n}_{\rm obs}$ is the line-of-sight unit vector pointing from the black hole to the observer. We assume a fixed azimuthal viewing angle without loss of generality, such that the observer lies in the $x$–$z$ plane. The line of sight is then specified by a single polar angle $\theta_{\rm obs}$ measured from the $z$-axis:

[TABLE]

The geometric delay is then given by:

[TABLE]

Therefore, the complete set of parameters in our model is:

[TABLE]

where the parameters $(q_{r,0},q_{z,0},q_{\phi,0},\phi_{\rm{disk},0})$ are choices of initial phase which do not carry interesting physical information. For the purposes of the simulated injections in this study, we fix these initial phases to zero. It is possible that the initial phases act as nuisance parameters affecting the measurability of the other physically interesting parameters, particularly for observing times short compared to the precession timescales, but our analysis is not sensitive to this effect. In the next section, we convert the semilatus rectum $p$ to the semi-major axis333The semi-major axis $a$ should not to be confused with $\chi_{\bullet}$, which is the spin of the massive black hole. The semi-major axis $a$ is a commonly used orbital parameter in QPE literature, whereas the EMRI gravitational wave community typically describes the secondary’s orbit using the semilatus rectum $p$. $a=p/(1-e^{2})$.

For a given window of continuous observations, the $i$th observed QPE is then matched against the timing of the $i$th computed disk crossing:

[TABLE]

where $t_{\text{model},i}(\vec{\xi})$ are the zeros of $D(t,\vec{\xi})$, modified by the Shapiro and geometric delays, for a given set of parameters. To account for observational uncertainties, we also introduce a vector of timing errors $\{\sigma_{i}\}$ corresponding to the measurement uncertainty on each of the QPE timings. For our simulated data, we assume a constant error of $\sigma_{i}=100\;\rm sec$ for all QPEs (though with constant errors, any choice would be equivalent, as it would uniformly scale the likelihoods). We then divide the squared residuals of Eq. 23 by the squared error to compute the $\chi^{2}$ statistic across each of $N_{\rm obs}$ total observing windows:

[TABLE]

The likelihood function is then expressed as:

[TABLE]

In this framework, the predicted disk-crossing times $t_{\text{model},i}(\vec{\xi})$ explicitly depend on the model parameters, via the trajectory of the orbiter. By minimizing the residuals between the observed and predicted timings, the likelihood assigns higher probabilities to parameter combinations that better match the data. In performing inference with this likelihood function, we assume flat priors such that the posterior is proportional to the likelihood via Bayes’ theorem. We note that while geometric considerations suggest uniform priors in $\cos\theta$ for the angular parameters $\theta_{\rm obs}$ and $\theta_{\rm disk}$, we adopt flat priors in $\theta$ for simplicity. This choice has minimal impact for $\theta_{\rm obs}$ given its weak effect on timings (Appendix C), though the appropriate prior for $\theta_{\rm disk}$ remains uncertain and depends on the unknown formation/evolution mechanisms of QPE-hosting disks.

We performed Bayesian inference using the affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler (Goodman & Weare, 2010) as implemented in the emcee Python package (Foreman-Mackey et al., 2013). These calculations were carried out within our inference framework, QPE-FIT, which features trajectory and likelihood computations that are parallelized and GPU-accelerated using the CuPy Python library (Okuta et al., 2017). CuPy provides a Python interface to CUDA, enabling high-throughput numerical computations on GPUs. Computations were performed on NVIDIA A100 GPUs available on the engaging cluster operated by the MIT Office of Research Computing and Data. We present the main results of our MCMC inference in Section 3. In addition to those results, we performed a large grid of simulations to explore the parameter space extensively, presented in Appendix E.

3 Results

Using QPE-FIT (Section 2), we explored a large grid of synthetic QPE timings across the parameter space of $\{a,e,i,M_{\bullet},\chi_{\bullet},\theta_{\rm disk},T_{\rm disk}\}$, also varying the observing strategy, to explore the ability of MCMC sampling to recover the true parameters. For each sampling run we inject a solution computed with the exact Kerr geodesic trajectory (Park & Nasipak, 2024), then perform inference with the less accurate, but more efficient, PN sampler. We begin with only the static disk model; we then generalize to a misaligned precessing disk to show its confounding effect on parameter inference; and finally apply our model to the particularly irregular QPE timing data of eRO-QPE1 to showcase its ability to describe observations of even the erratic sources, albeit with extremely multimodal posteriors.

3.1 The effects of varying semimajor axis, eccentricity, spin, and $M_{\bullet}$ with a non-precessing disk

We first consider the effect of varying EMRI/MBH parameters in the simplest-case model, i.e. a static planar accretion disk aligned with the MBH spin axis. All results from the MCMC sampling runs reported in Sections 3.1-3.2 were carried out with 1000 independent walkers for 100,000 steps, and used an uninterrupted 100 ks observing window. We use a large walkers-to-steps ratio because of the multimodality of the parameter space, with starting positions chosen from a uniform random distribution outlined in Table 2. We discarded the first 50% of steps as burn-in, then retained the 100,000 highest-likelihood points from the resulting distributions. We discuss the convergence of these MCMC chains in Appendix B. Given the broad, conservative priors, compatible with actual uncertainties on the underlying system properties given real data, it is remarkable that several of the parameters (e.g. $a$, $e$, $M_{\bullet}$) converge to the true values with posterior uncertainties of $\sim 10\%$.

In our first example, we inject a signal with semimajor axis $a=100R_{g}$ (where $R_{g}\equiv GM_{\bullet}/c^{2}$), detector-frame MBH mass $M_{\bullet}=10^{5}M_{\odot}$, and MBH spin $\chi_{\bullet}=0.5$. We perform inference for four different choices of eccentricity, $e\in\{0.01,0.1,0.3,0.5\}$. Our results are shown in the top panel of Fig. 3. We observe that very low-eccentricity EMRIs ($e\sim 0.01$) will not be as powerful for $M_{\bullet}$ mass constraints, as the timing modulation effects of apsidal precession are not as pronounced as in EMRIs with modest eccentricity ($e\sim 0.1-0.3$). On the other hand, at higher eccentricities ($e\sim 0.5$) the parameter estimation once again becomes worse, with wider errors on $a$, $e$, and $\log M_{\bullet}$. However, this is likely due to the breakdown of the PN eccentricity expansion used in our orbital trajectory calculation, rather than an intrinsically larger degeneracy. The MBH spin ($\chi_{\bullet}$) is unconstrained in all cases we explored, which are limited to relatively wide orbits ($\geq 100R_{g}$) and short observing baselines.

For the middle panel of Fig. 3, we inject a signal with $e=0.1$, $M_{\bullet}=10^{5}M_{\odot}$, $\chi_{\bullet}=0.5$, and three different choices of $a\in\{100,200,400\}R_{g}$. Remarkably, varying $a$ has relatively little effect on measurements of $M_{\bullet}$. While larger $a$ results in a slower apsidal precession (thus making $M_{\bullet}$ harder to measure), it also results in larger-amplitude timing modulations; we posit that the two effects offset each other to result in a similar overall mass constraint.

For the bottom panel of Fig. 3, we inject a signal with $a=100R_{g}$, $e=0.1$, $\chi_{\bullet}=0.5$, and two different choices of $\log(M_{\bullet}/M_{\odot})\in\{5,6\}$. The effect of $M_{\bullet}$ enters predominantly through the QPE recurrence timescale: for a fixed $a/R_{g}$, the orbital timescale increases linearly with $M_{\bullet}$. Increasing $M_{\bullet}$ by $100\times$ thus requires a (usually infeasible) corresponding $100\times$ longer observing baseline, making lower-mass MBHs better-suited for QPE discovery and parameter estimation.

In Fig. 4, we show a corner plot from a sampling run with injected parameters of $a=100R_{g}$, $e=0.1$, $M_{\bullet}=10^{5}M_{\odot}$, $\chi_{\bullet}=0.5$, $i=60^{\circ}$, and $\theta_{\rm obs}=\pi/4$. The orbital semimajor axis is significantly degenerate with $M_{\bullet}$ (as expected, see e.g. Fig. 4 of Kejriwal et al. 2025), and the main effect disentangling them is the timescale of Schwarzschild (apsidal) precession. Eccentricity is essentially well-constrained within just a few bursts, whereas the EMRI orbital inclination ($i$) and MBH spin ($\chi_{\bullet}$) are completely unconstrained. The former is due to the relatively large radii at which QPEs occur ($a\gtrsim 100R_{g}$): the Kerr metric is not spherically symmetric, so spacetime does vary with $i$, but only weakly except for strong-field orbits.

The dominant observable effect of the black hole spin $\chi_{\bullet}$ is Lense-Thirring (nodal) precession, but only over timescales significantly longer than the ones probed by our experiments. Long-term dissipative effects are also likely to set in over these baselines, but our analysis does not incorporate these; in follow-up work, we will examine parameter estimation over much longer baselines while accounting for secular evolution due to interaction with the environment. Other QPE timing studies, e.g. Xian et al. (2021) and Zhou et al. (2025a), have explored spin constraints by considering these longer baselines, and in some cases produce good constraints with uncertainties $<0.1$.

We also varied several of the other parameters, such as the observing angle $\theta_{\rm obs}$ and MBH spin ($\chi_{\bullet}$), and found that they had limited impact on the parameter recovery. See Appendix C for further discussion.

3.2 The effects of an inclined, rigidly-precessing disk

As discussed in Section 2.2 and in prior studies (Franchini et al., 2023; Linial & Metzger, 2024; Zhou et al., 2025b), the static in-plane disk model of Section 3.1 appears inconsistent with observations in some QPEs. Systems such as eRO-QPE1 (Chakraborty et al., 2024) and GSN 069 (Miniutti et al., 2025) show long-term modulations in QPE arrival times not predicted by static disk models, pointing instead to a further modulating timescale. One possible (but not unique) explanation is rapid disk precession, naturally expected for an isotropic distribution of TDEs with respect to the black hole spin axis. Given that the weak-field per-orbit apisdal motion is approximately $\delta\varphi\approx 6\pi(a/R_{g})^{-1}(1-e^{2})^{-1}$, the EMRI apsidal precession timescale ($\Delta\varphi=2\pi$) is given by:

[TABLE]

This expression occupies values between $\sim 25-133P_{\rm orb}$ for the range of eccentricities and semimajor axes explored by our trials. The predominant effect qualitatively impacting QPE timing behavior is how $T_{\rm aps}$ compares with $T_{\rm disk}$.

Before proceeding to injection/recovery trials, we first provide some illustrative examples to showcase these effects, for cases where $T_{\rm disk}\gg T_{\rm aps}$, $T_{\rm disk}\approx T_{\rm aps}$, and $T_{\rm disk}\ll T_{\rm aps}$. In Fig. 5 we show synthetic $O-C$ (observed minus calculated) diagrams, where the “calculated” timings assume a constant recurrence time and the “observed” timings are modified by the combined orbital and disk precessions. We construct $O-C$ diagrams considering all bursts together, as well as odd and even bursts separately, as this effectively isolates the orbital evolution of the EMRI subject to precession effects.

The two parameters affecting the deviation from the static, aligned-disk case are the disk inclination angle with respect to the MBH spin axis ($\theta_{\rm disk}$) and the disk precession period ($T_{\rm disk}$). In the top panel, we show the $O-C$ diagram expected from the static disk model (which is equivalent to taking $T_{\rm disk}\rightarrow\infty$), in which the odd and even QPE flares trace sinusoids in antiphase over the apsidal precession period (red/blue points) while also producing the characteristic “long-short” recurrence time alternation on shorter timescales (black points). For a slowly precessing disk, with $T_{\rm disk}=33P_{\rm orb}\approx T_{\rm aps}$, the apsidal precession is still measurable, and the characteristic long-short alternation is preserved (black points), but the odd and even bursts now vary together in phase (red/blue points). For sufficiently rapid disk precession, the long-short alternation is eliminated completely, with the disk precession becoming the dominant timing effect—in the bottom panel we show the case of $T_{\rm disk}=5P_{\rm orb}\ll T_{\rm aps}$, in which a sinusoidal trend on the disk precession timescale becomes the dominant effect. The middle panel is qualitatively similar to the $O-C$ diagram computed for GSN 069 in Miniutti et al. (2025), in the sense that there is a long/short alternation on timescales $\approx P_{\rm orb}$ and an in-phase modulation on timescales $\gg P_{\rm orb}$, whereas the bottom panel is more akin to eRO-QPE1 (Chakraborty et al., 2024), suggesting a range of disk dynamical properties across the QPE population. No known QPEs have thus far been shown to be compatible with the static-disk picture.

Following this toy example, we now proceed to injection/recovery runs including a misaligned rigidly precessing accretion disk. In the top panel of Fig. 6 we show, by analogy to Fig. 3, the eccentricity constraints with a precessing disk for a 100 ks observing window (about 50 consecutive QPEs). We do not duplicate the posterior distribution of $\chi_{\bullet}$ because it is equally unconstrained compared to the non-precessing case. In the top panel of Fig. 6 we inject solutions with $a=100R_{g}$, $M_{\bullet}=10^{5}M_{\odot}$, $\chi_{\bullet}=0.5$, $\theta_{\rm disk}=20^{\circ}$, $T_{\rm disk}=20P_{\rm orb}$, and four different eccentricities $e\sim\{0.01,0.1,0.3,0.5\}$. In the middle panel, we hold fixed $e=0.1$, use two different disk inclinations $\theta_{\rm disk}\in\{5^{\circ},20^{\circ}\}$, and leave all other parameters the same. In the bottom panel, we fix $\theta_{\rm disk}=20^{\circ}$, use three different disk precession periods $T_{\rm disk}\in\{5,20,50\}P_{\rm orb}$, and leave other parameters unchanged. In the middle/bottom panels we vary the disk inclination angle $\theta_{\rm disk}$ and precession period $T_{\rm disk}$, respectively.

Relative to the static-disk scenario, the parameter constraints degrade significantly. This is an expected consequence of introducing additional degrees of freedom without augmenting the data with new constraining features such as light curve morphology or spectra. The observed sequence of QPE arrival times cannot in general break the degeneracies between multiple parameter combinations that produce similar timing modulations, at least for the baselines we examined. This inherent limitation underscores a key challenge in using QPE timing data alone to constrain EMRI models with realistic disk physics. Using additional observables or longer baselines may help to lift these degeneracies and better constrain the underlying physical parameters.

In Fig. 7 we show a corner plot from a sampling run with $a=100R_{g}$, $e=0.1$, $M_{\bullet}=10^{5}M_{\odot}$, $\chi_{\bullet}=0.5$, $\theta_{\rm obs}=\pi/4$, $\theta_{\rm disk}=20^{\circ}$, and $T_{\rm disk}=20P_{\rm orb}$. There is a significant degeneracy between $i$ and $\theta_{\rm disk}$, and a slight correlation between $e$ and $\theta_{\rm disk}$. The increased parameter dimensionality significantly worsens constraints on the inference of most parameters, with only $e$, $\theta_{\rm disk}$, and $T_{\rm disk}$ appearing to converge. Given that disk precession (or some other modulating timescale) appears to be common, we advise longer observational baselines for parameter inference.

3.3 The effect of neglecting disk precession when performing parameter inference

We further consider the effect of fitting observations assuming a static disk, when the true underlying process is more complex. Omitting precession leads to systematic biases across all parameters, driven by the model’s failure to capture key time-dependent structure. To evaluate the impact of this effect on parameter estimation, we inject a solution with a precessing disk, then perform inference with 500 walkers$\times$5000 steps assuming the static-disk model. Figure 8 shows the resulting posterior distributions, while Table 3 summarizes maximum a posteriori (MAP) values, uncertainties, and deviations from the true parameters

The largest bias appears in the observing angle, with a MAP shift to $\theta_{\rm obs}=0.38$, nearly $4\sigma$ from the truth, reduced to $0.92\sigma$ when precession is included. Biases in $\log M_{\bullet}$, semi-major axis, and eccentricity are similarly reduced, with deviations falling below $0.2\sigma$ under the precessing model. Posterior widths also decrease, except for $\theta_{\rm obs}$, which remains weakly constrained due to its weak effect on the observed burst timings. These results underscore the importance of modeling disk precession to accurately recover system parameters in the presence of longer-period modulations in the signal. While acceptable fits can be sometimes be obtained with static-disk models over short baselines, longer baselines in particular require precession to be correctly modeled in order to avoid biased inferences.

3.4 The effects of varying observational cadence and baseline

We now consider how varying the observation strategy modifies parameter constraints. For a given set of parameters, observing more QPEs with higher density and longer baseline generally improves constraints, at the cost of more resources. To quantify the information gain of additional observing time, we performed sampling runs for three different observing strategies: $1\times 100$ kilosecond observation (chosen to approximately match the capability of XMM-Newton); $4\times 100$ kilosecond observations, with short gaps between snapshots of 200 kiloseconds (roughly mimicking the observing strategy of four contiguous XMM-Newton orbits); and $4\times 100$ kilosecond observations with long gaps of 1 megasecond between observations.

In Fig. 9, we show the effect of varying the invested observation time, first for an in-plane static disk (top panel). For an assumed system with $a=100R_{g}$, $e=0.1$, $M_{\bullet}=10^{6}M_{\odot}$ (note the higher-than-typical value chosen here), $\chi_{\bullet}=0.5$, we see that $1\times 100$ kilosecond observation—which observes 5 consecutive eruptions—provides weak constraints on the EMRI semimajor axis, eccentricity, and $M_{\bullet}$, due to the small number of bursts. Increasing to $4\times 100$ ks observations results in improved constraints, but for the same invested time, we find that a shorter baseline with smaller observing gaps is better able to constrain $a$, $e$, and $M_{\bullet}$. This is a surprising result, as this strategy covers a smaller fraction of the overall apsidal period. We suggest this may be due to the difficulty of uniquely attributing a burst number following a longer gap. This is even more pronounced for the case including disk precession (bottom panel of Fig. 9), which further adds the ambiguity of uniquely modeling the apsidal motion.

3.5 Application to data of eRO-QPE1

To showcase the ability of our model to describe QPE timing data over short baselines, particularly for the more erratic sources which are difficult/impossible to analyze within the standard static-disk model, we fit the burst timings of eRO-QPE1, a particularly irregular source in both its timing properties and eruption profiles Arcodia et al. (2022). Using the NICER data from Chakraborty et al. (2024), we fit each epoch separately due to the long (month- to year-long) gaps between observations and the non-inclusion of secular evolution effects in our timing code. We used uniform priors (Table 2) with $a\sim[100,500]R_{g}$, $e\sim[0,0.5]$, $i\sim[0^{\circ},90^{\circ}]$, $\theta_{\rm disk}\sim[0,60^{\circ}]$, and $T_{\rm disk}\sim[2,100]P_{\rm orb}$. We also added two additional parameters: $t_{0}\sim[-50000,50000]\;\rm sec$, which allows for an arbitrary initial time offset of the EMRI-disk system to match the observed timings, necessitated because the three initial phases of the Kerr geodesics were set to zero; and $\phi_{\rm disk,0}\sim[0,2\pi]$, which allows for an arbitrary initial phase offset between the accretion disk and EMRI orbit. We assumed a fixed $M_{\bullet}=10^{5.8}M_{\odot}$ (Wevers et al., 2022), as in this exercise we are not interested in independently inferring the MBH mass, and making some choice is necessary to break the $a-M_{\bullet}$ degeneracy. We also fix $\chi_{\bullet}=0.5$ and $\theta_{\rm obs}=\pi/2$, because they have little to no effect on the timings for these large orbital radii and short baselines. In total, the free parameters of this fit are {$a$, $e$, $i$, $\theta_{\rm disk}$, $T_{\rm disk}$, $\phi_{\rm disk,0}$, $t_{0}$}.

We find the posteriors are extremely multimodal, so we do not report the recovered parameters for further astrophysical inference. In Fig. 10 we show the predicted orbiter-disk collision timings compared to the real data, showcasing that as a proof-of-concept, rigid disk precession appears able to roughly capture the observed erratic burst timings for reasonable parameter choices. As discussed in Section 3.2, sources such as these make inferring $M_{\bullet}$ and $a$ difficult, and we find even $e$, $i$, $\theta_{\rm disk}$ and $P_{\rm disk}$ are not uniquely constrained here.

The fact that no single solution can fit all epochs together is almost certainly indicative of secular evolution in the EMRI-disk system, e.g. due to orbital decay, disk alignment, and changes in the disk density profile or precession speed. As we have not accounted for these effects in our short-timescale model, it cannot fully fit all epochs coherently in its current form. However, as QPE monitoring baselines grow longer, including secular effects in timing models is well-motivated, and will be the subject of future work (see Zhou et al. 2025b; Xian et al. 2025 for further discussion).

4 Discussion

QPEs provide a novel window into orbital dynamics in galactic nuclei. If powered by repeated interactions between the accretion disk of an MBH and an orbiting secondary, their timing encodes key orbital and system parameters, including the spin and mass of the central black hole, and the eccentricity and semi-major axis of the secondary. In Section 3 we have explored the constraining power of QPE timings alone in parameter inference of the MBH and EMRI properties. Two noteworthy shortcomings of our methods are:

MCMC-based inference, which is notoriously susceptible to stalling at local minima for multimodal posterior surfaces (e.g. Fig. 7). As a result, it is possible the errors we report are over-estimates of the best-case constraints that may be obtained by more sophisticated sampling methods such as nested sampling (see e.g. Zhou et al. 2025b for comparison). 2. 2.

Short observational baselines consisting of $\mathcal{O}(10-100\rm s)$ of QPEs, over which longer-term effects such as orbital decay and EMRI nodal precession driven by frame-dragging around a spinning MBH, do not set in. As a result, we do not produce any constraints on $\chi_{\bullet}$ and do not address prospects for measuring $\dot{P}$, which have also been addressed in other studies (Xian et al., 2021; Zhou et al., 2025b).

Nevertheless, our work shows that, at least for some EMRI orbital parameters and disk properties, QPE timing may be used to infer $a$, $e$, $M_{\bullet}$, $\theta_{\rm disk}$, and $T_{\rm disk}$, to varying degrees of precision. Once inferred, these parameters become powerful tools for probing a range of astrophysical questions: How do massive black holes grow? What seeds them? Which channels form EMRIs? And how do QPEs couple to their accretion environments? In this section, we discuss what astrophysics becomes accessible through inference on each parameter.

4.1 What can we learn from $M_{\bullet}$?

Constraining the mass of the central black hole in QPE-hosting systems offers a unique opportunity to probe the less well-understood low-mass end of the black hole mass function (BHMF), defined as

[TABLE]

which quantifies the comoving number density of black holes per logarithmic mass interval at a given redshift. While the high-mass end of the BHMF is relatively well constrained through stellar and gas dynamical measurements in massive galaxies (McConnell & Ma, 2013) and reverberation mapping of luminous AGN (Peterson et al., 2004), the low-mass end ($M_{\bullet}\lesssim 10^{6}\,M_{\odot}$) remains uncertain due to observational selection effects and the difficulty of detecting and characterizing faint or inactive nuclei in low-mass galaxies (Greene & Ho, 2007).

In this regime, even the black hole occupation fraction is not well known, and different MBH seeding models (e.g., light vs. heavy seeds) make divergent predictions. Light seed models predict a steep low-mass slope and high occupation fraction in dwarf galaxies, while heavy seed models yield a suppressed BHMF below $10^{5}-10^{6}M_{\odot}$ and a low occupation fraction in low-mass hosts (Volonteri, 2010; Greene et al., 2020; Ricarte & Natarajan, 2018). QPEs provide a new dynamical channel to infer MBH masses in otherwise quiescent or low-luminosity systems, potentially enabling direct constraints on the local BHMF and the demographics of the smallest massive black holes. Thus far, TDEs have been the dominant observational probe of low-mass MBHs; QPE timing may provide a promising method for refining the mass measurements in this population (see also Zhou et al. 2025a), but as shown in Section 3, large numbers of bursts may be required for reliable constraints.

4.2 What can we learn from $e$ and $a$?

The distribution of orbital semimajor axes ($a$) and eccentricities ($e$) of QPE-hosting EMRIs offer constraints on their formation channels, which is relevant for estimating the EMRI rate observable by LISA. EMRIs formed via the dry loss-cone channel or the Hills mechanism are expected to be highly eccentric ($e\gtrsim 0.9$), as they are injected onto plunging orbits through dynamical scatterings or binary disruptions in the nuclear stellar cluster (Hopman & Alexander, 2005; Preto & Amaro-Seoane, 2010; Broggi et al., 2022; Miller et al., 2005; Raveh & Perets, 2021). In contrast, EMRIs formed in the wet channel, via migration through an accretion disc, experience efficient eccentricity damping due to interaction with the gaseous environment, yielding nearly circular orbits with $e\sim h\ll 1$, where $h$ is the local disc aspect ratio (Lyu et al., 2024; Pan et al., 2021b; Pan & Yang, 2021). EMRIs formed via the Hills mechanism could retain high eccentricities even after extended GW-driven evolution (see e.g., Fig. 4 of Raveh & Perets (2021)). Therefore, measurements of eccentricity from QPE-hosting systems can help provide evidence in favor of a wet or dry formation scenario.

The semi-major axis distribution is sensitively dependent on the orbital decay mechanisms. The GW-driven inspiral timescale scales steeply with the semi-major axis, approximately as $t_{\rm GW}\propto a^{4}$ for high-eccentricity systems (Peters, 1964), allowing predictions for when the secondary will either merge with the central MBH or be tidally disrupted. Moreover, since the companion crosses the accretion disc at a radius $R_{\rm cross}\sim a(1-e)$, the local gas density $\rho(R)\propto R^{-s}$ and relative velocity $v_{\rm rel}\sim\sqrt{GM_{1}/R_{\rm cross}}$ directly influence the energy deposited during each crossing, and therefore the amount of orbital energy viscously dissipated to accelerate the inspiral.

4.3 What can we learn from $\theta_{\rm disk}$ and $T_{\rm disk}$?

The properties of accretion disks hosted in QPE systems—and more broadly, those around massive black holes—are a rich area of active study, as there remain longstanding uncertainties in their formation, stability, and structural properties (Kara & Garc´ıa, 2025). Most theoretical treatments of QPE accretion disks thus far assume either a standard steady-state $\alpha$-disk model (Linial & Metzger, 2023; Yao et al., 2025; Vurm et al., 2025; Tsz-Lok Lam et al., 2025) or do not explicitly model the disk (Xian et al., 2021; Franchini et al., 2023). Given the emergence of QPEs in the accretion flows formed by TDEs (Chakraborty et al., 2021, 2025a; Miniutti et al., 2023; Quintin et al., 2023; Nicholl et al., 2024), it is very likely the disks in these systems diverge significantly from this steady-state assumption. Recent work has explicitly taken steps toward treating the disk properties with the finite mass/energy/angular momentum reservoirs of compact TDE disks (Guolo et al., 2025a; Mummery, 2025) by fitting the disk SEDs to infer the surface density profiles and radial extents. They indeed appear compatible with compact TDE disks extending to the EMRI orbital radii (Mummery & Balbus, 2020; Mummery et al., 2024; Nicholl et al., 2024; Chakraborty et al., 2025a; Wevers et al., 2025; Guolo et al., 2025b), adding mild circumstantial evidence to the EMRI-disk collision paradigm, albeit with some standing questions. Preliminary time-resolved studies of QPE X-ray spectra are also beginning to explore the accretion disk environment (Kosec et al., 2025; Chakraborty et al., 2025b), probing the properties of possible ionized outflows ejected from the disk in the QPE-generating process.

QPE timings probe orthogonal properties to the SED/spectral models, via the inclination angle and precession period. It has long been suggested that Lense-Thirring precession should be an observable phenomenon in TDE disks—e.g. via quasi-periodic oscillations of the disk emission—provided the accretion rate is high enough that the precession period exceeds the sound-crossing time in the disk. This is equivalent to requiring thick-disk conditions, i.e. $H/R\gtrsim\alpha$, where $\alpha$ is the viscosity parameter, $H$ the scale height, and $R$ the precession radius (Stone & Loeb, 2012; Franchini et al., 2016). As the disk precession frequency depends jointly on the MBH spin, mass, and disk density profile, detection of precession and progressive alignment could also constrain the MBH properties. QPEs thus offer the possibility of studying the newly-formed, rapidly-evolving compact disks around MBHs via their dynamical imprints on the collision timings, adding a new method alongside the standard spectral fitting approach.

4.4 Comparisons with prior work

During the preparation of this work, two studies were posted (Zhou et al., 2025b, a) which explore a similar Bayesian inference paradigm for QPE timings with a rigidly precessing accretion disk. Our work was carried out independently, using different methods including the post-Newtonian formalism rather than exact Kerr geodesic trajectories, MCMC rather than nested sampling, and GPU acceleration, with the primary goal of computational efficiency. In particular, the costs of nested sampling and exact calculation of the Kerr geodesics are prohibitively expensive in many cases, requiring up to 640 CPU cores over 20 days ($\gtrsim 300,000$ CPU-hours) for a full inference run including the effects of disk precession and orbital decay (Zhou et al., 2025b). While we have not yet introduced orbital evolution into our code—and also therefore limited our focus to shorter observational baselines than considered in those works—our work provides a first step toward developing more inexpensive inference tools.

The speed of our inference tools allowed us, for the first time, to systematically explore the entire EMRI/disk parameter space broadly. The QPE timing problem warrants independent investigations with a range of different approaches, particularly as the observational data grows and its relevance to the broader astrophysical community becomes more compelling.

5 Conclusion

We have developed QPE-FIT, an efficient model to compute the timings of orbiter-disk collisions for a low-mass companion around a massive black hole (MBH) hosting an accretion disk. The model is relevant for the recently-identified Quasi-Periodic Eruptions (QPEs), which have been proposed to be powered by star- or BH-disk collisions in an extreme mass-ratio inspiral (EMRI). In principle, the burst timings encode information about the EMRI orbital trajectories, including the general relativistic precessions affecting both the orbiter and the accretion disk; thus, precise measurement of QPE timings may allow inference of the properties of the MBH (mass, spin), EMRI orbit (semi-major axis, eccentricity), and accretion disk (inclination angle, nodal precession timescale).

QPE-FIT computes the geodesic trajectory of the companion in the post-Newtonian formalism up to 3PN order, following the eccentricity expansion of Sago & Fujita (2015). These analytic approximations, along with the significant acceleration afforded by GPUs (Table 1), enable rapid computation of the trajectories, making them computationally amenable to Bayesian parameter inference. By carrying out an extensive grid of MCMC parameter inference simulations across a wide grid of injected EMRI/disk/MBH properties, we have quantified the characteristic errors on each inferred parameter. The results are presented and described in Section 3, Figures 3-9, and Table LABEL:tab:errors. We find the following:

• •

For the simple-case model of an EMRI colliding with a static, non-precessing accretion disk, QPE timings provide $\lesssim 10$% error-bars on EMRI semimajor axis/eccentricity and MBH mass within tens of orbital periods (Fig. 4). The best constraints are for EMRIs with modest eccentricities $0.1\lesssim e\lesssim 0.5$ (Fig. 3).

• •

However, given that observed QPE timings appear more complicated than this simple picture, it is important to consider further effects influencing the timing, e.g rigid-body precession of the accretion disk (Fig. 5). Introducing a precessing disk complicates parameter inference considerably, most significantly degrading MBH mass constraints (Fig. 6), although the disk properties themselves (e.g. inclination, precession rate) can be constrained to the $\sim 10-50\%$ level (Fig. 7, Table LABEL:tab:errors). Performing inference with the static-disk model when the true underlying process involves a precessing disk can systematically induce biases in the inferred parameters (Fig. 8).

• •

Better constraints typically arise from monitoring campaigns with shorter observational gaps (Fig. 9), motivating high-cadence observations for constraints on MBH mass/EMRI properties.

• •

To showcase the ability of QPE-FIT to describe real QPE timing data, we applied it to the irregular timings of eRO-QPE1 (Fig. 10), finding reasonable fits but with time-varying parameters. This motivates self-consistently modeling the effects of secular period/eccentricity/disk evolution alongside the EMRI trajectory calculation.

In this paper, we only explored relatively short observational baselines of $\mathcal{O}(10-100\rm s)$ of QPEs, meaning we have not considered timing effects from the secular evolution of the EMRI-disk system. These baselines also prohibit inference of MBH spin, as the effect of general relativistic frame-dragging only sets in on timescales $\gg P_{\rm orb}$. Both of these effects are inevitable consequences of the EMRI model, which predicts orbital decay driven by a combination of gravitational waves, gas drag, and/or dynamical friction (Kejriwal et al., 2024; Linial & Quataert, 2024; Zhou et al., 2025b; Xian et al., 2025). Long-term monitoring observations of some sources have even tentatively detected orbital period evolution in some sources (Arcodia et al., 2024b; Zhou et al., 2025b; Miniutti et al., 2025; Xian et al., 2025), though not yet with unambiguous solutions. Continued monitoring of the known QPEs is sure to place further constraints on models, and in future work, we will add secular evolution to the model presented here. Moreover, the posterior error values on some parameters are likely to be overestimated due to the inherent shortcomings of Markov Chain Monte Carlo (MCMC) inference, particularly for the multi-modal parameter spaces relevant here. Exploring more sophisticated inference schemes is thus well-motivated, and may provide significant improvements in parameter inference from QPE timings.

Acknowledgements

We thank Katerina Chatziioannou and Ethan Payne for useful discussions. We also thank the anonymous referee for their helpful comments. The authors acknowledge the MIT Office of Research Computing and Data for providing high performance computing resources that have contributed to the research results reported within this paper. Additionally, some of the calculations in this work were conducted in the Resnick High Performance Computing Center, a facility supported by Resnick Sustainability Institute at the California Institute of Technology. LVD is supported by the Sherman Fairchild Postdoctoral Fellowship at the California Institute of Technology. SAH is supported by NSF Grants PHY-2110384 and PHY-2409644; NSF Grant PHY-2110384 also supported LVD at MIT during a portion of this work. GM thanks the Spanish MICIU/AEI/10.13039/501100011033 and ERDF/EU grants n. PID2020-115325GB-C31 and n. PID2023-147338NB-C21 for support. MG is funded by Spanish MICIU/AEI/10.13039/501100011033 and ERDF/EU grant PID2023-147338NB-C21. MB acknowledges support from the Italian Ministry for Universities and Research (MUR) program “Dipartimenti di Eccellenza 2023-2027”, within the framework of the activities of the Centro Bicocca di Cosmologia Quantitativa (BiCoQ). This work makes use of the Black Hole Perturbation Toolkit (Black Hole Perturbation Toolkit Collaboration, 2025), in particular the KerrGeoPy package (Park & Nasipak, 2024).

Appendix A Impact of the Post-Newtonian approximation in the trajectory

Figure A.1 quantifies the deviation of 3PN orbital trajectories from exact Kerr geodesics for a range of orbital configurations. The deviation increases noticeably above eccentricities of 0.5, demonstrating the reduced accuracy of the PN model in describing highly eccentric motion. At fixed $e=0.3$, we also observe that deviation generally decreases with increasing $a$, consistent with improved PN convergence in the weak-field regime. The PN model performs best for low-eccentricity, large-radius orbits.

To examine the range over which the Post-Newtonian (PN) timing model remains accurate for parameter inference, we inject synthetic timing data using exact KerrGeoPy trajectories and then recover the eccentricity ($e$) and MBH spin ($\chi_{\bullet}$) with a 2D PN-based likelihood. All other EMRI parameters (mass, inclination, etc.) are held fixed at their true values, such that any deviations in the posterior directly reflect PN truncation errors in the recovered parameters. Figure A.2 presents the resulting posterior distributions. The two panels on the left show recovered $e$ and $\chi_{\bullet}$ for three injected eccentricities ($e=0.3,\,0.5,\,0.7$ at fixed $a=100\,R_{g}$), and the panels on the right show the same recovered quantities for three injected semi-major axes ($a=10,\,50,\,100\,R_{g}$ at fixed $e=0.3$).

In the top-left panel, both the width of the posterior distributions and the bias in the recovered eccentricities increase with the true value of $e$. At $e=0.5$, which marks the formal limit of post-Newtonian (PN) applicability, the posterior is noticeably broader and shifted, indicating a clear bias in the recovered eccentricity. At $e=0.7$, the PN expansion fails, resulting in a breakdown of the recovery. In the bottom-left panel, which shows the recovered semi-major axis ($a$) as a function of $e$, accurate recovery of $e$ (e.g., at $e=0.3$) leads to well-constrained posteriors on $\chi_{\bullet}$. However, for $e=0.5$ and $e=0.7$, where the eccentricity estimates are biased or highly uncertain, the posterior on $a$ becomes extremely broad and effectively unconstrained.

In the bottom right panel where the $\chi_{\bullet}$ value is recovered, spin inference shows a clear dependence on the injection radius. At large radii ($a=100\,R_{g}$), the timing signal becomes less sensitive to the effects of $\chi_{\bullet}$. As a result, the posterior widens and its median shifts modestly. At small radii ($a=10\,R_{g}$), deep in the strong-field regime, spin effects such as frame dragging are enhanced, but PN truncation errors dominate, leading to a broad and significantly biased posterior. In contrast, at $a=50\,R_{g}$, the PN model remains valid and the timing remains sufficiently sensitive to spin, yielding the most accurate and tightly constrained $\chi_{\bullet}$ posterior. This trend reflects competing influences: diminished spin sensitivity in the weak field versus poor PN convergence in the strong field, with $a\sim 50\,R_{g}$ emerging as an optimal regime for PN-based spin recovery. In the recovered-$e$ panels for varying $a$, eccentricity is accurately recovered at $a=100\,R_{g}$, with modest broadening at $a=50\,R_{g}$ as the PN approximation begins to deteriorate. At $a=10\,R_{g}$, the posterior becomes broad and biased, consistent with failure of the PN model. Unlike spin, the influence of eccentricity is not strongly relativistic and can still be inferred reasonably well in the weak-field, provided the PN series remains convergent.

These 2D posteriors define practical boundaries for PN-based inference; the PN timing model reliably recovers $(\chi_{\bullet},e)$ for $e\lesssim 0.4$ at $a\gtrsim 50\,R_{g}$, and should be replaced by a fully relativistic trajectory (e.g. KerrGeoPy (Park & Nasipak, 2024)) outside that regime. Note that because we performed a restricted two-parameter inference in this section, the posterior widths are inherently narrower, making biases appear more pronounced than they would in a full-parameter model.

Appendix B MCMC convergence

In Fig. B.1 we show two illustrative metrics addressing the convergence of the MCMC chains, for two sampling runs presented in Figs. 4 & 7 ($a=100R_{g}$, $e=0.1$, $M_{\bullet}=10^{5}M_{\odot}$, $\chi_{\bullet}=0.5$, and for the case with disk precession, $\theta_{\rm disk}=20^{\circ}$, $T_{\rm disk}=20P_{\rm orb}$). In the top panel, we show the integrated autocorrelation time (Goodman & Weare, 2010), which is a crude proxy for assessing whether the chains are properly converged. While there is no formally defined criteria for acceptable convergence, typically an autocorrelation time ($\tau$) smaller than $\tau\lesssim N/50$ (where $N$ is the number of MCMC steps) is taken as indication that the chain is reasonably converged (e.g. see the discussion in https://emcee.readthedocs.io/en/stable/tutorials/autocorr/). Our samplers reach this threshold roughly at 20,000 steps, and are well below it within $10^{5}$ steps, providing some indication that the results presented in Section 3 are reasonable.

The bottom panel of Fig. B.1 shows an even clearer indication that running the chains longer will not produce significant additional returns. We show the mean and standard deviation of the posterior likelihood distributions pictured in Figs. 4 & 7, which plateau within about $10^{4}$ steps and increase only modestly thereafter. We compare this to the “theoretical” maximum expected likelihood, which is obtained by evaluating the likelihood function for the true injected parameters. The posterior likelihood distribution falls noticeably short of this maximum, indicating the inherent shortcoming of MCMC for exploring this likelihood space within a reasonable number of steps.

Appendix C Impact of observing angle

We examine the impact of varying the observing angle $\theta_{\rm obs}$ on parameter recovery. Table C.1 summarizes the recovered MAP values and deviations from the true values for the semi-major axis, eccentricity, and black hole mass across $\theta_{\rm obs}=\pi/2$, $3\pi/4$, and $\pi/4$, while Figure C.1 shows the corresponding posterior distributions. We find that the inferred parameters are largely insensitive to the choice of $\theta_{\rm obs}$: deviations from the truth remain below $0.3\sigma$ in all cases. The semi-major axis exhibits the largest spread, with the MAP value decreasing from $103.00$ at $\theta_{\rm obs}=\pi/4$ to $95.38$ at $3\pi/4$, indicating a weak trend with observing angle. The greatest deviation in eccentricity occurs at $\theta_{\rm obs}=\pi/2$ (0.26$\sigma$), while the largest deviation in black hole mass appears at $3\pi/4$ (0.26$\sigma$).

These trends are modest, and overall changes in $\theta_{\rm obs}$ introduce only minor variations in parameter recovery. We conclude that uncertainties related to the observing angle are subdominant compared to other sources of modeling bias, such as whether disk precession is properly included.

Appendix D Corner plots of data fits to eRO-QPE1

In Fig. D.1, we show corner plots of each of the individual-epoch fits to the NICER data of eRO-QPE1, as outlined in Section . The posteriors are smoothed by a $2\sigma$ gaussian kernel. The blue horizontal/vertical lines indicate the locations of the maximum-likelihood parameter combination, which were used to produce the simulated timings presented in Fig. 10. In many cases, the maximum-likelihood parameters differ significantly from the mean/median of the posteriors, indicating a complex multimodality in the parameter space. While we have shown in principle that the EMRI+rigid disk precession model is capable of producing timings consistent with QPE observations, this suggests that significant caution must be taken when interpreting parameter estimates from QPE data.

Appendix E Results across all simulations

To explore as wide a region of parameter space as computationally feasible, we ran a large grid of MCMC sampling runs for different EMRI/disk parameter combinations with and without disk precession included. For the static aligned-disk model, we sampled over every combination spanning $a\in\{100,200,400\}R_{g}$, $e\in\{0.01,0.1,0.3,0.5\}$, $\chi_{\bullet}\in\{0.5,0.9\}$, $\log(M_{\bullet}/M_{\odot})\in\{5,6,7\}$, and observing strategies from {$1\times 100$ ks, $4\times 100$ ks with 200 ks gaps, and $4\times 100$ ks with 1 Ms gaps}, for a total of $3\times 4\times 2\times 3\times 3=216$ combinations. For the misaligned precessing disk model, we included all of the above parameters, plus $T_{\rm disk}=\{5,20,50\}P_{\rm orb}$ and $\theta_{\rm disk}=\{5^{\circ},20^{\circ}\}$, for a total of $216\times 3\times 2=1296$ samples. All MCMC sampling was performed on NVIDIA A100 GPUs, using 1000 walkers $\times$ 20000 steps with a 10-second time resolution for the trajectory computation. The sampling runs with a static disk spanned between 1-1.5 hr, with an average duration of 70.5 minutes, while the sampling runs for a precessing disk spanned between 1.5-8 hr, with an average duration of 280.8 minutes. In total, this amounts to 6319 GPU-hours, which we carried out across 16 NVIDIA A100s. In Table LABEL:tab:errors we report the average 16th/50th/84th percentile parameter recovery, averaged over all combinations, to show the typical error-bars expected from QPE timing inference. We used uniform priors of $a\sim\mathcal{U}(30,500)R_{g}$, $e\sim\mathcal{U}(0,0.9)$, $i\sim\mathcal{U}(0^{\circ},90^{\circ})$, $\chi_{\bullet}\sim\mathcal{U}(0,0.998)$, $\log(M_{\bullet}/M_{\odot})\sim\mathcal{U}(4,8)$, $\theta_{\rm disk}\sim\mathcal{U}(0^{\circ},60^{\circ})$, and $T_{\rm disk}\sim\mathcal{U}(2,100)P_{\rm orb}$. The restrictions of $i\sim\mathcal{U}(0^{\circ},90^{\circ})$ and $\theta_{\rm disk}\sim\mathcal{U}(0^{\circ},60^{\circ})$ are made for simplicity, but they exclude regions of parameter space which are formally possible—retrograde EMRIs and highly-inclined, rigidly precessing disks. We note, however, that EMRI inclination has minimal effect on QPE timings (Section 3), and it is furthermore unclear whether extremely inclined disks can undergo rigid precession as unmodeled non-linear effects such as disk tearing (Nixon & King, 2012; Kaaz et al., 2023) become potentially important.