Testing EMRI Models for Quasi-periodic Eruptions with 3.5 yr of Monitoring eRO-QPE1

Quasi-periodic eruptions (QPEs) are luminous X-ray outbursts recurring on hour timescales, observed from the nuclei of a growing handful of nearby low-mass galaxies. Their physical origin is still debated, and usually modeled as (a) accretion disk instabilities or (b) interaction of a supermassive black hole (SMBH) with a lower mass companion in an extreme mass-ratio inspiral (EMRI). EMRI models can be tested with several predictions related to the short- and long-term behavior of QPEs. In this study, we report on the ongoing 3.5 yr NICER and XMM-Newton monitoring campaign of eRO-QPE1, which is known to exhibit erratic QPEs that have been challenging for the simplest EMRI models to explain. We report (1) complex, non-monotonic evolution in the long-term trends of QPE energy output and inferred emitting area; (2) the disappearance of the QPEs (within NICER detectability) in 2023 October, and then the reappearance by 2024 January at a luminosity of ∼100× fainter (and temperature of ∼3× cooler) than the initial discovery; (3) radio non-detections with MeerKAT and Very Large Array observations partly contemporaneous with our NICER campaign (though not during outbursts); and (4) the presence of a possible ∼6 day modulation of the QPE timing residuals, which aligns with the expected nodal precession timescale of the underlying accretion disk. Our results tentatively support EMRI–disk collision models powering the QPEs, and we demonstrate that the timing modulation of QPEs may be used to jointly constrain the SMBH spin and disk density profile.


INTRODUCTION
Quasi-Periodic Eruptions (QPEs) are a rare class of extragalactic transient characterized by high-amplitude, repeating soft X-ray flares with recurrence times of hours.They were first discovered in the nucleus of GSN 069 (Miniutti et al. 2019), and since then there have been five more confirmed QPE hosts (Giustini et al. 2020;Arcodia et al. 2021Arcodia et al. , 2024) ) and two candidates (Chakraborty et al. 2021;Quintin et al. 2023).Related properties, including recurring large X-ray outbursts and super-soft thermal spectra, are also observed from a handful of other supermassive black holes (SMBHs) (Terashima et al. 2012;Tiengo et al. 2022;Evans et al. 2023;Guolo et al. 2023).
QPEs have been observed so far from the nuclei of nearby low-mass galaxies (M * ≈ 10 9−9.5 M ⊙ , Arcodia et al. (2021)), hosting black holes of M BH ≈ 10 5−6.6 M ⊙ (Wevers et al. 2022).Contrary to other kinds of Xray transients which show multi-wavelength counter-parts, the flares from QPEs have so far been observed only in the soft X-rays (though it is possible that other wavelengths are dominated by the host galaxy or the disk itself).Their spectra are consistent with the exponentially decaying Wien tail of a blackbody (kT ≈ 100−200 eV), with soft X-ray luminosities reaching up to 10 42 −10 43 erg s −1 (∼ 0.01−0.1LEdd ).In eruption, these sources increase their count rate by up to 10-150 times compared to the quiescence state in between flares.Quiescence, when detected, is ultra-soft (kT ≈ 50 − 80 eV) and it is interpreted as the emission from the accretion disk around the central massive black hole.
Currently, the origin of QPEs is still debated.Existing models separate broadly into (a) recurring limit-cycle instabilities within the SMBH accretion disk (Raj & Nixon 2021;Pan et al. 2022Pan et al. , 2023;;Kaur et al. 2023;Śniegowska et al. 2023); or (b) the interaction of the SMBH with a lower-mass orbiting companion, which allows many different prescriptions for precisely how the X-ray emission is produced (Xian et al. 2021;Suková et al. 2021;King 2022;Metzger et al. 2022;Krolik & Linial 2022;Zhao et al. 2022;Linial & Sari 2023;Franchini et al. 2023;Linial & Metzger 2023a;Tagawa & Haiman 2023;?).The latter class of models has generated much excitement, as they may reveal QPEs to be the firstobserved electromagnetic counterparts of extreme massratio inspirals (EMRIs) (Arcodia et al. 2021), a source of millihertz gravitational waves in the sensitivity band of upcoming space-based gravitational wave detectors (though the detectability of the thus-known population is limited, Chen et al. (2022)).
Finally, multiple lines of evidence have linked QPEs to Tidal Disruption Events (TDEs).The QPEs in GSN 069 and eRO-QPE3 are seen alongside a long-term decay of the quiescent flux level (Miniutti et al. 2019;Arcodia et al. 2024).Two QPE candidates were also discovered in TDEs, which were initially detected by XMM Slew (Chakraborty et al. 2021) and the Zwicky Transient Facility (Quintin et al. 2023).GSN 069 shows an abnormal C/N ratio (Sheng et al. 2021) and a compact nuclear [O III] region suggesting a young accretion system (Patra et al. 2023).Although optical spectra of QPE hosts show signatures for nuclear ionization, no QPE host nuclei are associated with broad optical/UV emission lines, which may suggest the disks around their SMBHs are too compact to support a mature broad line region (Wevers et al. 2022).Recent modelling efforts have begun to explicitly account for the coincidence with TDEs to explain the QPE phenomenon (Linial & Metzger 2023a;Franchini et al. 2023), as the rarity of TDEs themselves is an extremely strong constraint of the prerequisite conditions needed to produce QPEs.Contin-ued monitoring of known sources (Miniutti et al. 2023b) and late-time follow-up of TDE candidates to identify unusual X-ray variability will be a powerful tool in determining the relationship between TDEs and QPEs.
With the initial sample of QPE sources now growing to multi-year observational baselines, we are able to begin studying the secular evolution of these sources, which allows important tests of theoretical models.Miniutti et al. (2023b) and Miniutti et al. (2023a) studied the 12-year XMM dataset of GSN 069 (with QPEs appearing in the past 5 years), finding significant longterm variability in the quiescence luminosity.They also found coupling of the QPE properties to the quiescent accretion disk: as accretion rate increases, the QPE temperature and amplitude both decrease, with the QPE temperature asymptotically approaching the disk temperature.Giustini et al. (in prep) found modest variations of the quiescent luminosity of RX J1301.9+2747, and the constant presence of QPEs with a complex pattern of variability, over more than 20 years of X-ray observations.
Here, we report on the 3.5-year NICER and XMM-Newton dataset of eRO-QPE1, a source initially reported in Arcodia et al. (2021).The strikingly complex flare patterns and unpredictable recurrence times were immediately noteworthy, and the comparatively erratic QPE behavior has presented a challenge for the simplest EMRI-based models which predict roughly regular outbursts on the orbital timescale.The initial XMM observations of this source were studied in-depth by Arcodia et al. (2022), who found evidence for overlapping QPEs and a similar L − kT hysteresis pattern as seen in GSN 069.Our long-baseline NICER dataset allows us to probe secular evolution of the QPEs and their timing properties across three years and a total of 92 flares.In Section 2 we report on our analysis procedures, including a non-standard NICER data reduction process, and subsequent analysis methods.In Section 3 we report the results of time-resolved spectroscopy and timing, as well as implications for theoretical models.
2. METHODS 2.1.NICER light curves and spectra eRO-QPE1 was observed in 141 NICER observations (PI: Arcodia) for a total of 964.2 kiloseconds between August 2020-January 2024.The data were processed using HEAsoft v6.32.1 and NICERDAS v11a.The standard procedure of using nicerl3 to generate source light curves was insufficient here, due to the low signal/noise ratio and variable background; thus we devised a nonstandard procedure, detailed here.Reliably estimating light curves for faint sources like eRO-QPE1, in which the source count rate is generally comparable to (or less than) the background, presents a challenge for NICER.In most observations, the combination of background counts produced by optical loading (at low energies) and particle interactions with the Si detectors (at high energies) produces a higher count rate than the source itself.This background level can also vary slightly within observations, creating a degeneracy of the observed total variability between source and background components-because NICER is a nonimaging instrument, variability cannot be resolved directly into source and background components, instead requiring other metrics to estimate total background contribution.We approach the problem using timeresolved spectroscopy by splitting the data into many Good Time Intervals (GTIs), then spectroscopically estimating the source and background components separately in each GTI.
First we use nimaketime, with custom filtering choices of unrestricted undershoot (underonly range=*-*) and overshoot rates (overonly range=*-*), as well as per-FPM and per-MPU autoscreening disabled to prevent unnecessarily aggressive event filtering.We then split the intervals produced by nimaketime into GTIs of maximum 200 seconds to allow time-resolved estimation of the variable background.In each GTI, we manually discard focal plane modules (FPMs) with 0-0.2 or 5-15 keV count rates > 4σ higher than the average across all GTIs within the OBSID, or above an absolute threshold of 20 counts sec −1 .These cutoffs were chosen because eRO-QPE1 is super-soft (kT max ∼ 120 eV) and faint (∼2 counts sec −1 at peak).Thus the 5+ keV band is entirely background-dominated (so 20+ counts sec −1 signifies extremely high particle-induced background), whereas the 0-0.2 keV band is undershoot-dominated (so 20+ counts sec −1 signifies extremely severe light leak conditions).
After screening the event lists, we used the SCORPEON 1 template-based background model to estimate the contribution from astrophysical and non X-ray backgrounds separately for each GTI.We fit the entire broadband (0.2-15 keV) array counts with the PyXspec2 interface to XSPEC (Arnaud 1996).We leave the Solar Wind Charge Exchange (SWCX) oxygen emission line normalizations free to vary to account for partially ionized oxygen fluorescence from the solar wind (which is particularly im-portant during NICER dayside GTIs).Along with the SCORPEON background we fit each GTI with a source model represented by tbabs×zbbody (motivated by the pure thermal-like emission spectrum of the source, Arcodia et al. ( 2021)) with z = 0.0505 fixed to the host galaxy redshift and N H = 2.23 × 10 20 cm −2 fixed to the lineof-sight galactic column density, resulting in a distribution of broadband fit statistics (cstat/dof) with µ = 1.2, σ = 0.3.We sum only the counts contained within the zbbody model to create the background-subtracted light curves presented in Fig. 1.We consider a source "detection" as any GTI in which the blackbody normalization is >1σ inconsistent with zero, i.e. a non-background component is required by the fit at the 1σ level.See Appendix A for examples of robust and marginal detections, as well as further discussion of the source and background spectral components.
We note that NICER was not able to detect the quiescent disk emission in eRO-QPE1 at any epoch.As eRO-QPE1 is the most distant known QPE (z = 0.0505), the quiescence is detected only at 1σ by XMM (Arcodia et al. 2021) despite a similar L X to other QPE hosts.
We also perform time-resolved spectroscopy on relative-intensity bins to track the evolution of the QPE bolometric luminosity (L bol ), temperature (kT ), and the inferred blackbody radius (R bb ) over the different phases of each flare (Fig. 3).We followed the same filtering and background estimation procedures described above, then grouped all GTIs with a source detection into 5 relative-intensity bins for each QPE, based on its flux compared to peak: 1-80% F peak , 80-100%, 100-80%, 80-50%, and 50-1%.We choose 2 rise bins and 3 decay relative-intensity bins, because the faster rise means they are typically more poorly sampled.Not all QPEs are well-sampled across all 5 segments, because NICER did not always fully sample the rise and decay (Figs. 3, B.2).

XMM-Newton
Upon observing the possible re-emergence of QPEs in Nov 2023 with NICER (Fig. 1), we requested a 124 ks Director's Discretionary Time (DDT) observation with XMM-Newton, which was taken on Jan 5-6, 2024.The data were reduced with standard tools and prescriptions (XMM SAS v20.0.0 and HEAsoft v6.29).For eRO-QPE1, source products were extracted within a circle of 13.8" to avoid a nearby contaminant, while background was extracted from a source-free area within a circle of 39".We restrict our analysis to the EPIC-PN instrument due to its higher count rate compared to the MOS CCDs.Source and background rates were then extracted, and grouped into 2000-second bins to produce the background-subtracted light curve in Fig. 2. Phaseresolved spectra were extracted with evselect, represented by the differently-colored regions in Fig. 2 (quiescence comprises the beginning and ending segments of the observation).We analyzed the spectra using XSPEC v12.13.1, fitting with the Cash statistic and using data from 0.2-1.5 keV (as the source is backgrounddominated at higher energies).We fit spectra using the same model as the NICER data (tbabs×zbbody) discussed in Sec.2.1, with fitting results and comparison to previous XMM-Newton results (Arcodia et al. 2021) reported in Table 1.

Flare profiles
NICER has observed eRO-QPE1 for several observing epochs lasting 6-15 days over the past 3.5 years.QPEs are visible in all epochs except Oct. 2023 (but have since emerged in Nov 2023), with complex, non-monotonically varying properties including QPE luminosity, recurrence time, duration, temperature, and the intra-epoch scatter of these properties.To allow for precise characterization of these properties, we fit each observed QPE with the following exponential flare model motivated initially for long gamma-ray bursts (Norris et al. 2005) and adopted for eRO-QPE1 in Arcodia et al. (2022): where A is the flare amplitude; t peak is the time of peak flux; τ 1 , τ 2 are the e-folding times of rise and decay, respectively; λ = e √ τ1/τ2 is a normalization to join rise and decay; and t as = √ τ 1 τ 2 sets the asymptote time such that flux= 0 for t < t peak − t as .We compute integrated energy outputs of the flares (Section 3.2) by integrating from (t peak −3τ 1 ) to (t peak +3τ 2 ), i.e. the total energy contained within ±3 e-folds of the flare peak.
The observed flares are generally well-described by this exponential profile (Fig. 1), and we use these fits to produce most of the results below.Arcodia et al. (2022) found that one XMM observation of eRO-QPE1 (OBSID 0861910201) displayed multiple overlapping bursts, making it the only QPE source showing this behavior.One NICER flare clearly shows this behavior (MJD 59413 in Jul 2021).The presence of at least two overlapping QPEs in eRO-QPE1, between XMM 2020 and NICER 2021, presents an interesting constraint for EMRI models, perhaps pointing towards a high-eccentricity EMRI orbit or unusual disk geometry.It is worth noting that the double-flare on MJD 59413 occurred unusually long after the previous burst a MJD 59412-a difference of 1.28 days, which is slightly under twice the mean recurrence time of nearby epochs.Thus, a geometry in which one burst was missed, then appeared alongside the next QPE, is possible.The more poorly-sampled QPE at MJD 59900 may be another instance of this double-flare behavior associated with a longer recurrence time, though that case is not as clear.

Radio nondetections
In addition to the NICER campaign, we undertook radio monitoring of eRO-QPE1.Our campaign consisted of 2 × 1 hour observations with MeerKAT in Sep 2021, and 7 × 20 minute observations with the Karl G. Jansky Very Large Array (VLA; see §2.4.2 for details), partially concurrent with our Jan 2023 NICER campaign.No observations resulted in a radio detection, though we note that unfortunately none of our snapshots were taken during an eruption (Fig. 1).The 3σ upper limits are listed in Table 2.4.2.Our most constraining radio limit comes from the stacked image of all seven VLA observations; 3 − σ upper-limit of 15.6 µJy at 6 GHz.

MeerKAT observation details
The MeerKAT radio telescope is a 64-dish interferometer based in the Karoo Desert, South Africa.We obtained two observations through a Director's Discretionary Time proposal (PI: R. Arcodia, DDT-20210908-RA-01).The observations were made on the 11 th and 19 th September 2021, each lasting 2.33 hours.MeerKAT observes at a central frequency of 1.28 GHz with a total bandwidth of 0.86 GHz.For each observation, we started with five-minutes observing the flux and bandpass calibrator J0408-6545 followed by switching between the target for 40 minutes and the phase calibrator J0240-2309 for two minutes.
Both epochs were reduced using oxkat, a set of python scripts specifically designed for the reduction of MeerKAT observations (Heywood 2020).The scripts apply flags to the calibrator fields before calculating delay, bandpass and complex gain calibrations using both calibrator sources and applies them to the target, all of which is done in casa.The target field is split out and flagged then imaged using wsclean.For these observations, we also performed phase-only self-calibration.
We did not detect any radio emission at the position of eRO-QPE1 in either epoch.We provide 3σ upper limits in Table 2.4.2.

VLA observation details
To efficiently sample the X-ray QPE cycle, we ran simulations based on previous X-ray light curves of eRO-QPE1 to estimate the combination of observation num-ber/separation that gave us the highest probability of a VLA observation occurring during an X-ray eruption.The result was 7 observations separated by 5-15 hours.Therefore, we observed eRO-QPE1 with the VLA (Project Code: 23A-059) for 7 epochs in 2023 January.The array was in the B configuration at the time of all of our observations.We used the 8-bit samplers, observing in C (4 -8 GHz) band, comprised of 2 base-bands, with 8 spectral windows of 64 2-MHz channels each, giving a total bandwidth of 1.024 GHz per base-band.We reduced and imaged the data within the Common Astronomy Software Application (casa version 5.6), using standard procedures outlined in the casaGuides3 for VLA data reduction (i.e., a priori flagging, setting the flux density scale, initial phase calibration, solving for antenna-based delays, bandpass calibration, gain calibration, scaling the amplitude gains, and final target flagging).For all observations, we used 3C147 (J0542+498) as a flux/bandpass calibrator, and J0239-0234 as a phase calibrator.Imaging was performed with natural weighting to maximize sensitivity.

Short-term and secular evolution of the inferred emitting region
QPEs show differing evolution at different energies, with narrower profiles and earlier peaks in higher-energy bands (Miniutti et al. 2019;Giustini et al. 2020;Arcodia et al. 2021;Chakraborty et al. 2021;Quintin et al. 2023;Arcodia et al. 2024).Said another way, Miniutti et al. (2023b) found that QPEs in GSN 069 show L−kT hysteresis, i.e. the blackbody temperature peaks early and begins to decay before the luminosity peaks.This is physically consistent with the inferred blackbody radius  Miniutti et al. (2023b).Over the long-term, the peak inferred R bb appears to grow (from Aug 2020 to Feb 2022), then shrinks (by Oct-Nov 2022).This may be a geometric (viewing-angle) effect attributed to the EMRI apsidal precession (∼ 10s of days), or a change in the mutual EMRI-disk inclination due to the EMRI nodal precession (∼ 100s of days) resulting in changing ejecta properties (see discussion in Section 4.2).Bottom: Same plot for kT .
(assuming a spherical source): expanding over time (shown in their Fig.18).The QPEs in GSN 069 are consistent with an initial R bb ≈ R ⊙ , expanding by a factor ∼ 3 over the course of the flares (Miniutti et al. 2023b).Arcodia et al. (2022) found a similar hysteresis pattern is obeyed in the XMM outbursts of eRO-QPE1 (their Figs.7 & 8).The longer temporal coverage of NICER compared to XMM allows us to track the average evolution of the inferred blackbody emission over the course of many QPEs, as well as how this quantity evolves in the long-term-an important constraint for distinguishing between theoretical models.
The early NICER observations of eRO-QPE1 broadly agree with these trends, but this does not hold true at later epochs.In Aug 2020 and Aug 2021, the inferred blackbody radius R bb indeed appears to grow over the course of each QPE by a factor ∼ 2−3 while the emission temperature cools by ∼50% (Fig. 3, also Figs.B.2 & B.3).Then, there is long-term evolution in the typical value of R bb , which peaks around 6 × 10 10 cm in Aug 2020 but increases up to > 10 11 cm by Feb 2022, before returning to the previous value in subsequent epochs.The typical rate of change dR bb /dt correlates with the typical peak R bb in each epoch (though this cannot be directly interpreted as the physical expansion speed of the emission region; see Section 4.2).The Aug 2020 and Aug 2021 epochs, which observed 15 QPEs each, show remarkable consistency in the evolution of R bb (as well as the flare profiles themselves), but this grows more erratic in later epochs, with larger short-term variation in typical QPE profiles, luminosities, R bb , etc. epochs are also where R bb appears to invert its evolution (Fig. 3).The ∝ T 4 lines are plotted for visualization only, and are not fits to the data.

Secular evolution of other QPE properties
Miniutti et al. (2023a) noted the QPEs of GSN 069 roughly follow L bol ∝ T 4 at peak, which suggests that the X-rays comprise the majority of the emission of QPEs and that the peak R bb is roughly constant across flares.The peaks of eRO-QPE1 do not appear to obey this trend across all epochs (Fig. 4).At the lower luminosities typical of the Oct-Dec 2022 epochs, this relationship is less precisely followed, indicating the emission at lower L bol may not be consistent with pure black- body, or a significant fraction is emitted below soft X-ray energies.Moreover, as clearly visible in Fig. 3, the typical peak R bb is far less consistent in the Oct-Dec 2022 epochs, also contributing to the breaking of the L ∝ T 4 trend.
The increasingly erratic behavior of R bb (Fig. 3) is also accompanied by an overall decrease in the typical QPE energy output (Fig. 5), as well as growing irregularity in the flare recurrence times (Fig. 6).The typical peak R bb increases by a factor of 3 over these years, before returning to its initial level; and the expansion speed dR bb /dt grows with the peak R bb .We do not see any evolution in the typical decay times over the 3.5-year baseline (Fig. C.5), though the NICER sampling cannot constrain the rise times precisely enough to make any comment about their change.
The single QPE detected by XMM-Newton in Jan 2024 (Fig. 2) appears at a flux ∼10x lower than the most recent NICER-detected flares of 2023(Fig.4), and ∼100x fainter/2.8xcooler than the previous XMMdetected QPEs of July-August 2020 (Arcodia et al. 2021).The most recent XMM flare also exhibits an expanding emitting region, though its peak at 2.2 × 10 10 cm (Fig. 2) is a factor of a few smaller than the typical R bb of the NICER QPEs (Fig. 3).Due to the large uncertainties of the quiescence spectra, the quiescence L bol and R bb are 3σ consistent with the previous XMM detection in July-August 2020 (Table 1).This is in contrast to GSN 069, where a drop in the QPE luminosity was seen alongside a significant rise in the quiescence luminosity (Miniutti et al. 2023b).

3.3.
Evidence for a ∼6-day modulation in the QPE timing residuals

QPE timing analysis methods
After fitting the flares as described in Section 2.3, we use the computed t peak values to perform an O-C (Observed Minus Computed) analysis to search for any significant patterns in the timing behavior of the QPEs.The O-C technique, which is best known for its use in discovering the orbital period decay of the Hulse-Taylor binary pulsar (Taylor & Weisberg 1982), computes the timing deviation of a recurring signal with respect to a constant period.The residuals are computing by subtracting the expected (Computed) arrival time of the N th strictly periodic repetition (with period P const ) after some fixed starting time, from the Observed actual arrival time of that repetition.In the case of gravitational-wave driven orbital decay, these O-C values are fit with a quadratic because of Ṗ < 0; thus, the actual transit times arrive successively earlier and earlier than expected from naive extrapolation of P const .
In our case, we separately compute an O-C within each epoch because we cannot phase-connect across themi.e., we have no way to tell how many QPEs have occurred between epochs, and the average recurrence time also changes significantly (∼ 30%) between them (Fig. 6) for reasons not yet physically understood.Nevertheless, we are interested in probing the short-timescale behavior of the QPE timings independently of the secular evolution.Thus, on the timescale of individual NICER epochs, we assume a P const given by the average re-currence time within that epoch, then compute O-C as described above.
As the O-C timing residuals are roughly sinusoidal (Sec.3.3.2),we perform a least-squares fit to the residuals with a sine+constant model.We choose not to apply a Fourier-based technique to estimate the power spectral density (PSD) for two reasons: 1) our NICER observing windows of 6-15 days limit us to only 1-2 periods per epoch, and 2) as the QPEs "sample" the sinusoidal profile only on their recurrence time (∼1 day), we have too few data points to make a meaningful estimate of the PSD.Comparing our sine+constant fit to a constantonly fit allows us to make a rough estimate of the recurrence period while remaining in the time domain.We fit each epoch individually, then all epochs together with a fixed overall period (but amplitude, phase, and constant offset free to vary within each epoch) to assess improvement over a constant model hypothesis.Our fit results are presented in Section 3.3.2,and the implications are discussed in Section 4.3.

Results of O-C analysis
The bursts of eRO-QPE1 show significant scatter compared to other QPEs (Arcodia et al. 2022).To assess the degree of irregularity in the NICER data, we folded the QPEs within each epoch into segments of three consecutive bursts, then overplotted each of these segments beginning at the time of the first burst (Fig. 7).The QPE timings deviate by ∼ 50% from the average recurrence time computed in each epoch (compared to ∼ 10% in GSN 069 and eRO-QPE2).
The recurrence times do not correlate with any other properties of the flare (luminosity, rise time, and decay time), nor do those properties strongly correlate with each other (see Appendix C).Moreover, the average recurrence time varies across epochs, appearing to increase slightly from 0.76 days in Aug 2020 to 1.03 days in Jan 2023 (Fig. 6).
We also report a tentative finding in the timing structure of the bursts (based on the procedure described in Section 3.3.1):the timing residuals (i.e.bursts arriving early/late compared to a strict period) exhibit possible periodic structure within several NICER epochs, at an average period of 5.73 days (Fig. 7).The bestfit period of the timing residuals are fairly consistent for most epochs (ranging from 5.2-6.1 days).Our fits of a sine+constant model provide a significant improvement over a constant model hypothesis, with an average ∆χ 2 = 935 in the per-epoch fits (gray dashed lines in Fig. 7) and a total ∆χ 2 = 5881 in the global fixed-period fit (blue dashed lines).
The Aug 2020 and first Dec 2022 epochs are significant outliers, with best-fit super-periods of 9 days and 2.5 days respectively.Thus we cannot definitively claim a ≈ 6 day super-period, but it is noteworthy, and followup observations will provide more clarity.Interestingly, a super-orbital period of ∼6 days has been predicted in some EMRI precession models, which we will describe in the discussion Section 4.3.

DISCUSSION
Our key observational findings in Section 3 are: • The QPEs show secular evolution in R bb , total energy output, and recurrence times (Figs. 3, 5, 6) over the epochs from Aug 2020-Jan 2023.As seen in other sources, the QPEs are consistent with an emission region that grows in size and cools in temperature over the course of a flare.The characteristic peak R bb changes by a factor ∼2 over a timescale of years, and the emitting area evolution becomes increasingly erratic.
• QPEs are entirely undetected by NICER in the Oct 2023 epoch, but the peak of one flare is seen in Nov 2023 (Fig. 1).We confirmed the ongoing presence of QPEs with XMM-Newton in Jan 2024 (Fig. 2) at a luminosity ∼ 10x lower than recent NICER detections (Fig. 4), and ∼100x lower than the initial discovery (Arcodia et al. 2021) • The QPE recurrence times in most epochs are modulated on a ∼6-day super-period, and the long-short recurrence pattern of other QPEs is generally not seen with the exception of a brief window in early Dec 2022 (Fig. 7).

QPE models
There are a large number of models proposed to explain the origin of QPEs, fitting broadly into two categories.The first class is recurring limit-cycle instabilities within the SMBH accretion disk (Raj & Nixon 2021;Pan et al. 2022Pan et al. , 2023;;Kaur et al. 2023;Śniegowska et al. 2023), inspired by observations of recurring highamplitude outbursts (the "heartbeat" states) seen in black hole X-ray binaries GRS 1915+105 (Belloni et al. 2000) and IGR J17091-3624 (Altamirano et al. 2011).In fact, Wang et al. 2024 reported high-amplitude variability at 0.5 Hz, which scaled up for a 10 5 M ⊙ SMBH corresponds to ∼ 2 days, similar to the QPE timescales.However, the largest weakness of these models is they make fewer testable predictions due to our poorer understanding of accretion instabilities in the radiationpressure dominated regime.

Overplotted even/odd QPEs O-C timing residuals
Overplotted even/odd QPEs O-C timing residuals The second class is interaction of the SMBH with a lower-mass orbiting companion in an extreme mass-ratio inspiral (EMRI), whether by direct accretion/tidal stripping of the companion (King 2022;Metzger et al. 2022;Krolik & Linial 2022;Zhao et al. 2022;Linial & Sari 2023;Lu & Quataert 2023) or its interaction/collision with the SMBH accretion disk (Xian et al. 2021;Suková et al. 2021;Franchini et al. 2023;Tagawa & Haiman 2023;Linial & Metzger 2023a).EMRI-disk collision models have thus far been most successful at reproducing QPE observational properties within GSN 069 and eRO-QPE2, the two comparatively "well-behaved" QPE sources.Moreover, recent iterations of EMRI-disk models (Franchini et al. 2023;Linial & Sari 2023) suggest the underlying accretion disk is comprised of debris from the disrupted star following a TDE-thus explaining the suggestive coincidence of QPEs with TDE hosts (Chakraborty et al. 2021;Miniutti et al. 2023b;Quintin et al. 2023), as well as the compact emission radii and hot blackbody temperatures observed (Miniutti et al. 2023a).Thus, in the following discussion, we discuss our results in light of these EMRI-disk collision models.
GSN 069 and eRO-QPE2 obey a quasi-period to within 10% (Miniutti et al. 2019;Arcodia et al. 2021), and also a characteristic "long-short" recurrence time (Miniutti et al. 2023b), which provided motivation for models to consider an EMRI companion on a mildly eccentric orbit (Franchini et al. 2023;Linial & Metzger 2023a).Within this picture, the alternating trend can be conveniently produced by associating the long recurrence times with the orbiter's passage through apocenter, and the short recurrence time through pericenter.Each EMRI-disk collision then ejects part of the disk mass, which is shock-heated by the supersonically moving EMRI.The disk ejecta, which is radiationpressure dominated (as expected for the typical Ṁ of the quiescent disk), expands adiabatically over the QPE timescale, resulting in the observed declining kT and increasing R bb (Fig. 3).
eRO-QPE1, on the other hand, has proven difficult for these models to explain.In previous studies (Arcodia et al. 2021(Arcodia et al. , 2022)), this source has shown 1) significantly larger, apparently unpredictable scatter in burst recurrence times (up to 50%, Fig. 7); 2) apparently no long-short recurrence pattern; and 3) occasional overlapping double-flares (Section 2.3), raising concerns for the simple picture of a mildly eccentric EMRI generating QPEs twice per orbit.We may consider that the above models imply the QPE timing should be affected by several precession frequencies (as also discussed in Linial & Metzger (2023a) and Franchini et al. (2023)): the underlying accretion disk should precess nodally via gravitomagnetic frame dragging (Lense & Thirring 1918;Kato 1990;Merloni et al. 1999), and the EMRI companion itself should also precess apsidally (P ∼tens of days) and nodally (P ∼hundreds of days) (Franchini et al. 2023;Linial & Metzger 2023a).The precise observational im-plications of these precession effects has not yet been explored extensively, and we do so here.

Secular evolution
The long-term evolution in the flare properties L, kT , R bb , and dR bb dt (Figs. 3,4,5) can be explored within the framework of EMRI-disk interaction models, where the QPE emission is generated by shocked disk material from collisions with the companion and the disk.
The amount of ejected disk material M sh is set by the interaction cross-section of the companion and disk, which should depend on the geometric cross-section of the companion if it is a star (Suková et al. 2021;Linial & Metzger 2023a) or the Bondi-Hoyle radius if the companion is a compact object (Franchini et al. 2023).There should also be a dependence on the relative inclination of the companion and the disk (in the face-on impact scenario, the orbiter collides with the minimum amount of disk mass, whereas for more highly inclined impacts there is a longer impact duration, thus more ejected mass).The resulting ejecta will have a kinetic energy ∝ ⃗ v 2 sh , where ⃗ v sh is the relative velocity between the orbiter and the disk.⃗ v sh also depends on the relative EMRI-disk inclination: for face-on impacts, the (Keplerian) orbital velocity of the disk material is entirely perpendicular to the EMRI velocity, but in general the vector difference ⃗ v sh = ⃗ v disk − ⃗ v EM RI will have an inclination-dependent magnitude.The ejecta then expands adiabatically (driven by radiation pressure from the significant energy within the shock-heated disk material), cooling over time as its optical depth declines; this expansion of the shocked ejecta causes the observed increase in R bb (Fig. 3, Miniutti et al. (2023b)).
It is still an open question whether the emission itself is dominated by blackbody emission, bremsstrahlung, Compton up-scattering, or some other radiative process.Linial & Metzger (2023a) noted that the rapid timescale of QPEs may be too short for the ejecta to reach thermal equilibrium, with photon production instead dominated by free-free emission.It is difficult to clearly distinguish the two with only X-ray observations, but deep longerwavelength observations may eventually be able to identify the dominant mechanism.In either case, fitting a blackbody model to the X-ray data allows us to infer a rough scale of the physical emitting region (which, in the bremsstrahlung case, does not correspond to the actual physical extent of the photosphere).
From Fig. 3 we observe that the typical R bb , as well as expansion rate dR bb dt , appears to increase over the epochs, which appears to suggest an increasing ⃗ v sh (thus a changing relative inclination).We cannot exactly use this to constrain the relative EMRI/disk inclination an-gles, as there are a number of different configurations (with a prograde/retrograde disk/EMRI) that create the same trend.One natural way to explain a change in the relative velocity is by invoking nodal precession of the disk and nodal+apsidal precession of the companion.We discuss the observable effect of these timescales further in Section 4.3.
The most recent Jan 2024 XMM-Newton-detected QPE is a factor ∼100x fainter than initial discovery (Arcodia et al. 2021).Moreover, the peak temperature has decreased by a factor ∼3x (Table 1).It is worth noting, if this trend continues, QPEs may continue while ceasing to be detectable at X-ray wavelengths, instead appearing as an ultraviolet transient (Linial & Metzger 2023b).

Order in the recurrence times
One natural explanation for the scatter in recurrence times within the EMRI-disk interaction picture is (independent) precession of the EMRI and the disk.The GR effects at the hours-long periods of QPEs should significantly modulate the EMRI and disk orbits, via nodal precession (induced by gravitomagnetic frame dragging; Lense & Thirring (1918); Kato (1990); Merloni et al. (1999)) and apsidal precession on observable timescales.We first discuss the expected disk nodal precession period for a given set of system parameters, following the arguments of Franchini et al. (2016) and Franchini et al. (2023).Assuming the angular momentum of the disk to be misaligned with respect to the SMBH spin axis (which is generally expected for TDE disks, for an isotropic distribution of disrupted stars), frame dragging causes differential precession of the disc radii.If the disc viscosity is smaller than the disc thickness, the propagation of these disturbances occurs in the so-called 'bending waves' regime, allowing the disk to rigidly precess (Franchini et al. 2016) on a timescale that depends on the SMBH mass (M BH ) and spin (a), as well as the disk outer radius (R out ) and surface density profile.
Assuming the disk has a power-law surface density profile, with index p and normalization Σ 0 : and the nodal precession frequency at a given orbital radius is given by:

R Rg
3/2 + a the total disk precession frequency is Ω LT (R) integrated along the radial extent of the disk, weighted by the angular momentum where Ω K = (GM BH /R 3 ) 1/2 is the Keplerian orbital frequency: Observationally, the effect of the disk nodal precession is to change the relative inclination between the disk and EMRI companion, thus delaying or advancing the impact/interaction times resulting in QPEs.If we assume our reported QPE timing residuals (Fig. 7) are associated with the disk LT precession, for M BH ∼ 10 5.8 M ⊙ (Arcodia et al. 2021), M D = 2M ⊙ , and R out = 100R g , we can constrain (although degenerately) the SMBH spin a and disk surface density profile Σ(R) (Fig. 8).It is important to note that these constraints depend sensitively on M BH and R out , increasing either of which will increase the inferred disk precession period.
The timing residuals are not strictly sinusoidal (Fig. 7), as expected due to the presence of other relevant frequencies which also modulate the burst arrival times.For example, three possible effects are light travel-time delays from apsidal precession of the EMRI (P ∼tens of days), changing relative EMRIdisk inclination due to nodal precession of the EMRI (P ∼hundreds of days) (Linial & Metzger 2023a;Franchini et al. 2023), or variable disk surface density due to stochastic mass accretion rate fluctuations as proposed for X-ray binaries (Ingram & Done 2011).Still, the hierarchy of precession timescales is such that, over individual NICER epochs (∼10 days) the disk nodal precession should dominate, resulting in the observed ∼6-day period (largely consistent with the assumed 7.5-day nodal precession period used for the model of eRO-QPE1 in Franchini et al. (2023)).
The apparent modulation of the QPE arrival times on a timescale consistent with the predicted disk precession period supports QPEs being generated by some direct interaction with the accretion disk (Suková et al. 2021;Xian et al. 2021;Linial & Metzger 2023a;Franchini et al. 2023;Tagawa & Haiman 2023).Other models that propose that QPEs are fed by the companion itself (i.e.where the quiescent disk does not play a part in modulating the timing of the burst) require some other tuning to account for the erratic modulation of the QPE period.The precise radiative processes generating the emission following the EMRI-disk interactions are still uncertain, and will require further observations to constrain.
Another observational consequence of disk precession is modulation of the quiescent emission (it is this reason, for example, that QPOs in accreting stellar-mass black holes are sometimes interpreted in light of rigid precession).However, since the quiescent level of eRO-QPE1 is below NICER detectability, we are unable to verify this.The presence of a quiescent-level QPO in GSN 069, during epochs showing QPEs, is worth noting (Miniutti et al. 2023b).
One possible shortcoming of the disk precession interpretation is the alignment timescale of the disk with the SMBH spin (Franchini et al. 2016).As the disk cools (lower Ṁ ) and becomes thinner, it eventually violates the conditions necessary for warp propagation in the bending waves regime, therefore no longer allowing rigid precession.In addition, the intrinsic viscosity of the disk should also damp oscillations.The precession amplitude should thus decay over a timescale of ∼years, with the exact value dependent on the disk viscosity, SMBH spin, and initial misalignment angle.As these quantities are all fairly uncertain, matching the observed years-timescale precession decay would require some fine-tuning.One interpretation of this is that the accretion disk in eRO-QPE1 is dynamically younger than those in e.g.GSN 069 or eRO-QPE2, where the recurrence pattern is comparatively well-behaved and shows the long/short recurrence of a mildly eccentric orbit, indicative of a relatively stable disk configuration.Indeed, the QPEs in GSN 069 appeared up to 8 years after the appearance of a TDE-like increase in quiescent flux (Miniutti et al. 2019), suggesting the accretion system had time to evolve.No such constraints are available for eRO-QPE1, eRO-QPE2, or eRO-QPE3, which were first discovered during their QPE-exhibiting phases.In eRO-QPE4, the quiescence must have bright-ened with or after the first-detected QPEs, inidicating a short-lived nature/evolution of the accretion flow (Arcodia et al. 2024).
The Aug 2020 and first Dec 2022 epochs are significant outliers from the ∼6-day super-period (Fig. 7).The Aug 2020 timing residuals appear roughly sinusoidal, but with a significantly longer P = 8.96 days.It is unclear what causes this deviation, though one possible explanation is a changing disk surface density profile (a shallower decay increases the LT precession period for a fixed SMBH spin, Fig. 8), or viscous spreading of the disk resulting in a longer precession period.The first Dec 2022 epoch shows a brief long-short recurrence phase, seen in other sources but generally not in eRO-QPE1, possibly indicative of a lower eccentricity e ≈ 0, but with large uncertainties considering other deviations due to precession effects or secular evolution.It may be that this epoch occurred during a particular disk/EMRI configuration resulting in double-impacts, but was short-lived, thus returning to the ∼6-day period in successive Dec 2022 and Jan 2023 epochs.Still, we emphasize that these explanations are speculative, and the reason these epochs are outliers is very unclear.An alternative explanation to the disk/EMRI precession picture is the presence of a third body in the system, e.g. a second solar-mass companion in orbit around the SMBH which modulates the orbital frequency/configuration between epochs.

CONCLUSIONS
We analyzed data from the 3.5-year NICER campaign of eRO-QPE1, with the following key findings: • The QPEs show complex, non-monotonic secular evolution, including variations by a factor of several in the characteristic luminosities, temperatures, blackbody radii, total energy output, and recurrence times.In the Oct 2023 epoch, the QPEs disappeared within NICER detectability (Fig. 1), but we confirmed they are ongoing with XMM-Newton in Jan 2024 (Fig. 2), at a level ∼10x fainter than recent NICER detections and ∼100x fainter (and ∼3x cooler) than initial discovery (Arcodia et al. 2021).The increasing irregularity of QPE recurrence times (Fig. 6) and R bb evolution (Fig. 3) coincides with the decreasing energy output (Fig. 5).The large dynamic range of QPEs even within a single source is an important constraint for theoretical models.
• The QPEs show L − kT hysteresis (Fig. 3) as also seen in Arcodia et al. (2022), which is physically consistent with a slowly-cooling emitting region which grows by a factor ∼ 2 − 3 over the course of a flare.We infer R bb ∼ R ⊙ , which is similar to the size for GSN 069 (Miniutti et al. (2023b) Fig. 18).The long-term increase, then decrease of peak R bb may point to a long-term super-period caused by e.g. by the apisdal/nodal precession cycles of an EMRI companion, but any such periodicity will take multiple years of observation to confirm.
• Radio campaigns with MeerKAT and VLA did not detect the source at any epochs (Tab.2.4.2,Fig. 1).Though the VLA campaign was partly simultaneous with a NICER observing epoch (Jan 2023), none of the radio observations aligned with a burst; thus we are unable to conclude whether the QPEs are definitely undetected at radio wavelengths.The lack of detected radio emission during quiescence rules out the presence of strong AGN activity, thus the presence of a long-lived accretion disk.
• There is a possible ∼6-day modulation in the timing residuals of the bursts (Fig. 7) during 6 out of 8 epochs.Interestingly, this period aligns with the inferred accretion disk nodal precession period.We may explain the significant scatter in QPE recurrence times, which is an unusual feature of eRO-QPE1, by a periodically changing relative inclination between the EMRI companion and the quiescent disk altering the relative delay/advance of each EMRI-disk interaction.This also provides a convenient explanation for the general lack of a long-short recurrence pattern (though it appears briefly in the first Dec 2022 epoch alongside a highlow amplitude, Fig. 7), and the atypically long quasi-period of eRO-QPE1.
• Associating the ∼6-day modulation with the nodal precession frequency of the disk (due to gravitomagnetic frame dragging) allows us to make a (degenerate) constraint on the SMBH spin and disk radial density profile (Fig. 8).The result is sensitive to the choice of M BH and disk R out , and only weakly dependent on M D .
The decreasing luminosity/altogether disappearance of QPEs (Fig. 1) has been noted in other sources already (Chakraborty et al. 2021;Miniutti et al. 2023b;Arcodia et al. 2024).Notably, in GSN 069 the QPEs re-appeared after one year (Miniutti et al. 2023a), revealing a quiescent luminosity threshold for the appearance of QPEs and a dependence of the burst amplitude/temperature on the quiescent disk properties (tending asymptotically towards kT QP E /kT disk = 1 as the quiescent lu-minosity increases).We cannot say definitively whether QPEs disappeared altogether in Oct 2023 (as NICER is not sensitive enough to detect the quiescence), and indeed the presence of one peak in the Nov 2023 observations suggests they may have been ongoing at an undetectable faintness all along.The detection of a faint QPE in Jan 2024 with XMM-Newton seems to agree with this, though the lack of a significant change in quiescence luminosity alongside the faint QPE (in contrast to GSN 069) now presents an open puzzle, complicating the dependence on the disk Ṁ obeyed by QPEs in general.Continued monitoring with XMM and NICER will reveal the evolution of eruptions in this eRO-QPE1, and whether more unexpected behaviors will continue to arise.
Example robust detection (13σ, rate=1.6 count sec -1 ) Example marginal detection (1.2σ, rate=0.17count sec -1 ) Figure A.1.NICER spectra for a sample marginal detection (1.2σ, top) and robust detection (13σ, bottom).We show spectra fit first with only the default SCORPEON components (representing the diffuse X-ray background and non X-ray background), then with the addition of ionized oxygen lines from the solar wind charge exchange (SWCX), then with the addition of a tbabs×zbbody representing the source contribution.

C. QPE SAMPLE-WIDE CORRELATIONS AND ENERGY DISTRIBUTIONS
Here we present plots correlating the flare features across the population of QPEs across all epochs.Not all 92 QPEs are used in each plot, as some bursts are too poorly-sampled to determine e.g.rise/decay time, and thus total energy output.Fig. C.4 shows that there is no apparent correlation of the QPE total energy outputs with time before/after the burst.Fig. C.5 shows there is a slight anticorrelation of peak luminosity with decay time, though we cannot constrain a similar relationship for the rise times with the NICER sampling.
Figure3.Top: Inferred blackbody radius R bb (filled points) plotted over fitted double-exponential model luminosity (gray), folded within each epoch.The short-term evolution is consistent with a cooling emission region expanding by a factor ∼ 2 − 3 over each flare, in agreement withMiniutti et al. (2023b).Over the long-term, the peak inferred R bb appears to grow (from Aug 2020 to Feb 2022), then shrinks (by Oct-Nov 2022).This may be a geometric (viewing-angle) effect attributed to the EMRI apsidal precession (∼ 10s of days), or a change in the mutual EMRI-disk inclination due to the EMRI nodal precession (∼ 100s of days) resulting in changing ejecta properties (see discussion in Section 4.2).Bottom: Same plot for kT .

Figure 4 .
Figure4.Peak L bol vs. kT for each QPE.In general the trend appears consistent with L bol ∝ T 4 , with the exception of the Feb 2022 and Oct-Dec 2022 epochs.The Oct-Dec 2022 epochs are also where R bb appears to invert its evolution (Fig.3).The ∝ T 4 lines are plotted for visualization only, and are not fits to the data.

Figure 5 .
Figure5.Secular evolution of QPE energy output (luminosity integrated over +/-3 e-folds from peak).The most recent epochs appear to show an overall decreasing trend, and QPEs are not seen in Oct 2023.We compute the Oct 2023 upper limit of the total energy by assuming the duration of possible QPEs stays roughly consistent with previous epochs.The peak of one QPE is seen in Nov 2023, but we cannot estimate the total energy output without better sampling of the overall flare.

Figure 6 .
Figure6.The QPE recurrence time trec shows significant long-term evolution, varying between 0.76 and 1.10 days across 3 years (and doing so non-monotonically, possibly hinting at a longer-term precession period of ∼ 10 − 100s of days).

Figure 7 .
Figure 7. Overplotted even/odd burst sequences (demonstrating the high scatter and lack of any long/short recurrence pattern), and the corresponding O-C timing residuals (i.e.whether each burst arrives early/late compared to the average trec within the epoch, Section 3.3.1).As noted in Arcodia et al. (2022), the scatter in recurrence time (∼ 50%) is significantly higher than other QPE sources.The O-C plots are overplotted with best-fit periods per-epoch (gray dashed line) and overall (blue dashed).The Aug 2020 epoch (P peak = 8.96 days) and first Dec 2022 epoch (P peak = 2.5 days) are significant outliers from the 5.73 average period.The long-short recurrence pattern seen in other QPEs is not generally seen in eRO-QPE1, with the possible exception of a short-lived phase in early Dec 2022.

Figure C. 4 .Figure C. 6 .
Figure C.4.Total energy output of each QPE (y-axis) vs. time delay from previous burst (x-axis, left) or to next burst (x-axis, right).No clear correlation is apparent.

and the model is tbabs×zbbody. Compared to the
previous XMM -detected QPEs, the peak temperatures and luminosities are now significantly lower.On the other hand, due to the low count rate (hence large errors) of quiescence, the quiescence L bol and R bb are within 3σ consistent with the previous XMM detection in July-August 2020.
Figure8.Joint constraints on SMBH spin a, and the powerlaw index p of the disk radial density profile Σ(R) ∝ R −p .We assume MBH = 10 5.8 M⊙(Arcodia et al. 2021), and the accretion disk has outer radius Rout = 350Rg.These constraints depend strongly on the assumed MBH and Rout.