X-ray quasi-periodic eruptions from two previously quiescent galaxies

Quasi-periodic eruptions (QPEs) are very-high-amplitude bursts of X-ray radiation recurring every few hours and originating near the central supermassive black holes of galactic nuclei1,2. It is currently unknown what triggers these events, how long they last and how they are connected to the physical properties of the inner accretion flows. Previously, only two such sources were known, found either serendipitously or in archival data1,2, with emission lines in their optical spectra classifying their nuclei as hosting an actively accreting supermassive black hole3,4. Here we report observations of QPEs in two further galaxies, obtained with a blind and systematic search of half of the X-ray sky. The optical spectra of these galaxies show no signature of black hole activity, indicating that a pre-existing accretion flow that is typical of active galactic nuclei is not required to trigger these events. Indeed, the periods, amplitudes and profiles of the QPEs reported here are inconsistent with current models that invoke radiation-pressure-driven instabilities in the accretion disk5–9. Instead, QPEs might be driven by an orbiting compact object. Furthermore, their observed properties require the mass of the secondary object to be much smaller than that of the main body10, and future X-ray observations may constrain possible changes in their period owing to orbital evolution. This model could make QPEs a viable candidate for the electromagnetic counterparts of so-called extreme-mass-ratio inspirals11–13, with considerable implications for multi-messenger astrophysics and cosmology14,15.

The second new eROSITA QPE, hereafter eRO-QPE2, showed similar variability patterns and Xray spectra as eRO-QPE1 during the X-ray sky survey (Fig. 2a and 2b).We associated it with the galaxy 2MASX J02344872-4419325 and determined a spectroscopic redshift of z=0.0175 (see Methods, 'The host galaxies of QPEs').The related intrinsic 0.5-2 keV luminosity of the quiescent (1σ upper limit) and peak phase is <4.0x10 40 erg s -1 and ~1.0x10 42 erg s -1 , respectively.A follow-up observation with XMM-Newton revealed 9 eruptions in a single day, oscillating between ~1.2x10 41 erg s -1 and ~1.2x10 42 erg s -1 in the 0.5-2 keV band (Fig. 2c).In neither eRO-QPE1 nor eRO-QPE2 is there evidence of simultaneous optical/UV variability (see Fig. 1c and 2c), in agreement with the behaviour of GSN 069 1 .eRO-QPE1 shows a range of QPE rise-to-decay duration with a mean (dispersion) of ~7.6 hours (~1.0 hours) and peak-to-peak separations of ~18.5 hours (~2.7 hours), as derived from the NICER light curve (Fig. 1d).The duty-cycle (mean duration over mean separation) is ~41%.Conversely, eRO-QPE2 shows much narrower and more frequent eruptions (see Fig. 2c): the mean (dispersion) of the rise-to-decay duration is ~27 minutes (~3 minutes), with a peak-to-peak separation of ~2.4 hours (~5 minutes) and a duty-cycle of ~19%.Interestingly, compared to the two previously known QPEs 1,2 , eRO-QPE1 and eRO-QPE2 extend the parameter space of QPE widths and recurrence times towards longer and shorter timescales, respectively.We also note that eRO-QPE1 is the most luminous and the most distant QPE discovered to date, and the most extreme in terms of timescales.The outbursts duration and recurrence times in eRO-QPE1 are approximately an order of magnitude longer than in eRO-QPE2.This could simply be an effect of the timescales scaling with black hole mass 17 .We estimated the total stellar mass of the two host galaxies, which is 4-8 times higher in eRO-QPE1 compared to eRO-QPE2.Assuming a standard scaling of the black hole mass with stellar mass (see Methods, 'The host galaxies of QPEs'), this is broadly in agreement with their different X-ray timing properties.Furthermore, peak soft X-ray luminosities of ~ 2x10 43 erg s -1 and ~10 42 erg s -1 , for eRO-QPE1 and eRO-QPE2 respectively, exclude a stellar-mass black hole origin and their X-ray positions, within uncertainties, suggest a nuclear origin (Extended Data Fig. 1a and  2a).
The optical counterparts of eRO-QPE1 and eRO-QPE2 are local low-mass galaxies with no canonical AGN-like broad emission lines in the optical nor any infrared photometric excess indicating the presence of hot dust (the so-called torus 18 ).In this sense they are similar to GSN 069 3 and RX J1301.9+2747 4 , whose optical spectra, however, show narrow emission lines with clear AGN-driven ionisation 3,4 .Instead, the optical counterpart of eRO-QPE1 is easily classified as a passive galaxy from the absence of any significant emission line (Extended Data Fig. 1b) and in eRO-QPE2 the strong narrow emission lines observed classify it as a star forming galaxy (Extended Data Fig. 2b and Methods, 'The host galaxies of QPEs').This in turns suggests that the two newly discovered galaxies have not been active for at least the last ≈10 3 -10 4 years, assuming narrow-line region light-travel timescales 19 .While the number of known QPEs is too low to reach firm statistical conclusions, our blind search is inherently designed to sample the QPEs' host galaxies population without bias, as opposed to serendipitous or archival discoveries which relied on the source being previously active and known 1,2 .These results hint that perhaps the parent population of QPE hosts consists more of passive, than active galaxies.The QPEs X-ray spectra in quiescence are consistent with an almost featureless accretion disk model 1,2 (see Methods, 'X-ray spectral analysis'), although the inactive nature of the host galaxies of our sources argues against a preexisting AGN-like accretion system.
A few scenarios to explain the QPEs have been suggested 1,10 , some based on the presumed active nature of the QPEs host black holes.These include so-called limit-cycle radiation-pressure accretion instabilities (see Methods, 'On accretion flow instabilities'), proposed for GSN 069 1 based on the similarities between its observed properties and two extremely variable stellar-mass black holes, namely GRS 1915+105 20,21 and IGR J17091-3624 22 .However, the observed properties of our newly discovered QPEs, as well as those of RX J1301.9+2747 2 , are inconsistent with the theoretical predictions of this scenario [5][6][7][8][9] .In particular, the faster rise and slower decay of eRO-QPE1 would imply a thicker flow in the cold and stable phase than in the hot and unstable phase, contrary to theory 6 .Moreover, the theory predicts that once the period, the duty cycle and the luminosity amplitude are known, only specific values of black hole mass M BH and viscosity parameter α are allowed 8 : for eRO-QPE1 (eRO-QPE2) one solution is found for M BH ~ 4x10 6 M ⊙ and α ~ 5 (M BH ~ 3x10 6 M ⊙ and α ~ 3), therefore for the expected masses 1,2 an unphysically high viscosity parameter would be required; alternatively, more reasonable values of α ~ 0.1 and ~ 0.01 would yield very small M BH ~ 2.4x10 3 M ⊙ and M BH ~ 60 M ⊙ (M BH ~ 4.3x10 3 M ⊙ and M BH ~ 30 M ⊙ ) for eRO-QPE1 (eRO-QPE2).Even in this latter scenario and pushing α as high as ~ 0.2, the resulting thermal timescales for eRO-QPE1 (eRO-QPE2) are τ th ~ 20s (35s) at 20 r g , which is orders of magnitude smaller than the observed QPE timescales (more details in Methods, 'On accretion flow instabilities').
Extreme or sinusoidal quasi-periodic variability as seen in QPEs is also typically associated with compact objects binaries, a scenario which would not require the galactic nuclei to be previously active, as our new evidence suggests.Drawing a simplistic scenario, we assumed the mass of the main body to be in the range ~10 4 -10 7 M ⊙ for both eRO-QPE1 and eRO-QPE2 and computed the expected period decrease of a compact binary due to emission of gravitational waves.We inferred that an orbiting body with a similar mass, namely a supermassive black-hole binary with a mass ratio of order unity 23 , is unlikely given the properties of the observed optical, UV, infrared and Xray emission in QPEs and the lack of evident periodicity and/or strong period decrease therein.If QPEs are triggered by the presence of a secondary orbiting body, our data suggest its mass (M 2 ) to be much smaller than the main body.This is in agreement with at least one proposed scenario for the origin of GSN069, for which the average luminosity in a QPE cycle can be reproduced by a periodic mass-inflow rate from a white dwarf orbiting the black hole with a highly eccentric orbit 10 .Our current data for eRO-QPE1 only exclude M 2 larger than ~ 10 6 M ⊙ (~10 4 M ⊙ ) for zero (high, 0.9) eccentricity (as a function of the mass of the main body, Extended Data Fig. 7a); instead, for eRO-QPE2 we can already argue that only an orbiting M 2 lower than ~10 4 M ⊙ (~10 M ⊙ ) is allowed for zero (~0.9) eccentricity (Extended Data Fig. 7b).More details are reported in Methods ('On the presence of an orbiting body').
Future X-ray observations on longer temporal baselines (months or years) will help to constrain or rule out this scenario and to monitor the possible orbital evolution of the system.This picture is also reminiscent of a suggested formation channel of extreme mass ratio inspirals 24,25 and it could make QPEs their electromagnetic messenger 13,26 .Regardless of their origin, the QPEs seen so far seem to be found in relatively low-mass super-massive black-holes (~10 5 -10 7 M ⊙ ) and finding more will help us to understand how black holes are activated in low-mass galaxies, a poorly explored mass range so far in their co-evolution history 27,28 , which is however crucial for synergies with future LISA gravitational waves signals 29 .

Methods
Blind search for QPEs with eROSITA eROSITA 16 is the main instrument aboard the Spectrum-Roentgen-Gamma (SRG) mission (R. Sunyaev et al., in preparation), which was launched on 13 July 2019.On 13 December 2019 it started the first of eight all-sky surveys (eRASS1-8), each completed in six months, observing in the 0.2-8 keV band.In each survey, as the field of view moves every point of the sky is observed for ~40 seconds every ~4 hours with the number of times (typically six) varying with the location in the sky, increasing towards high ecliptic latitudes.Our search for QPE candidates starts with a systematic screening of all eROSITA light curves, produced for each detected source on a weekly basis by the eROSITA Science Analysis Software (eSASS; H. Brunner et al., in preparation).Light curves are binned to yield one data point for each 4-hour revolution (called an 'eROday').A light curve generated by the eSASS pipeline will trigger a 'QPE alert' if it shows two or more high-count states with (at least) one low-count state in between (see Fig. 1a and 2a as examples) in any of its standard energy bands (0.2-0.6, 0.6-2.3,2.3-5.0 keV).As thresholds, we fixed a relative count rate ratio (including uncertainties) of 5, if both high and low states are detected, or 3, if the low-count state is consistent with the background.Since neither the survey scans nor QPEs are strictly periodic, every eRASS can be treated as an independent sky to find new candidates.This produces a census of X-ray sources varying on hours timescales for each eRASS, albeit only for the specific intermittent pattern described above.Unsurprisingly, the vast majority of the automatically generated alerts are produced by Galactic sources (mainly flaring coronally active stars), but we can filter them out by finding the multi-wavelength counterpart associated to every X-ray source (M.Sa.et al., in preparation).Good QPE candidates are then selected screening the handful of alerts with a secure or possible extra-galactic counterpart.Thanks to this process, we identified the two best eROSITA QPE candidates which were worth immediate follow-up, promptly obtained with both XMM-Newton and, in one case, NICER.Given the success of our initial search over the first nine months of the survey, we are confident that we can detect up to ~3-4 good eROSITA QPE candidates every year.Therefore, by the end of the last eROSITA all-sky survey in December 2023 this search may provide a sample of up to ~10-15 QPEs.
Previous X-ray activity.eRO-QPE1 has not previously been detected in X-rays, although upper limits can be obtained from the XMM-Newton upper limits server for ROSAT 49 , both from the survey and a pointed observation (taken in 1991 and 1992, with ~270 and ~5300 seconds, respectively), and the XMM-Newton Slew Survey 50 (taken in 2004, 2007, 2008 and 2017, all between ~3-8 seconds of exposure).The ROSAT pointed observation puts a stringent upper limit at ≤3.8x10 -14 cgs in the 0.2-2.0keV band.However, given the very short exposures compared with the timescales of eRO-QPE1, we can not rule out that QPEs were already ongoing and that all previous missions caught eRO-QPE1in a faint state.Like eRO-QPE1, eRO-QPE2 has not been previously detected in X-rays.Upper limits were again computed for ROSAT (taken in 1990, ~480 seconds of exposure) and the XMM-Newton Slew survey (taken in 2004, 2008, 2012 and 2013, all between ~4-8 seconds).The most stringent upper limit, at ≤8.8x10 -14 cgs in the 0.2-2.0keV band, comes from ROSAT.It is slightly below the flux observed by XMM-Newton in quiescence in the same band (Extended Data Fig. 4), perhaps indicating that the QPE behaviour only started more recently.For both QPE sources however, the ROSAT and Slew exposures are too shorter than the evolving timescales (the QPE quasi-period and its dispersion), hence they do not provide meaningful constraints on the start of the QPE behavior.

Data reduction
In this Section we report details of the processing of the complete data-set.We show a summary of the observations in Extended Data Table 1.
eROSITA.Members of the German eROSITA consortium (eROSITA-DE) have full and immediate access to survey data at Galactic longitudes 180 • <l <360 • .These data were processed using eSASS v946 (H.Brunner et al. in preparation).For eRO-QPE1 (eRO-QPE2), photons were extracted choosing a circular aperture of radius 80'' (67''), while background counts were extracted from an annulus (off-centered circle) of inner and outer radii 178'' (382'') and 996'', respectively, excluding all the other sources detected within the area.eRO-QPE1 was detected with a detection likelihood of 440 and a total number of 119 counts in the 0.2-5.0keV band.eRO-QPE2 was detected with a detection likelihood of 125 and a total number of 48 counts in the 0.2-5.0keV band.XMM-Newton.XMM-Newton data from EPIC MOS1-2 33 and EPIC-PN 34 cameras and the Optical Monitor (OM 35 ) were processed using standard tools (SAS v. 18.0.0and HEAsoft v. 6.25) and procedures.Event files from EPIC cameras were filtered for flaring particle background.Source (background) regions were extracted within a circle of 38'' and 34'' in eRO-QPE1 and eRO-QPE2, respectively, centered on the source (in a source-free region).eRO-QPE1 was consecutively observed three times with the U filter, then seven times with UVW1 and nine (eight) times with the UVM2 in the first (second) XMM-Newton observation, each exposure ~4400 s long.The source was detected only in the U and UVW1 with mean magnitudes ~19.9 and ~20.3 in both XMM-Newton observations (OM light curves in Fig 1c).eRO-QPE2 was consecutively observed twice with the U filter, then ten times with UVW1, six with UVM2 and three with UVW2 with all exposures beding 4400s.It was almost always detected in all filters with mean magnitudes of ~17.4,~17.5, ~18.0, and ~18.1, for U, UVW1, UVM2 and UVW2 filter, respectively (OM light curves in Fig 2c).eRO-QPE2 was flagged as extended in the U, UVW1 and UVM2 filters, therefore the reported absolute magnitudes include at least some contamination from the host galaxy.
NICER.NICER's X-ray Timing Instrument (XTI 36,37 ) onboard the ISS observed eRO-QPE1 between 17 August 2020 and 31 August 2020.Beginning late on 19 August, high-cadence observations were performed during almost every ISS orbit, which is roughly 93 minutes.All the data were processed using the standard NICER Data Analysis Software (NICERDAS) task 'nicerl2'.Good time intervals (GTIs) were chosen with standard defaults, yielding ~186 ks of exposure time.We further divided the GTIs into intervals of 128 s, and on this basis we extracted the spectra and applied the '3C50' model (R.R. et al., submitted) to determine the background spectra.The light curve for eRO-QPE1 in soft X-rays (Fig. 1d) was determined by integrating the background-subtracted spectrum for each 128-s GTI over the range 0.3-1.0keV.More detailed spectral analyses of these data will be discussed in a follow-up paper.
SALT.Optical spectra of eRO-QPE1 and eRO-QPE2 were obtained using the Robert Stobie Spectrograph (RSS 38 ) on the Southern African Large Telescope (SALT 32 ) on the nights of 2020 September 24 and 8, respectively.The PG900 VPH grating was used to obtain pairs of exposures (900 s and 500 s, respectively) at different grating angles, allowing for a total wavelength coverage of 3500-7400Å.The spectra were reduced using the PySALT package, a PyRAF-based software package for SALT data reductions 39 , which includes gain and amplifier cross-talk corrections, bias subtraction, amplifier mosaicing, and cosmetic corrections.The individual spectra were then extracted using standard IRAF procedures, wavelength calibration (with a calibration lamp exposure taken immediately after the science spectra), background subtraction and extraction of 1D spectra.We could only obtain relative flux calibrations, from observing spectrophotometric standards in twilight, due to the SALT design, which has a time-varying, asymmetric and underfilled entrance pupil 40 .

X-ray spectral analysis
In this work, X-ray spectral analysis was performed using v3.4.2 of the Bayesian X-ray Analysis software (BXA 41 ), which connects a nested sampling algorithm (UltraNest 42 ; J.B. in preparation) with a fitting environment.For the latter, we used XSPEC v12.10.1 43 with its Python oriented interface pyXSPEC.eROSITA source plus background spectra were fit including a model component for the background, which was determined via a principal component analysis (PCA) from a large sample of eROSITA background spectra 44 (J.B. et al., in preparation).XMM-Newton EPIC-PN spectra were instead fit using wstat, namely XSPEC implementation of the Cash statistic 45 , given the good counts statistics in both source and background spectra.We quote, unless otherwise states, median values with the related 16th and 84th percentiles and upper limits at 1σ. Results are also reported in Extended Data Table 2 and 3. eRO-QPE1.For eRO-QPE1 both eROSITA and XMM-Newton EPIC-PN spectra were fit with a simple absorbed black-body (using the models tbabs 46 and zbbody) or accretion disk (tbabs and diskbb 47 ), with absorption frozen at the Galactic column density of N H ~2.23x10 20 cm -2 , as reported by the HI4PI Collaboration 48 .For eROSITA, we jointly extracted and analysed spectra of the faint states (red points in Fig1a) and, separately, of the two bright states observed in eRASS1 (orange points in Fig1a).In the eROSITA bright states the temperature, in terms of k B T in eV, is 138 131 146 eV and 180 168 195 eV, using zbbody and diskbb as source model, respectively.The related unabsorbed restframe 0.5-2.0keV fluxes are 1.6 1.4 1.8 x10 -12 cgs and 1.5 1.4 1.7 x10 -12 cgs, respectively.The eROSITA spectrum of the faint states combined is consistent with background, with the temperature and unabsorbed rest-frame 0.5-2.0keV flux constrained to be ≤124 eV (≤160 eV) and ≤3.5x10 -14 cgs (≤3.4x10 -14 cgs) for zbbody (diskbb).We also analysed the observations of eRO-QPE1 obtained six months later during eRASS2, which triggered our QPE search again: two bright states were observed separated by several faint ones, with fluxes consistent with eRASS1.
We performed time-resolved X-ray spectral analysis on XMM-Newton data, extracting a spectrum in each 500s time bin of the EPIC-PN light curve, with the exception of the quiescence spectrum, which was extracted and analysed combining all the related time bins of both observations (i.e.before t~26500s in eRO-QPE1-XMM1 and before t~35788s in eRO-QPE1-XMM2, with times as defined in Fig1c).Fit results obtained using XMM-Newton EPIC-PN spectra with diskbb as the source model component are shown in Extended Data Fig. 3. Furthermore, we show for visualization three EPIC-PN spectra and related best-fit models (Extended Data Fig. 5a) corresponding to the quiescence phase and the peak of both XMM-Newton observations.A more thorough X-ray spectral analysis with other models and additional components for the bright phase will be presented in future work.eRO-QPE2.For eRO-QPE2 eROSITA's faint and bright phases were also separately combined and analysed (Fig 2a,b).The faint state as observed by eROSITA is consistent with background.The temperature and normalization of the source can not be constrained, thus we only quote an upper limit for the unabsorbed rest-frame 0.5-2.0keV flux of ≤1.9x10 -14 cgs (≤5.7x10 -14 cgs) using zbbody (diskbb).The spectrum of the eROSITA bright states constrains the temperature to 164 149 182 eV and at 209 185 241 eV, using zbbody and diskbb as source model, respectively.The related unabsorbed restframe 0.5-2.0keV fluxes are 1.4 1.2 1.8 x10 -12 cgs and 1.5 1.2 1.8 x10 -12 cgs, respectively.The triggering eROSITA observation was obtained during eRASS2, although a single bright state (thus not satisfying our trigger criterion) was also detected in eRASS1 with the same flux level.For eRO-QPE2, in addition to the Galactic column density (N H ~1.66x10 20 cm -2 [48]) we included an absorption component at the redshift of the host galaxy (i.e. with the models tbabs, ztbabs, and zbbody or diskbb).This excess absorption was inferred to be present based on the XMM-Newton spectrum (see below).
For XMM-Newton, we performed time-resolved X-ray spectral analysis for each 150 s time bin of the EPIC-PN light curve.The absorption in addition to the Galactic value was first fitted in the XMM-Newton quiescence spectrum, which was extracted combining all the low states in the XMM-Newton light curve (Fig. 2c and Extended Data Fig. 4).The fit yielded a N H =0.35 0.30 0.40 x10 22 cm -2 .In all other observations, including all eROSITA spectra and the rises, peaks and decays in the XMM-Newton light curve, the additional N H was left free to vary between the 10th and 90th percentile of the fitted posterior distribution of the quiescent spectrum.Under the assumption that absorption did not vary throughout the observation, this ensures that no spurious effects is imprinted on the fit temperature and normalisations due to degeneracies with N H ; at the same time, in this way parameters are marginalised over a reasonable range in N H . Freezing the value instead would artificially narrow the uncertainties on the temperature and normalisations.Fit results obtained with diskbb as the source model are shown in Extended Data Fig. 4. Furthermore, we show for visualization the EPIC-PN spectra and best-fit models of the quiescence and peak phases (Extended Data Fig. 5b).Similar results are obtained using zbbody as the source model.

Timing analysis
In Fig. 3 we show the median (with related 16 th and 84 th percentile contours) light curve profiles obtained by folding the light curve at the eruptions peaks.First, a random representative burst is selected and cross-correlated with the whole light curve.The peaks of this cross-correlation identify the times when the phase is zero.Data are then folded at these phase zero times to obtain a template median profile, which is then used to repeat the same operation and yield Fig. 3.A phase bin of ~0.1 corresponds to ~6600s and ~820s for eRO-QPE1 and eRO-QPE2, respectively.Moreover, XMM-Newton and NICER light curve profiles were fit with UltraNest 42 .Motivated by the strong asymmetry in eRO-QPE1 (Fig 1c,d and Fig 3a), we adopted a model with Gaussian rise and an exponential decay, a generic model often adopted for transients 51 .eRO-QPE2, on the other hand, can be fit with a simple Gaussian profile (Fig. 3b), possibly due to the much shorter timescales.Here we simply highlight the most evident results of timing analysis, while a more in depth study of the variability properties of QPEs is deferred to future work.Here the modeling allows us to determine mean values for the QPEs duration and recurrence time, which were used for comparison with models of accretion instabilities (see 'On accretion flow instabilities') and compact object binaries (see 'On the presence of an orbiting body').The mean rise-to-decay duration for eRO-QPE1, as observed from the NICER light curve (Fig. 1d), is ~7.6 hours (dispersion of ~1.0 hours), while the mean peak-to-peak separation is ~18.5 hours (dispersion of ~2.7 hours).The related dutycycle (here computed simply as mean duration over mean separation) is ~41%.Conversely, eRO-QPE2 shows much narrower and more frequent eruptions (see Fig. 2c): the mean the rise-to-decay duration is ~27 minutes (dispersion of ~3 minutes), with a mean peak-to-peak separation of ~2.4 hours (dispersion of ~5 minutes) and a duty-cycle of ~19%.

The host galaxies of QPEs
Very little was known on both galaxies from published multi-wavelength catalogs, except for WISE infrared monitoring indicating a quite stable W1~W2 emission, typical of inactive galactic nuclei, for the last few years.Most of our knowledge is based on optical spectra taken with SALT after the X-ray discovery.The optical counterpart of eRO-QPE1 is classified as a passive galaxy from the absence of any significant emission line (Extended Data Fig. 1b 4000) is the ratio of the continuum level after and before the 4000Å break 53 .Using standard line diagnostics plots 54 we can confirm that indeed eRO-QPE2 can be classified as star forming.Spectroscopic classification of future QPEs will be crucial to confirm whether their host galaxies are indeed preferentially inactive, as our pilot study suggests, or not.A first census in a statistically significant sample may bring new insights as has been the case for other transients, such as Tidal Disruption Events (TDEs) [55][56][57][58] .
A preliminary analysis of the properties of QPEs' host galaxies was performed by fitting the optical spectra (Extended Data Fig. 1b and 2b) with Firefly 59,60 .We first re-normalized the flux of the optical spectra using the most recent g-band and r-band archival magnitudes, since SALT spectra are not calibrated to absolute values 40 .For eRO-QPE1, gri-band photometry (g=18.7±0.06,r=18.0±0.05,i=17.8±0.05mag) was taken on July 30 th 2020 with the Rapid Response Robotic Telescope at Fan Mountain Observatory, indicating that the source did not change significantly with respect to archival photometry 61 .The total stellar masses inferred with Firefly from the optical spectra are M * ~3.8 − 1.9 + 0.4 x10 9 M ⊙ and ~1.01 − 0.50 +0.01 x10 9 M ⊙ for eRO-QPE1 and eRO-QPE2, respectively.Systematic errors and degeneracy due to the use of different stellar population models 62 would push M * to higher values for eRO-QPE1 (~4.8x10 9 M ⊙ ) and lower values for eRO-QPE2 (~0.6x10 9 M ⊙ ), enhancing their relative difference.Firefly also yields an estimate of the age of the stellar population and the star formation rate (SFR), although for medium and low signal-to-noise ratios these estimates are more prone to biases 59 .For eRO-QPE2, the mean signal-to-noise ratio (~23) is high enough to yield a fairly reliable SFR~0.078− 0.066 +0.001 M ⊙ /yr, which is also consistent within uncertainties with the SFR that can be estimated from the [O II ] and Hα luminosities 63,64 .For eRO-QPE1 the mean signal-to-noise ratio (~8) is lower and no reliable estimate of the SFR was obtained.We therefore inferred an upper limit of ~0.01M ⊙ /yr from the absence of significant narrow emission lines 63,64 .We report in Extended Data Fig. 6 the M * -SFR plane with our two newly discovered QPEs, together with normal galaxies, and hosts of known TDEs 57 and CLAGN 65 , all taken below z<0.1 and within the Sloan Digital Sky Survey MPA-JHU DR7 catalog 66 .Evidence is mounting that both TDEs 57,67 and CLAGN 65 might be over represented in galaxies in the so-called 'green valley', perhaps indicating that they are activated in specific periods of galaxy co-evolution with their central black holes.For QPEs, a statistically meaningful sample still needs to be built before reaching any conclusion.
We have estimated that the host galaxy of eRO-QPE1 is more massive than that of eRO-QPE2.We here refrain from quoting absolute values for black hole masses using their scaling relations with the host galaxies properties, since our stellar masses are lower than the ones used to calibrate them 68 .However, it is worth mentioning that the relative ratio of ~4-8 in stellar masses, between eRO-QPE1 and eRO-QPE2, would propagate to a black hole mass ratio of the order of ≈10 [68].This is in line with the X-ray timing properties in eRO-QPE1 and eRO-QPE2, since their peak-topeak separation and rise-to-decay duration scale roughly by the same amount.Finally, X-ray emission from eRO-QPE1 and eRO-QPE2 was observed to be positionally consistent with the galaxy nucleus for both objects (Extended Data Fig. 1a and 2a; Methods, 'The two new eROSITA QPEs').If a future QPE is found in a more nearby galaxy we can aim to constrain more precisely the X-ray position with respect to the galactic nucleus.This will allow us to determine conclusively whether or not these phenomena are nuclear.

On accretion flow instabilities
Accretion disks 69 with an accretion rate such that radiation-pressure dominates in the inner flow are thought to be subject to thermal-viscous instabilities 70 .The net result of these instabilities is that the luminosity is predicted to oscillate [5][6][7][8] with timescales and amplitude proportional to the black-hole mass and bolometric luminosity 71,72 .The predicted light curves profiles show first a stable and slow rise in luminosity, as both temperature and surface density increase while matter is slowly accumulated.Thereafter a sharp luminosity burst is produced by a runaway increase (decrease) in temperature (surface density) propagating outwards within the unstable region.Finally, the inner flow, devoid of the matter supply, cools down rapidly and cycles back to the initial stable state with low temperature and density.Both heating and cooling fronts propagate following thermal timescales 6 , where τ th ~ α -1 (GM BH /R 3 ) -1/2 .These so-called limit-cycle or heartbeat instabilities have been successfully applied to a few accreting sources across all mass scales, for instance to the stellar-mass black holes GRS 1915+105 20,21,73 , IGR J17091-3624 22 and 4XMM J111816.0-324910 74nd to super-massive black holes in a statistical fashion 71,72 .The similarity of their timing properties with QPEs in GSN069 is tantalizing and naturally led to the proposed connection with limit-cycle instabilities for that object.In particular, the symmetry of the eruptions in GSN069 was compared to the fast heating and cooling phases of the instability 1 , both following similar τ th under the assumption of invariant α across the two phases 75 .The lack of a slow rise before the eruptions in QPEs, predicted by the instability models, could be due to our limited coverage of the full disk temperature profile in the soft X-ray band.
With the two newly discovered QPEs we can now argue against at least this type of accretion disk instability as the origin of QPEs.Specifically, the strong asymmetric nature of the eruptions in eRO-QPE1, which show a faster rise and a much slower decay (Fig. 3a), argues against this interpretation.Qualitatively, our data would suggest that QPEs are not related to τ th , since α is not expected to change between the hot and cold phases in AGN 75 .Moreover, if instead it is the front propagation timescales, following τ front ~ (H/R) -1 τ th , or viscous timescales, following τ visc ~ (H/R) -2 τ th , that regulates the rise (decay) in the cold (hot) phase, it would imply a thicker flow in the cold and stable phase than in the hot and unstable phase.This runs contrary to the theoretical expectation that unstable flows should be thicker 6 .The limit-cycle oscillation theory further predicts that once the period, duty cycle and luminosity amplitude are known and a viscosity prescription for the accretion flow is adopted, there are only specific values of M BH and α that unstable sources are allowed to follow 8 .Here we adopt for eRO-QPE1 (eRO-QPE2) a peak-to-peak period T pp =18.5 h (2.4 h), an amplitude A~294 (~11) and a duty-cycle D=41% (19%).The amplitude A is the ratio of the disk luminosity (computed within the 0.001-100 keV range) for peak and quiescence, taken as proxy of the maximum and minimum bolometric luminosity, while D is here defined as the ratio of the flare duration (rise-to-decay T rd ) and the period T pp .We begin by adopting a standard viscosity prescription, with the average stress between the annuli proportional to the gas plus radiation pressure P tot [69].The allowed M BH and α values for eRO-QPE1 (eRO-QPE2) are M BH ~ 4x10 6 M ⊙ and α ~ 5 (M BH ~ 3x10 6 M ⊙ and α ~ 3), therefore an unphysically high viscosity parameter would be required.Considering alternative viscosity prescriptions 5,8 , for eRO-QPE1 (eRO-QPE2) a more reasonable α ~ 0.1 or ~ 0.01 would correspond to allowed M BH ~ 2.4x10 3 M ⊙ or M BH ~ 60 M ⊙ (M BH ~ 4.3x10 3 M ⊙ or M BH ~ 30 M ⊙ ), respectively.The above combinations are either unphysical or very unlikely.Adopting α ~ 0.2 and alternative viscosity prescriptions eRO-QPE1 (eRO-QPE2) would yield an intermediate-mass black hole (IMBH) at M BH ~ 0.9x10 4 M ⊙ (M BH ~ 1.6x10 4 M ⊙ ) accreting at ~0.1 (~0.3)Eddington in quiescence and at ~30 (~3) Eddington at the peak.However, this IMBH scenario would not account for the opposite asymmetry shown by eRO-QPE1 compared to the theoretical predictions, nor would the resulting thermal timescales be self-consistent for either of the two: for eRO-QPE1 (eRO-QPE2) τ th ~ 20s (35s) at 20 r g adopting M BH ~ 0.9x10 4 M ⊙ (M BH ~ 1.6x10 4 M ⊙ ), which is orders of magnitude smaller than the observed QPE durations, and the rise-to-peak times would be only reconciled with τ th at ~ 780 r g (~ 250 r g ).Analogous results can be obtained using the observed properties of RX J1301.9+2747 2 , adopting T pp ~20 ks (or the second period T pp ~13 ks), D=6% (9%) and A~9.4, the latter obtained taking the ratio of the quoted quiescence and peak 0.3-2.0keV flux as proxy for a bolometric luminosity ratio: adopting 2 α ~ 0.15 the allowed black hole mass is ~ 2.2x10 4 M ⊙ (~ 1.5x10 4 M ⊙ ), much lower than the quoted 2,76 ~ 0.8-2.8x10 6M ⊙ .
When a given source is in a 'sweet-spot' regime in mass accretion rate , some more recent modified accretion disks viscosity prescriptions predict the presence of a narrow unstable zone placed within an inner inefficient advection-dominated flow and an outer standard geometrically-thin and stable flow 9 .This model would reduce the propagation timescales by a factor ~dR/R, where dR is the radial extent of the unstable zone at a distance R from the black hole, which may help reconcile the model with the dramatic and fast variability observed in CLAGN 77 .This would not, however, change the inconsistency with the asymmetric shape of the newly observed QPEs, nor was it successful in modeling all the observables in GSN 069 9 .In summary, our data for both our newlydetected QPEs are inconsistent with published models for radiation pressure instability [5][6][7][8][9] .The role of more complex, or exotic phenomenology 9 should be further explored.
We also note that a fast-rise exponential-decay profile, like the one in eRO-QPE1, can be naturally produced by ionization instability models which are used for some bursting stellar-mass accreting compact objects 78 .To our knowledge there is no evidence so far of such instabilities taking place in AGN 79 .In addition, the predicted timescales are many orders of magnitude longer than QPEs for both AGN [79][80][81] and IMBH masses 82 .Finally, we discuss disk warping and tearing induced by Lense-Thirring precession 83,84 , which has been recently qualitatively compared also to QPE sources 85 .In this work we presented new evidence of QPEs being observed in previously inactive galaxies, therefore the accretion flow in these systems should be young.Moreover, a key element of disk warping and tearing due to Lense-Thirring precession is that mass needs to flow in from large inclination with respect to the black hole spin.Both conditions are satisfied if the accretion flow is formed, for instance, by a fullystripped TDE.However, in this case the warped inner flow would be damped quite fast 86 , which would be in contrast with QPEs lasting at least months 1 (Figure 1 and 2) or even years 2 .A more quantitative comparison is beyond the reach of this work and of current disk warping and tearing simulations, but this is a promising scenario worth exploring in the future.

On the presence of an orbiting body
Periodic variability is also often associated to binary systems of compact objects 23 and the connection with the quasi-periodic nuclear emission observed in QPEs is tempting.We here assume the main body to be a super-massive black hole ranging between ~10 4 -10 7 M ⊙ and we first consider the presence of a second orbiting super-massive black hole with a similar mass.There are several reasons which, when combined, disfavor such a scenario.Firstly, simulations show that the accretion flow of such objects is composed by a circum-binary disk with two inner small minidisks [87][88][89] , which are thought to produce a quasi-sinusoidal modulated emission 90,91 .This signature can be detected in transient surveys 92,93 or in single sources 94 , with a well-known extreme case being OJ 287 95,96 .However, so far there is no evidence of such variability in optical and UV data 1 of QPEs (Fig. 1c and 2c), in particular in eRO-QPE1, which was covered in g-and r-band by the Zwicky Transient Facility DR3 97 until the end of 2019.Nor can this prediction be reconciled with the dramatic non-sinusoidal eruptions observed in X-rays, even in the case of binary self lensing 98 which can produce sharper bursts, albeit achromatic therefore in contrast with the energy dependence of QPEs 1,2 .Moreover, we do not observe peculiar single-or double-peaked emission lines [99][100][101] and this can not be reconciled by enhanced obscuration 102 , since infrared photometry in QPEs is not AGN-like (WISE observed W1~W2 for the past 6-7 years) and X-rays do not indicate the presence of strong absorption.Secondly, super-massive black hole binaries are expected to form mostly via galactic mergers 103,104 , but the host galaxies of the two newly discovered QPEs look unperturbed (Extended Data Fig. 1a and 2a).Perhaps most importantly, a binary of super-massive black holes observed with a periodicity of the order of hours, such as the four observed QPEs, would show a large period derivative, due to gravitational wave emission, and would be relatively close to merger.To have (at least) four such objects very close to merger within z~0.02-0.05 is very unlikely 105 and would imply that they are much more common in the local Universe than observations suggest 92,93 .
Under the simplified assumption that the orbital evolution is dominated by gravitational waves emission, Extended Data Fig. 7 shows the allowed parameter space in terms of Ṗ and M 2 for a range of M BH.1 ~10 4 -10 7 M ⊙ and zero or high orbit eccentricity (e O ~0.9), given the rest-frame period of both eRO-QPE1 and eRO-QPE2.We have additionally imposed M 2 ≤M BH,1 .For both sources we can draw a tentative line at the minimum period derivative that, if present, we would have measured already within the available observations: quite conservatively, we adopt a period decrease of one cycle over the 15 observed by NICER for eRO-QPE1 and the 9 observed by XMM-Newton for eRO-QPE2 (Fig 1d and 2c).Our constraint on Ṗ is not very stringent for eRO-QPE1 and only high M 2 and eccentricities are disfavored; instead, for eRO-QPE2 only an orbiting IMBH, or smaller, is allowed for zero eccentricity, while only a stellar-mass compact object is allowed for high eccentricity (e O ~0.9).Future observations of eRO-QPE1 and eRO-QPE2 in the next months may lead to tighter constraints on the mass and eccentricity of the putative orbiting body.
The preliminary conclusion of our analysis is that, if QPEs are driven by the presence of an orbiting body around a central black hole, it is more likely that this is a compact object with a mass significantly smaller than the ~10 4 -10 7 M ⊙ assumed for the main body.This scenario could make QPEs a viable candidate for the electromagnetic counterparts of the so-called extreme mass ratio inspirals (EMRI [11][12][13] ), with considerable implications for multi-messenger astrophysics and cosmology 14,15 .Interestingly, it has been recently suggested for GSN 069 that a stellar-mass compact object orbiting around a super-massive black hole could be the origin of QPEs: a white dwarf of ~0.2M ⊙ on a highly eccentric orbit (e O ~0.94) could reproduce the mass inflow rate needed to produce the observed X-ray luminosity averaged over a QPE cycle 10 .This is reminiscent of a suggested, albeit still observationally elusive, EMRI formation channel 13,[24][25][26] .For GSN 069, a possible explanation of the QPE-free X-ray bright and decaying phase could be given by an accretion flow expanding and intercepting the body at a later time 1 ; or if the orbiting body was originally a massive star and the stripped envelope produced the TDE-like behavior of the past decade 1 while the remaining core started interacting with the newly born or expanded accretion flow only at a later stage, which would also explain the relatively small mass required by the white dwarf calculations 10 .For the other QPEs which did not show evidence of a past X-ray bright and decaying phase, this scenario is not necessary and the interaction with a second stellar-mass (or more massive) compact object could qualitatively reproduce the periodic behavior (Extended Data Fig. 7b).Future X-ray observations of the known QPEs would help in further constraining the possible orbital evolution.It should be pointed out that these calculations so far only matched the average observed QPE luminosity with the mass inflow rate required to produce it 10 , but details on the exact physical mechanism that would produce these X-ray bursts are still elusive (see 'On accretion flow instabilities').
Finally, we address the lack of UV and optical variability in the scenario of an orbiting body.The X-ray plateau at minimum shows an almost featureless accretion disk thermal spectrum [1][2] (Extended Data Fig. 5), which could have been built up during the first orbiting cycles.This accretion flow should be unusually small due to the lack of a broad line region [3][4] (Extended Data Fig. 1b and 2b), which would respond in light-days and that, if present, should have been observed in the SALT spectra taken months after the X-ray QPEs.The lack of strong UV and optical variability might be then due to the fact that the accretion disk is not large enough to even emit strong enough UVoptical radiation to emerge above the galaxy emission, which we can assume to be most of the observed L~4.0x10 41 erg s -1 (L~4.3x10 41erg s -1 ) in the OM-UVW1 filter at 2910Å for eRO-QPE1 (eRO-QPE2).Using a simplified but physically motivated accretion disk model 106 for a spin zero black hole accreting at ~0.1 Eddington, we computed the distance at which the bulk of 2910Å radiation would be emitted, namely ~1100 and ~500 r g for a mass of 10 5 and 10 6 M ⊙ , respectively.This would shift to even larger radii for increasing accretion rate (e.g.~1850 and ~860 r g at ~0.5 Eddington), while even for high spinning sources (spin ~0.9) the peak of OM-UVW1 flux would still come from ~775 and ~360 r g .Furthermore, the predicted OM-UVW1 disk luminosity would be at least one or two orders of magnitude lower than the observed L~4.0x10 41 erg s -1 in the most luminous scenario.Therefore, even an UV-optical eruption 100 times brighter than the plateau would be barely detectable above the galaxy component.
Predicted numbers.Detailed self-consistency calculations on the predicted rate of such EMRI events, as compared to QPE rates, are required but are beyond the scope of this paper.Instead, we can provide here a rough model-independent estimate of the expected number of observed QPEs, regardless of their origin.Convolving the black hole mass function 107 between ~10 4.5 -10 6.5 M ⊙ up to z~0.03 with the eROSITA sensitivity yields N≈100, which is then reduced with some educated guesses on a number of unknowns: during what fraction of their X-ray bright phase such sources undergo QPE behavior (the biggest unknown; e.g.>20% for GSN 069); how many such sources are obscured and missed (≈2/3); how many times we detect ongoing QPEs given the eROSITA sampling (depends on the period and the burst duration; e.g.≈20% for GSN 069).This results in a (extremely uncertain) number of order unity per eRASS scan in the eROSITA-DE hemisphere, which is however in agreement with our pilot study of the first few months of eROSITA operations.Thus, the low observed numbers do not necessarily imply that these events are a rare phenomenon intrinsically and they can actually be a fairly common product of the black holes co-evolution with their host galaxies 28 .With a statistically meaningful sample of QPEs, inverting this calculation may allow us to constrain the black hole mass function in a poorly known mass regime 27 .Extended Data Fig. 6 -The properties of QPEs' host galaxies.Stellar mass M * and star formation rate (SFR) for eRO-QPE1 (blue) and eRO-QPE2 (red), with related 1σ uncertainties; for eRO-QPE1 SFR is largely unconstrained (see Methods, 'The host galaxies of QPEs').For a comparison, normal galaxies 66 , TDEs 57 and CLAGN 65 , all below z<0.1, are also shown.

Extended Data
Extended Data Fig. 7 -Constraints on a secondary orbiting body.a, Allowed parameter space in terms of period derivative and secondary mass M 2 for a range of primary mass M BH,1 ~10 4 -10 7 M ⊙ and zero (solid lines) or high orbit eccentricity (e O ~0.9, dotted lines), in which can reproduce the rest-frame period of eRO-QPE1.We have additionally imposed M 2 ≤M BH,1 .We have drawn an approximate threshold at the minimum period derivative that, if present, we would have measured already within the available observations, corresponding to a period decrease of one QPE cycle over the 15 observed by NICER (Fig. 1d).The excluded region is shaded in red.b, same as a but for eRO-QPE2 and adopting as tentative minimum Ṗ a period decrease of one cycle over the 9 observed with XMM-Newton (Fig. 2c).The median value and related 16 th and 84 th percentile values are reported for every quantity; for unconstrained values 1σ upper limits are quoted using the 84 th percentile value of the parameter posterior distribution and are denoted with ↓.

Extended
Reported results are obtained with the model tbabs x diskbb, with Galactic NH frozen at 2.23x10 20 cm -2 , as reported by the HI4PI Collaboration 48 .Fluxes and luminosities are unabsorbed and rest-frame.The two eROSITA states are shown in Fig. 1a, whilst the three XMM-Newton observations in the table correspond to the three spectra in Extended Data Fig. 5a.Fdisk is computed within 0.001 and 100 keV.Same as Extended Data Table 1, but for eRO-QPE2.Reported results are obtained with the model tbabs x ztbabs x diskbb, with Galactic NH frozen at 1.66 x10 20 cm -2 , as reported by the HI4PI Collaboration 48 ; absorption in excess was estimated from 'XMM quiescence' and was allowed to vary within its 10 th and 90 th percentiles for all the other observations.The two eROSITA states are shown in Fig. 2a and model parameters in the low state are unconstrained; the two XMM-Newton observations in the table correspond to the spectra in Extended Data Fig. 5b.

Fig. 1 -
Fig. 1 -The first eROSITA QPE.a, eROSITA light curve in the 0.2-0.6 keV and 0.6-2.3keV energy bands (circles and squares, respectively), with red and orange highlighting faint and bright observations, respectively.The start of the light curve is t eRO,0 is MJD~58864.843.b, eROSITA Xray spectra of the bright and faint states in orange and red as in a. c, background subtracted XMM-Newton X-ray light curves with 500 s bins for EPIC PN (dark gray), MOS1 (green) and MOS2 (red) in the energy band shown in the legend.The beginning of both observations was contaminated by flares in the background and excluded; the dark grey solid line and contours show the underlying ≤1 keV EPIC-PN light curve to give a zeroth-order extrapolation of the rate, excluding the presence of obvious soft X-ray eruptions.t XMM,0 corresponds to the start of the cleaned MOS2 exposure in the first observation, namely MJD~59057.805.XMM-Newton optical and UV fluxes are shown in the lower sub-panels (units of erg cm -2 s -1 ), with non-detections shown as upper limits.d, background

Fig. 2 -
Fig. 2 -The second eROSITA QPE.a,b same as Fig1a,b but for eRO-QPE2.The start of the eROSITA light curve is MJD~59023.191.c, same as the Fig1c but for the XMM-Newton observation of eRO-QPE2.t XMM,1 corresponds to the start of the cleaned MOS1 exposure, namely MJD~59067.846.The mean (and related dispersion) of the rise-to-decay duration is ~27 minutes (~3 minutes), with a peak-to-peak separation of ~2.4 hours (~5 minutes).In all panels 1σ uncertainties are shown, as error bars or shaded regions.

Fig. 1 -
eRO-QPE1 position and identification.a, Legacy DR8 image cut-out around the optical counterpart of eRO-QPE1.Red and green circles represent the astrometrycorrected eROSITA and XMM-PN positions, respectively, with 1σ positional uncertainties.The EPIC-PN position was corrected excluding the target (blue cross) to ensure an unbiased estimate of the possible positional offset.b, SALT spectra of eRO-QPE1 shown in black and blue with related 1σ errors as shaded regions.The cyan spectrum represents a re-normalized sky spectrum to guide the eye for the residual sky feature around 5577Å.Extended Data Fig. 2 -eRO-QPE2 position and identification.Same as in Extended Data Fig.1, but for eRO-QPE2.Extended Data Fig. 3 -eRO-QPE1 spectral fit results.XMM-Newton PN light-curve (top panel) and time-resolved spectroscopy fit results for spectra extracted in the 500 s time bins (bottom panels) of the two XMM-Newton observations of eRO-QPE1 using an accretion disk model (diskbb): in particular, the evolution of the peak accretion disk temperature and the normalization, which is proportional to the inner radius once distance and inclination are known.The quiescence level is fit combining the first part of both XMM-Newton observations.It is shown with a dashed line because due to low counts the fit is more uncertain (see Extended Data Fig. 5a).Median fit values and fluxes of the high and low eROSITA states are reported with orange and red arrows pointing left (upper limits are denoted with diagonal arrows).1σ uncertainties on the fit results are shown with shaded regions around the median.Extended Data Fig. 4 -eRO-QPE2 spectral fit results.Same as Extended Data Fig. 3, but for eRO-QPE2.Extended Data Fig. 5 -eRO-QPE1 and eRO-QPE2 spectra.a, XMM-Newton EPIC-PN source plus background spectra for eRO-QPE1.Red, orange and green data correspond to quiescence and to the peak of the second and first XMM-Newton observation, respectively, with error bars showing 1σ uncertainties.The related solid lines show the unabsorbed source model obtained with diskbb, just for visualization.The grey line represents the background spectrum alone.The plateau is shown with a dotted line because due to low counts the fit is more uncertain.b, same as a but for eRO-QPE2.Here green data represent one of the peaks and the additional dashed lines indicate the absorbed source model.