© The Authors 2024

1. Introduction

The advent of wide-field cadenced surveys across the electromagnetic spectrum during the last decade opened a new window to the realm of extra-galactic transients, with the discovery of a few classes of repeaters. Some galactic nuclei have been shown to undergo outbursts recurring on timescales of months to decades in optical and/or X-rays (Payne et al. 2021; Wevers et al. 2023; Liu et al. 2023; Malyali et al. 2023), others over tens of days (Evans et al. 2023; Guolo et al. 2024) and even down to surprisingly fast recurrences of a few hours, such as the X-ray quasi-periodic eruptions (QPEs; Miniutti et al. 2019; Giustini et al. 2020; Arcodia et al. 2021, hereafter M19; G20; A21). Whether all these flavors are connected within a single recipe remains to be seen, although the leading models of each subclass appear to relate somehow to the dynamics of stellar objects in galactic nuclei (e.g., Linial & Quataert 2024).

In particular, QPEs are sharp X-ray bursts that last a few hours and repeat in a quasi-periodic manner every several hours to about a day. So far, these peculiar bursts have been observed from massive black holes (MBHs) of M_BH ≈ 10^{5 − 6.7} M_⊙ (Shu et al. 2017; M19; G20; A21; Wevers et al. 2022) residing in the nuclei of nearby low-mass galaxies (e.g., stellar masses of M_* ≈ 1 − 3 × 10⁹ M_⊙; A21). To date, no associated optical, UV, IR or radio flares have been observed, although most of the current multiwavelength photometry is likely contaminated or dominated by either the host galaxy emission or that of the accretion disk. In fact, when in quiescence in between QPEs, most sources are still detected in X-rays, with a soft spectrum reminiscent of the Wien tail of a radiatively efficient accretion disk (with a peak temperature kT ∼ 40 − 70 eV; M19; G20; A21). Eruptions reach a soft X-ray luminosity of ∼10⁴² − 10⁴³ erg s⁻¹ at their peak, ∼10 − 100 times brighter than the quiescence level, with a spectrum that is hotter when brighter and remains soft and thermal in shape (with a peak temperature kT ∼ 100 − 200 eV; M19; G20; A21). While the timing properties show significant diversity in how regular the eruption arrival times are (M19; G20; A21), all eruptions from all the different QPE sources seem to follow the same spectral evolution: bursts progress with a harder rise than decay at the same X-ray count rate or, modeling the emission as a black body, with a hotter rise than decay at the same luminosity (Arcodia et al. 2022; Miniutti et al. 2023a; Giustini et al., in prep.).

Here, we consider secure QPE sources GSN 069 (M19), RX J1301.9+2747 (Sun et al. 2013; G20), eRO-QPE1 and eRO-QPE2 (A21). In addition, the QPE candidate XMMSL1 J024916.6−041244 (Chakraborty et al. 2021) showed remarkably similar spectral properties but only 1.5 eruptions, impeding its secure classification to date. Recently, the start of a flare consistent, in terms of spectral and timing properties, with those of eRO-QPE1 was seen in an optical tidal disruption event (TDE) “Tormund” (Quintin et al. 2023); however, the lack of repetition to date makes its classification still ambiguous. Further, the source Swift J023017.0+283603 shows repeating soft X-ray outbursts every ∼21 d (Evans et al. 2023; Guolo et al. 2024), although its opposite asymmetry in the burst shape and opposite spectral evolution during the bursts (with a harder decay than rise) does not allow a secure association with QPEs for the time being.

The origin of QPEs is still actively debated, with some models proposing some kind of accretion disk instability (Raj & Nixon 2021; Śniegowska et al. 2023; Kaur et al. 2023; Pan et al. 2023), whilst most suggest that QPEs are triggered by a binary system including a central MBH and a much smaller body orbiting it (King 2020, 2022; Suková et al. 2021; Xian et al. 2021; Zhao et al. 2022; Wang et al. 2022; Metzger et al. 2022; Krolik & Linial 2022; Linial & Sari 2023; Lu & Quataert 2023; Franchini et al. 2023; Tagawa & Haiman 2023; Linial & Metzger 2023). For instance, most of the observational properties can be qualitatively reproduced by shocks caused by the interaction between the smaller orbiter and the accretion flow around the MBH (Xian et al. 2021; Franchini et al. 2023; Tagawa & Haiman 2023; Linial & Metzger 2023). In this framework, the disk is pierced by the orbiter once or twice per orbit, and shocks are caused by an initially optically thick cloud of gas ejected by the collisions (Linial & Metzger 2023). Orbital precession and Lense-Thirring precession of the disk would provide the required departure from exact periodicity (Franchini et al. 2023). A requirement, if not a prediction, for the collisions model is that the accretion flow around these MBHs is not an extended flow typical of active galactic nuclei (AGN), but rather a more compact density distribution. Otherwise the orbiter would sink into the disk plane and/or be ablated by the collision in a more extended disk. This led the latest models to suggest that the disk originates from a TDE, whose role in the QPE emission is solely to throw gas at the preexisting extreme mass ratio inspiral (EMRI) orbiters (Franchini et al. 2023; Linial & Metzger 2023). Interestingly, this theoretical connection between QPEs and a preexisting behavior akin to that of TDEs was previously suggested based on observational properties alone: both GSN 069 (Shu et al. 2018; Sheng et al. 2021; Miniutti et al. 2023a) and the two candidate QPEs (Chakraborty et al. 2021; Quintin et al. 2023) show evidence of past TDEs prior to the onset of QPE (or candidate QPE) behavior.

Clearly, at this stage it is fundamental to find more QPEs to draw any significant conclusions on the population and its physical origin. After the serendipitous discovery of the first QPE emitting source (M19), QPEs were either found through dedicated searches in the X-ray archives (G20; Chakraborty et al. 2021; Quintin et al. 2023) or with a blind search within the live data stream of the eROSITA X-ray telescope (A21). As much as the former method is important to make sure no interesting source was overlooked, wide-area surveys in X-rays provide the only way to systematically detect new QPE sources as they happen in the sky. In A21, we reported results from the first two eROSITA all-sky surveys (Merloni et al. 2024). Here, we report on two new discoveries based on the subsequent eROSITA surveys.

2. Data processing and analysis

2.1. X-ray spectral analysis

Spectral analysis is performed with the Bayesian X-ray Analysis software (BXA) version 4.0.7 (Buchner et al. 2014), which connects the nested sampling algorithm UltraNest (Buchner 2019, 2021) with the fitting environment XSPEC version 12.13.0c (Arnaud 1996), in its Python version PyXspec¹. The continuum model adopted is absorbed by Galactic column density from HI4PI (HI4PI Collaboration 2016) and redshifted to rest-frame using spectroscopic redshifts. These are z = 0.024 for eRO-QPE3 and z = 0.044 for eRO-QPE4 (more details in Sect. 2.3). We quote median and 16th and 84th percentiles (∼1σ) uncertainties from fit posteriors, unless otherwise stated, for fit parameters, flux and luminosity. The bolometric luminosity (labeled with “bol”, e.g., L_disk, bol) is obtained integrating the adopted source model between 0.001 and 100 keV. For non-detections, we quote ∼1σ (∼3σ) upper limits using the 84th (99th) percentiles of the fit posteriors, unless otherwise stated. eROSITA source plus background spectra (Sect. 2.2) were fit including a model component for the background, which was determined via a principal component analysis (e.g., Simmonds et al. 2018) from a large sample of eROSITA background spectra (e.g., Liu et al. 2022). XMM-Newton spectra (Sect. 2.4) were instead fit using wstat, namely XSPEC implementation of the Cash statistic (Cash 1979), given the good counts statistics in both source and background spectra. In plots (e.g., Appendix A.1), data are rebinned, with uncertainty on the summed counts (CTS_tot) computed as $1+\sqrt{{\mathrm{CTS}}_{\mathrm{tot}}+0.75}$, only for visualization purposes.

As discussed in Miniutti et al. (2023a,b), we interpret the QPE emission as thermal-like and we use a simple model black body (zbbody in XSPEC). This is a rather model-independent interpretation given the observed spectral shape, with consequences only on the QPE’s bolometric emission. The soft X-ray emission in quiescence, namely in between the eruptions, is often detected in QPE-sources and interpreted as the inner regions of a radiatively efficient accretion disk. Therefore, we model the quiescence, if detected, with a diskbb in XSPEC². This component is held fixed whilst fitting of the eruption, therefore QPE fluxes are integrated only under the QPE model. For spectra with good counts statistics (e.g., XMM-Newton spectra) we allow the quiescence parameters free to vary within the 10th and 90th percentiles of the posterior distributions of the quiescence-only spectral fits. In this way, the spectral fit parameters and fluxes of the eruptions are marginalized over the uncertainties of the quiescence model. For spectra with lower counts statistics (e.g., eRASS spectra) we freeze the quiescence model in the QPE fits. Additional spectral components are added, if needed to model residuals, and discussed individually in the following Sections. The best-fit model among the ones adopted is selected by inspection of the residuals and by comparing the logarithmic Bayesian evidence (Z), adopting the model with the highest value. In general, these hard residuals, if present, are here interpreted as thermal Comptonization of the disk emission. We use either a simple phenomenological powerlaw (zpowerlw in XSPEC) or a more physically motivated model (e.g., nthComp in XSPEC; Zdziarski et al. 1996; Życki et al. 1999), depending on the counts statistics in the spectra. We note that alternative models may also be used in future work, although we make the simple assumption that disk residuals are due to Comptonization of the disk photons. This is for the sake of simplicity and it is motivated by its nearly ubiquitous presence in accretion flows around black holes.

Fig. 1. eRASS1-5 light curves of eRO-QPE3 in the 0.2 − 2.0 keV band, each scaled by a reference starting time (MJD ∼58874.111, ∼59057.576, ∼59233.903, ∼59422.701 and ∼59604.611 from teRO, 1 to teRO, 5, respectively). Orange data points highlight the putative detection of eruptions, whilst the gray line connects the background level across the eROdays in each eRASS.

2.2. eROSITA

eROSITA (Predehl et al. 2021) is the soft X-ray instrument aboard the Spectrum-Roentgen-Gamma (SRG) mission (Sunyaev et al. 2021). On 13 December 2019 it started observing in survey mode and to this date it has completed four (and started the fifth) of the foreseen eight all-sky surveys (eRASS1-8), each of which is completed in six months. In each survey, every source in the sky is observed for ∼40 s every ∼4 h (i.e., a so-called eROday) for a total number of times within a single eRASS that depends on its location in the sky: around six times on the ecliptic plane and increasing toward higher (ecliptic) latitudes. We have developed an algorithm to look for significant and repeated high-amplitude variability in the eROSITA sources. Light curves are systematically extracted by the eROSITA team with the srctool task of the eROSITA Science Analysis Software System (eSASS; Brunner et al. 2022) from event files version 020. In a nutshell, we select sources showing high-amplitude variability within eROdays of the single eRASS, or across multiple eRASS. After excluding secure Galactic objects using Gaia DR3 (Gaia Collaboration 2016, 2021), we visually inspect their X-ray and multiwavelength properties, although only the former (truly alternating variability pattern and a soft spectrum) are used to consider a source as QPE candidate. A first version of this method was described in A21 and we present more details in a companion paper (Arcodia et al. 2024) presenting the calculation of intrinsic volumetric rates based on the eROSITA QPE search method.

2.2.1. eRO-QPE3

The source eRASSt J140053–284557 (hereafter eRO-QPE3) is located at the astrometrically corrected X-ray position of (RA J2000, Dec J2000) = (210.2222, −28.7665), with a total 1σ positional uncertainty (including a systematic error of ∼1.5″) of 1.6″, based on the latest-available internal catalog from the first four surveys (eRASS:4). The source is detected in eRASS1 with the identifier 1eRASS J140053.3−284558 (Merloni et al. 2024). It was also independently found in a search for TDEs, as a new, relative to the archival upper limit, soft eRASS1 source (Grotova et al., in prep.). During eRASS1 the source was detected in each eROday with a compatible flux level, with the exception of one eROday showing count rate higher by a factor ∼2 − 3 (top left of Fig. 1). However, this intra-eROday variability was not deemed significant by the flare searching code. eRO-QPE3 then triggered an alert from the QPE-candidates searcher in eRASS2, eRASS3 and eRASS4. In eRASS2, the source was fainter compared to eRASS1 (in terms of average count rate), while still detected at all eROdays with an overall constant flux with the exception of one (top middle panel in Fig. 1). In eRASS3-4-5, nearly all the counts were detected in a single eROday, with most of (or all, depending on the eRASS) the remaining eROdays consistent with background (see Fig. 1).

Details on the spectral fits performed on eROSITA data of eRO-QPE3 are reported in Table A.1. The eRASS1 spectrum in quiescence was modeled with a diskbb plus zpowerlw (top panel of Fig. A.3). The median (and 16th, 84th percentiles) peak temperature of the disk model is $k{T}_{\mathrm{in}}=k{T}_{\mathrm{disk}}={100}_{-20}^{+18}$ eV, with a rest-frame flux of $F_{0.2-2.0\mathrm{keV}}=1.4_{-0.6}^{+0.4}\times {10}^{-12}$ erg s⁻¹ cm⁻². Here, we consider the orange point in the top-left panel of Fig. 1 as a bright state part of an eruption. However, we do not exclude the possibility that in eRASS1 no QPEs were present, as the total eRASS1 spectrum can be also modeled with the quiescence model alone (see Fig. A.2). The quiescence model was held fixed during the putative eruption and a black body component was added (bottom panel of Fig. A.3). The median (and 16th, 84th percentiles) peak temperature of the QPE component is $kT=k{T}_{\mathrm{QPE}}={122}_{-15}^{+19}$ eV and its flux is F_{0.2 − 2.0 keV} = (3.0 ± 0.7)×10⁻¹² erg s⁻¹ cm⁻². The median temperature is harder than the quiescence disk’s peak temperature, but compatible within uncertainties. This confirms that the contrast between the putative quiescence and the QPE, if present, is low during eRASS1. The total flux during the bright eROday, quiescence included, in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=4.9_{-0.6}^{+0.7}\times {10}^{-12}$ erg s⁻¹ cm⁻². The eRASS2 spectrum in quiescence is fainter than in eRASS1 (see top panel of Figs. A.5 and 1). The median (and related 16th, 84th percentiles) temperature of the disk is $k{T}_{\mathrm{disk}}={89}_{-60}^{+42}$ eV, with a flux $F_{0.2-2.0\mathrm{keV}}=3.{12}_{-3.11}^{+1.51}\times {10}^{-13}$ erg s⁻¹ cm⁻². The total flux in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=5.8_{-2.9}^{+2.7}\times {10}^{-13}$ erg s⁻¹ cm⁻². As done for eRASS1, this quiescence model is held fixed during the brighter eROday (the orange point in the top-medium panel of Fig. 1) and a black body component is added (bottom panel of Fig. A.5). The median (and 16th, 84th percentiles) peak temperature of the QPE component is $k{T}_{\mathrm{QPE}}={82}_{-9}^{+10}$ eV and its flux is $F_{0.2-2.0\mathrm{keV}}=3.2_{-0.7}^{+1.0}\times {10}^{-12}$ erg s⁻¹ cm⁻². The median temperature is colder compared to the bright eROday in eRASS1, despite the flux being compatible. However, given the short ∼40 s exposure of an eROday compared to the typical QPE duration (0.5 − 7 h, A21), eRASS data catch the eruption at different phases and they must be considered as lower limits of the true QPE peaks. The total flux during the bright eROday, quiescence included, in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=3.9_{-0.8}^{+1.0}\times {10}^{-12}$ erg s⁻¹ cm⁻². Therefore, despite the decreasing quiescence level, the flux in the bright eROday is compatible between eRASS1 and eRASS2. In eRASS3, eRASS4 and eRASS5, most of the signal is observed during a single eROday, therefore no quiescent state is detected. We adopt a disk model to compute upper limits (Table A.1) and we fit the bright state with a black body QPE component only. In eRASS5, the fit is compatible with background (within 3σ, Fig. A.8) even in the putative bright state. We report the fit values in Table A.1 for completeness.

Fig. 2. Same as Fig. 1, but for eRASS4 data of eRO-QPE4. The reference starting time is MJD ∼59457.807. eRO-QPE4 was undetected in the previous eRASS1-3 survey.

Fig. 3. Left panels: 45″ × 45″ cutout of the DESI Legacy Imaging Surveys Data Release 10 (Legacy Surveys/D. Lang (Perimeter Institute)) with the X-ray 1σ position circles in red (eROSITA) and green (XMM-Newton), for eRO-QPE3 (top) and eRO-QPE4 (bottom). Right panels: SALT optical spectra of eRO-QPE3 (top) and eRO-QPE4 (bottom), with inferred spec-z of ∼0.024 and ∼0.044, respectively. The sky spectrum of the red spectral exposure is shown in cyan in arbitrary units in both subpanels.

2.2.2. eRO-QPE4

The source eRASSt J044534–101201 (hereafter eRO-QPE4) is located at the astrometrically corrected X-ray position of (RA J2000, Dec J2000) = (71.3921, −10.2005), with a total 1σ positional uncertainty (including a systematic error of ∼1.5″) of 4.9″, based on the latest-available internal catalog from the first four surveys (eRASS:4). eRO-QPE4 was undetected in the first three eRASS, with a 1σ (3σ) upper limit of F_{0.2 − 2.0 keV} <  3.8 × 10⁻¹⁵ erg s⁻¹ cm⁻² (< 2.8 × 10⁻¹⁴ erg s⁻¹ cm⁻²). In eRASS4, most of the X-ray signal was observed in a single eROday (see Fig. 2). Adopting a black body model for the QPE emission (a behavior in the source confirmed with follow-up observations), we obtain $k{T}_{\mathrm{QPE}}={93}_{-14}^{+18}$ eV and a flux $F_{0.2-2.0\mathrm{keV}}=3.4_{-1.1}^{+1.7}\times {10}^{-12}$ erg s⁻¹ cm⁻². Excluding this eROday, the source is undetected with a 1σ upper limit of F_{0.2 − 2.0 keV} <  6.9 × 10⁻¹⁴ erg s⁻¹ cm⁻², using a diskbb model. The position of eRO-QPE4 was not observed in the partial eRASS5 survey.

2.3. Optical identification

2.3.1. eRO-QPE3

eRO-QPE3 is associated with the galaxy 2MASS 14005331−2846012 at (RA, Dec) = (14:00:53.315, −28:46:01.26). We show the Legacy Survey DR10 optical image in the top left panel of Fig. 3, centered at the X-ray coordinates. The eROSITA position and 1σ accuracy (based on the cumulative eRASS:4 image) are shown with a red cross and circle, respectively, whilst in green we show the analog from the deeper, hence more accurate, XMM-Newton observation. Similarly to other QPE sources, eRO-QPE3 is consistent with the nucleus position, given the available positional accuracy (M19; G20; A21). 2MASS 14005331−2846012 was observed in the grizJHK bands with the seven-channel imager GROND (Gamma-Ray Burst Optical/Near-Infrared Detector; Greiner et al. 2008) at the MPG 2.2 m telescope at La Silla Observatory (Chile) on 2020 February 27th. The data were reduced with the standard GROND pipeline (Krühler et al. 2008), which performs the bias and flat-field corrections, image stacking and astrometric calibration. The photometric calibration was achieved against PanSTARRS (griz) and 2MASS (JHK; Skrutskie et al. 2006). The observed Kron magnitudes in the AB system are 18.37 ± 0.06, 17.60 ± 0.05, 17.29 ± 0.07, 17.05 ± 0.05, 16.82 ± 0.11, 16.60 ± 0.13 and 16.95 ± 0.14 mag, respectively.

We instigated optical spectroscopic follow-up with the Robert Stobie Spectrograph (RSS, Burgh et al. 2003) on the Southern African Large Telescope (SALT, Buckley et al. 2006) on the night of April 24, 2022 (top right panel of Fig. 3). Data were reduced following the procedure outlined in A21. The host galaxy shows a seemingly inactive spectrum, with z = 0.024 based on Calcium H and K absorption lines (top right panel of Fig. 3). We fit the SALT optical spectrum with Firefly, which is a fitting code for deriving the stellar population properties of galaxy spectra (Wilkinson et al. 2017; Maraston et al. 2020). SALT spectra are not calibrated to absolute values (Buckley et al. 2018), therefore we first renormalized the spectrum flux using optical photometry from the GROND data using the Specutil Python package (Earl et al. 2023). Firefly was run twice using two different stellar population models (Maraston & Strömbäck 2011). The run with the population model ELODIE (Maraston & Strömbäck 2011) was adopted to estimate the mean values for stellar-mass (M_*) and star formation rate (SFR). We added in quadrature to the measured statistical errors a systematic error to account for the choice between the two stellar population models, which was computed from the difference between the two mean values. The stellar mass inferred with Firefly for eRO-QPE3 is $M_{\ast }=2.{56}_{-1.40}^{+0.24}\times {10}^9{M}_{\odot }$, whist SFR is $0.{20}_{-0.14}^{+0.02}{M}_{\odot }$ yr⁻¹. The low M_* value is in line with the very compact nature of the host galaxy (Fig. 3, top), and is in line with those inferred for eRO-QPE1 and eRO-QPE2 (A21), which is remarkable given the blind nature of our search. Using scaling relations between black hole and total stellar mass of the galaxy (Reines & Volonteri 2015), we obtain $M_{\mathrm{BH}}=5.3_{-3.5}^{+0.7}\times {10}^6{M}_{\odot }$. Although these scaling relations are not necessarily well-calibrated at low masses, the relation with bulge stellar mass has been shown to hold sufficiently well at lower masses too (Schutte et al. 2019). We performed SED photometry fitting to have an independent estimate of M_*. We collated Legacy Survey DR10 photometry with RainbowLasso ³ and fit the SED with Genuine Retrieval of AGN Host Stellar Population (GRAHSP; Buchner et al., in prep.). The model used in GRAHSP includes AGN components (continuum, emission lines, torus) and galaxy components, both attenuated by dust and redshifted. We obtain $M_{\ast }=0.{74}_{-0.29}^{+0.41}\times {10}^9{M}_{\odot }$, which is slightly lower than the Firefly estimate. This would correspond (Reines & Volonteri 2015) to a black hole mass of $M_{\mathrm{BH}}=0.{93}_{-0.47}^{+0.79}\times {10}^6{M}_{\odot }$.

No significant emission lines are apparent in the spectrum (top right panel of Fig. 3), although subtracting the continuum model from Firefly the presence of faint narrow Hα, [O III] λ5007 and O Iλ6302 lines, and very faint Hβ can be inferred. Since we are mainly interested in computing line flux ratios for narrow-lines diagnostics, we computed them from the continuum-subtracted spectrum. We used the lmfit package (Newville et al. 2014) and adopted a polynomial model plus Gaussians for emission lines. The algorithm fits a narrow line to Hβ, both [O III] lines, O Iλ6302 and Hα, whilst no narrow [N II] lines are visible. We obtain from this log([O III]/Hβ) ∼ 0.42 and log([N II]/Hα) ≲ −0.87, after tentatively computing an upper limit on the [N II] λ6549 line by imposing the presence of a narrow line. These values would place eRO-QPE3 in the star-forming region of narrow lines classifications, whilst using the log([O I]/Hα) = −0.45 the galaxy is instead classified as a LINER (Kewley et al. 2006). Given the limited resolution and tentative detection of some narrow lines, we defer a conclusive classification to future work. However, this result confirms that for these nuclear transients narrow lines classifications are ambiguous and have to be interpreted with care, particularly so since they trace past nuclear activity potentially unrelated to the current transient activity. We obtained an independent estimate of the black hole mass by fitting the SALT spectrum with the Penalized PiXel-Fitting method (pPXF; Cappellari & Emsellem 2004; Cappellari 2017, 2023). The fit velocity dispersion is ∼204 km s⁻¹. Estimating an instrumental broadening of 6 Å from the narrowest arc lines for the grating angle settings adopted, the instrumental velocity dispersion computed around the median wavelength ∼150 km s⁻¹. Subtracting this dispersion in quadrature and summing in quadrature the template dispersion of ∼83 km s⁻¹, we obtain a corrected velocity dispersion of ∼152 km s⁻¹. Using the scaling relation from Gültekin et al. (2009), we infer M_BH ∼ 4 × 10⁷ M_⊙, which is significantly larger than the estimate obtained from the stellar mass. Given the good agreement of the logM_* ∼ 9 estimate from both spectroscopy and photometry fitting, we adopt the M_BH estimate from stellar mass as reference for eRO-QPE3 at this stage (thus in the range M_BH ∼ (0.9 − 5.3)×10⁶ M_⊙). Further analysis with spectra obtained at higher resolution is needed for a more robust narrow-line classification and black hole mass estimate from velocity dispersion.

2.3.2. eRO-QPE4

eRO-QPE4 was associated with the galaxy 2MASS 04453380−1012047 at (RA, Dec) = (04:45:33.80, −10:12:04.74). We show the 1σ eROSITA positional accuracy (based on the cumulative eRASS:4 image) in the bottom left panel of Fig. 3 with a red circle, whilst in green the more accurate XMM-Newton circle. We took an optical spectrum with RSS on SALT on February 23, 2023 (bottom right panel of Fig. 3), data were reduced as reported in A21. Similarly to the other eRO-QPEs, the galaxy does not show prominent emission lines and a redshift of 0.0437 is estimated from Calcium absorption lines (bottom right panel of Fig. 3). We fit the salt spectrum with Firefly, after renormalizing it using optical photometry from VEXAS DR2 catalogs (Khramtsov et al. 2021), which are also consistent with SkyMapper petrosian magnitudes (Wolf et al. 2018). We refer to the eRO-QPE3 fit for more details on the stellar population models and the related systematic uncertainties. The stellar mass inferred with Firefly for eRO-QPE4 is $M_{\ast }=1.6_{-0.6}^{+0.7}\times {10}^{10}{M}_{\odot }$, whist SFR is $2.{26}_{-2.21}^{+1.06}{M}_{\odot }$ yr⁻¹. Using Reines & Volonteri (2015), we obtain $M_{\mathrm{BH}}=6.8_{-3.2}^{+4.8}\times {10}^7{M}_{\odot }$. This is the highest M_* estimate inferred from any of the eROSITA QPEs (A21) and the highest black hole mass ever inferred for a QPE source (Wevers et al. 2022). With the GRAHSP SED fitting we obtain $M_{\ast }=0.{60}_{-0.12}^{+0.29}\times {10}^{10}{M}_{\odot }$, which is slightly lower than the Firefly estimate, even including uncertainties, but still much higher than what is estimated for eRO-QPE3 with the same method. This corresponds (Reines & Volonteri 2015) to $M_{\mathrm{BH}}=1.{74}_{-0.47}^{+1.28}\times {10}^7{M}_{\odot }$. From the SALT spectrum we can only identify faint narrow Hα, [O II] and [O III] λ5007 in emission, while no Hβ nor [N II] lines are visible, even from the residual spectrum provided by Firefly after subtracting the continuum. Hence, we are unable to securely classify the galaxy at this stage, although in similarity with other eROSITA QPEs a strong preexisting AGN is disfavored. We obtained a corrected velocity dispersion of ∼133 km s⁻¹ with pPXF. Using the scaling relation from Gültekin et al. (2009), we infer M_BH ∼ 2.4 × 10⁷ M_⊙, which is reasonably consistent with that obtained via the M_* estimate. Hence, for eRO-QPE4 all diagnostics used point to M_BH in the range ∼(1.7 − 6.8)×10⁷ M_⊙.

Fig. 4. XMM-Newton 0.2 − 10.0 keV EPICpn light curve of eRO-QPE3, of the first (eRO3-XMM1, left) and second observation (eRO3-XMM2, right). The two eruptions in eRO3-XMM1 are separated by ∼20.4 h. Here, teRO3-XMM1, 0 corresponds to MJD ∼59779.799 and teRO3-XMM2, 0 to ∼59799.745.

2.4. XMM-Newton

Two XMM-Newton observations were performed on eRO-QPE3, namely ObsID 0883770101 starting on 19 July 2022 (hereafter eRO3-XMM1) and 0883770701 starting on 8 August 2022 (hereafter eRO3-XMM2). A single XMM-Newton observation was taken on eRO-QPE4, namely ObsID 0883770401 starting on 10 March 2023 (hereafter eRO4-XMM). We reduced XMM-Newton data of EPIC-MOS1 and 2 (Turner et al. 2001) and EPIC-PN (Strüder et al. 2001) cameras and the Optical Monitor (OM; Mason et al. 2001) using standard tools and prescriptions (SAS v. 20.0.0 and HEAsoft v. 6.29). Event files from EPIC cameras were screened for flaring particle background. For eRO-QPE3, source (background) regions were extracted within a circle of 30″ centered on the source (in a nearby source-free region). For eRO-QPE4, apertures of 40″ were used. The X-ray position obtained from eROSITA during the survey was refined using XMM-Newton data and the task eposcorr. We cross-correlated the X-ray sources in the EPIC-PN image with external optical and infrared catalogs (available through the Processing Pipeline Subsystem). The counterparts of eRO-QPE3 and eRO-QPE4 were excluded to obtain a more unbiased estimate of the possible offset from the nucleus. The resulting 1σ positional circle from XMM-Newton is shown in green in Fig. 3 for both eRO-QPE3 (accuracy of 1.5″) and eRO-QPE4 (accuracy of 4.9″). For both sources, the refined X-ray position is consistent with being nuclear within the current uncertainties.

2.4.1. eRO-QPE3

In eRO3-XMM1 two eruptions are observed, with faster rise than decay, separated by ∼20 h, while in eRO3-XMM2 only one burst was detected (Fig. 4). We use the burst model described in Arcodia et al. (2022), since it has shown itself to be successful in parametrizing asymmetric QPE light curves:

$$\begin{array}{c}\hfill \{\begin{array}{cc}A\lambda e^{{\tau }_1/(t_{\mathrm{peak}}-t_{\mathrm{as}}-t)}\hfill & \mathrm{if}t<t_{\mathrm{peak}}\hfill \\ A{e}^{-(t-t_{\mathrm{peak}})/{\tau }_2}\hfill & \mathrm{if}t\geqq t_{\mathrm{peak}}\hfill \end{array}\end{array}$$(1)

which is evaluated at zero for times smaller than the asymptote at t_peak − t_as, where $t_{\mathrm{as}}=\sqrt{{\tau }_1{\tau }_2}$. τ₁ and τ₂ are the characteristic timescales, although only the latter is directly related to the decay timescale. A is the amplitude at the peak and λ = e^t_λ a normalization, where $t_{\lambda }=\sqrt{{\tau }_1/{\tau }_2}$. Rise and decay times can be defined as a function of 1/eⁿ factors with respect to the peak flux, with n being an integer. Rise and decay are obtained from the width and asymmetry factors (w and k), where w(n)=τ₂ n − τ₁/(t_λ n)+t_as and k(n)=[τ₂ n + τ₁/(t_λ n)−t_as]/w. We define rise and decay timescales as τ_rise = w (1 − k)/2 and τ_decay = w (1 + k)/2 using n = 1. We quote median values and related 16th and 84th percentiles. Referring to the three eruptions observed by XMM-Newton with time, the rise is ${2908}_{-435}^{+95}$ s, ${2522}_{-192}^{+290}$ s, ${2903}_{-175}^{+219}$ s, respectively, and the decay is ${5896}_{-308}^{+681}$ s, ${5194}_{-357}^{+672}$ s, ${5015}_{-512}^{+395}$ s, respectively. The two eruptions in eRO3-XMM1 are separated by a median recurrence time (and related 16th, 84th percentiles) $20.{40}_{-0.11}^{+0.25}$ h. We show an example of this model and its application to QPEs in Fig. A.9, for the first burst of observation eRO3-XMM1. Given the relatively low signal-to-noise in eRO-QPE3, we test the QPEs’ energy dependence dividing the light curve in two energy bins only, namely 0.2 − 0.6 keV and 0.8 − 2.0 keV. Adopting the first eruption of eRO3-XMM1 as representative, we obtain a rise (decay) of ${3661}_{-244}^{+202}$ s (${6020}_{-342}^{+445}$ s) in the low-energy bin and of ${2237}_{-398}^{+895}$ s (${4946}_{-1234}^{+1288}$ s) in the high-energy bin. Furthermore, we obtain that the peak time of the low-energy light curve is (13.3 ± 0.3) ks after the start of eRO3-XMM1, while it is $12.0_{-0.5}^{+0.6}$ ks for the high-energy light curve. Therefore, the eruptions of eRO-QPE3 follow the known energy-dependence (M19; G20; A21; Chakraborty et al. 2021; Arcodia et al. 2022), in that they are wider and start later at lower energies.

Fig. 5. Source plus background spectrum of eRO-QPE3 in quiescence (dark gray, background dominated) and at the peak of the first burst in the observation eRO3-XMM1 (orange, following the color-coding of Fig. A.9). The light-blue line and contour represent the background spectrum and its uncertainties. Source plus background data points are shown to put the source-only models (lines and contours with the same colors) into context. Darker model contours are 1σ percentiles, lighter ones are 3σ.

The quiescent state of eRO-QPE3 is too faint to be securely detected above background. We adopt a disk spectrum to provide a flux upper limit of F_{0.2 − 2.0 keV} <  1.1 × 10⁻¹⁴ erg s⁻¹ cm⁻² (Fig. 5). The fit posterior on the flux, or disk normalization, is unconstrained, while the disk temperature is loosely constrained to $k{T}_{\mathrm{disk}}={46}_{-18}^{+21}$ eV. Despite the weak constraint on the temperature, we conservatively consider the quiescence of eRO-QPE3 to be undetected at the XMM-Newton epoch with a soft X-ray luminosity upper limit of L_{0.2 − 2.0 keV} <  1.6 × 10⁴⁰ erg s⁻¹. We divided the first eruption of eRO3-XMM1, as example, in five epochs, namely rise1, rise2, peak, decay1 and decay2 (see Fig. A.9). We fit each epoch with a black body (see more details in Appendix A.2) and report best-fit parameters in Table A.1. We present more details on the spectral evolution during the eruptions in Sect. 3. Here, we quote results of the peak spectrum (see Fig. 5): we obtained $k{T}_{\mathrm{QPE}}={111}_{-5}^{+6}$ eV and $F_{0.2-2.0\mathrm{keV}}=3.1_{-0.2}^{+0.3}\times {10}^{-13}$ erg s⁻¹ cm⁻² (L_{0.2 − 2.0 keV} = (4.2 ± 0.4)×10⁴¹ erg s⁻¹), thus ≳25 − 30 times brighter than the quiescence flux in the soft X-ray band. The bolometric luminosity of the black body component is L_QPE, bol = (4.9 ± 0.5)×10⁴¹ erg s⁻¹.

Fig. 6. XMM-Newton 0.2 − 10.0 keV EPICpn light curve of eRO-QPE4. The three eruptions are separated by ∼9.82 h and ∼14.70 h. Here, teRO4-XMM, 0 corresponds to MJD ∼60013.127.

The OM UVW1 (U) filter was used throughout eRO3-XMM1 (eRO3-XMM2) with a series of ∼4400 s-long exposures. The subpanels of Fig. 4 show the OM light curves in the respective filters. Similarly to other QPE sources (M19; G20; A21), no simultaneous variability is observed, with the caveat that these data are likely contaminated by the host galaxy and possibly nuclear star cluster. The mean flux is λF_{λ, 344 nm} = (2.2 ± 0.5)×10⁻¹³ erg s⁻¹ cm⁻² and λF_{λ, 291 nm} = (1.0 ± 0.5)×10⁻¹³ erg s⁻¹ cm⁻² in the U (344 nm) and UVW1 (291 nm) filter, respectively, with magnitudes of 19.44 ± 0.09 and 19.87 ± 0.11, respectively.

2.4.2. eRO-QPE4

In eRO4-XMM, three eruptions are observed (Fig. 6). Their duration appears shorter compared to eRO-QPE3 and they appear more symmetric. We fit each eruption in eRO4-XMM with both a Gaussian model and the asymmetric model of Eq. (1). The latter model performs significantly worse (in terms of residuals and goodness of fit) and we adopt the Gaussian model for the eruptions in eRO-QPE4. We show an example in Fig. A.9, for the second burst of observation eRO4-XMM. The median (with related 16th, 84th percentile values) of the Gaussian width is ${1678}_{-90}^{+83}$ s, ${1765}_{-23}^{+24}$ s and ${2081}_{-29}^{+30}$ s, for the three eruptions, respectively. The peak-to-peak separation between the QPEs is $9.{82}_{-0.02}^{+0.03}$ h and $14.{70}_{-0.01}^{+0.02}$ h. We divided the light curve in energy bins of 0.2 − 0.4 keV, 0.4 − 0.6 keV, 0.6 − 0.8 keV, 0.8 − 1.0 keV and 1.0 − 2.0 keV (hereafter E₁ to E₅, respectively). Adopting the second eruption in eRO4-XMM as example, given its higher signal-to-noise, we obtain that the median (with related 16th, 84th percentile values) of the Gaussian width is (1759 ± 22) s, ${1630}_{-27}^{+25}$ s, ${1417}_{-32}^{+33}$ s, ${1292}_{-52}^{+51}$ s and ${1299}_{-116}^{+120}$ s, respectively, hence decreasing from E₁ to E₅. Adopting the E₅ Gaussian peak as reference (which occurs $10.{28}_{-0.03}^{+0.04}$ h after the start of eRO4-XMM), we obtain delays of ${603}_{-30}^{+31}$ s, ${464}_{-30}^{+35}$ s, ${165}_{-40}^{+43}$ s and $-{30}_{-70}^{+67}$ s, for E₁, E₂, E₃ and E₄, respectively. Apart from E₄ which is consistent with the properties of E₅ within uncertainties, the eruptions of eRO-QPE4 follow the known energy-dependence (M19; G20; A21; Chakraborty et al. 2021; Arcodia et al. 2022), in that they are wider and start later at lower energies.

Fig. 7. Spectrum of eRO-QPE4 in quiescence (dark gray) and at the peak of the second burst in the observation eRO4-XMM (light green, following the color-coding of Fig. A.9). The dark-red line and contour represent the background spectrum and its uncertainties. Source plus background data points are shown to put the source-only models (lines and contours with the same colors) into context. Darker model contours are 1σ percentiles, lighter ones are 3σ. Individual model components are labeled.

The quiescent state of eRO-QPE4 is detected with a soft spectrum (Fig. 7) similar to that of other QPE sources in quiescence (M19; G20; A21; Chakraborty et al. 2021). We fit the quiescence with a disk model and subsequently add a Comptonization component (Fig. 7), which yields a disk temperature of kT_disk = (43 ± 2) eV with a flux F_{0.2 − 2.0 keV} = (3.9 ± 0.4)×10⁻¹³ erg s⁻¹ cm⁻². This corresponds to L_{0.2 − 2.0 keV} = (1.7 ± 0.3)×10⁴² erg s⁻¹ for the disk component alone, or to a bolometric luminosity of $L_{\mathrm{disk},\mathrm{bol}}=2.0_{-0.3}^{+0.6}\times {10}^{43}$ erg s⁻¹. The Comptonization component is much fainter, with $F_{0.2-2.0\mathrm{keV}}=1.4_{-0.4}^{+0.8}\times {10}^{-14}$ erg s⁻¹ cm⁻², and its slope is unconstrained. More details are presented in Table A.2. As done for eRO-QPE3, the QPE flares of eRO-QPE4 were separated in five epochs, namely rise1, rise2, peak, decay1 and decay2 (bottom panel of Fig. A.9 for the second burst of the observation eRO4-XMM). The quiescence model was held fixed by imposing its parameters to vary only within the 10th–90th percentile interval of the posteriors obtained during the quiescence fit alone. The fit results for the different burst phases are reported in Table A.2. For the epochs rise1, rise2, decay1 and decay2 the best-fit model adopted for the burst is that of the quiescence plus a black body component. For the peak spectrum we add a Comptonization component (see Fig. 7). More details are presented in Table A.2. This model yields kT_QPE = (123 ± 2) eV and F_{0.2 − 2.0 keV} = (2.9 ± 0.1)×10⁻¹² erg s⁻¹ cm⁻² (L_{0.2 − 2.0 keV} = (1.27 ± 0.03)×10⁴³ erg s⁻¹). The bolometric luminosity of the black body component is L_QPE, bol = (1.44 ± 0.04)×10⁴³ erg s⁻¹. Similarly to eRO-QPE3, the OM UVW1 filter was used throughout and no simultaneous variability is observed. The mean flux is λF_{λ, 291 nm} = (2.1 ± 0.5)×10⁻¹³ erg s⁻¹ cm⁻² and the mean magnitude 18.96 ± 0.15.

2.5. NICER

NICER data were processed using HEAsoft v6.32.1, using NICERDAS v11a. We grouped the data into 200-s Good-Time Intervals (GTIs) and adopted custom filtering choices of unrestricted undershoot (underonly_range=*-*) and overshoot rates (overonly_range=*-*), with per-FPM and per-MPU autoscreening disabled to prevent aggressive event filtering. We manually discarded focal plane modules (FPMs) with 0 − 0.2 keV rates or 5 − 18 keV count rates > 3σ higher than average, or above an absolute threshold of 5 counts s⁻¹. As both eRO-QPE3 and eRO-QPE4 are super-soft and faint, the emission above 5 keV is entirely background-dominated, whereas the 0 − 0.2 keV band is undershoot-dominated, making these bands effective proxies for assessing the background conditions.

Following the methodology of Chakraborty et al. (2024) we create the light curves in Figs. 8 and 9 with time-resolved spectroscopy across the GTIs. After screening the event lists, we used the SCORPEON ⁴ background model to estimate the contribution from the diffuse X-ray background and non X-ray background. We fit the entire broadband (0.2 − 18 keV) array counts with PyXspec⁵. Along with the SCORPEON background we fit each GTI with a source model represented by tbabs × zbbody. We consider a source detection as any GTI in which the blackbody normalization is > 1σ inconsistent with zero, namely a non-background component is required by the fit at the 1σ level. The source count rates are the summed counts contained only in the blackbody component. We consider all other cases as non detections and we adopt as upper limit the 3σ upper error uncertainty on the flux of the blackbody model that did not make the significance cut.

Fig. 8. NICER background-subtracted 0.2 − 2.0 keV light curve of eRO-QPE3, in the epoch eRO3-NICER1 (top, starting at MJD 59697.644) and eRO3-NICER2 (bottom, starting at MJD 59822.782). Black points represent secure detections of a source component on top of background emission, contrary to the epochs shown in red which are consistent with background-only. Vertical dashed lines represent evenly spaced recurrences of 20 h, as a guide for the eye.

Fig. 9. As in Fig. 8 but for eRO-QPE4, starting at MJD 59971.637. Vertical dotted lines represent approximate locations for the eruptions’ peak (thus they are not evenly spaced), as a guide for the eye.

2.6. eRO-QPE3

eRO-QPE3 was observed for a total of ∼26.2 ks between 28 April 2022 and 3 May 2022, and again for a total of ∼40.3 ks between 31 August 2022 and 11 September 2022, namely a few months before and about a month after eRO3-XMM1 (Fig. 8). These observation sets are hereafter named eRO3-NICER1 and eRO3-NICER2, respectively. In both sets, the source is significantly detected in a transient fashion, although in most cases the full eruption is not resolved, at times due to nonoptimal background conditions preventing a secure detection of a source component.

Only two clear bright eruptions are seen in eRO3-NICER1 (Fig. 8, top), separated by roughly ∼17 h, which is significantly shorter than the ∼20.4 h recurrence seen in eRO3-XMM1 (Fig. 4). eRO-QPE3 is then significantly detected above background for a few short exposures in eRO3-NICER1. As these transient detections are roughly separated by ∼20 h (see vertical dashed lines in Fig. 8, as guide for the eye), we are confident that they all correspond to eruption phases. Although current data do not allow for a more precise constraint on the exact recurrence times, we can conclude that in eRO-QPE3 QPEs have a typical recurrence in the range ∼17 − 20 h, with significant scatter as seen in other sources (e.g., eRO-QPE1; A21). The brightest eruptions in eRO3-NICER1 reach a maximum flux of $F_{0.2-2.0\mathrm{keV}}=9.8_{-1.3}^{+1.2}\times {10}^{-13}$ erg s⁻¹ cm⁻² and a temperature $kT={133}_{-8}^{+10}$ eV, whilst the average of all bursts including the other transient detections is $F_{0.2-2.0\mathrm{keV}}=6.5_{-3.1}^{+2.3}\times {10}^{-13}$ erg s⁻¹ cm⁻² and $kT={107}_{-20}^{+24}$ eV. The maximum values are larger than those observed in eRO3-XMM1, perhaps indicating a long-term decay of peak flux and temperature. In eRO3-NICER2, only part of a bright eruption is observed (Fig. 8, bottom), with a peak flux and temperature roughly consistent with those of the maximum values observed in eRO3-NICER1. At several other epochs the source is found with a flux of F_{0.2 − 2.0 keV} ∼ 5 × 10⁻¹³ erg s⁻¹ cm⁻², similar to the peak flux observed in eRO3-XMM1. This would suggest that most QPEs peak at a level consistent with eRO3-XMM1, but occasional brighter eruptions reach a brighter flux, as that seen at the maximum of eRO3-NICER1 and eRO3-NICER2. The repetition is unclear and not as close to ∼20 h as in eRO3-NICER1, although in most snapshots we were not able to securely assess the presence of a source component. However, these lower-flux detections close to the background count rate are much harder to be securely disentangled from the latter, compared to brighter eruptions. Hence, we refrain from further interpretations of eRO3-NICER2 at this stage and defer more thorough analysis to future work.

2.7. eRO-QPE4

eRO-QPE4 was observed for a total of ∼66 ks between 27 January 2023 and 2 February 2023 (Fig. 9), hereafter named eRO4-NICER. The source is significantly detected in a transient fashion, although in most cases the full eruption is not resolved due to the short duration of eruptions (≈30 min of FWHM, based on XMM-Newton, Fig. 6) compared to the typical gap in NICER’s monitoring. In Fig. 9, vertical dotted lines represent approximate locations for the eruptions, as a guide for the eye. Inferred separations are in the range ∼(11 − 15.5) h and appear somewhat regular (considering that an eruption was likely lost in the gap around MJD_{eRO4 − NICER} + 1.5 in Fig. 9). However, they do not come in a clear alternating pattern such as in GSN 069 and eRO-QPE2 (M19; A21; Arcodia et al. 2022), nor do they show very large scatter such as the eruptions in eRO-QPE1 (A21; Chakraborty et al. 2024) or RX J1301.9+2747 (G20; Giustini et al., in prep.), but rather show an intermediate behavior. The peak flux and temperature range from $F_{0.2-2.0\mathrm{keV}}=1.6_{-0.1}^{+0.2}\times {10}^{-12}$ erg s⁻¹ cm⁻² and $kT={121}_{-7}^{+6}$ eV for the lowest peaks, to F_{0.2 − 2.0 keV} = (2.4 ± 0.2)×10⁻¹² erg s⁻¹ cm⁻² and kT = (118 ± 5) eV. This is reasonably consistent with XMM-Newton (Fig. 6 and Table A.2), given that not all NICER eruptions are resolved and some peaks have been missed.

2.8. Swift

Two serendipitous Swift-XRT observations were taken with eRO-QPE3 in the field of view, on 2020-08-25 and 2020-12-26, for an exposure of ∼584 s and ∼767 s, respectively. Aperture photometry was performed at the location of eRO-QPE3 with the Python package photutils. The source aperture adopted was 40″, while background counts were extracted in an annulus with inner and outer radius of 40″ and 40″, respectively. Count rates were converted to fluxes using WebPIMMS and adopting the best-fit models of the closest eROSITA observation. eRO-QPE3 was detected during the first observation, close to eRASS2, at about the same flux level of the eRASS2 bright state, perhaps indicative of a serendipitous QPE observation (see Sect. 4.1). The second observation did not detect the source, at the same depth of the eRASS3 upper limit of the faint phase. eRO-QPE3 was in the UVOT field of view only in the second observations, for a total exposure of ∼390 s with the UVW1 filter at 268 nm. Performing aperture photometry with the task uvotsource, eRO-QPE3 is only detected at 2.6σ, hence it is formally undetected. The corresponding flux upper limit is < 1.3 × 10⁻¹³ erg s⁻¹ cm⁻², so consistent with the OM-UVW1 detection of ∼1.0 × 10⁻¹³ erg s⁻¹ cm⁻² at 291 nm.

We triggered Swift ToO follow up observations of eRO-QPE4 after the eRASS4 detection. eRO-QPE4 was also serendipitously observed with Swift before the eRASS4 discovery. We used the XRT online data analysis tool⁶ (Evans et al. 2009) to check whether eRO-QPE4 was detected for each observation. The tool was also used to extract the X-ray spectra for observations in which it was detected, and to estimate the 3σ count rate upper limits for non-detections. The X-ray spectra were rebinned to have at least one count in each bin. An absorbed accretion disk model (diskbb) model was used to calculate the flux and upper limits. Most observations are non detections, shallower than the eROSITA and XMM-Newton fluxes, hence uninformative. In between the eRASS4 detection and eRO4-XMM, there are a few Swift-XRT detections in the range F_{0.2 − 2.0 keV} ∼ (1 − 4)×10⁻¹² erg s⁻¹ cm⁻², which is compatible with the eruptions fluxes observed in eRASS4 and with XMM-Newton, thus supporting a serendipitous detection during an eruption. eRO-QPE4 was not in the UVOT FoV of the archival observations. We thus only analyzed the UVOT data for the ToO observations taken after eRASS4. We used the UVOT analysis pipeline provided in HEASOFT software (version 6.31) with UVOT calibration version 20201215 to analyze the UVOT data. Source counts were extracted from a circular region with a radius of 5″ centered at the source position. A 20″ radius circle from a source-free region close to the position of eRO-QPE4 was chosen as the background region. The task uvotsource was used extract the photometry. The UVOT/UVW1 (UVOT/UVW2) flux is fairly consistent among the several snapshots, with an average flux at ∼2.8 × 10⁻¹³ erg s⁻¹ cm⁻² (∼2.4 × 10⁻¹³ erg s⁻¹ cm⁻²).

2.9. Radio observations

We observed the coordinates of eRO-QPE3 and eRO-QPE4 on two occasions each with the Australia Telescope Compact Array (ATCA) between October 2022 and August 2023 (proposal ID C3513/C3527). All radio data were reduced using standard procedures in the Common Astronomy Software Application (CASA v5.6.3, CASA Team 2022) including flux and bandpass calibration with PKS 1934−638 and phase calibration with PKS 1406−267 (eRO-QPE3), PKS 0458−020 (eRO-QPE4 4 cm), and PKS 0420−014 (eRO-QPE4 16 cm). For eRO-QPE3 we observed the target on two occasions separated by approximately 8 months with the ATCA dual 4 cm receiver and the CABB correlator producing 2 × 2 GHz of bandwidth split into 2048 × 1 MHz channels centered on 5.5 GHz and 9 GHz. For eRO-QPE4 we initially observed at 5.5 and 9 GHz with the dual 4 cm receiver and then followed up 3 months later with an ATCA observation with the 16 cm receiver and the CABB correlator producing 2 GHz of bandwidth split into 2048 × 1 MHz channels centered on 2.1 GHz. Images of the target field were made with the CASA task tclean. No radio emission was detected at the coordinates of either QPE source in any of the observations. A summary of the ATCA observations including the 3σ flux density upper limits is given in Table 1. The radio non-detections of eRO-QPE3 and eRO-QPE4 indicate it is unlikely the host galaxies have strong AGN activity in the nuclei, in agreement with the optical proxies.

Fig. 10. Spectral evolution of QPEs in eRO-QPE3 (left) and eRO-QPE4 (right). The top-right (bottom) subpanels show the luminosity (temperature) evolution with time of the QPE component. Darker to lighter colors represent the evolution from start to end of the eruptions. The first eruptions of the eRO3-XMM1 observation and the second of the eRO4-XMM observation are taken as reference for the QPEs in eRO-QPE3 and eRO-QPE4, respectively. The top-left subpanels show the luminosity-temperature coevolution, analogous to the rate-hardness evolution in Arcodia et al. (2022) which first reported the hysteresis cycle.

3. Spectral evolution of QPEs

Given the growing number of repeating X-ray transients, it is useful to identify common properties that indicate a common emission process, or suggest otherwise. For instance, Arcodia et al. (2022) first reported the specific spectral evolution of QPEs in eRO-QPE1 as a soft X-ray spectrum which, for the same count rate or flux, is harder during the rise than during the decay. Moreover, the peak hardness of the spectrum is reached before the peak flux of the eruption. The same behavior was identified in GSN 069 (Miniutti et al. 2023a) and RX J1301.9+2747 (Giustini et al., in prep.). Identifying this hysteresis cycle in a repeating X-ray transient is the smoking gun of the same emission process. This is particularly important since the timing properties of the known QPEs are quite diverse (M19; G20; A21; Arcodia et al. 2022), and can even significantly change in the same source (Miniutti et al. 2023b,a).

We show the spectral evolution of eRO-QPE3 and eRO-QPE4 in Fig. 10. We choose model-dependent quantities such as luminosity and temperature, compared to pure data-driven quantities such as hardness-ratio and count rate. This is because, at least in eRO-QPE4, the quiescence component needs to be decomposed (e.g., as done for GSN 069, Miniutti et al. 2023a). This was not required for eRO-QPE1 (Arcodia et al. 2022) as the source is not detected between the eruptions. Therefore, Fig. 10 assumes that QPEs can be modeled by a black body. Results would be identical with any model with an exponential decay due to a characteristic temperature, for instance Bremsstrahlung or an accretion disk. In Fig. 10, the first eruption of eRO3-XMM1 and the second of eRO4-XMM are used as representative, although other eruptions are shown in Fig. A.11. The eruption evolves from darker to lighter colors of the respective color maps, as shown in Fig. A.9. The eruptions were divided in 5 epochs with the goal to compare rise and decay at similar count rate or luminosity. As there is not a one-to-one relation between count rate and luminosity, this is not necessarily achieved. As a matter of fact, as much as the L_bol and kT values would change slightly depending on the definition of the epochs, the hysteresis behavior in the L_bol − kT plane (or count rate versus hardness ratio) would remain evident. This is confirmed by the same pattern in other eruptions of both eRO-QPE3 and eRO-QPE4 (Fig. A.11) and other QPE sources (Arcodia et al. 2022; Miniutti et al. 2023a). We also note that, similarly to eRO-QPE1 (Arcodia et al. 2022) and GSN 069 (Miniutti et al. 2023a), the peak temperature is reached before the peak luminosity.

As done in Miniutti et al. (2023a), we estimate the size of this emitting region, under the assumption that it is indeed black body emission. We assume that the observer sees half of a radiating spherical surface. For eRO-QPE3 (Fig. 10, left) we obtain radii of ∼1.2 × 10¹⁰ cm², ∼1.2 × 10¹⁰ cm², ∼1.7 × 10¹⁰ cm², ∼1.6 × 10¹⁰ cm² and ∼1.8 × 10¹⁰ cm² for rise1, rise2, peak, decay1 and decay2, respectively. Therefore, we infer an increase in size of a factor ∼1.5 from start to end of the eruptions, lower than what was inferred for GSN 069 (Miniutti et al. 2023a). For eRO-QPE4 (Fig. 10, right), we obtain radii of ∼2.4 × 10¹⁰ cm², ∼7.0 × 10¹⁰ cm², ∼9.5 × 10¹⁰ cm², ∼17.3 × 10¹⁰ cm² and ∼11.6 × 10¹⁰ cm² for rise1, rise2, peak, decay1 and decay2, respectively. Therefore, we infer an increase in size of a factor ∼7 from rise1 to decay1 and of a factor ∼5 from rise1 to decay2, which is larger than what was inferred for GSN 069 (Miniutti et al. 2023a). Regardless of the fine details which are deferred to future homogeneous work extended to all QPE sources, the common thread is that during the L_bol − kT hysteresis the size of the emitting region increases, which is qualitatively consistent with the latest models of disk-orbiter collisions (Franchini et al. 2023; Linial & Metzger 2023; Tagawa & Haiman 2023). The absolute numbers are, of course, model-dependent and any departure from this spectral model would likely change the absolute values of the size, but not the qualitative evolution during the eruptions, provided the X-ray spectrum is indeed the exponential tail of a thermal spectrum.

Fig. 11. Long-term evolution of eRO-QPE3, separating flux states between quiescence and eruption (which includes the former). eROSITA points are shown as circles, squares for XMM-Netwon, diamonds for NICER (the brightest eruptions in the light curves of Fig. 8) and stars for Swift-XRT. All uncertainties are 3σ. Red contours on orange points highlight observations in which the identification of an eruption is less secure. White symbols for eRASS1 and eRASS2 correspond to the flux of the full observation without separating different states. Horizontal lines highlight ROSAT and XMM-Newton archival upper limits, as stated in the legend.

4. Long-term evolution

4.1. eRO-QPE3: Quiescent flux decays until disappearance

In Fig. 11, we show the long-term evolution of eRO-QPE3 inferred with eROSITA, XMM-Newton, NICER and Swift-XRT. We separate flux states between quiescence (dark red data points) and eruption (orange). This separation is straight-forward in XMM-Newton and NICER data and at least in eRASS3 and eRASS4 data (see Fig. 1). In eRASS1 and eRASS2, since the quiescence is detected and given eROSITA’s sampling, the identification of the eruption state is informed by the presence of QPEs at later times. However, formally one cannot exclude high-amplitude variability of the quiescent state with no eruptions. This potential no-QPE flux state is shown in white with dark red contours and the related possible eruption states are shown with red contours on orange points. For eRASS5, we consider the eruption state as ambiguous since the variability is not as significant and the flux chain has a tail extending to faint fluxes, therefore it is formally not well constrained at 3σ (hence the large uncertainty in Fig. 11). The Swift-XRT detection is also ambiguous: its flux is as bright as the eruption state in eRASS2, but no variability can be inferred from the few hundred seconds of exposure to unambiguously identify it as an eruption. Despite these caveats, throughout the rest of the discussion we assume that all the orange points caught eRO-QPE3 in eruption, and the dark red points represent the flux of the intra-QPEs’ quiescence.

Quite interestingly, eRO-QPE3 showed a disappearing quiescence component (Fig. 1 and Table A.1), being detected in eRASS1 and eRASS2 and never since, including Swift-XRT and eRASS3 observations taken within 6 months. This is consistent with what has been observed for GSN 069 (Shu et al. 2018; Miniutti et al. 2023a), albeit over a decade and with a much better data coverage before the QPEs’ discovery. Furthermore, the two QPE candidates presented in Chakraborty et al. (2021) and Quintin et al. (2023) also showed evidence of eruptions on top of a decaying baseline flux. In general, a possible interpretation of the decaying quiescence in GSN 069 and the two candidates is that the accretion flow was induced by a TDE. The presence of a precursor TDE is secure in AT2019vcb (“Tormund”; Quintin et al. 2023), but since only the possible start of an X-ray eruption was observed, the association between QPEs and TDEs is still open. Were this connection to be unambiguously proven, one can interpret the long-term evolution of eRO-QPE3 in a similar way. However, the available data (only two epochs are detected in eRO-QPE3) do not allow us to achieve a precise constraint or prediction on this long-term evolution. We can test whether the spectral evolution between the disk component in eRASS1 and eRASS2 is consistent with a cooling thin disk. The disk temperature decreased from ∼100 eV to ∼89 eV from eRASS1 to eRASS2 (Table A.1). Approximating the disk emission as a black body with constant area with the temperature of the inner radius, the eRASS2 bolometric emission would be ∼(89/100)⁴ ∼ 0.6 times that of eRASS1. The disk bolometric emission of eRASS1 and eRASS2, estimated from the spectral fits, is $L_{\mathrm{bol}}=4.1_{-2.3}^{+1.6}\times {10}^{42}$ erg s⁻¹ and $L_{\mathrm{bol}}\sim 8.4_{-8.3}^{+4.7}\times {10}^{41}$ erg s⁻¹, respectively. This corresponds to a decrease to ≈20% of the eRASS1 luminosity, although given the large 1σ uncertainties of the eRASS2 value, a 60% decrease is still compatible. Regardless of the exact type and nature of the observed decay, eRO-QPE3 is the first QPE source in which the intra-QPE quiescence is detected at some earlier phases and never again. Other secure QPE sources have either always shown it, or never (M19; G20; A21). Hence, eRO-QPE3 would bridge the gap between QPE sources that always showed a (possibly time-evolving) quiescence spectrum in between QPEs (e.g., GSN 069, Miniutti et al. 2023a) and QPE sources which never showed a quiescence spectrum (e.g., eRO-QPE1, A21). This suggests that the absence of intra-QPE quiescence in some QPE sources is merely observational, as it might have faded compared to an earlier epoch.

Furthermore, we also show in Fig. 11 3σ archival upper limits from ROSAT (from 1991) and the XMM-Newton Slew Survey (from 2000 and 2013; Saxton et al. 2008). Comparing these with the eRASS1 flux (both if the quiescence includes or not the putative eruption), we notice that ROSAT would have detected the source if it were as bright as eRASS1. Therefore, the accretion flow of eRO-QPE3 likely became brighter (e.g., perhaps radiatively efficient) in the soft X-rays some time between 1991 and 2020.

Fig. 12. Long-term evolution of eRO-QPE4, separating flux states between quiescence and eruption. eROSITA points are shown as circles, squares for XMM-Netwon, diamonds for NICER and stars for Swift-XRT. Swift-XRT and NICER eruptions points are considered such based on their flux level similar to the range shown by XMM-Netwon. All uncertainties are 3σ. The lower-panel shows Swift-UVOT data with the UVW1 (light blue) and UVW2 filter (azure).

Despite XMM-Newton data and most NICER data (unless there are large gaps in the monitoring) allow the identification of the peak, the eROSITA and Swift-XRT data points would catch the source in an unknown part of the eruptions. Therefore, all data points except XMM-Newton and NICER are lower limits to the peak. Without NICER monitoring covering several burst cycles, one might have concluded that QPEs must have significantly faded over only a few months. However, NICER data taken both before and after eRO3-XMM1 and eRO3-XMM2 (which showed F_{0.2 − 2.0 keV} ∼ 3.1 × 10⁻¹³ erg s⁻¹ cm⁻²) unveil the presence of some eruptions significantly brighter (F_{0.2 − 2.0 keV} ∼ 1.0 × 10⁻¹² erg s⁻¹ cm⁻²) than the others, which are instead compatible with the eRO3-XMM1 peak flux. This indicates that whilst on average the peak flux might still be overall decreasing, current data also unveil the presence of significant diversity in amplitudes, as seen in other QPE sources (G20; A21). We show the flux of these brightest eruptions in Fig. 8. Finally, we note that no significant optical variability is apparent within the available public datasets from optical sky monitors (e.g., ASAS-SN, Shappee et al. 2014) covering the epochs shown in Fig. 11.

4.2. eRO-QPE4: Quiescent flux appears or brightens after the eruptions

Figure 12 shows the long-term evolution of the X-ray and UV flux of eRO-QPE4. The UV flux is overall constant within the available exposures. On the contrary, the X-ray quiescence flux (dark gray points) is only constrained by XMM-Newton with a detection. However, it is interesting to note that both the archival ROSAT (dotted line) and eROSITA 3σ upper limits suggest that the quiescent accretion disk must have been much fainter or absent at the ROSAT/eROSITA epochs. Current data surely rule out the presence of a radiatively efficient accretion flow, as bright as that seen by XMM-Newton, prior to the eRASS4 QPE discovery. This would be in agreement with many other QPE sources with no evidence of bright X-ray sources much earlier than the onset of QPEs. However, we are not able to constrain the start of the QPE behavior. Regarding the peaks of the eruptions, the eROSITA, Swift-XRT and NICER fluxes are compatible with the range spanned by the faintest and brightest eruption in eRO4-XMM (Fig. 6). Finally, we note that similarly to eRO-QPE3 no significant optical variability is observed (e.g., ASAS-SN and ZTF; Shappee et al. 2014; Bellm et al. 2019) during the epochs shown in Fig. 12.

5. Discussion

5.1. QPEs in preexisting AGN?

The X-ray properties of all QPE sources suggest that the nucleus is currently active, as the soft X-ray spectrum observed in quiescence is indicative of radiatively efficient accretion (M19; G20; A21; Chakraborty et al. 2021). In eRO-QPE4 the quiescence emission was a few times fainter, or nonexistent, prior to the QPE discovery (Fig. 12). This would support a transient nature for the quiescent accretion disk and disfavor a long-lived radiatively efficient accretion flow. However, in other QPE sources it is much harder to say whether or not this was the case. In fact, most QPE sources do not have constraining archival X-ray observations prior to the current activity phase, and results from the radio and optical bands are currently ambiguous. For instance, eRO-QPE2 is detected as a faint radio source, although its radio emission is consistent with that predicted from star-formation alone (Arcodia et al., in prep.). All other eROSITA QPE sources are undetected in radio (e.g., this work), which also rules out a preexisting powerful AGN in these nuclei and agrees with the lack of a canonical hard and bright X-ray power-law component in QPE sources. From the optical band, given the lack of broad lines detected, we do not have information on the recent nuclear activity. Narrow-line diagnostics place the select number of QPE sources in different places of the typical diagnostic diagrams including regions indicating the need of a nuclear ionizing source (e.g., an AGN, Wevers et al. 2022). However, these diagnostics cannot securely determine whether these lines were excited by the current activity epoch, while spatially resolved observations will help solving this ambiguity. As a matter of fact, GSN 069 was found to be evidently powered by an AGN using narrow-line diagnostics (Wevers et al. 2022). However, recent work by Patra et al. (2023) making use of archival HST data of GSN 069 found evidence of a compact nuclear [OIII] region, which indicates that the accretion activity is only ∼10 − 100 yr old. More radially extended [OIII] regions, which are not indicative of the current accretion activity, would instead contaminate typical spectra which populate the narrow-line diagnostic diagrams. In absence of a comparable systematic study with spatially resolved observations on all QPE sources, simplistic narrow-line diagnostics have to be interpreted with caution. Significant progress can be made in the near future using integral-field unit data, and high angular resolution spectroscopic and photometric observations, to understand the spatial distribution of the gas in these galaxies and how it relates to the current QPE activity.

Fig. 13. Relation between recurrence and peak luminosity (top) and recurrence and duration timescales (bottom) in secure QPE sources, shown as labeled. The extent of shaded areas represents the standard deviation, when available, or the span within observed values in a given source. eRO-QPE3 is a clear outlier of the top panel, ruling out the presence of a predictive correlation between timescales and peak luminosity.

5.2. What is new on QPEs based on these discoveries?

The discovery of QPEs on top of a decaying X-ray continuum in eRO-QPE3 (Fig. 11) offers a new piece of evidence that strengthens the connection between QPEs and previous TDEs in the same nucleus (Shu et al. 2018; Chakraborty et al. 2021; Sheng et al. 2021; Miniutti et al. 2023a; Quintin et al. 2023). Moreover, eRO-QPE3 showed the faintest peak luminosity (average of $L_{0.5-2.0\mathrm{keV}}^{\mathrm{peak}}\sim 9.7\times {10}^{41}$ erg s⁻¹ in the eRASS1-4 epoch, down to ∼4.2 × 10⁴¹ erg s⁻¹ in the eRO3-XMM1 epoch) and the highest recurrence times (≈20 h) among secure QPEs, therefore invalidating any strict relation between peak luminosity and recurrence time across the QPE population. In fact, eRO-QPE1 shows similar recurrence times but a luminosity up to 10−100 times higher (A21). We show this in the top panel of Fig. 13, where the extent of shaded areas represents the standard deviation, when available, or the span within observed values in a given source (M19; G20; A21; Arcodia et al. 2022; Miniutti et al. 2023a,b; Chakraborty et al. 2024; Giustini et al., in prep.). This lack of correlation disfavor accretion-based QPE emission models, in which both peak luminosity and accretion timescales would depend on the black hole mass. Instead, for EMRI models, the somewhat more diverse orbital parameters (semi-major axis, eccentricity) and the orbiter’s mass and nature (i.e., stars of different kind, degenerate stars and/or black holes) would more easily accommodate scenarios in which the peak luminosity and average recurrence are decoupled. Conversely, as shown in previous work (Chakraborty et al. 2021; Guolo et al. 2024) recurrence and duration timescales do appear to be somewhat correlated (bottom panel of Fig. 13), although we refrain from attempting to fit any scaling due to the low number of sources and the lesson learned from the top panel of Fig. 13. Nonetheless, for the sake of designing a somewhat informed follow-up monitoring based on a single observed eruption (e.g., Quintin et al. 2023), this tentative relation may be used based on the current QPE sources.

eRO-QPE4 shows evidence of a harder component appearing close to the QPE peak (Fig. 7), although at this stage it is unclear what its origin might be (Appendix A.2). Modeling it as Comptonization of the QPE thermal continuum, we infer a flux ratio of ∼1 : 94 in the 0.2 − 2.0 keV band compared to the latter. However, this is merely a possible interpretation, as after all the inferred slope is rather unconstrained and its signal very weak. Nonetheless, we note that until now, a harder component was only seen in quiescence (and not ubiquitously, M19; G20; A21; Chakraborty et al. 2021), but never during the eruptions, let alone in a transient fashion.

In eRO-QPE4 we see clearly, for the first time, that the quiescence emission as detected during the QPE epoch must have been either absent or much fainter previous to the QPE behavior (Fig. 12). What is particularly intriguing is that in eRASS4, when the first QPE was observed, the quiescence emission was also much fainter than the eRO4-XMM epoch. Given the lower sensitivity with Swift-XRT and NICER, their upper limits are shallower than the eRO4-XMM level therefore we are not able to constrain how fast the brightening was.

The long-term emission of both eRO-QPE3 and eRO-QPE4 can only be constrained qualitatively. In the former case the fainter peak flux, close to NICER’s background level, prevents us from securely resolving the fainter eruptions, which appear to be more frequent than the bright ones around the NICER and XMM-Newton epochs (Figs. 4 and 8). In eRO-QPE4 eruptions have a FWHM of ≈30 min, which is on the order of the typical separation in a NICER monitoring given its ∼90-min orbit. Nonetheless, based on NICER data we are able to identify the timing properties of eRO-QPE4 as somewhat intermediate compared to the other known QPEs. The eruptions are not strictly regular or in alternation, as from Fig. 9 we roughly constrain the last five uninterrupted recurrences to be ∼12.1 h, ∼14.2 h, ∼14.0 h, ∼12.4 h and ∼15.5 h, combined with the ∼9.8 h and ∼14.7 h separations seen by XMM-Newton (Fig. 6). They are somewhat similar to the irregular recurrences in eRO-QPE1 (A21; Chakraborty et al. 2024) or RX J1301.9+2747 (G20; Giustini et al., in prep.), although in eRO-QPE4 they do not show burst superposition such as in eRO-QPE1 (Arcodia et al. 2022), nor closely spaced eruptions such as RX J1301.9+2747 (Giustini et al., in prep.). Hence, eRO-QPE4 shows intermediate properties (between regular and irregular arrival times) confirming what the latest orbital models suggest (e.g., Linial & Metzger 2023; Franchini et al. 2023), namely that precession would induce different timing behavior in a continuum from regular (no disk precession, lower eccentricity) and irregular (disk and orbital precession, higher eccentricity).

We also note that in eRO-QPE3, in which the bursts recur every ≈20 h, the asymmetry (with a fast rise and slow decay) appears evident (Fig. A.9, top), whereas in eRO-QPE4 bursts are consistent with being symmetric within the current statistics (Fig. A.9, bottom). Interestingly, the QPE sources with the clearest asymmetry (eRO-QPE1 and eRO-QPE3) are also the ones with undetected quiescence.

The two discoveries in this work confirm the trend of finding QPEs from low-mass galaxies (logM_* ≈ 9 − 10) and low-mass black holes (logM_BH ≈ 5 − 7), where the high-mass end has been now pushed to larger values by eRO-QPE4. We remind readers that the eROSITA search applied to discover the first four eROSITA QPEs is blind in terms of the host galaxy nature, as we only select extragalactic sources with Gaia without selecting specific galaxy counterparts. This might be a selection effect due to observational and model-dependent reasons. An observational bias could be that for more massive galaxies and black holes, the Wien tail of the accretion flow would be shifted out of band. In this case, if the temperature of the QPE component is somewhat dependent on the temperature of the accretion flow, the QPE emission would also be softer and currently unobservable. So far, the observed disk and QPE temperatures have spanned a fairly consistent range within all the known QPE sources, although the absence of harder eruptions is remarkable, despite being easily observable. Larger sample statistics on the QPE population is required for more conclusive statements. In addition, for partial TDEs and stellar EMRI models the detectability of QPEs is limited to smaller galaxies and MBHs for which the tidal radius is not within the innermost stable circular orbit (e.g., Hills 1975; van Velzen 2018; King 2023; Linial & Metzger 2023).

6. Summary and prospects

We report on the discovery of two more galaxies showing QPEs in the eROSITA all-sky survey data, eRO-QPE3 and eRO-QPE4. eRO-QPE3 was found to flare once in multiple eROSITA surveys (Fig. 1) and its nature as a QPE emitter was later confirmed by XMM-Newton (Fig. 4) and NICER (Fig. 8). Eruptions are observed to last ∼2.1 − 2.4 h and to recur every ∼17 − 20 h, although the exact recurrence pattern is currently poorly constrained. eRO-QPE4 was instead detected only once in a transient fashion across all eROSITA surveys (Fig. 2) and both XMM-Newton (Fig. 6) and NICER (Fig. 9) confirmed the repeating nature of the source. Eruptions in eRO-QPE4 are observed to last ≈0.5 h and recur every ∼9.8 − 15.5 h.

Many of their properties are in agreement with those of known secure QPEs (M19; G20; A21), for instance:

• Their X-ray spectra during the eruptions are soft and thermal in shape, with typical peak temperatures of $kT={111}_{-5}^{+6}$ eV and kT = (123 ± 4) eV for eRO-QPE3 and eRO-QPE4, respectively, as observed by XMM-Newton.

• The emission in quiescence, when detected, is consistent with the Wien tail of a radiatively efficient accretion disk, with inner temperatures kT_in ∼ (90 − 100) eV and kT_in ∼ 50 eV for eRO-QPE3 and eRO-QPE4, respectively (all spectral parameters are reported in Tables A.1 and A.2).

• Both sources show a clear energy dependence during the eruptions, with a harder rise than decay at the same flux level (Figs. 10 and A.11), which was first reported for eRO-QPE1 in Arcodia et al. (2022) and then recovered in both GSN 069 (Miniutti et al. 2023a) and RX J1301.9+2747 (Giustini et al., in prep.). At this stage, we consider it the only consistent property across all secure QPE sources, and therefore it is a key requirement to name an X-ray repeater as such.

• Their optical counterparts are local galaxies with a seemingly inactive spectrum (Fig. 3) and even if a secure classification is compromised by the current data quality (e.g., see Wevers et al. 2022), a preexisting powerful AGN is ruled out. This is in agreement with the lack of a bright hard power-law spectrum indicative of a corona, the lack of broad lines, infrared photometry suggestive of an AGN-like torus obscurer and with the deep radio non-detections (Table 1).

• Both eRO-QPE3 and eRO-QPE4 are found in low-mass galaxies. In eRO-QPE3 M_* is in the range (0.7 − 2.6)×10⁹ M_⊙, which in turn implies M_BH in the range (0.9 − 5.3)×10⁶ M_⊙ (Reines & Volonteri 2015). Using the stellar velocity dispersion, the inferred M_BH (Gültekin et al. 2009) is however ∼4 × 10⁷ M_⊙, which only confirms how these scaling relations are poorly calibrated at low masses. In eRO-QPE4, all estimates are more consistent and point toward a slightly more massive galaxy, with M_* being in the range (0.6 − 1.6)×10¹⁰ M_⊙ and black hole, with M_BH in the range (1.7 − 6.8)×10⁷ M_⊙, with the velocity dispersion also yielding ∼2.7 × 10⁷ M_⊙. eRO-QPE4 therefore extends the known QPE population to slightly higher masses.

Furthermore, there are many novelties brought by the new discoveries, which only further highlight the importance to find more QPE sources to build meaningful statistics for a population study. For instance, we report the following:

• eRO-QPE3 showed QPEs on top of a decaying quiescence emission (Fig. 11), in support of other evidence that connects QPEs to a previous TDE (Chakraborty et al. 2021; Sheng et al. 2021; Miniutti et al. 2023a; Quintin et al. 2023).

• eRO-QPE3 shows the longest recurrence time (≈20 h) and the lowest peak luminosity ($L_{0.5-2.0\mathrm{keV}}^{\mathrm{peak}}\sim \mathrm{few}\times {10}^{41}$ erg s⁻¹) among known QPEs (Fig. 13, top), thus further strengthening that there is no correlation between these two quantities. Such timescale-luminosity correlations are instead predicted in many accretion-instability models.

• eRO-QPE4 shows evidence, for the first time, of a harder component appearing close to the QPE peak (Fig. 7), although it is very faint and of an unclear origin (Appendix A.2).

• In eRO-QPE4 we see clearly, for the first time, that the quiescence emission detected during the QPE epoch must have been either absent or much fainter previous to the QPE behavior (Fig. 12), supporting a short-lived nature of the radiatively efficient accretion flow seen in-between QPEs.

• The long-term recurrence pattern of both eRO-QPE3 and eRO-QPE4 cannot be accurately constrained based on current data, although in eRO-QPE4 (Fig. 9) it appears somewhat intermediate between regular alternating and irregular, perhaps bridging the gap between the apparent dichotomy that would have been inferred from the only four sources available (e.g., Arcodia et al. 2022); a continuum of variability patterns would be expected in the case of orbital phenomena with orbital and/or disk precession in play.

The detection of eRO-QPE4 as a non-repeating bright and soft extragalactic flare, compared to the repeated nature of eRO-QPE1, eRO-QPE2 (A21) and eRO-QPE3 (this work), opens up a new way to select more candidates within the available eROSITA surveys. Furthermore, the wealth of eROSITA data can still be used for searches for QPEs performed in a more informed way, by targeting galaxies strikingly similar to those of the known QPE sources (e.g., Wevers et al. 2022). Looking at the future, wide-area X-ray telescopes that are sensitive in the soft X-rays (< 2 keV) are needed to efficiently detect QPEs. In the meantime, thorough archival searches are still very useful (Webbe & Young 2023) and have already proved successful in providing QPE candidates (Chakraborty et al. 2021; Quintin et al. 2023). Finally, we note that the homogeneous nature of our search with the eROSITA telescope allows us to estimate the intrinsic abundance rate of QPEs, which will be presented in a companion paper (Arcodia et al. 2024).

Acknowledgments

We thank the teams of both XMM-Newton and NICER for their prompt efforts in the scheduling of our observations. R.A. thanks G. Miniutti and M. Giustini for the many lively discussions about QPEs. R.A. received support for this work by NASA through the NASA Einstein Fellowship grant No. HF2-51499 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. R.A. thanks T. Nowa and B. Banks and their sets for inspiring the writing. G.P. 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). This work was supported by the Australian government through the Australian Research Council’s Discovery Projects funding scheme (DP200102471). 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: https://doi.org/10.26093/cds/vizier). The original description of the VizieR service was published in Ochsenbein et al. (2000). Part of the funding for GROND (both hardware as well as personnel) was generously granted from the Leibniz-Prize to Prof. G. Hasinger (DFG grant HA 1850/28-1). M.K. acknowledges support from DLR grant FKZ 50 OR 2307. We acknowledge the use of the matplotlib package (Hunter 2007).

References

Appendix A: X-ray spectral analysis

A.1. eROSITA data of eRO-QPE3

The eRASS1 spectrum in quiescence, if modeled with a simple accretion disk, shows marginal residual at high-energies (top panel of Fig. A.1) and a median peak temperature of kT_disk ∼ 133 eV. This temperature would be unusually high compared to other QPE-sources in quiescence, in which the presence of a fainter harder component is often observed (M19; G20; A21; Chakraborty et al. 2021). Therefore, we add a power-law component to account for the residuals, mimicking thermal Comptonization of the accretion disk photons (top panel of Fig. A.3). Compared to a disk model alone, logZ improves by 2.4, however we note that spectral simulations are needed to calibrate the δlogZ for this model comparison. This more complex model is mostly chosen as best-fit model due to the comparison with other quiescence spectra of QPE sources, obtained with higher signal-to-noise spectra (M19; G20; A21; Chakraborty et al. 2021). With the addition of the power-law component at higher energies, the median (and 16th, 84th percentiles) peak temperature of the disk model becomes $k{T}_{\mathrm{disk}}={100}_{-20}^{+18}$ eV. The slope of the power-law component Γ_X is unconstrained, with a slight preference toward softer values⁷. The median (and related 16th, 84th percentiles) flux of the disk component is $F_{0.2-2.0\mathrm{keV}}=1.4_{-0.6}^{+0.4}\times {10}^{-12}$ erg s⁻¹ cm⁻², while that of the power-law component is $F_{0.2-2.0\mathrm{keV}}=2.7_{-1.7}^{+8.4}\times {10}^{-13}$ erg s⁻¹ cm⁻². The total flux in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=1.9_{-0.6}^{+0.7}\times {10}^{-12}$ erg s⁻¹ cm⁻². This model is held fixed during the eruptions and a black body component is added (bottom panel of Fig. A.3). Since the eRASS1 light curve (top-left panel of Fig. 1) does not show variability of the same amplitude compared to eRASS2-3-4, we explore the possibility that no QPEs were present during eRASS1 and that the quiescence continuum was simply highly variable. For instance, in GSN 069, the only other secure QPE source with a clearly decaying quiescence emission over months-years as in eRO-QPE3, QPEs were not detected when the quiescence emission was the brigthest (Miniutti et al. 2023a). Therefore we also fit the quiescence model (diskbb and zpowerlw) to the full eRASS1 observation, without separating faint eROdays from the bright one (top-left panel of Fig. 1). In this work, we do not aim to reach a conclusive statement on the presence of QPEs in eRASS1 or not. Model comparison cannot be performed due to the use of different good-time intervals (GTIs) for the spectra in the two scenarios. We show this fit in Fig. A.2 and we report the fit parameters in Table A.1 together with other eRASS1 fit results.

Fig. A.1. Spectral fit for the eRASS1 (top) and eRASS2 (bottom) spectra in quiescence modeled with a diskbb. Black points are source plus background data, empty gray points show the background alone. The green line and related light green (dark green) shaded regions are the source plus background model median and 1st-99th (16th-84th) percentiles, respectively. The orange dashed lines shows the background model alone. In the lower panel, the data-model ratio is shown. The individual absorption-corrected source model component (here a diskbb) is shown with a red line, as indicated by the legend.

Fig. A.2. As in Fig. A.1, but for the quiescence model (diskbb and zpowerlw) applied to the eRASS1 spectrum of the full observations, without separating faint eROdays from the bright (top-left panel of Fig. 1).

Fig. A.3. As in Fig. A.1, but for the chosen best-fit models for both eRASS1 in quiescence (top, diskbb and zpowerlw) and in eruption (bottom, adding a zbbody).

The eRASS2 spectrum in quiescence appears fainter than in eRASS1 (see Fig. A.1 and Fig. 1) and can be adequately fit by either a diskbb or a power-law continuum. Models provide equally good fits due to the low counts statistics (32 counts within 0.2 − 8.0 keV versus, e.g., the 95 in eRASS1). Adopting the disk model, based on results of other QPE sources, we obtain kT_disk ∼ 151 eV, which is higher than the disk temperature in eRASS1 (Table A.1). We consider unlikely that the quiescence flux is lower and the temperature higher (e.g., Miniutti et al. 2023b,a, for the long-term evolution in GSN 069). Therefore, we attribute this result on the eRASS2 quiescence spectrum to the use of an incomplete source model. Based on eRASS1 and results of other QPEs in quiescence, we adopt the diskbb plus zpowerlw model as reference, despite not being statistically required compared to a single component. If left free to vary, the power-law component dominates the fit and the disk component is fit at background level (Fig. A.4), with a highly-uncertain temperature ($k{T}_{\mathrm{disk}}={80}_{-59}^{+85}$ eV). The median (and related 16th, 84th percentiles) flux of the power-law component is $F_{0.2-2.0\mathrm{keV}}=4.2_{-4.0}^{+1.7}\times {10}^{-13}$ erg s⁻¹ cm⁻² while the disk component is an upper limit at < 4.2 × 10⁻¹³ erg s⁻¹ cm⁻², with a median value of the flux posterior at ∼1 × 10⁻¹⁴ erg s⁻¹ cm⁻². Based on eRASS1 data of eRO-QPE3 and results obtained for other QPE sources with higher signal-to-noise spectra (M19; G20; A21; Chakraborty et al. 2021), we assume that this fit result is likely due to the greater flexibility of the power-law model and model degeneracies. A more realistic result with a diskbb plus zpowerlw would find the former brighter than the latter. We therefore impose that the normalization of the power-law is tied to that of the disk with the same ratio observed in the eRASS1 spectrum. The assumption is that both components have faded. We note that, for the scope of this work, the main aim is to obtain a reliable quiescence model to be held fixed in the QPE spectrum, to isolate this component. We verified a posteriori that the use of either free or tied power-law components does not affect the fit of the QPE spectrum in the bright state. We therefore adopt this tied diskbb plus zpowerlw model as reference for the quiescence (top panel of Fig. A.5). The median (and related 16th, 84th percentiles) temperature of the disk becomes $k{T}_{\mathrm{disk}}={89}_{-60}^{+42}$ eV, thus slightly more constrained and still compatible with the former value. Its flux becomes $F_{0.2-2.0\mathrm{keV}}=3.{12}_{-3.11}^{+1.52}\times {10}^{-13}$ erg s⁻¹ cm⁻², therefore still compatible with background. This confirms that we did not artificially imposed the presence of a disk, but merely allowed the posterior of its flux to be more reasonably brighter. The flux of the power-law component is $F_{0.2-2.0\mathrm{keV}}=2.4_{-1.1}^{+0.9}\times {10}^{-13}$ erg s⁻¹ cm⁻². The total flux in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=5.8_{-2.9}^{+1.7}\times {10}^{-13}$ erg s⁻¹ cm⁻². As done for eRASS1, this quiescence model is held fixed during the brighter eROday (the orange point in the top-medium panel of Fig. 1) and a black body component is added (see bottom panel of Fig. A.5). The median (and 16th, 84th percentiles) peak temperature of the QPE component is $k{T}_{\mathrm{QPE}}={82}_{-9}^{+10}$ eV and its flux is $F_{0.2-2.0\mathrm{keV}}=3.2_{-0.7}^{+1.0}\times {10}^{-12}$ erg s⁻¹ cm⁻². The median temperature is colder compared to the bright eROday in eRASS1, despite the flux being compatible. However, given the short ∼40 s exposure of an eROday compared to the typical QPE duration (0.5 − 7 h, A21), eRASS data catch the eruption at different phases and they are lower limits of the QPE peaks. The total flux during the bright eROday, quiescence included, in the observed X-ray band is $F_{0.2-8.0\mathrm{keV}}=3.9_{-0.8}^{+1.0}\times {10}^{-12}$ erg s⁻¹ cm⁻². Therefore, despite the decreasing quiescence level, the flux in the bright eROday is compatible between eRASS1 and eRASS2.

Fig. A.4. As in Fig. A.1, but for the quiescence spectrum of eRASS2, where both disk and power-law components are free to vary.

Fig. A.5. As in Fig. A.3, but showing the best-fit models for eRASS2, in quiescence (top) and eruption (bottom).

Fig. A.6. As in Fig. A.1, but for the QPE spectrum of eRASS3.

Fig. A.7. As in Fig. A.1, but for the QPE spectrum of eRASS4.

Fig. A.8. As in Fig. A.1, but for the QPE spectrum of eRASS5. As the fit model is compatible with background within the 3σ contours, eRASS5 fit results are to be interpreted with caution.

A.2. XMM-Newton spectral analysis

A.2.1. eRO-QPE3

The quiescent state of eRO-QPE3 was selected from GTIs excluding the bursts from both eRO3-XMM1/2 observations: namely taking GTIs between times 774645000 − 774651000 s and 774675000 − 774724000 s of eRO3-XMM1 and GTIs < 776410000 s and > 776440000 s of eRO3-XMM2. The resulting quiescence spectrum is too faint to be securely detected above background (dark-gray points in Fig. 5). We adopt a disk spectrum to provide a flux upper limit of F_{0.2 − 2.0 keV} <  1.1 × 10⁻¹⁴ erg s⁻¹ cm⁻². The fit posterior on the flux, or disk normalization, is unconstrained (dark-gray line and contours in Fig. 5), while the disk temperature is loosely constrained to $k{T}_{\mathrm{disk}}={46}_{-18}^{+21}$ eV. We conservatively consider the quiescence of eRO-QPE3 to be undetected at the XMM-Newton epoch, although the weak constraints of temperature and luminosity (L_disk, bol <  2.3 × 10⁴¹ erg s⁻¹, with a median value of the unconstrained posterior at L_disk, bol ∼ 3.5 × 10⁴⁰ erg s⁻¹) are consistent with quiescent states of other QPEs (M19; G20; A21). Therefore, we suggest that the quiescent state of eRO-QPE3 is compatible with background (hence not formally detected), but not significantly fainter.

Fig. A.9. Top: a chunk of the eRO3-XMM1 observation (see Fig. 4). The best-fitting model (Eq. 1) is superimposed as a green line (median) and related percentile contours (equivalent to 1σ and 3σ). The vertical shaded areas represent the phases of the burst (namely rise1, rise2, peak, decay1, decay2, with the same color-coding as Fig. 10 for eRO-QPE3). Bottom: same as the top panel, but for the eRO4-XMM observation. The best-fitting model is here shown in red, and the phases of the burst follow the color-coding of eRO-QPE4 in Fig. 10.

The QPE flares are separated in two rise states, one peak and two decays. We show an example of this in Fig. A.9 for the first burst of the observation eRO3-XMM1. The color-coding, darker to lighter going from rise1 to decay2, follows that of Fig. 10. Separations are defined in terms of total count rate, so that, for instance, the separation between rise1 and rise2 occurs roughly at the same count rate as that between decay1 and decay2. Naturally, they are somewhat arbitrary and, as much as temperature or flux values would slightly change with different definitions, no major conclusion of this work would, not the hysteresis shown in Fig. 10. For the example in Fig. A.9, the separating times are 774651000, 774654500, 774657500, 774661000, 774664000 and 774668500 s, respectively. Since the quiescent state is formally undetected, we fit all phases of all bursts with a black body model, representative of the QPE. The individual temperature and luminosities are shown in Fig. 10 and Table A.1. For the peak spectrum (orange in Fig. 5), we obtain $k{T}_{\mathrm{QPE}}={111}_{-5}^{+6}$ eV and $F_{0.2-2.0\mathrm{keV}}=3.1_{-0.2}^{+0.3}\times {10}^{-13}$ erg s⁻¹ cm⁻² (L_{0.2 − 2.0 keV} = (4.2 ± 0.4)×10⁴¹ erg s⁻¹), thus ≳25 − 30 times brighter than the quiescence flux upper limit in the soft X-ray band. The bolometric luminosity of the black body component is L_QPE, bol = (4.9 ± 0.5)×10⁴¹ erg s⁻¹. Based on the eRASS1 and eRASS2 results (and XMM-Newton analysis of eRO-QPE4; see Sect. 2.4.2 and see below), we test the addition a Comptonization component to the QPE phases. However, during all phases the flux of this component is compatible with background (thus, with zero flux) and the log Z of the fit is orders of magnitude worse. Therefore, this component is not statistically required given the current data quality. We adopt the simple black body as reference model for the eruptions of eRO-QPE3 observed by XMM-Newton.

A.2.2. eRO-QPE4

The quiescent state of eRO-QPE4 was selected from GTIs excluding the bursts from the eRO4-XMM observation: namely taking GTIs between times 794814000 − 794837000 s and 794850500 − 794889800 s. The source is detected, with a soft spectrum (Fig. 7) reminiscent of that of other QPE sources in quiescence. We initially fit the quiescence with a disk model (Fig. A.10) and subsequently add a Comptonization component (Fig. 7). As it can be noted in Fig. A.10, a thermal disk model alone fails to account residual signal between ∼0.6 − 1.5 keV. This signal is significantly above background (Fig. 7), as it is comparable to it after background subtraction (Fig. A.10). Adding a Comptonization component vastly improves the logZ of the spectral fit by 92. Therefore, we adopt the disk plus Comptonization model as reference for the quiescent state of eRO-QPE4, which yields a disk temperature of kT_disk = (43 ± 2) eV with a flux F_{0.2 − 2.0 keV} = (3.9 ± 0.4)×10⁻¹³ erg s⁻¹ cm⁻². This corresponds to $L_{0.2-2.0\mathrm{keV}}=1.7_{-0.1}^{+0.2}\times {10}^{42}$ erg s⁻¹ for the disk component alone, or to a bolometric luminosity of $L_{disk,bol}=2.0_{-0.3}^{+0.6}\times {10}^{43}$ erg s⁻¹. The Comptonization component is much fainter, with $F_{0.2-2.0\mathrm{keV}}=1.4_{-0.4}^{+0.8}\times {10}^{-14}$ erg s⁻¹ cm⁻², and its slope is unconstrained with preference for the softest values allowed (i.e., Γ_X = 3.5). As discussed for eRASS data of eRO-QPE3, Comptonization is merely one possible interpretation for the excess signal observed. In this work, the scope of the quiescence fit is to obtain a good-enough model in order to better isolate the QPE component. We report both models in Table A.2. We have also tested a Comptonization component alone to fit the observed X-ray photons, with the assumption that a colder disk is present. Residuals are much poorer and ΔlogZ ∼ 156 compared to the disk plus Comptonization quiescence model, hence the latter is preferred.

Fig. A.10. Same as Fig. 7, but showing background-subtracted spectra. The models shown are a simple disk spectrum for quiescence and a black body in addition for the peak, both showing significant residuals at higher energies. These residuals are comparable with background or even brighter. The background spectrum (subtracted to yield the data shown) is represented in in red for visualization.

The QPE flares of eRO-QPE4 were separated in five epochs, namely rise1, rise2, peak, decay1 and decay2 (bottom panel of Fig. A.9 for the second burst of the observation eRO4-XMM). The color-coding, darker to lighter going from rise1 to decay2, follows that of Fig. 10. For the example in Fig. A.9, the separating times are 794837000, 794839000, 794840500, 794843500, 794845500 and 794850500 s, respectively. The quiescence model was held fixed by imposing its parameters to vary only within the 10th-90th percentile interval of the posteriors obtained during the quiescence fit alone. Therefore, the bursts parameters are marginalized over the uncertainties of the quiescence model. The fit results for the different burst phases are reported in Table A.2. We fit each phase with a black body alone and also with a Comptonization component of the QPE spectrum in addition. For all phases excluding the peak, the additional Comptonization component is not required by the fit, as established by comparing the logZ and visualizing the normalization or flux posteriors of the component. Therefore, for these epochs (rise1, rise2, decay1 and decay2) the best-fit model adopted for the burst is that of the quiescence plus a black body component. For the peak spectrum the logZ of the more complex (black body plus Comptonization) is compatible with that of the simpler model (black body) within typical uncertainties of 0.3 on logZ, thus the former is not formally required. However, the flux and normalization posteriors of the Comptonization component are constrained and the marginal higher-energy residuals accounted for (comparing Fig. 7 and A.10), all at the expense of a marginal change in the black body properties: as it can be noted in Table A.2, both temperature and flux of the much brighter black body QPE component are unaffected by the much fainter harder component. Therefore, we adopt the more complex model as reference for the QPE peak of eRO-QPE4 at the XMM-Newton epoch (see Fig. 7). This model yields kT_QPE = (123 ± 2) eV and F_{0.2 − 2.0 keV} = (2.9 ± 0.1)×10⁻¹² erg s⁻¹ cm⁻² (L_{0.2 − 2.0 keV} = (1.27 ± 0.03)×10⁴³ erg s⁻¹), thus ∼7 − 8 times brighter than the quiescence luminosity in the same soft X-ray band. The bolometric luminosity is L_QPE, bol = (1.44 ± ,0.04)×10⁴³ erg s⁻¹. In this case, the additional Comptonizaton component requires further discussion, since it is not ubiquitous to QPE sources, not throughout the burst of eRO-QPE4 itself. Again, for the scope of this paper this is irrelevant since the QPE black body properties do not change. Future work should establish whether this signal is significant and physical, or whether, for instance, it could be pile-up effect due to the bright and soft QPE spectrum, since it is only observed at the peak.

Fig. A.11. Same as Fig. 10, but for the eruption of the eRO3-XMM2 observation (top) and the third of the eRO4-XMM observation (bottom).

All Tables

All Figures

	Fig. 1. eRASS1-5 light curves of eRO-QPE3 in the 0.2 − 2.0 keV band, each scaled by a reference starting time (MJD ∼58874.111, ∼59057.576, ∼59233.903, ∼59422.701 and ∼59604.611 from teRO, 1 to teRO, 5, respectively). Orange data points highlight the putative detection of eruptions, whilst the gray line connects the background level across the eROdays in each eRASS.
In the text

	Fig. 2. Same as Fig. 1, but for eRASS4 data of eRO-QPE4. The reference starting time is MJD ∼59457.807. eRO-QPE4 was undetected in the previous eRASS1-3 survey.
In the text

Fig. 3. Left panels: 45″ × 45″ cutout of the DESI Legacy Imaging Surveys Data Release 10 (Legacy Surveys/D. Lang (Perimeter Institute)) with the X-ray 1σ position circles in red (eROSITA) and green (XMM-Newton), for eRO-QPE3 (top) and eRO-QPE4 (bottom). Right panels: SALT optical spectra of eRO-QPE3 (top) and eRO-QPE4 (bottom), with inferred spec-z of ∼0.024 and ∼0.044, respectively. The sky spectrum of the red spectral exposure is shown in cyan in arbitrary units in both subpanels.

In the text

	Fig. 4. XMM-Newton 0.2 − 10.0 keV EPICpn light curve of eRO-QPE3, of the first (eRO3-XMM1, left) and second observation (eRO3-XMM2, right). The two eruptions in eRO3-XMM1 are separated by ∼20.4 h. Here, teRO3-XMM1, 0 corresponds to MJD ∼59779.799 and teRO3-XMM2, 0 to ∼59799.745.
In the text

Fig. 5. Source plus background spectrum of eRO-QPE3 in quiescence (dark gray, background dominated) and at the peak of the first burst in the observation eRO3-XMM1 (orange, following the color-coding of Fig. A.9). The light-blue line and contour represent the background spectrum and its uncertainties. Source plus background data points are shown to put the source-only models (lines and contours with the same colors) into context. Darker model contours are 1σ percentiles, lighter ones are 3σ.

In the text

	Fig. 6. XMM-Newton 0.2 − 10.0 keV EPICpn light curve of eRO-QPE4. The three eruptions are separated by ∼9.82 h and ∼14.70 h. Here, teRO4-XMM, 0 corresponds to MJD ∼60013.127.
In the text

Fig. 7. Spectrum of eRO-QPE4 in quiescence (dark gray) and at the peak of the second burst in the observation eRO4-XMM (light green, following the color-coding of Fig. A.9). The dark-red line and contour represent the background spectrum and its uncertainties. Source plus background data points are shown to put the source-only models (lines and contours with the same colors) into context. Darker model contours are 1σ percentiles, lighter ones are 3σ. Individual model components are labeled.

In the text

Fig. 8. NICER background-subtracted 0.2 − 2.0 keV light curve of eRO-QPE3, in the epoch eRO3-NICER1 (top, starting at MJD 59697.644) and eRO3-NICER2 (bottom, starting at MJD 59822.782). Black points represent secure detections of a source component on top of background emission, contrary to the epochs shown in red which are consistent with background-only. Vertical dashed lines represent evenly spaced recurrences of 20 h, as a guide for the eye.

In the text

	Fig. 9. As in Fig. 8 but for eRO-QPE4, starting at MJD 59971.637. Vertical dotted lines represent approximate locations for the eruptions’ peak (thus they are not evenly spaced), as a guide for the eye.
In the text

Fig. 10. Spectral evolution of QPEs in eRO-QPE3 (left) and eRO-QPE4 (right). The top-right (bottom) subpanels show the luminosity (temperature) evolution with time of the QPE component. Darker to lighter colors represent the evolution from start to end of the eruptions. The first eruptions of the eRO3-XMM1 observation and the second of the eRO4-XMM observation are taken as reference for the QPEs in eRO-QPE3 and eRO-QPE4, respectively. The top-left subpanels show the luminosity-temperature coevolution, analogous to the rate-hardness evolution in Arcodia et al. (2022) which first reported the hysteresis cycle.

In the text

Fig. 11. Long-term evolution of eRO-QPE3, separating flux states between quiescence and eruption (which includes the former). eROSITA points are shown as circles, squares for XMM-Netwon, diamonds for NICER (the brightest eruptions in the light curves of Fig. 8) and stars for Swift-XRT. All uncertainties are 3σ. Red contours on orange points highlight observations in which the identification of an eruption is less secure. White symbols for eRASS1 and eRASS2 correspond to the flux of the full observation without separating different states. Horizontal lines highlight ROSAT and XMM-Newton archival upper limits, as stated in the legend.

In the text

Fig. 12. Long-term evolution of eRO-QPE4, separating flux states between quiescence and eruption. eROSITA points are shown as circles, squares for XMM-Netwon, diamonds for NICER and stars for Swift-XRT. Swift-XRT and NICER eruptions points are considered such based on their flux level similar to the range shown by XMM-Netwon. All uncertainties are 3σ. The lower-panel shows Swift-UVOT data with the UVW1 (light blue) and UVW2 filter (azure).

In the text

Fig. 13. Relation between recurrence and peak luminosity (top) and recurrence and duration timescales (bottom) in secure QPE sources, shown as labeled. The extent of shaded areas represents the standard deviation, when available, or the span within observed values in a given source. eRO-QPE3 is a clear outlier of the top panel, ruling out the presence of a predictive correlation between timescales and peak luminosity.

In the text

Fig. A.1. Spectral fit for the eRASS1 (top) and eRASS2 (bottom) spectra in quiescence modeled with a diskbb. Black points are source plus background data, empty gray points show the background alone. The green line and related light green (dark green) shaded regions are the source plus background model median and 1st-99th (16th-84th) percentiles, respectively. The orange dashed lines shows the background model alone. In the lower panel, the data-model ratio is shown. The individual absorption-corrected source model component (here a diskbb) is shown with a red line, as indicated by the legend.

In the text

	Fig. A.2. As in Fig. A.1, but for the quiescence model (diskbb and zpowerlw) applied to the eRASS1 spectrum of the full observations, without separating faint eROdays from the bright (top-left panel of Fig. 1).
In the text

	Fig. A.3. As in Fig. A.1, but for the chosen best-fit models for both eRASS1 in quiescence (top, diskbb and zpowerlw) and in eruption (bottom, adding a zbbody).
In the text

	Fig. A.4. As in Fig. A.1, but for the quiescence spectrum of eRASS2, where both disk and power-law components are free to vary.
In the text

	Fig. A.5. As in Fig. A.3, but showing the best-fit models for eRASS2, in quiescence (top) and eruption (bottom).
In the text

	Fig. A.6. As in Fig. A.1, but for the QPE spectrum of eRASS3.
In the text

	Fig. A.7. As in Fig. A.1, but for the QPE spectrum of eRASS4.
In the text

	Fig. A.8. As in Fig. A.1, but for the QPE spectrum of eRASS5. As the fit model is compatible with background within the 3σ contours, eRASS5 fit results are to be interpreted with caution.
In the text

Fig. A.9. Top: a chunk of the eRO3-XMM1 observation (see Fig. 4). The best-fitting model (Eq. 1) is superimposed as a green line (median) and related percentile contours (equivalent to 1σ and 3σ). The vertical shaded areas represent the phases of the burst (namely rise1, rise2, peak, decay1, decay2, with the same color-coding as Fig. 10 for eRO-QPE3). Bottom: same as the top panel, but for the eRO4-XMM observation. The best-fitting model is here shown in red, and the phases of the burst follow the color-coding of eRO-QPE4 in Fig. 10.

In the text

	Fig. A.10. Same as Fig. 7, but showing background-subtracted spectra. The models shown are a simple disk spectrum for quiescence and a black body in addition for the peak, both showing significant residuals at higher energies. These residuals are comparable with background or even brighter. The background spectrum (subtracted to yield the data shown) is represented in in red for visualization.
In the text

	Fig. A.11. Same as Fig. 10, but for the eruption of the eRO3-XMM2 observation (top) and the third of the eRO4-XMM observation (bottom).
In the text