© R. Arcodia et al. 2022

1. Introduction

A very peculiar and novel entry in the catalog of X-ray emitting sources are the so-called X-ray quasi-periodic eruptions (QPEs), so far comprising a handful of sources: GSN 069 (Miniutti et al. 2019, hereafter M19), RX J1301.9+2747 (Giustini et al. 2020, hereafter G20), eRO-QPE1 and eRO-QPE2 (Arcodia et al. 2021, hereafter A21), and possibly¹ XMMSL1 J024916.6−041244 (Chakraborty et al. 2021). These soft X-ray eruptions repeat quasi-periodically every few hours and show a count rate increase of more than one order of magnitude over a quiescence level. Despite some scatter in the arrival times, QPEs have so far appeared as orderly pulses (see the example in Fig. 1) separated by a steady plateau, usually detected at ≈0.3 − 1.6 × 10⁴¹ erg s⁻¹ in the soft X-ray band (M19; G20; A21; Chakraborty et al. 2021).

Fig. 1. X-ray eruptions of the sources eRO-QPE2 (top, see A21) and GSN 069 (middle, see M19), compared to those of eRO-QPE1 (bottom, see A21). The first two sources show alternating longer and shorter recurrence times and stronger and weaker bursts. In these plots (top and middle), eruptions of each source are overlapped in time-space by scaling each strong-weak-strong cycle (shown with a different color) with the related peak time of the first strong burst. The dashed vertical lines show the time corresponding to half of the separation between two strong bursts (i.e., between consecutive bursts aligned to zero), highlighting that the intermediate weak bursts occur systematically later in time (and that two recurrence times alternate). Bottom panel: the full light curve of eRO-QPE1 reported in A21 is here separated into chunks of two consecutive cycles and scaled in time at the peak of the odd eruptions to highlight a more complex distribution of recurrence and burst duration times.

Quasi-periodic eruptions have so far been identified in nearby low-mass galaxies (with stellar masses of M_star ≈ 1 − 3 × 10⁹ M_⊙; A21). These galaxies host black holes of M_BH ≈ 10^5 − 7 M_⊙, as inferred (with the usual large uncertainties) from methods that are both independent (Shu et al. 2017; M19; G20; Wevers et al. 2022) and dependent on the stellar mass (Strotjohann et al. 2016; A21). Given their peak X-ray luminosity (L_{0.5 − 2.0 keV} ≈ 10⁴² − 10⁴³ erg s⁻¹), the above black hole mass estimates suggest that QPEs peaks are close to being Eddington limited. Studies of the QPEs’ host galaxies have found them to be both active (Esquej et al. 2007; M19; G20) and, quite surprisingly, seemingly inactive galaxies (A21). However, recent studies performing a careful subtraction of the host galaxy’s stellar population have found evidence of an ionizing source in addition to star formation for the narrow lines in all five QPEs (Wevers et al. 2022). Nevertheless, the absence of broad lines in the optical spectra and of infrared photometry excess associated with the dusty environment typical of active galactic nuclei (AGN), intriguingly suggests that a pre-existing canonical AGN-like accretion flow is not a pre-requisite for QPEs. These fascinating objects might therefore provide a new channel to study how massive black holes are activated in the nuclei of low-mass galaxies through transient accretion events, which is to date a poorly studied regime in the black hole–galaxy co-evolution history (e.g., Kormendy & Ho 2013; Heckman & Best 2014; Reines & Volonteri 2015).

The origin of QPEs is still unclear. Proposed scenarios include some types of instability in the accretion disk (M19; Sniegowska et al. 2020; Raj & Nixon 2021) and gravitational lensing in a black hole binary of mass-ratio close to unity (Ingram et al. 2021), though both are currently disfavored (Ingram et al. 2021; A21). Lately, models based on a two- or three-body system with a massive black hole and at least one stellar-mass companion have gained significant attention (King 2020; A21; Suková et al. 2021; Zhao et al. 2022; Xian et al. 2021; Metzger et al. 2022). Such systems could also emit gravitational wave signals detectable by future detectors like LISA and Tianqin (e.g., Amaro-Seoane et al. 2007; Babak et al. 2017; Zhao et al. 2022, but see Chen et al. 2022) and represent an electromagnetic counterpart of the so-called extreme mass-ratio inspirals (EMRIs; e.g., Hils & Bender 1995).

The EMRI-related origin scenarios for QPEs are consistent overall with their multi-wavelength observational properties. The presence of EMRIs is, similarly to tidal disruption events (TDEs), only observable for relatively low black hole masses (e.g., M_BH ⪅ few × 10⁷ M_⊙, Hills 1975; Wevers et al. 2017; Stone et al. 2020), as the stellar-mass orbiting bodies can be otherwise silently swallowed. This mass range is consistent with QPE hosts (Wevers et al. 2022), and is also where nuclear star clusters are nearly ubiquitous (see, e.g., Neumayer et al. 2020, for a review), which was confirmed for at least GSN069 (M19; Sheng et al. 2021) and RX J1301.9+2747 (Shu et al. 2017), and which enhances interactions between the nuclear black hole and stellar-mass objects (Rauch 1999; Babak et al. 2017; Neumayer et al. 2020). At least some of the QPE sources have been associated with young post-starburst galaxies (Caldwell et al. 1999; Wevers et al. 2022), a type in which an enhanced rate of TDEs is observed (e.g., Arcavi et al. 2014; French et al. 2016) due to higher stellar density (Stone & van Velzen 2016; Stone et al. 2018; Law-Smith et al. 2017; Hammerstein et al. 2021) compared to other much more common types of quiescent galaxies. GSN 069 shows UV line-ratios consistent with a past TDE (Sheng et al. 2021), which is also supported by the long-term X-ray emission previous to the discovery of QPEs (Shu et al. 2018; M19). Moreover, a new QPE candidate has recently been found in the XMM-Newton archives (Chakraborty et al. 2021), which shows previous long-term X-ray emission indicative of a past TDE, further strengthening the EMRI interpretation.

However, even if the macro-scenario might have potentially been pinpointed, there are a lot of differences among the various EMRI-related models proposed (e.g., King 2020; Suková et al. 2021; Metzger et al. 2022; Zhao et al. 2022; Xian et al. 2021) and only qualitative comparisons with the QPEs spectral and timing properties have been made so far. With this work we aim to improve on this by exploring further the available XMM-Newton timing data of eRO-QPE1. We describe the data reduction and analysis performed in Sect. 2, we present our main results in Sects. 3–5, and we discuss their implications on our current understanding of QPEs in Sect. 6.

2. Data analysis

Two XMM-Newton observations of eRO-QPE1, also known as the z = 0.505 galaxy 2MASS 02314715−1020112, were analyzed, namely Obs. ID 0861910201 taken on 27 July 2020 (hereafter Obs1) and Obs. ID 0861910301 taken on 4 August 2020 (hereafter Obs2). Details on the processing are described in A21. We performed timing analysis on the light curve extracted in the full energy range (0.2 − 10.0 keV), although we note that QPEs are very soft events and most of the counts are below ∼2 keV even in the bright phase (M19; G20; A21). We show and discuss the results of the QPEs full-band light curve in Fig. 2 and Sect. 3. We also performed timing analysis divided into small energy bins: 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, represented throughout the paper with a color map from purple (0.2 − 0.4 keV) to yellow (1.0 − 2.0 keV). We show and discuss related results in Figs. 3–7 and Sect. 4.

Fig. 2. XMM-Newton background-subtracted light curve of the two observations of eRO-QPE1, named Obs1 and Obs2 (left and right panels, respectively). The full energy band was used, even though eRO-QPE1 shows little signal above ≳1.5 − 2.0 keV, even in the bright phase (A21). The two observations are separated by roughly a week; the times shown are scaled by the start time of each observation (with modified Julian date tObs1, 0 ∼ 59 057.834 and tObs2, 0 ∼ 59 065.988). Top panels: total model fitted (green line and 1σ shaded area), which is composed of one (Obs2) or five (Obs1) individual bursts (blue), each with a profile described by Eq. (1). A constant value for the quiescence is also added to the model. The inset in the top left panel shows the peak of the last eruption in Obs1, left out from the linear scale, which was adopted for clearer visualization of the fainter bursts. Bottom panels: ratio of the data to the model, where the shaded contours show the uncertainty in the model.

Fig. 3. XMM-Newton light curve of Obs2 extracted in small energy bins, represented with different colors (see legend). Each light curve is vertically shifted by a constant factor for visualization. The subpanels show the data-model residuals in each energy bin, color-coded accordingly; vertical dashed lines (with related 1σ uncertainties) indicate when the eruptions start in the given energy bin, while vertical solid lines identify the peak. The model adopted is described by Eq. (1), and it is shown in all panels with 1σ uncertainty as shaded intervals.

Fig. 4. Energy dependence of QPEs in the observation Obs2 of eRO-QPE1 (related to the light curve shown in Fig. 3). The times tstart and tpeak are shown in hours elapsed from tObs2, 0. QPEs start earlier at lower energies (top panel) and evolve more slowly in terms of rise and decay (middle panel), and peak at later times (bottom panel) at low energies.

Fig. 5. Same as Fig. 3, but for Obs1. Individual burst profiles are shown in black, with different line styles. Five bursts are required for the 0.2 − 0.4 keV and 0.4 − 0.6 keV bins, while four are needed for 0.6 − 0.8 keV and three in the remaining energy bins (see Table A.1). Lower panels: data/model ratios. Vertical lines indicate tstart of each burst at all energies, with the line style matching the one used for the burst profiles in the upper panels with the light curves. Shaded intervals show 1σ uncertainties.

Fig. 6. Same as Fig. 4, but for the first two bursts of Obs1 (see, e.g. Fig. 5). Squares are used in all panels for the first burst, stars for the second. The energy dependence of τrise, τdecay, and tpeak (from second to bottom panel) is consistent with Obs1 and other QPEs (M19; G20; Chakraborty et al. 2021).

Fig. 7. Plots showing results for Obs1 (resp. Obs2) of eRO-QPE1, on the left (right), represented with colored circles (triangles), for which darker to lighter colors with a red-to-yellow (blue-to-green) color map indicate the increase in time of the observation. Top central panels: count rate (CR) light curve as in Fig. 2, but here in logarithmic units and excluding the quiescence. Bottom panels: evolution of hardness ratio (HR = H/(H + S) = CR0.6 − 2.0 keV/CR0.2 − 2.0 keV) over time. External panels: HR vs. CR plot. In the panels related to Obs1, the central data points are highlighted with black squares for ease of comparison across the panels. The range of Obs2 in HR is shown as black dashed lines across the different panels. The HR reached at the peak of the eruptions is count rate dependent, and the eruptions undergo a counterclockwise hysteresis cycle in the HR vs. CR plot (top right).

The modeled profile shape adopted to fit the eruptions was chosen using Obs2 since it shows a single isolated burst (see right panel of Fig. 2). It is loosely based on the model proposed in Norris et al. (2005) for gamma-ray bursts, and defined as

$$\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 before the asymptote at t_peak − t_as, where $t_{\mathrm{as}}=\sqrt{{\tau }_1{\tau }_2}$, whereas A is the amplitude at the peak and λ = e^t_λ a normalization to join rise and decay, where $t_{\lambda }=\sqrt{{\tau }_1/{\tau }_2}$. We note that τ₁ and τ₂ are two characteristic timescales, although only τ₂ is directly related to the decay timescale. 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. Following again Norris et al. (2005), we define rises and decays 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 respectively as τ_rise = w (1 − k)/2 and τ_decay = w (1 + k)/2 using n = 1. We also define the start of a burst computing the time at which the flux is 1/e³ of the peak value, hence t_start = t_peak − w (1 − k)/2 with w and k evaluated at n = 3. We implemented the model and derived posterior probability distributions, Bayesian evidence, and Akaike information criterion (AIC) values (Akaike 1974) using the nested sampling Monte Carlo algorithm MLFriends (Buchner 2014, 2019) using the UltraNest² package (Buchner 2021). In this work the relative goodness of fit between two models is inferred by comparing AIC values, namely $\mathrm{AIC}=2N_{\mathrm{p}}-2log\widehat{L}$, where N_p is the number of parameters in the model and $log\widehat{L}$ is the maximum likelihood of the fit. The fit is considered improved with a more complex model (e.g., with more burst profiles) if its AIC value is smaller, specifically if there is a negative difference between the respective AIC values (quantity referred to as ΔAIC_m2, m1 = AIC_m2 − AIC_m1, where AIC_m2 is that of the more complex model). For each model comparison, we also checked the difference in logarithmic Bayesian evidence and find results consistent with the ΔAIC. We note that other burst profiles were tested, for instance the model as defined in Norris et al. (2005), which is smooth at the peak, or a model with a Gaussian rise and exponential decay (van Velzen et al. 2019), but they obtained worse fits on Obs2, and were therefore discarded. Throughout the paper, we quote median values of the posterior chains, with related 16th and 84th percentile values, unless otherwise stated.

For the comparison with eRO-QPE1 in Fig. 1, we also analyzed and presented XMM-Newton data of eRO-QPE2 (A21) and GSN 069 (M19) and NICER data of eRO-QPE1 (A21). For eRO-QPE2, data processing for Obs. ID 0872390101 (August 2020) was performed and explained in detail in A21. Similarly, details of the processing of NICER data (Obs. ID 3201730103) of eRO-QPE1 are reported in A21. For GSN 069, we processed Obs. ID 0831790701 (January 2019) and 0851180401 (May 2019) using standard procedures (SAS v. 18.0.0 and HEAsoft v. 6.25), and the source (background) region was extracted within a circle of 50″ centered on the source (in a source-free region). Photons were extracted between 0.2 and 10.0 keV.

3. Presence of multiple overlapping bursts

Although not exactly periodic, eruptions observed in a given QPE-emitting source have so far been estimated as single isolated bursts with a somewhat regular recurrence pattern (M19; G20; A21). For instance, eruptions in GSN 069 (M19) and eRO-QPE2 (A21) show a repeating alternation of longer and shorter recurrence times after alternating stronger and weaker bursts, respectively (M19; Xian et al. 2021; our Fig. 1).

In Fig. 1 we show the light curves of both eRO-QPE2 and GSN 069 (top and middle panel, respectively) divided in strong-weak-strong QPE cycles, each cycle defined as the separation between two strong QPEs and overlapping in time, by scaling every data point by the t_peak of the first eruption of each cycle. Dashed vertical lines show the time corresponding to half of the separation between two strong bursts (i.e., between the bursts aligned to zero in Fig. 1), highlighting that weak bursts occur systematically later in time than half of the strong-strong recurrences. Furthermore, all superimposed strong (or weak) eruptions appear to be overall compatible with one another.

In general, eRO-QPE1 immediately stands out with respect to the other known QPEs (see, e.g., Chakraborty et al. 2021), due to its peak X-ray luminosity (≳10⁴³ erg s⁻¹) and evolving timescales (mean recurrence of ∼0.8 d), which are larger by an order of magnitude (A21). Contrary to eRO-QPE2 and GSN 069, a quick look at a similar visualization for eRO-QPE1 (bottom panel of Fig. 1) also argues for a more complex evolution and shape of the eruptions. To investigate this, we show in Fig. 2 the two deep XMM-Newton observations of eRO-QPE1, Obs1 (left) and Obs2 (right), with the related best fit model and data/model ratios. It is evident from the data points alone that the Obs2 (right panel) shows evidence of a single burst, similar to what is observed for eRO-QPE2 and GSN 069 (in Fig. 1), albeit on much longer timescales and with an enhanced asymmetry. Instead, Obs1 (left panel) clearly shows a much more complex event. We first fit to the much simpler Obs2 the model in Eq. (1) with a constant plateau in addition, obtaining a very good fit and residuals throughout rise, peak, and decay. This flare template was adopted as a reference shape, which we then applied to Obs1, with the assumption that we could model its more complex behavior as a linear combination of the simple template. We started with a sum of three different bursts and increased their number by comparing the related AIC values and by visualizing the fit residuals. A model with four bursts improves the fit with a ΔAIC ∼ −167.8. Residuals around t − t_Obs1, 0 ∼ 20 h (see Fig. A.1) were significantly improved by adding a fifth burst, which improved the four-bursts model with ΔAIC ∼ −91.8. Therefore, we obtained that five independent profiles are needed (see Appendix A for more details). Since the model adopted (Eq. (1)) fits the simpler Obs2 well, we consider it unlikely that Obs1 needs a completely different and more complicated shape for the eruptions, but rather a combination of the simple burst template used for the simpler Obs2. A further relevant question is whether the complex behavior seen in Obs1 was a one-time or rare event, or whether it recurs throughout the long-term evolution of eRO-QPE1. We note that studying the latter with XMM-Newton is observationally expensive, and this will be the focus of an upcoming work (Arcodia et al., in prep.).

As a sanity check, we tested whether any statistical evidence of more complex structure is needed for Obs2 of eRO-QPE1. We obtained that Obs2 is best described with one single eruption (see Appendix A and Fig. A.2 for more details), which validates the use of its template to fit Obs1. Moreover, we performed a similar test on the seemingly much simpler QPE source eRO-QPE2 (see Appendix A). If a small degree of asymmetry is accounted for, no superposition of individual bursts is required to model the eruptions in eRO-QPE2 (see Appendix A and Fig. A.3). Therefore, no significant substructures are present in Obs2 of eRO-QPE1, nor in the bursts of eRO-QPE2.

4. Energy dependence of QPEs

Since the discovery of QPEs, it was soon noticed that once the light curve is decomposed into small energy bins, eruptions appear to peak at later times at lower energies and are broader compared to higher energies (M19; G20). However, it was not investigated in QPE observations whether there is an energy dependence at the start of the rise in QPEs. This is mostly due to the lower signal-to-noise ratio in the energy-resolved light curves and shorter evolution timescales in the first QPE sources (M19; G20). However, finding or excluding an energy dependence of the QPE start is important to understand whether QPEs are consistent with an emission component with a given spectrum simply getting brighter in luminosity, or if instead the spectrum evolves in temperature and/or energy over time during the start, as it seems to do during the rise (M19). Some sort of compact accretion flow should be present in QPEs as they show a stable spectrum during quiescence, which is as bright as ≈0.3 − 1.6 × 10⁴¹ erg s⁻¹ in the soft X-ray band (M19; G20; A21; Chakraborty et al. 2021). Furthermore, eRO-QPE1 is ideal to test the hypothesis of an energy dependence of the QPEs start times thanks to its luminosity and timescales which are one order of magnitude larger than all the other QPEs found so far (M19; G20; A21; Chakraborty et al. 2021).

Hence, we extracted light curves in small energy bins for both Obs1 and Obs2 (Sect. 2). We start from Obs2 due to its simplicity. Figure 3 shows Obs2 divided in energy bins, with the related best fit models obtained using Eq. (1) and data/model ratios in the lower panels. Shaded intervals indicate 1σ uncertainties on the model. We also highlight the median value (with shaded 1σ uncertainties) of the fit posteriors of t_start and t_peak (see Sect. 2) with a vertical dashed and solid line, respectively. Moreover, we show in Fig. 4t_start, τ_rise, τ_decay, and t_peak as a function of energy, as obtained from our light curve fit. The overall picture is that we confirm the previously observed trend for rise, peak, and decay times, namely that QPEs evolve faster and peak earlier at higher energies (M19; G20). Furthermore, for the first time we can resolve the start of the eruptions and find that QPEs start earlier at lower energies (see top panel of Fig. 4).

Figure 5 is analogous to Fig. 3, but for Obs1, confirming the trend found in Obs2 even during a more complex duty cycle. We fit all light curves in small energy bins (see Sect. 2) with a model including three, four, and five bursts. We recursively selected the best fit model in each energy bin by comparing their AIC values. We report the corresponding ΔAIC values in Table A.1. Intriguingly, we obtained that burst number three and four of the full-energy light curve (see Fig. 2, left) not only have the lowest amplitude, but are also much colder than the brighter eruptions. In particular, little or no significant signal is detected above ∼0.6 − 0.8 keV (Fig. 5). Therefore, good enough fits are obtained with four bursts in the 0.6 − 0.8 keV bin and three at the highest energy bins (e.g., ≳0.8 keV). The lower panels of Fig. 5 show the data/model ratios, while vertical lines show the fit t_start for each burst (with the line style matching that used in the upper panels for the burst profile). For bursts starting with some overlap with the previous burst, t_start is obviously less constrained. However, t_start in the first and the last (which is also the brightest) burst of Obs1 does show an increase with energy, consistently with Obs2. We also note that in the energy bins in the range 0.6 − 1.0 keV the difference between the two peaks around t − t_Obs1, 0 ∼ 10 h is enhanced (Fig. 5), confirming the need of a combination of burst profiles rather than a more complex single-burst model. We show in Fig. 6 the fit t_start, τ_rise, τ_decay, and t_peak of these two bursts in Obs1 as a function of energy (shown with different symbols in the legend). The energy dependence of τ_rise, τ_decay, and t_peak (from second to bottom panel) is consistent with Obs1 and other QPEs (M19; G20; Chakraborty et al. 2021). The time t_start of the first burst is also consistent with Obs2 (Fig. 4), in that it increases with energy. Instead, in the second burst t_start seems to be compatible within uncertainties across energies. Given the likely degeneracy with the decay of the previous burst, we refrain from overinterpreting this dependency.

Hence, from both the simpler Obs2 and the more complex Obs1, we inferred an earlier t_start of the eruptions at lower energies. We note that this is unlikely to be a model-dependent artifact as it is clearly observable by eye in the data points of Fig. 3 (comparing the lowest and highest energy bins). Finally, this result cannot be a spurious effect due to the instrument sensitivity, as XMM-Newton is much more sensitive around ∼1 keV than at ∼0.3 keV; and a possible incorrect determination of the counts redistribution in the current XMM-Newton calibration³ would be relevant for the opposite case, namely if higher energy photons were detected first.

Finally, we note that in other QPE sources an energy dependence of the burst amplitude is observed (M19; G20; Chakraborty et al. 2021). It is usually defined with respect to a well-detected and unabsorbed quiescence emission and leads to larger amplitudes in higher energy bins. Instead, in eRO-QPE1 the quiescence is not robustly detected and the amplitude would be computed against the background spectrum of XMM-Newton; therefore, comparison could be misleading. We propose in the next section an alternative way to compare the spectral evolution across different QPE sources.

5. Hardness ratio of the eruptions

We further tested the energy dependence of QPEs by computing the hardness ratio of Obs1 and Obs2, which confirms the results above in a model-independent way. It was computed as HR = H/(H + S) = CR_{0.6 − 2.0 keV}/CR_{0.2 − 2.0 keV}, where CR is the count rate in the related energy bands. We report in the left (right) half of Fig. 7 the results for Obs1 (Obs2), represented with circles (triangles), for which darker to lighter colors indicate the increase in time: the CR light curve is shown in the top middle panel, while the HR evolution over time at the bottom and the HR versus CR plot in the top external panel.

For simplicity, we start with Obs2. The evolution of HR over time indicates a gradual spectral hardening during the rise and softening during decay, as expected from previous results (M19; G20; A21; Chakraborty et al. 2021). However, the spectrum does not harden throughout the whole rise, in fact HR remains constant and then slightly decreases during the second part of the CR rise. We interpret this as the effect of the energy dependence of QPEs, as eruptions are known to peak later in time at lower energies (as we show in Figs. 3–6 and, e.g., M19; G20; Chakraborty et al. 2021), therefore enhancing the softening. Notably, it is the first time that QPEs are shown to undergo a counterclockwise hysteresis in a HR versus CR plot (top right for Obs2); specifically, there is not a unique spectral shape for a given count rate during both rise and decay. The source is softer during decay with respect to the rise for a compatible total count rate. Therefore, QPEs are not only asymmetric in their timing profiles, but also spectrally.

From a comparison with Obs1 (left half of Fig. 7), we see that the maximum hardness reached by an eruption increases with its peak count rate (bottom left and top left panels). The first part of Obs1 (redder colors) seems to show a counterclockwise hysteresis similarly to the eruption in Obs2, and the same is apparent for the final section of Obs1 (yellower colors), although the observation stopped just before the decay. Nonetheless, from the end of both the CR light curve and the HR evolution of Obs1 (top center-left and bottom left panel of Fig. 7), a flattening and possible start of decrease can be seen. Therefore, this suggests that the brightest eruption started a compatible counterclockwise hysteresis, albeit shifted to higher total count rates. Data points related to the weaker state in the middle of Obs1 (highlighted with black squares) seem to confirm the exceptional softness and weakness of the source (though still brighter than the quiescence), as suggested in Fig. 5 as well. We also note that the substructures in the first part of Obs1 are less remarkable in Fig. 7, although they remain quite evident in Fig. 5. As a sanity check for our model-dependent analysis (Sects. 3 and 4), we compared in Fig. 8 the HR versus CR of the modeled profiles of the single eruptions of Obs1 (see Fig. 5; but only for the bursts clearly detected in all energy bands used for the definition of HR) with that of Obs2 (black dashed line). For the latter, the comparison with the plot showing pure data (top right of Fig. 7) is straightforward. Intriguingly, a clear counterclockwise hysteresis is also evident for each of the bursts in Obs1 separately. This indicates the robustness of our decomposition of Obs1 in multiple bursts, as we managed to retrieve for them a spectral behavior compatible with that of the isolated burst of Obs2.

Fig. 8. Hardness ratio vs. count rate computed from the modeled profiles of the bursts in Obs1 (color-coded as in the legend) and Obs2 (dashed black line). For a visual representation of the model profiles see Fig. 2. For Obs1 the profiles are used of bursts that are clearly detected in all the energy bands used for the definition of HR. The thick solid lines represent the median of the modeled tracks, while the thinner lines highlight the span within the model posteriors. A counterclockwise hysteresis (highlighted with arrows on the tracks) is found for all modeled eruptions, which is also what the pure data indicate (Fig. 7).

6. Discussion

6.1. Two possible QPE populations

A first look at QPE observations suggested that in a given source, eruptions would occur more or less regularly, modulo some scatter in the quasi-period, with single seemingly isolated bursts (M19; G20; A21). This seems to indeed be the case for at least two of the four secure QPEs, namely GSN 069 (M19) and eRO-QPE2 (A21). In these two sources, eruptions occur with an alternation between long and short recurrence times, happening after relatively stronger and weaker bursts, respectively (see Fig. 1; M19; Xian et al. 2021). Moreover, strong (weak) eruptions have so far appeared to be overall comparable in their profiles with other strong (weak) ones, lacking any significant and observable signature of overlapping bursts (Fig. 1, Appendix A). Instead, in this work we show that the timing properties in the source eRO-QPE1 can be far more complex and irregular. In particular, we observed the presence of cycles with single isolated bursts (Fig. 2 right, which is similar to what we expected based on other QPEs), as well as other cycles in which we observed multiple overlapping bursts with different amplitudes (Fig. 2, left). This complex timing behavior may also be present in RX J1301.9+2747 (G20), though to a less dramatic extent, with the presence of two bursts with much lower amplitude than all the others (Giustini et al., in prep.). This might hint that two subclasses of QPEs exist, the first subgroup being GSN 069 (M19) and eRO-QPE2 (A21), while the second populated by RX J1301.9+2747 (G20) and eRO-QPE1 (A21). Only by finding more can we infer whether this is the case, or whether instead QPE timing properties are a continuum that spans from rather simple and regular, like GSN 069 and eRO-QPE2, to more complex and irregular, like eRO-QPE1, passing through RX J1301.9+2747.

Given the many similarities in the energy-dependent properties of all QPEs (M19; G20; A21; Chakraborty et al. 2021), we can consider the evidence, reported in this paper (Sect. 4, Figs. 3 and 4), of QPEs starting earlier at lower energies in eRO-QPE1, to be valid for all QPEs. Conversely, more complex timing behavior has so far only been observed to happen for eRO-QPE1 and RX J1301.9+2747 (to a less dramatic extent; Giustini et al., in prep.). Based on the current theoretical knowledge, it is possible that there are two different physical mechanisms in place for the two putative QPE subgroups. However, given the many common observed properties within all known QPEs, we can speculate that there could be one main mechanism, be it star–disk collisions⁴ (e.g., Suková et al. 2021; Xian et al. 2021), Roche-lobe overflow (RLOF) from one or two stars (e.g., Zhao et al. 2022; Metzger et al. 2022), or an alternative model yet to be proposed.

6.2. Implications for the origin of QPEs

Based on our findings, any model that aims to explain QPEs must be able to produce both regularly spaced eruptions, which may appear with alternating shorter and longer cycles, and a more erratic behavior with less clear periodicity and different amplitudes (during some of the cycles at least; see, e.g., Fig. 2). Even if an in-depth modeling of our new results is beyond the scope of this work, we speculate here on their impact on models that involve star–disk collisions (e.g., Xian et al. 2021) and RLOF from one or more stellar companions (e.g., Zhao et al. 2022; Metzger et al. 2022) that were proposed as QPEs. Current models of a pure lensing scenario (e.g., Ingram et al. 2021) now seem disfavored because of the known energy dependence of QPEs, the exact shape and amplitude of the bursts (Ingram et al. 2021), and also considering our new results (Sect. 3), which would need a very complex time-varying gravitational geometry.

In the framework of QPEs from star–disk collisions (e.g., Suková et al. 2021; Xian et al. 2021), a burst is produced at every crossing, each of which is assumed to leave the companion mostly unaffected over the time span studied (Dai et al. 2010; Pihajoki 2016), although long-term perturbations of the orbit are relevant in the case of a normal star (e.g., Syer et al. 1991; Vokrouhlicky & Karas 1993; MacLeod & Lin 2020). One possible outcome of these collisions is the production of a shock at the impact region within the disk (Pihajoki 2016; Suková et al. 2021). This shock in turn produces two fountains of expanding gas, which eventually radiate as they become optically thin (Lehto & Valtonen 1996; Ivanov et al. 1998; Pihajoki 2016), or produces episodic ejections of blobs from the accretion flow (Suková et al. 2021). Only generic predictions have been made for the expected spectral and timing evolution of the resulting emission component (e.g., Semerák et al. 1999; Dai et al. 2010), although it has been suggested that Bremsstrahlung or thermal radiation would be expected in the fountains scenario (e.g., Lehto & Valtonen 1996; Pihajoki 2016), which are both an overall good description for the QPEs spectra (M19; G20; A21). We note that the QPE energy dependence (e.g., Figs. 3 and 4) could be qualitatively reproduced if the collisions induce heating fronts propagating in the disk (e.g., Meyer 1984), although these shock waves are shown to dramatically destroy the inner disk (Chan et al. 2019), which is in tension with the stable quiescence detected in QPEs right after the eruptions (M19; G20; A21). Moreover, this model inherently assumes that throughout several collisions the system is in roughly the same state before each interaction (two collisions per orbit) in terms of physical conditions in both the accretion flow and the star (e.g., Ivanov et al. 1998; Dai et al. 2010; Pihajoki 2016). This is in tension with our results on eRO-QPE1, which sometimes shows cycles with single isolated bursts (Fig. 2, right) and at other times more complex behavior (left panel). This is particularly puzzling since the two observations were taken one week apart and, assuming a periodicity on the order of less than a day (A21), just a handful of cycles have occurred in between. Therefore, despite being a promising scenario (see, e.g., Xian et al. 2021, for GSN 069) more quantitative predictions from these models are needed to conclusively state whether they fail to reproduce the more complex timing observations of eRO-QPE1 (see Fig. 2) and RX J1301.9+2747.

In the model proposed by Metzger et al. (2022), two stellar EMRIs occur in a co-planar orbit around the central black hole, and the observed QPEs periodicity implies that they are counter-rotating. At least one star undergoes RLOF at each flyby and mass is accreted toward the black hole interacting with, and replenishing, an accretion disk, the presence of which is supported by the stable quiescence usually observed in between QPEs (M19; G20; A21). Possible differences in the mass transferred to the black hole can occur from cycle to cycle due to changes in the separation between the two stars, influenced by residual eccentricity in the system (Metzger et al. 2022). This could indeed provide the required diversity in the eruptions. In addition, it is possible that during a given cycle, or most of the cycles, a single star undergoes RLOF, but occasionally both stars do, with a slightly delayed mass transfer rate. Even though this is an extremely complex interaction to model, it is qualitatively in agreement with our observations of eRO-QPE1, which show isolated bursts (Fig. 2, right) as well as more complex cycles (Fig. 2, left). We speculate that this is also a promising model for RX J1301.9+2747 (G20). Interestingly, if one of the two EMRIs is a star and the second is a compact object, only the star can undergo RLOF; this scenario is worth exploring for GSN 069 and eRO-QPE2, which do not show evidence of multiple overlapping bursts. Other models with a single EMRI (King 2020; Zhao et al. 2022) would appear similar to the latter double-EMRI model with one star and one compact object, while they would be harder to reconcile with eRO-QPE1 and RX J1301.9+2747. Moreover, Metzger et al. (2022) show that the post-peak energy dependence of QPEs is reproduced using one-dimensional spreading disk equations (see, e.g., Pringle 1981) to model the evolution of the matter flowing in toward the black hole. However, more dedicated models or simulations are needed for comparison with Figs. 3–6 and to model the early phases, when most likely any steady-state assumption for the accretion flow is violated.

7. Conclusions

Quasi-periodic eruptions are the new frontier of variable accretion onto massive black holes (M19; G20; A21; Chakraborty et al. 2021). They are able to light up the nuclei of otherwise faint and unnoticed low-mass galaxies (A21), therefore providing a new channel to activate transients accretion events onto the black holes at their center. However, these peculiar repeating X-ray blasts still remain unexplained. Lately, models involving one or two stellar-mass companions around the central black hole have gathered significant attention (King 2020; A21; Suková et al. 2021; Metzger et al. 2022; Zhao et al. 2022; Xian et al. 2021). If this is indeed the correct origin scenario, QPEs could also emit low-frequency gravitational waves as EMRIs (e.g., Zhao et al. 2022; but see Chen et al. 2022) and could therefore revolutionize the future of multi-messenger astronomy with Athena (Nandra et al. 2013) and LISA (Amaro-Seoane et al. 2017).

We reported here two new observational results found by taking a closer look at the QPEs in eRO-QPE1 (A21):

• (i)

  At times, eruptions in eRO-QPE1 occur as single isolated bursts (Fig. 2, right), while they can also manifest as a complex mixture of multiple overlapping bursts with very diverse amplitudes (Fig. 2, left); this is in contrast to other known QPEs, namely GSN 069 and eRO-QPE2, for which so far only evidence of somewhat orderly bursts with a more regular recurrence pattern have been found (Fig. 1).

• (ii)

  Studying eruptions in small energy bins, we confirm previous works that found QPEs to peak at later times and to be broader at lower energies (M19; G20; Chakraborty et al. 2021), while for the first time we find that QPEs start earlier in time at lower energies (Figs. 3–6). Furthermore, eruptions appear to undergo a counterclockwise hysteresis cycle in a plane with hardness ratio versus total count rate, implying that the decay of each eruption is softer than its rise at a compatible total count rate (Figs. 7 and 8).

The first result implies that, if we demand a single trigger mechanism for all QPEs, this should be able to produce both regular and complex behaviors in some sources (e.g., in eRO-QPE1 and RX J1301.9+2747), while in others only the more regularly spaced eruptions (e.g., GSN 069 and eRO-QPE2). The second result implies that the X-ray emitting component is not achromatic, namely the resulting spectrum does not simply brighten and fade in luminosity and it has instead a specific energy dependence over time. If there is indeed an accretion flow around these black holes (with any kind of radiative efficiency and thickness) emitting the quiescence signal, the inferred energy dependence might imply the presence of inward radial propagation during the QPEs rise. Current and future models proposed for QPEs should reproduce our new results, as well as the other multi-wavelength properties of these fascinating sources.

Acknowledgments

We thank the referee for the positive and constructive review of our work. We thank N.C. Stone and B.D. Metzger for very insightful discussions on the possible theoretical interpretations of our results. R.A. thanks K. Dennerl for discussions on the EPICpn response. We acknowledge the use of the matplotlib package (Hunter 2007). 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]). G.M. acknowledges funding from Project No. MDM-2017-0737 Unidad de Excelencia “María de Maeztu” – Centro de Astrobiología (CSIC-INTA) by MCIN/AEI/10.13039/501100011033. MG is supported by the “Programa de Atracción de Talento” of the Comunidad de Madrid, grant number 2018-T1/TIC-11733.

References

Appendix A: More details on the light curve fitting

Obs1 requires multiple bursts to be fit. In Fig. A.1 we show the fit with three and four bursts. The latter improves on the former with a difference in AIC values of ΔAIC_4b, 3b ∼ −167.8. Residuals around t − t_Obs1, 0 ∼ 20 h (Fig. A.1 bottom) are significantly improved by adding a fifth burst (ΔAIC_5b, 4b ∼ −91.8 with respect to the four-burst model). Therefore, we adopted the five-burst model as reference for Obs1 in the full energy band. We adopted the same procedure for Obs1 light curves extracted in small energy bins (see Fig. 5) to determine the number of bursts required. We list the comparisons between the AIC values in Table A.1.

Fig. A.1. As in the left panel of Fig. 2, but fitting a model with three (top) and four (bottom) bursts.

To further validate our results, we tested whether a single-burst profile is indeed the preferred model for Obs2. We performed a fit with two overlapping bursts (see Fig. A.2). Despite the good residuals, it is immediately clear that one of the two bursts is relegated to fit a much smaller number of counts, and that Obs2 cannot be fit with two bursts with compatible amplitude, which is instead the case for the first part of Obs1 (left panel of Fig. 2). Finally, comparing the AIC values of the single-burst and double-burst fits (ΔAIC_2b, 1b ∼ 5.5), the former remains the preferred statistical description of Obs2 data. Hence, Obs2 is indeed best described with a single burst and since Obs1 requires many, the two observations of eRO-QPE1 show a fundamental difference in timing properties, which occurred within just a few days and a handful of QPE cycles.

Fig. A.2. As in the right panel of Fig. 2, but fitting a model with two bursts (shown in cyan and light blue).

We performed a similar test on the seemingly much simpler QPE source eRO-QPE2 (see Fig. 1, top). We fit the full light curve with one Gaussian and a superposition of two Gaussian profiles (obtaining AIC_2g and AIC_1g, respectively). The simpler model is still preferred (ΔAIC_2g, 1g ∼ 10). We also performed the same procedure for each single burst separately. In four bursts the AIC value of the two-Gaussian model is lower, while in the four remaining the single Gaussian is a good enough description of the data (ΔAIC_2g, 1g values of -4.6, -4.4, -3.4, 2.6, 1.5, -0.3, 3.9, 8.8, in order of time of the eruptions). However, as shown in Fig. A.3 with one example of a burst with an improved fit, the analysis with two overlapping Gaussian profiles leads to multi-modal features covering the full parameter space of the burst (highlighted by the cyan lines), which is surely not physical. Motivated by the clear asymmetry in eRO-QPE1 we also tested the asymmetric model shown in Eq. 1 on eRO-QPE2, obtaining for the full light curve an improvement in ΔAIC∼ − 22, with respect to the single-Gaussian model. The asymmetric model fit on the single bursts also yielded a compatible improvement, as did the fit with the two-Gaussian model when it was present in a single eruption. Clearly, dedicated analysis is needed for eRO-QPE2, which is beyond the scope of this paper; it will be tested in future work. In contrast, we note that in eRO-QPE1 both the asymmetry and the addition of multiple bursts are not only statistically significant, but also clearly evident in the data.

Fig. A.3. Eruption number four in the XMM-Newton light curve of eRO-QPE2 (see also Fig. 1). The red curve represents the fit with one Gaussian component, the blue curve the superposition of two Gaussian profiles is shown, and the cyan lines show the two individual Gaussian components (the median is shown with a thicker line).

All Tables

All Figures

Fig. 1. X-ray eruptions of the sources eRO-QPE2 (top, see A21) and GSN 069 (middle, see M19), compared to those of eRO-QPE1 (bottom, see A21). The first two sources show alternating longer and shorter recurrence times and stronger and weaker bursts. In these plots (top and middle), eruptions of each source are overlapped in time-space by scaling each strong-weak-strong cycle (shown with a different color) with the related peak time of the first strong burst. The dashed vertical lines show the time corresponding to half of the separation between two strong bursts (i.e., between consecutive bursts aligned to zero), highlighting that the intermediate weak bursts occur systematically later in time (and that two recurrence times alternate). Bottom panel: the full light curve of eRO-QPE1 reported in A21 is here separated into chunks of two consecutive cycles and scaled in time at the peak of the odd eruptions to highlight a more complex distribution of recurrence and burst duration times.

In the text

Fig. 2. XMM-Newton background-subtracted light curve of the two observations of eRO-QPE1, named Obs1 and Obs2 (left and right panels, respectively). The full energy band was used, even though eRO-QPE1 shows little signal above ≳1.5 − 2.0 keV, even in the bright phase (A21). The two observations are separated by roughly a week; the times shown are scaled by the start time of each observation (with modified Julian date tObs1, 0 ∼ 59 057.834 and tObs2, 0 ∼ 59 065.988). Top panels: total model fitted (green line and 1σ shaded area), which is composed of one (Obs2) or five (Obs1) individual bursts (blue), each with a profile described by Eq. (1). A constant value for the quiescence is also added to the model. The inset in the top left panel shows the peak of the last eruption in Obs1, left out from the linear scale, which was adopted for clearer visualization of the fainter bursts. Bottom panels: ratio of the data to the model, where the shaded contours show the uncertainty in the model.

In the text

Fig. 3. XMM-Newton light curve of Obs2 extracted in small energy bins, represented with different colors (see legend). Each light curve is vertically shifted by a constant factor for visualization. The subpanels show the data-model residuals in each energy bin, color-coded accordingly; vertical dashed lines (with related 1σ uncertainties) indicate when the eruptions start in the given energy bin, while vertical solid lines identify the peak. The model adopted is described by Eq. (1), and it is shown in all panels with 1σ uncertainty as shaded intervals.

In the text

	Fig. 4. Energy dependence of QPEs in the observation Obs2 of eRO-QPE1 (related to the light curve shown in Fig. 3). The times tstart and tpeak are shown in hours elapsed from tObs2, 0. QPEs start earlier at lower energies (top panel) and evolve more slowly in terms of rise and decay (middle panel), and peak at later times (bottom panel) at low energies.
In the text

Fig. 5. Same as Fig. 3, but for Obs1. Individual burst profiles are shown in black, with different line styles. Five bursts are required for the 0.2 − 0.4 keV and 0.4 − 0.6 keV bins, while four are needed for 0.6 − 0.8 keV and three in the remaining energy bins (see Table A.1). Lower panels: data/model ratios. Vertical lines indicate tstart of each burst at all energies, with the line style matching the one used for the burst profiles in the upper panels with the light curves. Shaded intervals show 1σ uncertainties.

In the text

	Fig. 6. Same as Fig. 4, but for the first two bursts of Obs1 (see, e.g. Fig. 5). Squares are used in all panels for the first burst, stars for the second. The energy dependence of τrise, τdecay, and tpeak (from second to bottom panel) is consistent with Obs1 and other QPEs (M19; G20; Chakraborty et al. 2021).
In the text

Fig. 7. Plots showing results for Obs1 (resp. Obs2) of eRO-QPE1, on the left (right), represented with colored circles (triangles), for which darker to lighter colors with a red-to-yellow (blue-to-green) color map indicate the increase in time of the observation. Top central panels: count rate (CR) light curve as in Fig. 2, but here in logarithmic units and excluding the quiescence. Bottom panels: evolution of hardness ratio (HR = H/(H + S) = CR0.6 − 2.0 keV/CR0.2 − 2.0 keV) over time. External panels: HR vs. CR plot. In the panels related to Obs1, the central data points are highlighted with black squares for ease of comparison across the panels. The range of Obs2 in HR is shown as black dashed lines across the different panels. The HR reached at the peak of the eruptions is count rate dependent, and the eruptions undergo a counterclockwise hysteresis cycle in the HR vs. CR plot (top right).

In the text

Fig. 8. Hardness ratio vs. count rate computed from the modeled profiles of the bursts in Obs1 (color-coded as in the legend) and Obs2 (dashed black line). For a visual representation of the model profiles see Fig. 2. For Obs1 the profiles are used of bursts that are clearly detected in all the energy bands used for the definition of HR. The thick solid lines represent the median of the modeled tracks, while the thinner lines highlight the span within the model posteriors. A counterclockwise hysteresis (highlighted with arrows on the tracks) is found for all modeled eruptions, which is also what the pure data indicate (Fig. 7).

In the text

	Fig. A.1. As in the left panel of Fig. 2, but fitting a model with three (top) and four (bottom) bursts.
In the text

	Fig. A.2. As in the right panel of Fig. 2, but fitting a model with two bursts (shown in cyan and light blue).
In the text

	Fig. A.3. Eruption number four in the XMM-Newton light curve of eRO-QPE2 (see also Fig. 1). The red curve represents the fit with one Gaussian component, the blue curve the superposition of two Gaussian profiles is shown, and the cyan lines show the two individual Gaussian components (the median is shown with a thicker line).
In the text