© The Authors 2024

1. Introduction

The interest in multiwavelength extragalactic nuclear transients has recently increased with the advent of sensitive wide-area surveys, which started discovering these transients in significant numbers. Particular attention is being dedicated to transients that are detectable in a multi-messenger fashion, via electromagnetic radiation and gravitational waves (GWs) and/or neutrino detectors. The main protagonists are tidal disruption events (TDEs) with the recent addition of quasi-periodic eruptions (QPEs). For TDEs, high-cadence optical surveys, such as the Zwicky Transient Facility, took a major step forward in detection numbers (van Velzen et al. 2021; Yao et al. 2023) after the first discoveries by ROSAT (Bade et al. 1996; Grupe et al. 1999), although X-ray discoveries are catching up (Sazonov et al. 2021; Grotova et al., in prep.) thanks to the extended ROentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021) on board the Spectrum-Roentgen-Gamma observatory (SRG; Sunyaev et al. 2021). For QPEs, the field is still in its early days, with the first discovery reported in 2019 (Miniutti et al. 2019, hereafter M19) and only a handful of secure QPE sources known to date (Giustini et al. 2020; Arcodia et al. 2021, 2024; hereafter G20; A21; A24), plus some candidates (Chakraborty et al. 2021; Quintin et al. 2023; Evans et al. 2023; Guolo et al. 2024). QPEs were discovered for the most part (i.e., four out of six) by the eROSITA X-ray telescope (A21; A24). As the search method was purposely designed for the discovery of QPEs (A21; A24) and the eROSITA surveys are relatively homogeneous over the entire sky (Merloni et al. 2024), we are able to provide in this work, for the first time, an intrinsic volumetric rate of QPEs.

Finally, most of the latest theoretical efforts interpret QPEs as triggered when the radiatively efficient accretion disk of a massive black hole is repeatedly pierced by an orbiter of much lower mass (e.g., Xian et al. 2021; Franchini et al. 2023; Linial & Metzger 2023; Tagawa & Haiman 2023; Zhou et al. 2024)¹. Intriguingly, this interpretation would make QPEs the electromagnetic counterparts of extreme mass ratio inspirals (EMRIs), which are detectable by future-generation GW detectors (e.g., Amaro-Seoane 2018). Under the assumption that QPEs are indeed EMRIs, a remarkable implication is that the QPE volumetric rate would be the first ever observationally driven constraint to the EMRI rate.

2. The search for QPEs with eROSITA

On 13 December 2019, eROSITA started observing in survey mode and completed four of the eight foreseen all-sky surveys (eRASS1 to eRASS8), and started the fifth. In each survey, lasting six months, every source in the sky was observed for ∼40 s (hereafter called a visit) every ∼4 h (i.e. a eROday) for a total number of times within a single eRASS depending on its location in the sky. This number is typically around six on the ecliptic plane and increases toward higher ecliptic latitudes (Merloni et al. 2024). Most of the area in the sky has thus around ≈6 − 10 data points separated by 4 h, for a total baseline of ∼1 − 2 d in the single eRASS (Merloni et al. 2024) and about ∼2 − 2.5 yr for the eRASS1–eRASS4 (or eRASS1–eRASS5) baseline. Leading to the first QPE discoveries with eROSITA (Arcodia et al. 2021, 2024), our team developed an algorithm to look for significant and repeated high-amplitude variability in the eROSITA sources. QPE candidates can be found within eROdays of the single eRASS or combining multiple eRASS. A first version of this method was briefly described in A21. In brief, the algorithm searching for variable sources is primarily based on the X-ray light curves with no pre-screening based on the counterpart’s nature. Then, more steps are taken in post-processing to classify each given source, such as the exclusion of secure Galactic objects using Gaia DR3 (Gaia Collaboration 2016, 2021) and some further visual inspection of their X-ray and multiwavelength properties for the handful of interesting cases per eRASS. Given eROSITA’s cadence, X-ray follow-up with other X-ray telescopes is required for secure confirmation (A21; A24).

Here, we present the QPE search method in more detail with the goal of providing the first estimate of the volumetric rate of QPEs. Light curves are extracted with the srctool task of the eROSITA Science Analysis Software System (eSASS; Brunner et al. 2022) from event files processing version 020, for all sources detected above detection likelihood DET_LIKE_0 = 20 (Merloni et al. 2024). We exclude an area with a radius of three degrees around the South Ecliptic Pole, as it requires separate data processing and analysis (e.g., Bogensberger et al. 2024). Light curves are provided with ten-second bins between 0.2 − 0.6 keV, 0.6 − 2.3 keV and 2.3 − 5.0 keV. Light curves are thus rebinned to have a single count rate (CR) estimate per eROday (A21; A24). This step is carried out by summing, within the eROday visit, the source and background counts in the unbinned light curve produced by eSASS. Then, we numerically compute the Poisson probability mass function (PPMF, with the scipy Python package Virtanen et al. 2020) for the count rate given the detected source and background counts, normalized by the areas of the extraction regions. We use the median value of the PPMF and 1σ percentiles as asymmetric upper (subscript +) and lower (subscript −) error bars. We compute the PPMFs for background alone and the total rate from source plus background. In the light curve extracted from the source area, visits compatible within the uncertainties with background alone are considered undetected.

An X-ray source is added to a pool of interesting alerts if it shows two (or more) visits that are significantly brighter than one (or more) fainter visit(s) in between. Our significance criteria are (i) a ratio > 2 of the lowest 1σ percentile of the PPMF of both bright eROday visits (dubbed HI; CRHI − CRHI₋) to the upper 1σ percentile of the fainter one in between (dubbed LO; CRLO + CRLO₊), and (ii) maximum amplitude significance > 2, computed as AMP/ERR where AMP = CRHI − CRHI₋ − (CRLO + CRLO₊) (e.g., Boller et al. 2022) and . On the order of ≈550 − 650 sources per eRASS satisfy these criteria, which are purposely loose as we aim to be complete, rather than pure, in finding QPEs from extreme variability from galaxies.

For each eRASS, this preliminary sample is cross-matched with proper motion and parallax estimates to exclude Galactic objects. We match the X-ray coordinates with Gaia DR3 (Gaia Collaboration 2016, 2021) using the Python package astroquery (Ginsburg et al. 2019) and store information of all Gaia sources within 10″. Gaia proper motion (GPM) significance (S_GPM) is computed as S_GPM = GPM/σ_GPM, where GPM is computed from the RA and Dec proper motions as and its uncertainty . Alerts associated with a single Gaia source with proper motion significance above 5σ are excluded. The remaining alerts are typically ≈140 − 150 per eRASS. Their X-ray and archival multiwavelength properties are visually checked to exclude possible spurious Galactic contaminants that may remain after the first automated screening², and others of extragalactic nature (e.g. known quasars and blazars). In addition, we also check alerts with no proper motion or parallax estimate in Gaia (≈70 − 100 per eRASS). The vast majority of these consist of stars masked out from Gaia, which can be easily excluded visualizing their optical images, but some local faint galaxies were also found to be missing a Gaia estimate. After this pre-screening, the remaining alerts represent a subsample of likely extragalactic X-ray sources. As mentioned above, our variability significance criteria are purposely loose, to favor detection completeness over purity. Hence, the algorithm also triggers on stochastic variability of bright sources, which make the significance cut due to their smaller uncertainties. Apart from the extreme light curves like those of the eROSITA discovered QPEs (A21; A24), the lower priority tiers most likely contain QPE sources with lower amplitude variability within eRASS, which cannot be easily distinguished from quasar-like stochastic variability around the eRASS detection limit. We note that if a strikingly interesting light curve originates from a known active galactic nucleus, this is not an exclusion criterion for QPE searches per se. However, typically these alerts are associated with fairly well-studied low-redshift sources, and thus a closer look at their X-ray products from other (previous or following) eRASS and archival X-ray and multiwavelength data is sufficient to exclude them as top-tier QPE candidates. A trivial case in which other eRASS can be used to assess the candidates is that of a given eRASS light curve showing seemingly high-amplitude variability (potentially indicating a QPE source), while other eRASS show a significant detection in most eROdays, but with insignificant variability. Cases where the source is only detected in a single eRASS are instead kept in the pool of potential QPE candidates. Furthermore, a trivial case where other X-ray or multiwavelength archival data can be used to assess the candidates is that of well-studied local Seyfert galaxies with many archival observations. We note that this step of visual inspection naturally imprints some degree of subjectivity in the choice of candidates. However, our best efforts were put into quantifying this as much as possible in our rate calculations (Sect. 3.2) and supporting simulations (Appendix A). After this screening, only a handful of highly significant variable sources remained. Their archival multiwavelength properties and photometry (e.g. via VizieR; Ochsenbein et al. 2000) were explored for further analysis and follow-up. For instance, the variability criterion described here (namely, that of significant high-amplitude variability repeated within the single eRASS) is the one that found eRO-QPE1 and eRO-QPE2 (A21).

In addition, we considered the case where significant high-amplitude variability is repeated across different eRASS, for instance with a single significant high state, but is present in more than one eRASS. For this, the pool of alerts was obtained from sources that show one or more visits significantly brighter than two (preceding and following) fainter ones. This single-flare search (in contrast with the repeated-flare version described above) yielded ≈3500 alerts per survey, of which ≈450 remained after high significance proper motion was excluded (as flaring coronally emitting stars heavily dominate among the triggers, Boller et al. 2022), and another ≈450 sources with no proper motion estimate. Then, we took advantage of the presence of multiple eRASS and we selected a subsample of sources which triggered the single-flare search in more than one eRASS (i.e. spreading the repeating variable behaviour across the different surveys rather than within the single one). This reduced the alerts to a handful of sources for further inspection. This is instead the method that found eRO-QPE3 (A24).

Furthermore, after the discovery of the first three QPEs we attempted a third set of criteria. Another subset of these ≈3500 alerts, obtained with a single flare per eRASS, selected alerts in which all the faint visits are consistent with background and the signal is detected in one or two consecutive visits only, regardless of whether it significantly occurred in more than one eRASS or not. Effectively, this would consist in relaxing the repeating nature of the flare, but limiting to the highest significance of single flares per eRASS. This yielded ≈200 alerts per survey, of which ≈10 remained after high significance proper motion was excluded, and another ≈30 of sources with no proper motion estimate. The same visualization procedure mentioned above was applied to these alerts as well. This is the method that found eRO-QPE4 (A24).

Finally, we note that all these searches have some obvious biases. For instance, they would likely miss potential QPE sources for which eruption last several days (i.e., longer than the baseline of the single eRASS) and repeat on timescales of tens of days, so that the source is alternatively detected or not in consecutive eRASS, but it is constant in each (e.g., Evans et al. 2023; Guolo et al. 2024). However, we note that the implementation of the single flare algorithm, cross-matched across the various eRASS, would allow eROSITA to detect short-duration eruptions (< 12 h, namely the span of three eROdays) that repeat on the timescales of several weeks or months. However, eROSITA has not yet found any such QPE sources (e.g., later confirmed to have recurrences over days–weeks and still durations over hours). Perhaps, this suggests that the observed relation between average recurrence and average durations (e.g., Chakraborty et al. 2021; Fig. 13 of A24) is intrinsic to the QPE population.

3. Intrinsic rates of QPEs

In this work, our aim is to provide a first estimate of the intrinsic volumetric rate of QPEs. We adopted the four sources discovered by eROSITA (A21,A24) as representative of the intrinsic population and computed their detection efficiency (Sect. 3.1) and resulting volumetric rates (Sect. 3.2). We discuss our assumptions and their validity in Sect. 3.3 and Appendix A.

3.1. Calculating the detection efficiency

Estimating the detection efficiency of a given survey is a function of the peak luminosity distribution. Since QPEs may show appreciably different peak luminosity values within the time span of a survey (M19; G20; A21; A24; Chakraborty et al. 2021), we estimated their averaged values from available data. For eRO-QPE1, we used the average between the main peaks of the two XMM-Newton observations (A21), which are representative of most of the time evolution shown by the source (Chakraborty et al. 2024), namely erg s⁻¹. For eRO-QPE2, the average peak luminosity in the XMM-Newton observation presented in A21: erg s⁻¹. For eRO-QPE3, since the peak luminosity decayed over time (A24), we adopted the average peak luminosity between the eRASS1 to eRASS4 epoch, erg s⁻¹, when the source was actually detected by our algorithm. For eRO-QPE4, we used the average between the three XMM-Newton peaks, namely erg s⁻¹ (A24).

Regarding the detection efficiency, the key element is adopting the correct sensitivity of the searching algorithm. In searches for one-off transients that aim to simply detect a given source, the instrument sensitivity can be used. For our QPE search in particular (Sect. 2), we required the presence of a flux state somewhat brighter than the eRASS sensitivity in the single eROday and a specific pattern in the light curve, both in a single eRASS, or across eRASS (see A21; A24). Furthermore, visual inspection was necessary to confirm the presence of purely alternating variability patterns adding a subjective element that is not straightforward to quantify. Nonetheless, we identified eRO-QPE2 (A21) as the QPE source with the lowest variability amplitude discovered so far with our search. In particular, we adopted as the QPE sensitivity threshold of a given visit the faintest count rate observed for eRO-QPE2 among the brightest that triggered the variability code (see Fig. 2 of A21). This corresponds to a count rate of 0.43 ± 0.12 in the 0.6 − 2.3 keV band (A21). This count rate can be interpreted as the lowest eRASS count rate that a QPE bright state can have to be considered significantly variable according to our search method. Based on our experience, any variable source with a putative bright state fainter than this threshold, would not have been considered a reliable QPE candidate and would not have been followed up. Thus, we convert this eROSITA sensitivity observed flux to the corresponding count rate that would be observed by XMM-Newton (for eRO-QPE2, eRO-QPE3 and eRO-QPE4) and NICER (for eRO-QPE1), using WebPIMMS³. We note that this count rate threshold, once converted to a flux with the spectral shape observed for each QPE source, also defines the maximum observability distance and the discovery volume of a given QPE source, given its observed average peak luminosity.

For each of the discovered eROSITA QPEs, we adopted a mock eRASS sampling composed of as many eROdays as in two eRASS surveys of the actual light curves (see A21; A24). This filter, which mimics eROSITA’s observational sampling during the survey, was then used to scan the XMM-Newton and NICER light curves of each confirmed QPE source. As an eROday visit exposure is ∼40 s long, we shift this filter by ∼40 s for enough times to cover at least one or a few full cycles of each QPE source. This time shift is adopted to be ∼4.5 d for eRO-QPE1 (given the high scatter observed; A21; Chakraborty et al. 2024), ∼5 h for eRO-QPE2, ∼20.5 h for eRO-QPE3 and ∼20 h for eRO-QPE4. If needed, the observed light curves are looped, ensuring that the observed average quasi-periodicity is maintained within the observed scatter. Then, we computed the detection efficiency as follows. For each scan with the eRASS mock filter, we consider it a QPE detection if at least two mock visits are found to be brighter than the count rate threshold (described above), with at least one visit between them that is below the threshold. Otherwise, the scan is considered unsuccessful. This directly compares with our actual QPE search (Sect. 2). The inferred efficiency, computed for each source separately, is therefore the ratio of the number of QPE detections to the total number of scans attempted with the filter. This efficiency can be interpreted as the probability that a given existing QPE source, with properties like those discovered, would trigger our search with a baseline of two eRASS. Since QPEs outlast the entire baseline given by all eRASS, we treat the two eRASS as consecutive for simplicity. We repeat the efficiency calculation 10 000 times, drawing the count rate threshold from a Gaussian distribution given by the observed mean and standard deviation (Fig. 1). We recall that this threshold is converted from the observed eRASS count rate of eRO-QPE2 with its 1σ uncertainty (A21). For a given QPE source (QPE_n), each draw i provides a different efficiency (ξ_i, Fig. 1) and a different maximum detectability luminosity distance (d_Lmax, i), volume (V_max, i, corrected by δA = 0.5 since eROSITA-DE contains half of the entire sky), and redshift (z_max, i).

Fig. 1. Histograms of the distribution of detection efficiency (ξ) of the QPE sources discovered by eROSITA. The sensitivity threshold required to significantly detect an eruption is drawn 10 000 times from a Gaussian distribution. Each iteration of ξi is then computed by downsampling the observed QPEs light curves through the real eRASS cadence and it represents the probability that our QPE search algorithm would trigger on the light curve of that given source.

3.2. Calculating volumetric rates

We adopted three bins in : 41.7 − 42.3, 42.3 − 42.8 and 42.8 − 43.3. The lowest bound of 41.7 was chosen as it is the lowest peak flux shown in the eRASS light curves by any eROSITA-discovered QPE source. This choice is also supported by the simulations reported in Appendix A, since the eRASS detection efficiency is zero below this luminosity even in the low-z Universe. These bins contain two, one, and one QPE source, respectively. For each luminosity bin Δlog L_j we computed the volumetric rate ℛ_vol, j by maximizing the mean Poisson likelihood obtained from the distribution of 10 000 iteration, each giving (given the rate and k detections) , where λ_i = ℛ_vol, j ξ_i V_max, i δA. This calculation approximates the distribution of QPE sources over luminosity and redshift with the known detected eROSITA QPEs, the total volume with V_max and that the nuisance parameters (e.g., ξ_i and V_max, i) are distributed according to the count rate measurement uncertainties used as sensitivity. The 1σ uncertainties on the mean Poisson distribution are taken at the points where the likelihood drops by a factor of e^−0.5. We plot the resulting X-ray luminosity function of QPEs in the bottom panel of Fig. 2. The top panel of Fig. 2 shows instead the redshift-luminosity plane, where the luminosity uncertainties reflect the diversity in peak luminosity observed rather than the statistical uncertainty on an individual peak estimate. Compared with the luminosity function of TDEs (Yao et al. 2023) we note a similar decrease of a factor of ≈100 per luminosity decade. We fitted a power-law relation through the luminosity function of Fig. 2, with uncertainties on both axes, using Ultranest (Buchner 2021), with , where Φ = dℛ_vol/d log L. We obtained a slope with median (and 16th–84th percentiles) of , with a normalization . We show the median (and 16th–84th percentile contours) power-law function in green in Fig. 2. Integrating this function within its uncertainties above we obtain a median (and related 16th–84th percentiles) volumetric abundance rate of Mpc⁻³.

Fig. 2. QPEs volumetric rates. Top: Distribution of the four eROSITA QPEs (A21; A24) in the luminosity-redshift plane. Bottom: X-ray luminosity function of eROSITA QPEs (green), expressed both in Mpc−3 and per galaxy, as a function of the average peak 0.5 − 2.0 keV luminosity. The three bins adopted in are 41.7 − 42.3, 42.3 − 42.8 and 42.8 − 43.3, respectively containing two, one and one QPE. The green line and related contours represent the median and 16th–84th percentiles of the fit power-law model, with a slope ≈ − 2. The two gray points represent the rates from eRO-QPE1 and eRO-QPE2 alone (A21) in two bins of 1 dex, to highlight the small effect of adding new discoveries.

We assume a volumetric density of galaxies with a stellar mass within 10^{8.5 − 10.5}, which includes that of the known eROSITA QPEs (A21,A24), of ∼1.65 × 10⁻² Mpc⁻³ (Baldry et al. 2012). Therefore, our volumetric abundance rate corresponds to a per galaxy abundance rate of gal⁻¹, above an intrinsic average . This implies that at any time, roughly one in 10 000 galaxies in that stellar mass range is erupting quasi-periodically. Relating this to a formation rate relies on assumptions on the typical QPE lifetime. Each secure QPE source has been emitting eruptions since discovery (even if perhaps not continuously, Miniutti et al. 2023a) for up to more than 20 yr (G20). Therefore, we assume a typical lifetime of τ_life ≈ 10 yr as an example, but leave it as a free parameter. This yields a formation rate ℛ_gal/τ_life of approximately a ≈0.4 × 10⁻⁵(τ_life/10 yr)⁻¹ per galaxy per year or ℛ_vol/τ_life ≈ 0.6 × 10⁻⁷(τ_life/10 yr)⁻¹ per Mpc³ per year.

3.3. Assumptions and approximations

Since the overall QPE activity, excluding potential short phases with disappearance (Miniutti et al. 2023a), appears to last from years to decades (M19; G20; A21; A24; Miniutti et al. 2023a), QPEs outlast a given survey. Therefore, the observed rates can be directly compared to abundance volumetric rates, while inferring formation rates strongly relies on the unknown typical lifetime of QPEs. The calculation of intrinsic abundance rates itself is also affected by two major factors. One is astrophysical, namely the knowledge of the intrinsic QPE population, its intrinsic peak luminosity, duration and recurrence distributions and if and how they relate. The other is observational, namely the detection efficiency of a given instrument. Estimating the latter is possible with eROSITA, given its relatively homogeneous survey strategy (Merloni et al. 2024) and our detection algorithm purposely designed to search for QPE-like flares (see Sect. 3.1). The former is by nature unknown and it is particularly difficult for QPEs given the wide diversity in timing properties (i.e. duration and recurrence of eruptions and their evolution over time) within only a handful of known sources (M19; G20; A21; A24; Chakraborty et al. 2021, 2024; Arcodia et al. 2022; Miniutti et al. 2023a,b; Giustini et al., in prep.).

In this work we adopt the properties of the four known eROSITA QPEs as representative of the population. The arguments presented in the previous sections had the objective of providing a qualitative validation of this approach. For instance, eRO-QPE2 (A21) and eRO-QPE3 (A24) have compatible average peak X-ray luminosities, but widely different burst durations (∼0.5 h and ∼3 h, respectively) and recurrences (∼2.4 h and ∼20 h, respectively). Most importantly, they reside in the most intrinsically populated luminosity bin. Therefore, as long as the burst properties of eRO-QPE2 and eRO-QPE3 span the bulk of the intrinsic population, we do not expect to dramatically underestimate or overestimate the intrinsic rates provided here. Furthermore, none of the QPE sources discovered so far seems to have a severely skewed efficiency distribution toward high values (Fig. 1), and thus eROSITA did not select only the sources that are the easiest to find with its sampling. More quantitatively, we verified these assumptions a posteriori by simulating QPE light curves with various values of peak luminosity, recurrence and duration of the eruptions. We computed the eRASS detection efficiency of these mock light curves and concluded, based on the current knowledge on QPEs, that the known eROSITA-discovered QPE sources are not a significantly biased draw from the intrinsic population. Hence, they can indeed be used for a meaningful first estimate of their intrinsic rate. We report more details of the simulations in Appendix A.

Nevertheless, several complications remain, that may affect the calculations presented in this work. In particular, the discovery rate of QPEs in the eROSITA data does, to some extent, rely on subjective criteria based on visual inspection, as we discussed at length in Sect. 2. However, we note that eRO-QPE1 and eRO-QPE2 were the only two QPEs found with a single consistent method (A21) using only eRASS1 and eRASS2 data. However, performing the same calculation above for eRO-QPE1 and eRO-QPE2 with two luminosity bins of one dex (41.5–42.5 and 42.5–43.5) seems to yield consistent results (gray points in Fig. 2). Finally, we note that the inferred rates are not corrected for X-ray absorption. Known QPE sources are rather unabsorbed (M19; G20; A21; A24; Chakraborty et al. 2021), with only eRO-QPE2 being absorbed by a moderate column density (N_H ∼ 3 × 10²¹ cm⁻²; A21). As a rough estimate, we can consider that the current rates underestimate the intrinsic QPE population by a factor ⪅2, following results of the fraction of obscured nuclear super-massive black holes (e.g., Buchner et al. 2015; Carroll et al. 2023).

4. Summary and future prospects

X-ray quasi-periodic eruptions are the most recent addition to the group of extragalactic nuclear transients. In this work, we report the first volumetric abundance rate of QPEs, computed with the first four eROSITA discoveries (A21,A24). We obtained Mpc⁻³ above an intrinsic average , or gal⁻¹. This yields a formation rate ℛ_gal/τ_life ≈ 0.4 × 10⁻⁵(τ_life/10 yr)⁻¹ gal⁻¹ yr⁻¹ or ℛ_vol/τ_life ≈ 0.6 − 7(τ_life/10 yr)⁻¹ Mpc⁻³ yr⁻¹, which depends, however, on the unknown QPE lifetime τ_life.

This formation rates agrees to zeroth-order with that inferred by some of the proposed models on QPEs triggered by secondary orbiters (e.g., Zhao et al. 2022; Lu & Quataert 2023; Linial & Metzger 2023), albeit within the large uncertainties and variables in play. However, the qualitative predictions on the emission mechanism and resulting luminosity of these modeled systems prevent us from making more quantitative comparisons with our observed rates, computed above . Instead, here we attempt a more agnostic comparison between QPE rates and those of other extragalactic nuclear transients like TDEs. Given the apparent dichotomy (e.g., Malyali et al. 2023; Guolo et al. 2023) between TDEs discovered in the optical (e.g., Yao et al. 2023) and X-rays (e.g., Sazonov et al. 2021), now with the latest discoveries via infrared selection (Masterson et al. 2024), this comparison requires particular care. The most direct and thus reliable comparison can be performed with the X-ray TDE rate from Sazonov et al. (2021). Thus, we integrate our QPE luminosity function (Fig. 2) starting from the same value () and obtain a formation rate ℛ_formation = ℛ_vol/τ_life ≈ 2 × 10⁻⁹(τ_life/10 yr)⁻¹ Mpc⁻³ yr⁻¹. We plot ℛ_formation as a function of τ_life in Fig. 3 together with the X-ray TDEs value from Sazonov et al. (2021). Based on Fig. 3, and assuming the QPE = TDE + EMRI scenario (e.g., Franchini et al. 2023; Linial & Metzger 2023), this would imply that on the order of ≈0.01 (τ_life/10 yr)⁻¹ galaxies with (X-ray bright) TDEs could show signatures of QPEs later on. Integrated formation rates of optically selected (Yao et al. 2023) and infrared-selected (Masterson et al. 2024) TDEs agree with the X-ray TDEs estimate, although they should be compared with caution given the complex time-evolving and inhomogeneous multiwavelength properties of TDEs (Guolo et al. 2023; Masterson et al. 2024). This prediction can be tested as future larger samples of confirmed X-ray TDEs and QPEs are unveiled, provided systematic and homogeneous follow-up campaigns are performed at late times.

Fig. 3. QPE formation rate ℛformation as a function of τlife for median and 1σ range (green line and contour, respectively). It is computed above for direct comparison with the X-ray TDE formation rate from Sazonov et al. (2021), shown in dark red. TDE rates from optical (Yao et al. 2023) and infrared (Masterson et al. 2024) searches are also shown (see legend). Within the QPE = TDE + EMRI scenario (e.g., Franchini et al. 2023; Linial & Metzger 2023), on the order of ≈0.01 (τlife/10 y)−1 galaxies with TDEs could show signatures of QPEs later on.

4.1. Predictions for LISA

In the case that the EMRI scenario is unambiguously confirmed as the origin for QPEs, the volumetric rate estimate presented in this work will be the first ever observational constraint on the EMRI rate, which is of paramount importance for future-generation GW detectors such as LISA and Tianqin (Amaro-Seoane 2018). Since not all EMRIs are expected to emit QPEs, the QPE rate will be a lower limit. Even if statistically every EMRI is likely to experience an independent TDE that throws gas on it with the potential of revealing it (Linial & Metzger 2023), the unknown modulo is the fraction of EMRIs for which the secondary is of the nature (a star, and of which kind, or a black hole), mass and orbital configuration required to produce observable X-ray eruptions. Currently, even within the proposed EMRI models it is still debated whether the orbiter is a star (Xian et al. 2021; Linial & Metzger 2023; Tagawa & Haiman 2023) or a black hole (Franchini et al. 2023), or even whether it can be both. Hence, we only attempt to provide some useful values under the assumption that QPEs are EMRIs, and we start by assuming they are black hole EMRIs for simplicity. We start by providing the abundance rate of QPEs per massive black hole (MBH, ℛ_MBH) integrating the mass function from Merloni & Heinz (2008) in the range log M_BH = 5.0 − 7.5, which yields MBH⁻¹. Hence, there are at least as many EMRIs per MBH, of which LISA is expected to detect a fraction of ∼(10 − 40)% (Babak et al. 2017). In particular, the EMRI population models predict a cosmic rate of EMRIs in units of ∼(10 − 2 × 10⁴) yr⁻¹ (Babak et al. 2017) Thus, we integrate the QPE abundance rate ℛ_vol up to z = 1 (where most EMRI population models are complete in terms of LISA detection rates; Babak et al. 2017) using the cosmology from Hinshaw et al. (2013), and obtain a 1σ interquantile range of (3 − 85)×10³ (τ_life/10 yr)⁻¹ yr⁻¹. We show this lower limit prediction in green in Fig. 4. Conversely, if QPEs are stellar-EMRIs, we can attempt to scale their rate to that of black hole EMRIs. Using the framework of Linial & Sari (2022, 2023) and assuming a 10% binary fraction at the influence radius and that the relative abundance of stellar-mass black holes and stars in the nuclear star cluster is 10⁻³, the rate of the former is ≈(1 − 10)% of the latter (Linial, I., priv. comm.). This implies that the 1σ interquantile range of the lower limit to the rate of EMRIs detectable by LISA is in the range ≈(30 − 8500) (τ_life/10 yr)⁻¹ yr⁻¹. We show this lower limit prediction in orange in Fig. 4. Both estimates of stellar and black hole EMRIs disfavor the pessimistic population models (Babak et al. 2017), shown as a gray hatched area in Fig. 4. In particular, the QPE rate lower limit would be quite stringent if QPEs are the electromagnetic counterparts of black hole EMRIs.

Fig. 4. Event rates of BH-EMRIs within z = 1. Predicted theoretical EMRI rates (hatched gray area) from Babak et al. (2017) are compared to those inferred from the QPE formation rate, which have to be interpreted as lower limits as not all EMRIs, even if QPEs are such, are expected to emit eruptions. In the case that QPEs are BH-EMRI, a direct comparison can be made (green 1σ interquantile range); instead, if QPEs are stellar-EMRIs the constraint for BH-EMRIs (orange 1σ interquantile range) requires further extrapolations. Nonetheless, and particularly in the former case, pessimistic EMRI population models would be disfavored, with huge implications for LISA detections.

4.2. Predictions for future-generation X-ray observatories

As a by-product of the inferred volumetric rates, we can provide predictions for future-generation X-ray telescopes for the first time. So far, known QPEs have flared in the soft X-ray band only⁴ (with most signal below 2 keV) and, as we have found in this work, their abundance rates are relatively low. Therefore, a soft X-ray telescope with a large effective area and a large field of view is desirable to systematically discover QPEs. In practice, as we have shown with our eROSITA search (Sect. 2), in order to significantly resolve an eruption and recognize it as such within the data stream of a given telescope, the peak flux has to be somewhat brighter than the instrument’s sensitivity. For simplicity, we adopt this eruption-detection threshold as three times the nominal sensitivity of a given instrument and assume a volumetric rate of ∼0.6 × 10⁻⁶ Mpc⁻³ at L_x, peak ∼ 10⁴² erg s⁻¹.

For instance, the wide-field X-ray telescope (WXT) on the upcoming Einstein Probe (EP) mission (Yuan et al. 2022) boasts a wide instantaneous field of view, although its sensitivity is only tens of times shallower than the brightest QPE peak observed so far. Therefore, simulations of stacked images are required to make further predictions. Instead, we find that the Wide Field Imager (Meidinger et al. 2020) on board Athena (Nandra et al. 2013) may detect ≈2 QPE sources per 15 deg² survey. Whether these sources are immediately discovered as QPE emitters or not depends on the exposure of the survey (t_exp) compared to the QPE duty cycle or recurrence time (t_recur), such as a factor η = t_exp/t_recur. With a typical recurrence time of ≈10 h = 36 ks, either a single moderately deep 10 ks 15 deg² survey or a wider (100 deg²) and shallower (1.5 ks) survey would discover of order unity QPE sources. Here, we adopted a WFI sensitivity of (three times) ∼10⁻¹⁵ erg s⁻¹ cm⁻² in 1 ks from Piro et al. (2022).

Finally, we note that for an efficient QPE discovery machine, repeating the same survey in the same area for a few times is desirable, as it would compensate the factor η after a number N_s ∼ 1/η of surveys. However, performing more than N_s would not be effective as eruptions from the same QPE source would be detected, rather than new QPE sources. Thus, an ideal QPE discovery machine would move to a different sky area after N_s ∼ 1/η surveys and complete as much of the entire sky as possible.

References

Appendix A: Simulating light curves of an empiric population of QPE sources

As discussed in Section 3.3, the accuracy of the intrinsic rates provided in this work relies on the assumption that the known QPE sources discovered by eROSITA so far (A21,A24) are fairly representative of the intrinsic observable QPE population. In particular, the assumption is that the eROSITA detection efficiency probed by the known QPEs is not skewed toward larger values compared to that of the intrinsic population. Here we verified that this assumption is indeed reasonable, by simulating QPE light curves and measuring their detection efficiency. The efficiency of any survey depends on the peak luminosity, and that of eROSITA (given its sampling; Merloni et al. 2024), depends on the duration and recurrence of eruptions as well. In this work, we chose to simulate an empirical QPE population based on the data currently available. In particular, fitting data from A24 we obtained a lack of correlation between average peak luminosity and average recurrence time. The fitted log-linear relation (top panel of Fig. A.1) shows a slope consistent with zero (0.3 ± 0.7). Hence, at least for luminosity and recurrence we adopt an agnostic 2D grid with no correlation between the two quantities. However, we took advantage from the fact that there seem to be a relation between the average recurrence and duration across QPE sources (see the bottom panel of Fig. A.1 and A24). We fitted a linear relation between the logarithm of average recurrence and duration, using the values reported in A24. We obtain a normalization of , a slope of , and intrinsic scatter (see the bottom panel of Fig. A.1). Hence, we generated the QPE population drawing from a grid in peak luminosity () and average recurrence (⟨t_recur⟩), the former in the range 41 − 44 log(erg s⁻¹) and the latter in the range 1 − 30 h. These values bracket those of known QPE sources. A scatter of 10% is assumed on both quantities, which is a representative value based on all the secure QPE sources (M19,G20,A21,A24). This scatter is not only observed in the sources discovered by eROSITA, and thus it is not a bias due to eROSITA’s cadence. Then, for each value of ⟨t_recur⟩, we extracted a value for the average QPE duration (⟨t_dur⟩) from the relation in the bottom panel of Fig. A.1, within its 1σ uncertainties.

Fig. A.1. Relation between average recurrence time (trecur) and average peak X-ray luminosity (top panel) or average eruption duration (tdur, bottom panel) in QPE sources. The data is from Fig. 13 in A24. In the top panel the slope is consistent with zero (red dotted line), and thus with quantities being uncorrelated given the available data. In the bottom panel the slope is consistent with being linear in the log-log plane.

We built 100 light curves for each grid value of and ⟨t_recur⟩ (and the inferred ⟨t_dur⟩) and for two redshift bins (z = 0.02 and 0.05, representative of the values spanned by known eROSITA QPEs, Fig. 2). From we obtained an average peak count rate using WebPIMMS, which we converted to a quiescent state assuming an amplitude of 100. Since the sensitivity within the single visit is typically much shallower than this value, the assumed amplitude is generally unimportant as the faint state will be eROSITA’s background. In addition to this quiescence, we generated Gaussian eruptions by drawing and ⟨t_recur⟩ within their 10% scatter, with FWHM converted from one third of the rise-to-decay duration inferred from the bottom panel of Fig. A.1. An example of generated light curve is shown in Fig. A.2, for and ⟨t_recur⟩∼10 h within their 10% scatter, and ⟨t_dur⟩∼1 h (from Fig. A.1).

Fig. A.2. Example of simulated light curve, drawn from and ⟨trecur⟩∼10 h within their 10% scatter, and ⟨tdur⟩∼1 h (from the relation in the bottom panel of Fig. A.1). An example of a successful (green points) and an unsuccessful (red) mock eROSITA scan is shown. The red line shows the eRASS sensitivity for bright states, the gray shows the threshold for faint states used in the simulations.

Each light curve has an associated efficiency value ξ, which can be seen as the probability that a given light curve will be seen by eROSITA, given its sampling, as QPE-like. This probability is computed as the fraction of successful scans out of the total possible, which are performed shifting a mock eRASS sampling through the light curve for enough times to cover a full QPE cycle. A scan is considered successful if at least two visits are found to be in a bright state, with at least a faint state in between, and combinations thereof (e.g., green points in Fig. A.2). A scan is otherwise considered unsuccessful (e.g., red points in Fig. A.2). A visit is considered in a high state if its count rate is higher than the sensitivity value of ∼0.42 c/s adopted in Sect. 2; a visit is considered in a low state if the count rate is lower than a tenth of this value. For the brightest luminosity bins in the simulation, the quiescent state might be brighter than both of these values, so we adopted a low state value as 1/50 of the peak and a high state value as ten times the low value. This is in qualitative analogy with the real QPE search presented in Sect. 2. We note that the choice of these threshold numbers have some impact on the single efficiency draw, but very little impact on the distribution after the 100 draws. Regarding the length of the mock eRASS baseline, we define a mock eRASS sampling as 18 data points, each with an exposure of exactly 40 s and separated by 4 h (in comparison with the real analysis in Sect. 2). The latter is the typical revolution time of the field of view around the ecliptic plane before coming back to the same point in the sky (Merloni et al. 2024). Instead, the number of visits simulated corresponds to two eRASS (i.e. nine per eRASS) and it is based on the typical length of the baseline for the four eROSITA QPEs (A21,A24). This preference for ∼10−eROday baselines is likely due to a sweet spot between having a large enough area on the sky and a long enough baseline. Longer baselines than these are obtained at the ecliptic poles where, however, the area is small (Merloni et al. 2024; Bogensberger et al. 2024). Since we correct for the efficiency given by the sampling, we do not expect that assuming this cadence imprints any significant bias on our calculations. Future work will be done in simulating a QPE population on the full sky and running mock source detection and light curve generation depending on the sky location, but it is beyond the scope of this work.

After this procedure each light curve has an associated ξ, and each bin in the simulations has 100 light curves generated, thus 100 ξ values. We show the outcome of these simulations in the two panels of Fig. A.3, one per redshift bin adopted (z = 0.02 and 0.05). The average ξ per bin is mapped in a 2D grid of luminosity and recurrence. The first takeaway is that in each redshift bin ξ is fairly homogeneous above the sensitivity and that this sensitivity does not significantly depend on the timing properties (⟨t_recur⟩, ⟨t_dur⟩). There are some obvious exceptions, corresponding to periods resonating with the 4 h cadence in the eRASS sampling. Further insights on the fact that eROSITA discoveries are likely not severely biased come from a comparison with their individual efficiencies (shown in Fig. 1) and the simulated ones. We plot the known eROSITA QPEs as squares in Fig. A.3, color-coded in the same way. Their location is not in the regions were the efficiency is the highest, showing that the eROSITA discoveries have not selected the easiest QPEs to be found. An obvious exception is eRO-QPE3 which appears as an outlier in the top panel, surrounded by a region with ξ ∼ 0. This is probably due to the fact that eRO-QPE3 showed time-varying amplitude over the eRASS surveys (A24), so our simplistic simulations have not grasped that sources can move across the grid. We also note that the average recurrence is only estimated crudely for eRO-QPE3 (A24). Figure A.4 shows the 2D maps collapsed on the luminosity axis, for both redshift bins. It shows, even more quantitatively, that the QPEs discovered by eROSITA (shown as black points) are representative of the efficiency values within their luminosity bins. We conclude the known QPEs are, given our current knowledge on QPEs, a fair draw of their intrinsic population. Hence, this justifies their use to infer the volumetric rates presented in this work.

Fig. A.3. Maps of the simulated eROSITA detection efficiency ξ for QPE sources in a 2D grid of average peak luminosity () and average recurrence time (⟨trecur⟩), for the case of z = 0.02 (top) and z = 0.05 (bottom). Squares with a red contour highlight the known QPE sources in the respective redshift bins. Apart from the grid points around the periods resonating with eROSITA’s 4 h cadence, ξ is fairly homogeneous above the cut provided by the sensitivity of each redshift bin.

Fig. A.4. Simulated detection efficiency ξ as a function of peak luminosity in the redshift bin z = 0.02 (blue) and z = 0.05 (orange). The plot is a version of Fig. A.3 collapsed on the y-axis. The black points indicate the known eROSITA QPEs (A21,A24). They are representative of the ξ values of the simulated population at their respective luminosity.

Appendix B: Acknowledgments

We thank the anonymous referee for their report, which improved the look of this manuscript. R.A. is grateful to I. Linial for the calculations on the stellar- versus BH-EMRI rates, and to A. Sesana, B. Metzger and N. Stone for insightful discussions on the volumetric rates and to E. Kara and R. Mushotzky for useful discussions on future X-ray missions. R.A. received support for this work by NASA through the NASA Einstein Fellowship grant No HF2-51499 awarded oby the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. GP acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 865637), support from Bando per il Finanziamento della Ricerca Fondamentale 2022 dell’Istituto Nazionale di Astrofisica (INAF): GO Large program and from the Framework per l’Attrazione e il Rafforzamento delle Eccellenze (FARE) per la ricerca in Italia (R20L5S39T9). We acknowledge the use of the matplotlib package (Hunter 2007).

This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tuebingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universitaet Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23.

The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. Pipeline processing and analyses of the data were supported by NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation.

NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. LBNL is managed by the Regents of the University of California under contract to the U.S. Department of Energy.

This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.

BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” Grant #XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant #114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant #12120101003, #11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

All Figures

Fig. 1. Histograms of the distribution of detection efficiency (ξ) of the QPE sources discovered by eROSITA. The sensitivity threshold required to significantly detect an eruption is drawn 10 000 times from a Gaussian distribution. Each iteration of ξi is then computed by downsampling the observed QPEs light curves through the real eRASS cadence and it represents the probability that our QPE search algorithm would trigger on the light curve of that given source.

In the text

Fig. 2. QPEs volumetric rates. Top: Distribution of the four eROSITA QPEs (A21; A24) in the luminosity-redshift plane. Bottom: X-ray luminosity function of eROSITA QPEs (green), expressed both in Mpc−3 and per galaxy, as a function of the average peak 0.5 − 2.0 keV luminosity. The three bins adopted in are 41.7 − 42.3, 42.3 − 42.8 and 42.8 − 43.3, respectively containing two, one and one QPE. The green line and related contours represent the median and 16th–84th percentiles of the fit power-law model, with a slope ≈ − 2. The two gray points represent the rates from eRO-QPE1 and eRO-QPE2 alone (A21) in two bins of 1 dex, to highlight the small effect of adding new discoveries.

In the text

Fig. 3. QPE formation rate ℛformation as a function of τlife for median and 1σ range (green line and contour, respectively). It is computed above for direct comparison with the X-ray TDE formation rate from Sazonov et al. (2021), shown in dark red. TDE rates from optical (Yao et al. 2023) and infrared (Masterson et al. 2024) searches are also shown (see legend). Within the QPE = TDE + EMRI scenario (e.g., Franchini et al. 2023; Linial & Metzger 2023), on the order of ≈0.01 (τlife/10 y)−1 galaxies with TDEs could show signatures of QPEs later on.

In the text

Fig. 4. Event rates of BH-EMRIs within z = 1. Predicted theoretical EMRI rates (hatched gray area) from Babak et al. (2017) are compared to those inferred from the QPE formation rate, which have to be interpreted as lower limits as not all EMRIs, even if QPEs are such, are expected to emit eruptions. In the case that QPEs are BH-EMRI, a direct comparison can be made (green 1σ interquantile range); instead, if QPEs are stellar-EMRIs the constraint for BH-EMRIs (orange 1σ interquantile range) requires further extrapolations. Nonetheless, and particularly in the former case, pessimistic EMRI population models would be disfavored, with huge implications for LISA detections.

In the text

Fig. A.1. Relation between average recurrence time (trecur) and average peak X-ray luminosity (top panel) or average eruption duration (tdur, bottom panel) in QPE sources. The data is from Fig. 13 in A24. In the top panel the slope is consistent with zero (red dotted line), and thus with quantities being uncorrelated given the available data. In the bottom panel the slope is consistent with being linear in the log-log plane.

In the text

	Fig. A.2. Example of simulated light curve, drawn from and ⟨trecur⟩∼10 h within their 10% scatter, and ⟨tdur⟩∼1 h (from the relation in the bottom panel of Fig. A.1). An example of a successful (green points) and an unsuccessful (red) mock eROSITA scan is shown. The red line shows the eRASS sensitivity for bright states, the gray shows the threshold for faint states used in the simulations.
In the text

Fig. A.3. Maps of the simulated eROSITA detection efficiency ξ for QPE sources in a 2D grid of average peak luminosity () and average recurrence time (⟨trecur⟩), for the case of z = 0.02 (top) and z = 0.05 (bottom). Squares with a red contour highlight the known QPE sources in the respective redshift bins. Apart from the grid points around the periods resonating with eROSITA’s 4 h cadence, ξ is fairly homogeneous above the cut provided by the sensitivity of each redshift bin.

In the text

	Fig. A.4. Simulated detection efficiency ξ as a function of peak luminosity in the redshift bin z = 0.02 (blue) and z = 0.05 (orange). The plot is a version of Fig. A.3 collapsed on the y-axis. The black points indicate the known eROSITA QPEs (A21,A24). They are representative of the ξ values of the simulated population at their respective luminosity.
In the text