© The Authors 2024

1. Introduction

Quasi-periodic eruptions (QPEs) are repeating soft X-ray bursts originating from galactic nuclei that last and recur on a timescale of hours (Miniutti et al. 2019; Giustini et al. 2020; Chakraborty et al. 2021; Arcodia et al. 2021, 2024a). Their unique characteristic is the soft X-ray spectrum during the eruptions, and its evolution with a harder rise than decay (Arcodia et al. 2022; Miniutti et al. 2023; Arcodia et al. 2024a). This peculiar behavior helps to distinguish their physical emission mechanism from that of other repeated nuclear transients, which are growing in number and flavor. To date, most models for the origin of QPEs relate to galactic nuclei stellar dynamics, although this is not the sole interpretation (e.g., Arcodia et al. 2024a, for a recent overview of the models). In particular, a popular scenario (e.g., Linial & Metzger 2023; Franchini et al. 2023; Tagawa & Haiman 2023; Zhou et al. 2024a,b) that we aim to test in this work involves a primary massive black hole (BH) with a compact and short-lived accretion flow (e.g., Patra et al. 2024), which in some cases may be fed by a tidal disruption event (Linial & Metzger 2023; Franchini et al. 2023, and references therein). This recent accretion event reveals a separate, preexisting, much smaller object (e.g., a star or a BH) on a low-eccentricity orbit, repeatedly ploughing through the disk twice per orbit and producing the observed QPEs. This interpretation has drawn particular interest since these systems are the so-called extreme mass ratio inspirals (EMRIs; e.g., Amaro-Seoane 2018) that would be detectable by future-generation gravitational-wave detectors such as the Laser Interferometer Space Antenna (LISA) and μAres (Sesana et al. 2021; Colpi et al. 2024). If this were the correct model, QPE volumetric rates would be the first ever observational constraint for EMRIs (Arcodia et al. 2024b).

In current datasets, the sparse light curve epochs of a given QPE source have gaps that are too long to maintain knowledge of both short-term and long-term evolution and, in particular, the number of eruptions occurring within the gaps, but some attempts have been made (e.g., Xian et al. 2021; Franchini et al. 2023; Chakraborty et al. 2024; Pasham et al. 2024; Zhou et al. 2024a,b). Here, we focus on the long-term (∼3.3 yr) evolution of eRO-QPE2, which has thus far demonstrated regular flaring patterns and overall stable emission properties. We find that the recurrence time decreases over time, albeit in a nonuniform way, while the long-short alternation seen in the first epoch (Arcodia et al. 2021) has not been significantly observed since. At the same time, the flux of both the eruptions and the quiescence emission has remained constant within uncertainties. We discuss possible mechanisms through which an EMRI system can reproduce these observations and the pitfalls. Our suggested scenario has clear predictions that future data will be able to verify.

2. Observations and analysis

Here, we discuss the main results. We report the timing and spectral analysis in detail in Appendix A, using four XMM-Newton datasets taken in August 2020, February 2022, June 2022, and December 2023, respectively. These are shown in Fig. 1.

Fig. 1. XMM-Newton 0.2–2.0 keV light curves of the four epochs, scaled by the arrival time of the first eruption in each observation. Time intervals contaminated by the flaring background are shown in grey. In the top panel, we additionally show the rise, peak, and decay phases selected for August 2020, as an example. In the second panel from the top, the first eruption is only partially resolved in EPIC-pn data and renormalized MOS1 data are shown with empty symbols. As the x-axis range is shared, we note the shortening of recurrence times over the ∼3.3 yr baseline, with vertical dashed lines showing the average arrival time interval of the previous epoch to guide the eye.

2.1. Decrease in the quasi-periodic eruption recurrence time with constant flux

The recurrence time (t_recur), defined as the peak-to-peak separation from light curve fits, was remarkable in the discovery epoch of August 2020 in its alternating pattern between consecutively longer and shorter recurrences with differences on the level of a few percent (Arcodia et al. 2021, 2022). Similar alternating patterns were noted in GSN 069 and RXJ 1301.9+2747 (Miniutti et al. 2023; Giustini et al. 2020). We show this in the leftmost panel in Fig. 2, which reports the evolution of t_recur. This alternating behavior has disappeared, within uncertainties, over the rest of ∼3.3 yr baseline probed by the data in this work (Fig. 2). At the same time, the mean recurrence time during a given epoch has decreased. The mean value (with the associated standard error of the mean) decreased from 8721.0 s (116.4 s) in August 2020 to 8468.4 s (149.3 s) in February 2022 and 8180.8 s (54.7 s) in June 2022. After that, it only slightly decreased further to 8128.6 s (59.5 s) in December 2023, which is still compatible within uncertainties with the previous epoch. Using the first and last epoch, we infer an apparent Δt_recur/Δt ≈ −6 × 10⁻⁶ s/s, or a ∼7% decrease compared to the August 2020 recurrence time. However, this decay was not uniform across all epochs and does not reproduce all of them. For instance, this recurrence decrease rate would overpredict the recurrence during June 2022 by ∼207 s, which is a difference around four times larger than the June 2022 uncertainty. As a matter of fact, there was a ∼6% decrease already between August 2020 and June 2022 (∼1.8 yr) and ≲1% (compatible with no decrease within uncertainties) in the following ∼1.5 yr (see, e.g., Fig. 2).

Fig. 2. Evolution of the recurrence time between eruptions over the ∼3.3 yr baseline. The start time t0 represents the start of the August 2020 observation. Data points with an orange contour contain an additional systematic uncertainty, as is described in the text. The mean value in each epoch, with an associated standard error of the mean, is shown with a dashed line and shaded contours.

Remarkably, all other properties of both flares and quiescent emission have remained constant within uncertainties and the available observations. In particular, the eruptions’ duration remain consistent within the standard deviation of the values of each epoch: the mean and standard deviations are 1666 ± 143 s, 1603 ± 36 s, 1741 ± 145 s, and 1736 ± 151 s, for the different epochs, respectively. Here, we defined the rise-to-decay duration of eruptions from the difference between the times at which the Gaussian model adopted for the fit (Appendix A) reaches a value of 1/e³ compared to the peak. Similarly, the flux and temperature of both the eruptions and the quiescent components have remained constant within 3σ uncertainties. We show detailed fit values in Table A.1 and the fit posteriors of all components in Fig. 3, for flux and temperature in the left and right panels, respectively. It is perhaps noteworthy that the QPE peak temperature marginally (within 1σ) increases from the first to the last epoch, although we stress the consistency of the measurements at 2σ. Future data on much longer baselines will perhaps confirm or rule out this tentative trend. We computed the energy emitted by the eruptions (e.g., as in Giustini et al. 2024) as the integral of the Gaussian model used for the light curve fit; namely, $\sqrt{2\pi }{L}_{\mathrm{bol},\mathrm{QPE}}{\sigma }_{\mathrm{QPE}}$, where σ_QPE is the best-fit Gaussian standard deviation (i.e., the characteristic timescale of the QPE duration). The resulting median (and related 16th–84th percentile) QPE energy is $1.2_{-0.2}^{+0.2}\times {10}^{45}$ erg, $1.0_{-0.3}^{+0.3}\times {10}^{45}$ erg, $1.1_{-0.2}^{+0.2}\times {10}^{45}$ erg, and $1.2_{-0.2}^{+0.3}\times {10}^{45}$ erg, respectively from August 2020–December 2023. Thus, it is consistent across the epochs of our ∼3.3 yr baseline.

Intriguingly, in the August 2020 epoch an additional component, harder than the accretion disk, but still rather soft, is required to account for residuals (Appendix A). This component is only marginally constrained, albeit still statistically required, in Feb. 2020, but it is not statistically required after that (Fig. 3). The nature of this component, occasionally seen in the quiescence of QPE sources (e.g., Giustini et al. 2020; Chakraborty et al. 2021; Arcodia et al. 2024a), is still unclear. Here, we modeled it as Comptonization of accretion disk photons, but this interpretation is assumed and not driven by data. In this work, given the low signal-to-noise (S/N), we do not attempt further interpretation and instead defer this to future work on all the QPE sources.

Fig. 3. Evolution of the emission properties. Left: Evolution of the 0.2–8.0 keV luminosity of individual spectral components (as is labeled in the top panel). Histograms represent the fitted flux chains, converted to luminosity and corrected for absorption. Vertical lines and shaded area (following the same coding) highlight the median and 16th–84th interquintile range. There is no significant evolution either of the accretion disk (dashed) or the eruptions (dash-dotted), while the additional spectral component required for the quiescence spectrum in the first epoch becomes less significant and ultimately undetected at later epochs. Right: Same as the left panels, but showing posterior chains of the fitted temperature values. There is no significant evolution of the accretion disk temperature (solid) over the ∼3.3 yr baseline. For the eruption temperature (dash-dotted), despite the median increasing over time and the first and last epochs not being compatible at 1σ, all values are consistent within 2σ.

2.2. The energy dependence of the eruptions

Arcodia et al. (2022) reported the evolution of the spectral shape of eruptions in eRO-QPE1, which shows a harder spectrum during the rise than during the decay. Since then, this behavior has been found in the other QPE sources (Miniutti et al. 2023; Giustini et al. 2024; Arcodia et al. 2024a), which is particularly useful to identify a common emission process among the growing number of repeating nuclear transients. Here, we test the same energy dependence in eRO-QPE2. Due to the lower S/N of eRO-QPE2 compared to other sources, we extracted a single rise and decay phase in addition to quiescence and peak, and we combine all eruptions within a single epoch. Fig. 1 shows how the different phases were selected during the August 2020 epoch, as an example, and Fig. A.2 shows the related spectra of the different phases. An analogous methodology was followed for the other epochs. After spectral analysis (Appendix A) we could isolate the eruption component and investigate the spectral evolution of the thermal component during the eruptions. Here, we confirm that the same spectral evolution is observed in eRO-QPE2 as well, consistently during each epoch (Fig. 4). In the top right (medium right) panel of Fig. 4 we show the QPE bolometric luminosity L_bol, QPE (the temperature kT_QPE) evolution with respect to the eruption phase. We note that, comparably to the other known QPE sources, the peak temperature is reached before peak luminosity, after which the spectrum cools down and becomes the coolest during decay. The top left panel of Fig. 4 shows the luminosity-temperature coevolution during the various epochs. Darker to lighter colors represent the evolution from start to end of the eruptions, color-coded as in the legend as a function of the epoch. Indeed, the rise is always harder than the decay and, albeit with fewer data points, the hysteresis cycle is confirmed for eRO-QPE2 as well. Effectively, a thermal spectrum transiting from a hotter state to a colder one at roughly comparable brightness could be attributed to an expansion of the associated blackbody radius (Miniutti et al. 2023; Chakraborty et al. 2024). In eRO-QPE2, we show this trend with time – R_em ∝ (L_bol, QPE/T_QPE⁴)^1/2 – in the bottom right panel of Fig. 4. Intriguingly, despite the increase in R_em from rise to decay being seen consistently at all epochs, the rate of increase is larger in later epochs. In particular, the median rise-to-decay expansion of R_em is a factor of ∼2.4, ∼3.3, ∼4.3, and ∼7.0, going from August 2020–December 2023. Since the phase-folded light curve profiles appear to be consistent between August 2020 and December 2023 (Appendix A and Fig. A.1), this is unlikely to be due to an evolution in the rise and decay phases over time, or biases due to how these phases are defined.

Fig. 4. Spectral evolution of QPEs. Darker to lighter colors represent the evolution from the start to end of the eruptions, color-coded as in the legend as a function of the epoch. The top right (medium right) panel shows the QPE bolometric luminosity (temperature) evolution with respect to the eruptions phase. The top left panel shows their coevolution. The bottom right panel shows the evolution of the emitting radius, assuming a thermal spectrum and spherical geometry.

3. Astrophysical interpretation

The two main results from the long-term timing analysis of eRO-QPE2 are (1) the detection of a nonconstant decrease in the QPE recurrence time t_recur and (2) an apparent disappearance of the long-short alternation pattern over the ∼3.3 yr observational baseline (Fig. 2).

Here, we attempt to put these results into the context of orbital models. We note that orbital models have already been tested on eRO-QPE2. For instance, Franchini et al. (2023) analysed the August 2020 epoch of eRO-QPE2 as a test case for their model and Zhou et al. (2024b) tested the 2020 and 2022 epochs independently, and, similarly to Linial & Metzger (2023), suggested gas drag as the energy loss mechanism driving the apparent period decrease. However, the December 2023 epoch (which highlights that a constant period decrease cannot reproduce the data, Fig. 2) was not included in their work; thus, their conclusions do not grasp the full behavior shown by eRO-QPE2. A further complication is that without model-dependent assumptions timing fits are only possible for individual epochs of current datasets, the main obstacle being the lack of knowledge of the number of eruptions that occurred within the gaps. This is why we are unable to distinguish whether Fig. 2 represents an overall continuous decreasing trend or whether the recurrence time and long-short alternation appearance and disappearance have also significantly oscillated over time. Thus, we discuss both options in Sect. 3.2 and 3.3. First, we outline in Sect. 3.1 the main model assumptions and equations used in this work.

3.1. Orbital model assumptions and main equations

The interpretation we aim to test is that QPEs are triggered by a secondary orbiter of much smaller mass with respect to the primary BH, which has M_• ∼ 10⁵ M_⊙ (Wevers et al. 2022); thus, we take M_•, 5 = M_•/10⁵M_⊙ as a reference. The secondary is on a low-eccentricity orbit and is considered to be a either a solar-mass star, à la Linial & Metzger (2023), or an ∼40–100 M_⊙ BH, à la Franchini et al. (2023). For simplicity, we refer to the orbital periodicity, P_orb, and its possible decrease, Ṗ_orb, to explain the behavior shown by the observed quantity, t_recur. Even assuming two collisions per orbit, we note that P_orb is not exactly t_recur, n + t_{recur, n + 1}. It would require, for instance, an idealized system with a fixed disk, the absence of precession, and that flares being produced promptly at disk passage. Moreover, we refer to orbital eccentricity, e_orb, and its possible decrease, ė_orb, to explain the initial presence (August 2020) and subsequent disappearance (June 2022 and December 2023), of the long-short alternation. However, the observed alternation actually probes a quantity that is an approximation of e_orb; namely, e ≡ δt_recur/⟨t_recur⟩, the difference between two consecutive recurrence times with respect to the average. For the scope of this work, we do approximate that trends in t_recur reflect trends in ∼1/2 P_orb and that e traces e_orb. Hence, with this in mind Fig. 2 unveils an apparent decrease for both P_orb and e_orb.

3.1.1. The relevant timescales

The relevant timescales at play are the apsidal precession time of the secondary:

$$\begin{array}{c}\hfill {\tau }_{\star ,\mathrm{ap}}\approx \frac{P_{\mathrm{orb}}}{3}\left(\frac{r_0}{R_{\mathrm{g}}}\right)\approx 22\mathrm{d}{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^{5/3}{M}_{\bullet ,5}^{-2/3},\end{array}$$(1)

and, for nonzero spin of the primary BH, the nodal precession timescale due to the Lense-Thirring (Lense & Thirring 1918) effect of the spinning BH on a misaligned orbit:

$$\begin{array}{c}\hfill {\tau }_{\star ,\mathrm{LT}}\approx \frac{1}{2a_{\bullet }}{P}_{\mathrm{orb}}{\left(\frac{r_0}{R_{\mathrm{g}}}\right)}^{3/2}\approx 1.6\mathrm{yr}{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^2{M}_{\bullet ,5}^{-1}{a}_{\bullet }^{-1},\end{array}$$(2)

where a_• is the massive BH’s dimensionless spin, P_orb = 2⟨t_recur⟩ is an approximation of the orbital period, r₀ is the orbit’s semi-major axis,

$$\begin{array}{c}\hfill r_0\approx 4.8\times {10}^{12}\mathrm{cm}{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^{2/3}{M}_{\bullet ,5}^{1/3}\approx 323R_{\mathrm{g}}{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^{2/3}{M}_{\bullet ,5}^{-2/3},\end{array}$$(3)

and R_g = GM_•/c² is the gravitational radius. An additional timescale could arise due to disk precession, although its effect depends on the exact model assumptions about its radial extent and mass distribution. Franchini et al. (2016) have shown that a compact, radiation-pressure-dominated accretion disk can rigidly precess around a spinning BH on timescales as short as days, shorter than the apsidal precession of the companion. For instance, in Franchini et al. (2023) the compact disk (e.g. as a result of a tidal disruption event, TDE) around a primary with a spin of a_• = 0.5 is misaligned by ι_d = 10° with respect to the x − y plane and precesses rigidly on a period of ∼5.7 days. However, other more complex geometries may exist (e.g., Raj & Nixon 2021; Ivanov & Zhuravlev 2024); for instance, with warps or tilts around the collision radii, which would complicate the dynamics even further. Motivated by the relative regularity of QPEs in eRO-QPE2, we adopt a simple scaling between the nodal precession of the secondary (τ_⋆, LT, Eq. 2) and that of the disk as τ_d, LT = α_LTτ_⋆, LT, where α_LT is a proportionality constant.

3.1.2. Mechanisms for orbital energy loss

For a BH, neglecting the interaction with the disk, a low-eccentricity orbit inspiraling due to the emission of gravitational waves (GWs) would yield a period derivative of the order

$$\begin{array}{cc}\hfill {\dot{P}}_{\mathrm{QPE}}{|}_{\mathrm{GW}}& \approx -\frac{96{\pi }^{8/3}}{5}\frac{m_{\star }}{M_{\bullet }}{(G{M}_{\bullet }/⟨t_{\mathrm{recur}}⟩>)}^{5/3}{c}^{-5}\hfill \\ \hfill & \approx -3\times {10}^{-8}{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^{-5/3}\left(\frac{m_{\star }}{100{\mathrm{M}}_{\odot }}\right)M_{\bullet ,5}^{2/3}.\hfill \end{array}$$(4)

For a fixed period (or semi-major axis), an eccentric orbit would undergo more rapid inspiraling (${\dot{P}}_{\mathrm{GW}}\propto {(1-e)}^{-7/2})$). However, the relatively mild eccentricity inferred, e ≈ 0.06 at most during August 2020, would only increase Ṗ by less than a factor of two relative to a circular orbit. In the case of a star, Ṗ_QPE is linearly lower with the mass.

For a star, energy losses due to hydrodynamic gas drag are more relevant (e.g., Linial & Metzger 2023; Linial & Quataert 2024). A star of mass m_⋆ and radius R_⋆, impacting a disk of surface density Σ twice per orbit on a highly inclined, nearly circular orbit, would experience an orbital decay at an approximate rate of

$$\begin{array}{c}\hfill ⟨\dot{P}⟩\approx -\frac{3\pi R_{\star }^2\mathrm{\Sigma }}{m_{\star }}\approx -2\times {10}^{-6}{\left(\frac{m_{\star }}{{\mathrm{M}}_{\odot }}\right)}^{-1}{\left(\frac{R_{\star }}{{\mathrm{R}}_{\odot }}\right)}^2\left(\frac{\mathrm{\Sigma }}{{10}^5\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2}}\right).\end{array}$$(5)

Naturally, the prediction from Eq. 5 is sensitive to the largely unknown disk structure and surface density. Since the observed decrease in recurrence time was not uniform over the observed baseline (Fig. 2), disfavoring a constant Ṗ throughout the system’s evolution, we further discuss the effect of inclination changes. As a toy model, we consider a star on a nearly circular orbit around a massive BH in the presence of a rigidly precessing accretion disk, both inclined with respect to the SMBH spin, with inclination angles ι_⋆ and ι_d, respectively. The relative inclination between the star and the disk thus varies periodically over the beat cycle, τ_LT, ⋆/|1 − α_LT|, in the range ι_rel ∈ (|ι_⋆ − ι_d|,ι_⋆ + ι_d). At a given time, when the relative inclination between the star and the disk is ι_rel, the orbital period decays at a rate of

$$\begin{array}{c}\hfill \dot{P}({\iota }_{\mathrm{rel}})\approx -\frac{6\pi \mathrm{\Sigma }{R}_{\star }^2}{m_{\star }}\frac{{sin}^2({\iota }_{\mathrm{rel}}/2)}{sin{\iota }_{\mathrm{rel}}}.\end{array}$$(6)

Here we assumed that the collision cross section is set by the star’s physical size, and neglected any enhancement due to gravitational focusing of the impacted disk material (e.g., Rein 2012; Arzamasskiy et al. 2018)¹. In the limit of small inclinations (10⁻² ≲ ι_rel ≪ 1), $\dot{P}\propto {\iota }_{\mathrm{rel}}$. We further note that the orbital decay scales linearly with the dimensionless combination η ≡ ΣR_⋆²/m_⋆, which we normalize to η ≡ 10⁻⁶ η₋₆.

3.1.3. Light travel time effects

Since the August 2020 epoch showed significant long-short alternation, but following epochs did not, here we discuss the effect of light travel time delays. For instance, in Franchini et al. (2023) these effects were included and considered small and negligible for most inclinations. However, the magnitude of this effect is maximal when the line connecting the two star-disk collision sites is near alignment with the observer’s line of sight (requiring that both the disk and orbital plane are viewed nearly edge-on). Variations in the recurrence time in this configuration for a circular stellar orbit, are of the approximate magnitude

$$\begin{array}{c}\hfill e=\frac{\delta t_{\mathrm{recur}}}{⟨t_{\mathrm{recur}}⟩}\approx \frac{2r_0/c}{t_{\mathrm{recur}}}\approx 0.044{\left(\frac{⟨t_{\mathrm{recur}}⟩}{2.5\mathrm{hr}}\right)}^{-1/3}{M}_{\bullet ,5}^{1/3},\end{array}$$(7)

similar to the observed e ≈ 0.05–0.06 in August 2020.

3.2. On whether the behavior shown in Fig. 2 is dominated by random sampling of the orbital phase

Here, we discuss the scenario in which the evolution shown in Fig. 2 for both ${\dot{P}}_{\mathrm{orb}}$ (decrease in t_recurr) and ${\dot{e}}_{\mathrm{orb}}$ (decrease in alternation, or its disappearance) is dominated by modulations and random sampling of different orbital phases. These could be imprinted by changes in the relative inclinations between the orbit, the disk and the observer due to the precession frequencies at play, and light travel time effects, as is shown in Sect. 3.1.

We remind the reader that this cannot be inferred by data in a model-independent fashion due to the gaps between observations. Thus, we discuss here all the possible impacts of these effects on the observations (Fig. 2). Even in the case of a non-spinning primary, apsidal precession of the companion object would modulate the periodicity; that is, naturally, if the orbit were not circular to start with. For now, we assume that August 2020 yields an apparent eccentricity close to the orbital one, and thus on the order of ≈0.05–0.06, and that the following epochs show a much lower variation (≲0.01) due to sampling effects. Considering the apsidal precession timescale on the order of ≈20 d (Eq. 1), alternations of ≲0.01 are expected to occur for roughly 0.5 d, twice per precession cycle, corresponding to a probability of ∼5% if the source is randomly sampled (given by 2 × (0.01/0.06)/(2π)≈0.05, the fraction of time in which the orbital eccentricity vector is near alignment with the plane of the disk). To have both June 2022 and Dec. 2022 as a relatively rare aligned disk-orbit phase is even rarer (∼0.25%), but formally possible. In this case, a strong modulation would not change the average P_orb from epoch to epoch; thus, orbital decay has to be invoked (see Sect. 3.3).

If the primary BH were spinning and the orbit and the disk were misaligned, both would evolve through nodal precession. As is shown in detail in Sect. 3.1, we estimate the former to be τ_⋆, LT ≈ 1.6 yr and the latter to be τ_d, LT = α_LTτ_⋆, LT, where α_LT is a proportionality constant smaller than one. The effect of τ_⋆, LT, since it is much longer and comparable to the baseline, is discussed in Sect. 3.3. Here, we discuss the possible sampling effect due to the shorter τ_d, LT. We first investigate whether disk alignment over the probed baseline of ∼3.3 yr could induce the observed decrease or disappearance of long-short alternation. Using the framework from Franchini et al. (2023), we find that the alignment process does not lead to a significant change in the recurrence pattern². We point out that the interplay between apsidal precession of the secondary and disk precession may give rise to more prolonged phases of low apparent eccentricity (compared to what we estimated above for when disk precession is absent), enhancing the probability of observing those phases by chance due to random sampling (Miniutti et al. in prep.). For instance, we evolve the solution for eRO-QPE2 from Franchini et al. (2023) for ∼3 months, in order to get enough cycles, with the goal being to estimate the likelihood of catching phases with low apparent eccentricity. Adopting ${\tau }_{\star ,\mathrm{ap}}\sim 20$ d, τ_d, LT ∼ 6 d, P_orb ∼ 5 h, and e_orb ∼ 0.05 (Franchini et al. 2023), we obtain that ≈12% of the time the system is found at an apparent eccentricity ≲0.02 and that ≈6% of the time it is found at one lower than ≲0.01. Similarly, but following the opposite reasoning, the August 2020 epoch could be the one in which we caught the system in a rare phase of the orbit with high apparent eccentricity, which is otherwise low. Following the same calculation, we estimate that ≈18% of the time the system is found at an apparent eccentricity ≳0.10. Thus, one can imagine that an intermediate e_orb between the values estimated for eRO-QPE2 could serendipitously give rise to samples with both higher and lower apparent eccentricities, such as the ones observed. In the same way, we estimate the impact of sampling this evolved model (Franchini et al. 2023) on the observed quantity t_recur, n + t_{recur, n + 1}, which only seemingly traces the real P_orb. Assuming that P_orb = 5 h, we note that the model predicts that ∼25% of the time t_recur, n + t_{recur, n + 1} is found in the range of 5.0–5.4 h, and that ∼25% of the time it is found in the range of 4.9–5.0 h. Therefore, in principle sampling effect of the orbital phase of a BH-EMRI system such that in Franchini et al. (2023) would be found at apparent periods significantly different than the true P_orb. However, at least in the phases with a higher apparent period, consecutive recurrence times would appear to monotonically increase or decrease within a typical XMM-Newton observation, which is instead not observed here. Furthermore, we note that if the disk were rigidly precessing (e.g., Franchini et al. 2023), we should see modulation in the flux of the quiescence component, which is assumed to come from the innermost radii of the disk. However, we do not, as the quiescence flux remains remarkably constant over the ∼3.3 yr probed (Fig. 3). On the other hand, no significant flux change would be expected if the inner disk (where the quiescence comes from) were aligned and the outer (where the collisions happen) were not, or if there were large differences in how the inner and outer disk precess.

Finally, as the relative inclination between orbit, disk, and observer changes due to the precession of both disk and orbit, light travel effects change. As we show in Sect. 3.1, the effect is in general small, although for edge-on inclinations and short observed periods like in eRO-QPE2 the effect could be non-negligible (∼4% of the QPE recurrence). Therefore, during the August 2020 epoch we could have found the system in temporary edge-on alignment (i.e., light travel effects of a few percent) and the following ones at higher inclinations (i.e., negligible light travel effects). As this evolution is expected to occur over a timescale of ∼ years (comparable to τ_⋆,LT), the observed trends from August 2020–December 2023 appear to be in overall agreement with this proposed interpretation.

We conclude that all of these sampling effects are reasonable, but hard to estimate with the available data, thus any explanation would inherently be fine-tuned at this stage. However, we note that these sampling effects are required to a larger extent if the orbiter is a BH, compared to a star. As is discussed in Sect. 3.1 a star may lose more orbital energy than a BH due to interactions with the disk, while for a BH the leading energy loss mechanism would be GW emission, which is much weaker than the former and unable to explain current data alone. Quite crucially, this scenario in which sampling dominates the observed long-term evolution inherently predicts that we shall likely have observations with a larger inferred recurrence time compared to December 2023 (i.e. positive apparent Ṗ_orb) and/or with large long-short alternation again (i.e. positive apparent ė). This can be tested indirectly with future observations.

3.3. On whether the behavior shown in Fig. 2 reflects a long-term trend

Here, we discuss the alternative scenario in which the evolution shown in Fig. 2 for both Ṗ_orb (decrease in t_recur) and ė_orb (decrease or disappearance of alternation) is dominated by a long-term trend. As much as some sampling effects due to precessions and inclination are expected (see Sect. 3.2) they would be second-order effects in this scenario. We observed a decrease in P_orb from ∼4.85 h to ∼4.53 h from August 2020–December 2023 (summing two consecutive t_recur and averaging per epoch), corresponding to r₀ = 320 R_g and r₀ = 307 R_g, respectively (e.g., Franchini et al. 2023), and e_orb decreased from ∼0.05–0.06 to ≲0.01. It is easy to picture these two separate solutions separately, although their evolution from one to the other should be modeled within the same framework. Taking the apparent period decrease from August 2020–December 2023 at face value, we need a mechanism that induces Ṗ_orb ≈ −6 × 10⁻⁶. As is shown in detail in Sect. 3.1, a BH orbiter would be disfavored since there is no mechanism to reproduce the observed decrease. Instead, for a star a modest orbital decay has been predicted to occur due to hydrodynamic drag; for instance, as proposed in Linial & Metzger (2023) and Zhou et al. (2024b). Predictions do not exactly reproduce the observed value, but they are sensitive to the largely unknown disk structure and surface density (Sect. 3.1).

What is instead harder to reproduce with this interpretation is the fact that the period decrease was nonconstant. As a matter of fact, there was a decrease of ∼6% between August 2020 and June 2022 and of ≲1% in the following ∼1.5 yr (Sect. 2, Fig. 2). In Sect. 3.1, we extended the hydro-drag interpretation with the effect of changing the relative inclination between the orbit and the disk due to τ_⋆, LT and τ_d, LT. In this framework, less orbital energy is lost when collisions are more face on, and viceversa. Here, we took Eq. (5) and for a given set of parameters {t₀,P₀,η,M_•,a_•,ι_d,ι_⋆,α_LT}, we calculated the orbital period $P(t)=P_0+{\displaystyle {\int }_{t_0}^t}\dot{P}(t\prime )\mathrm{d}t\prime$. Fig. 5 shows possible recurrence time evolution trends for a few different sets of parameters, as is described in the legend, alongside the mean recurrence times at the four observed epochs. We fixed M_∙ = 10⁵ M_⊙, ι_d = 0.8 and ι_⋆ = 0.78 (both approximately 45°) and manually adjusted the precession phase and P₀ to roughly align with the observations. Variations in the (negative) slope of the calculated curves correspond to variations in the relative inclination, such that the plateaus (low Ṗ) correspond to when the star and disk are near alignment. The duration of these plateaus (i.e., the time spent around the minimal ι_rel) is proportional to the beat cycle of the star and disk nodal precession, which scales ∝a_•⁻¹|1 − α_LT|⁻¹. All sets of parameters lead to a monotonic decay in the orbital period, and predict shorter recurrence times in the future. The transition between phases of rapid and slow decline could thus be attributed to the combined nodal precession of the star and the disk, coupled with hydrodynamic drag acting upon the star as it crosses the disk. We highlight that we do not show a fit to the data in Fig. 5, but that we do qualitatively reproduce the presence of nonconstant Ṗ within reasonable ranges of the parameters above. We also note that other effects are in place when the inclination changes (Sect. 3.2).

Fig. 5. Observed recurrence time evolution, overplotting some predicted evolution tracks due to gas drag and slowly evolving disk precession. Different lines are shown for different sets of parameters (as is shown in the legend), such as the spin (a•), the proportionality constant between disk and star nodal precessions (αLT), and the parameter η ≡ ΣR⋆2/m⋆. Variations in the slope of the calculated curves correspond to variations in the relative star-disk inclination, such that the plateaus (low Ṗ) correspond to when the star and disk are near alignment.

Even if this inclination effect is significant enough to explain the nonuniform Ṗ evolution, the related rate of eccentricity evolution is inconsistent with the relatively modest period evolution, if both orbital decay and circularization are due to the same hydrodynamic drag. As a matter of fact, we would expect in the low e limit that ė/e ∼ Ṗ/P (e.g., Appendix B of Linial & Metzger 2023). Hence, unless we are missing out a fundamental efficient circularization process, some sampling effects such as temporary alignment of (only) the major axis with the disk orbital plane are required to increase the apparent eccentricity in August 2020 compared to a true low e_orb, as is discussed in Sect. 3.2.

Finally, we note that a concrete prediction in case the apparent evolution of Fig. 2 is mostly an intrinsic period decrease trend, is that the mean QPE recurrence time should continue to decrease. If this is due to the above-mentioned drag picture, the recurrence time will decrease following alternating steeper and flatter epochs, following the trend of the relative star-disk inclination as shown in Fig. 5. Moreover, if the orbit has indeed circularized (regardless, for now, of the mechanism through which this has happened), no alternation will ever reappear. Conversely, if the orbit has not circularized and it is only an apparent transient alignment effect, long-short alternation phases may continue to reappear (and disappear) in future observations.

3.4. On the constant flux of quiescence and eruptions

Another key finding of this study is the striking consistency between the flux of both the quiescence and eruption components across different epochs (Fig. 3). The lack of significant disk evolution over this period is seemingly puzzling, if the origin of the disk is the accretion flow that follows the TDE of another star in the same galactic nucleus (e.g., Linial & Metzger 2023; Franchini et al. 2023), as the late time X-ray emission from TDEs typically falls off on timescale of years (e.g., Mummery & Balbus 2020). However, as was recently demonstrated by Linial & Metzger (2024), mass (and thermal energy) added to the disk through ablation of the star’s outer layers could result in a steady-state, equilibrium evolution of the disk, with its accretion rate matching the star’s stripping rate, even if the initial disk was seeded by a TDE. Considering a low-mass star of m_⋆ = 0.5 M_⊙, ⟨t_recur⟩ = 2.5 hr, M_∙ = 10⁵ M_⊙ and applying results from Section 3.2 of Linial & Metzger (2024), we find that the self-sustained disk should settle to a steady-state accretion rate of roughly ∼ 0.1 Ṁ_Edd and a surface density of about Σ_d ≈ 10⁴ g cm⁻², and have an expected lifetime on the order of 10² yr. Hence, the lack of flux evolution in eRO-QPE2 over only ∼3.3 yr is not a concern for this interpretation.

Next, we address the lack of substantial evolution in the flare temperature and luminosity over the observed baseline. If indeed the stability of the quiescent emission implies that the disk conditions remained nearly constant, an emission model that invokes collisions with the disk as the origin of the flares (e.g., Linial & Metzger 2023; Franchini et al. 2023) would indeed conform with stable flare properties (i.e., the star encounters similar disk conditions in different epochs). However, as variations in the relative star-disk inclination may unavoidably be invoked to explain the system’s nontrivial timing behavior (Sect. 3.3 and Sect. 3.1), the flare emission properties should be fairly insensitive to relative inclination, providing an important constraint on QPE emission models. Specifically, collisions of both low and high inclinations should produce flares of similar properties.

Finally, in Sect. 2 we obtained that the rate of increase in the blackbody radius from the rise to decay of the eruptions is apparently larger going from August 2020–December 2023 (i.e. from a factor of ∼2.4 to one of ∼7.0; Fig. 4). This might be expected if the ejecta are not spherical like the blackbody R_em calculation assumes, and if the system has evolved through changes in the inclination between the orbit and disk (hence collision and ejecta vector) and the observer. However, testing this effect quantitatively is beyond the scope of this work and the capabilities of the current emission models. We also stress that R_em does not necessarily represent the physical size of the X-ray emitting gas producing the flares. Specifically, in the context of the emission model discussed in Linial & Metzger (2023), in which the flares are generated by an expanding optically thick, shock-heated disk material, R_em does not trace the actual size of the emitting ejecta cloud (Vurm et al. in prep.).

3.5. On the detectability of quasi-periodic eruptions with gravitational wave detectors

QPEs have so far been identified only in the electromagnetic domain through strong X-ray flares. The conjecture of QPEs being caused by EMRI systems, however, opens the possibility of multimessenger observations of these sources. Indeed, EMRI systems are primary targets of low-frequency GW interferometers, like the recently adopted LISA mission (Colpi et al. 2024), planned to come online in the mid-2030s and covering the mHz band, as well as other mission concepts targeting even lower frequencies, such as μAres surveying the μHz band (Sesana et al. 2021). EMRI systems that could explain the detected QPEs are systems that generally are from hundreds to several thousands of years away from the coalescence. EMRI sources producing QPEs with a period of ∼2–20 hours emit GWs at frequencies between 10⁻⁵ and 10⁻⁴ Hz. This frequency range represents the lower end of LISA’s window, at which the sensitivity of the detector has undergone significant degradation. Nevertheless, the proximity of QPE sources, yielding a strain that can reach the level of ∼10⁻¹⁸–10⁻¹⁷ in the case of a BH companion (∼40–100 M_⊙), can compensate the LISA sensitivity loss. We can therefore expect S/N of order unity in the most optimistic cases. The situation could firmly improve with a detector targeting μHz frequencies, like μAres. There, the S/N of QPE sources is expected to be much higher and on the order of ∼10–100, which would make them potentially detectable even in the case of a stellar companion, depending on the exact mass of the central BH and distance to the source.

4. Summary

In this work, we have reported the long-term evolution of the QPE source eRO-QPE2 using four epochs over the 3.3 yr baseline since its discovery (Arcodia et al. 2021). The main observational results are:

• The average QPE recurrence time per epoch has decreased over time, albeit not uniformly: there was a decrease of ∼6% during the first ∼1.8 yr and of ≲1% (compatible with no decrease within uncertainties) in the following ∼1.5 yr (Fig. 2 and Sect. 2.1).

• The long-short alternating behavior found at discovery has not been observed since. In particular, it was only minor ∼0.8 yr after discovery (although only inferred with three consecutive eruptions) and it has been consistent with being absent since then (Fig. 2 and Sect. 2.1).

• The QPE duration has remained constant within uncertainties across all epochs (Sect. 2.1).

• The spectral properties namely, the flux and temperature of both eruptions and quiescence components have remained surprisingly consistent within 3σ uncertainties throughout all epochs (Fig. 3 and Sect. 2.1). The peak QPE temperature might have slightly increased going from the first to last epoch, although only at 1σ.

• The energy dependence during eruptions follows the known QPE trend (Arcodia et al. 2022; Miniutti et al. 2023; Arcodia et al. 2024a; Giustini et al. 2024), which shows a harder rise than decay. Modeling eruptions with a thermal spectrum, this leads to an increase in the fitted blackbody radius from rise to decay (Fig. 4). The rate of this increase appears to have changed over time, toward larger increases at later epochs (Sect. 2.2).

We attempted to interpret these results within orbital model prescriptions (Sect. 3.1), in particular the scenario of a secondary object (e.g., a star or a BH with a much smaller mass than the primary) on a low-eccentricity orbit that repeatedly pierces through the disk twice per orbit and produces the observed QPEs (e.g., Linial & Metzger 2023; Franchini et al. 2023). In this framework, the QPE recurrence decrease is interpreted as a (real or apparent) period decrease, Ṗ, and the decrease or disappearance of long-short alternation as a (real or apparent) eccentricity decrease, ė. However, since gaps between observations are too long we are not able to distinguish between an evolution dominated by just a decreasing trend and one dominated by modulations (e.g. due to the orbit and disk precession frequencies at play). We therefore discuss these two scenarios separately (Sects. 3.2 and 3.3). If the current observations are dominated by unconstrained modulations no conclusive remarks are possible on the orbital evolution and the nature of the orbiter (Sect. 3.2). If instead the observed evolution (Fig. 2) is dominated by an intrinsic decreasing trend in Ṗ and ė (Sect. 3.3), observations can be qualitatively reproduced by a star losing orbital energy due to hydrodynamic gas drag from disk collisions. The time-dependent apparent Ṗ would be induced by a change in relative inclination between the orbit and the disk due to the orbiter’s nodal precession, which for eRO-QPE2 is in the ballpark probed by our baseline (Fig. 5). However, the eccentricity decrease would be too fast to be only due to the same hydrodynamic gas drag regardless of the nature of the companion; thus, some sampling effects (Sect. 3.2) would need to be invoked. If the behavior shown in Fig. 2 is indeed dominated by a trend, a BH would not be able to lose enough orbital energy (through either GW emission or gas drag) to explain the observed decrease. A BH orbiter interpretation does require the presence of modulations in both apparent period and eccentricity due to precession frequencies, which are however possible and which remain untested with current data.

Fundamentally, these two cases (i.e., observation dominated by an intrinsic trend versus the effect of sampling different orbital phases, Sects. 3.2 and 3.3) come with very distinctive predictions for future X-ray data. In particular, in the case of a dominating decreasing trend we expect all future observations to show a shorter recurrence time than the latest epoch. Conversely, in the case of strong modulations dominating the observed evolution, we would expect at least some future observations to be found with a longer recurrence time, and hence an apparent temporary period increase compared to the last epoch shown in this work. Future X-ray observations over the next months and years will be able to indirectly discern between these two scenarios.

Acknowledgments

We thank the anonymous referee for their positive report. 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. IL acknowledges support from a Rothschild Fellowship and The Gruber Foundation and from Simons Investigator grant number 827103. GM was supported by grant PID2020-115325GB-C31 funded by MICIN/AEI/10.13039/501100011033. MG is supported by the “Programa de Atracción de Talento” of the Comunidad de Madrid, grant numbers 2018-T1/TIC-11733 and 2022-5A/TIC-24235. MB acknowledges support provided by MUR under grant “PNRR – Missione 4 Istruzione e Ricerca – Componente 2 Dalla Ricerca all’Impresa – Investimento 1.2 Finanziamento di progetti presentati da giovani ricercatori ID:SOE_0163” and by University of Milano-Bicocca under grant “2022-NAZ-0482/B”. AS acknowledges the financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691). GP acknowledges financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program HotMilk (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 research benefited from interactions at the Kavli Institute for Theoretical Physics, during the KITP TDE24 program, supported by NSF PHY-2309135.

References

Appendix A: Data processing and analysis

Four proprietary or publicly available observations taken with XMM-Newton were analyzed in this paper. The observations IDs are 0872390101, 0893810501, 0883770201, and 0931791301, and are referred to as ‘Aug. 2020’, ‘Feb. 2022’, ‘Jun. 2022’, and ‘Dec. 2023’, respectively, throughout this work. Data were reduced using SAS v. 20.0.0 and HEAsoft v. 6.31. Products were extracted from source (background) regions selecting a circle of ∼30 − 40" centered on eRO-QPE2 (in a nearby source-free region). For all epochs, we extracted and fitted light curves in the ∼0.2 − 2.0 keV range (see Fig. 1). Timing analysis was performed using UltraNest (Buchner 2019, 2021) and assuming a Gaussian profile for the eruptions, which works generally well for eRO-QPE2 (Arcodia et al. 2021, 2022). Asymmetric models (e.g., Arcodia et al. 2022) will be systematically tested in future work, to account for the occasional residuals found during the decay in some eruptions. Here, we simply account for possible systematic uncertainties on the peak eruptions times by adding in quadrature the median difference between peak times obtained with the Gaussian and asymmetric models (which is ∼57 s). For the August 2020 dataset (top panel of Fig. 1), an additional eruption (with respect to the light curve shown in the original discovery paper; Arcodia et al. 2021) was included in the EPIC-pn data by not performing background flaring screening. This is why the first data point of the August 2020 dataset in Fig. 2 is highlighted with an orange contour, and why it is shown in gray in the top panel of Fig. 1. For the February 2022 dataset (second panel from the top in Fig. 1), the first eruption is only partial in the EPIC-pn data, therefore we used EPIC-MOS1 data to estimate the first recurrence time. Thus, we added in quadrature to the uncertainties on February 2022 arrival times the average difference between EPIC-MOS1 and EPIC-pn arrival times from the full light curve fit (∼58 s). For the June 2022 dataset (third panel form the top in Fig. 1), we extract the first three eruptions from background-contaminated intervals. For December 2023, the same applies to the fourth and sixth eruption (bottom panel of Fig. 1). Thus, for these eruptions and recurrence times we proceeded in the same way as for the first eruption of August 2020. All peak times shown with orange contours in Fig. 2 (i.e., estimated in epochs contaminated by background flaring) have an additional error added in quadrature, estimated by the difference between arrival times or eruptions in light curves with or without background flaring filtering (∼110 s). Similarly to Arcodia et al. (2021), we phase-fold light curve profiles at the eruption peaks for each epoch. In Fig. A.1 we show the median profile (with related 16th and 84th percentile contours) for the August 2020 and December 2023 epochs, in comparison. Despite a slight difference in the median profile toward longer or wider decays, the profiles are consistent at 1σ level.

Fig. A.1. Median light curve profile (with the associated 16th and 84th percentile contours) for the August 2020 (blue) and December 2023 (red, dashed) epochs, folded at the eruption peaks.

Fig. A.2. XMM-Newton spectra for different phases, namely quiescence (circles), rise (up triangles), peak (stars), and decay (down triangles), extracted following the same colors and symbols as in the top panel of Fig. 1. The median models from the best-fit posterior are shown with a dashed (disk component, here diskbb), dotted (additional component, here nthComp) and dash-dotted (eruption component, here zbbody) line.

Event files from EPIC cameras were screened for flaring particle background for the purpose of X-ray spectral analysis. This was 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³. XMM-Newton EPIC pn spectra were fit in the 0.2 − 8.0 keV band using wstat, namely XSPEC implementation of the Cash statistic (Cash 1979). We adopted a Galactic column density of N_H = 1.66 × 10²⁰ cm⁻² from HI4PI (HI4PI Collaboration 2016) and redshifted the source model to rest-frame using the spectroscopic redshift of z = 0.018 (Arcodia et al. 2021), corresponding to a luminosity distance of ∼77 Mpc. Uncertainties are quoted using 16th and 84th percentiles (∼1σ) vlues from fit posteriors, unless otherwise stated. The bolometric luminosity (labeled with “bol”, e.g., L_bol, QPE in Fig. 4 or $F_{\mathit{bol}}^{\mathit{disk}}$ in Table A.1) is obtained integrating the adopted source model between 0.001 and 100 keV. For non-detections, we may quote ∼1σ (∼3σ) upper limits using the 84th (99th) percentiles of the fit posteriors.

We extracted phase-folded spectra for quiescence, rise, peak, and decay selecting the related time intervals. We show this for the August 2020 epoch in the top panel of Fig. 1, as an example, and follow a similar procedure for all epochs. We show the related spectra in Fig. A.2 for the August 2020 epoch, with circles, up triangles, stars and down triangles for quiescence, rise, peak, and decay, respectively. Following the recent QPE literature, we model the quiescence as accretion disk emission (here diskbb) and possible residuals with an additional spectral component (here Comptonization with nthComp). Given the obvious presence of additional absorption (Arcodia et al. 2021), we included in the fit a component at the redshift of the host galaxy (ztbabs): the median (and related 16th-84th percentiles) is (0.47 ± 0.07)×10²² cm⁻², (0.44 ± 0.12)×10²² cm⁻², $0.{43}_{-0.14}^{+0.18}\times {10}^{22}$cm⁻², and (0.51 ± 0.14)×10²² cm⁻², for epochs in chronological order, respectively. Thus, there is no apparent variability in absorption column across the epochs. We modeled the eruptions component with a thermal model (here zbbody). The quiescence spectrum, including the additional host galaxy absorption, is held fixed during the QPE epochs by letting its parameters free to vary only within the 10th-90th percentiles of the posteriors of the quiescence fit alone. We show the median models from the best-fit posterior of the August 2020 epoch in Fig. A.2 with a dashed (disk component, here diskbb), dotted (additional component, here nthComp) and dash-dotted (eruption component, here zbbody) line. The unabsorbed luminosity posteriors for all epochs are shown in Fig. 3, and the line style follows the same coding. We show in Table A.1 the fit results and parameters for all epochs. We note that the additional harder component in quiescence is only statistically required for the August 2020 at high significance, with an improvement in logarithmic Bayesian evidence of a factor Δlog Z ∼ 1100. However, while the flux chain of this additional component is only marginally constrained in February 2022 (Fig. 3), the fit significantly improves (Δlog Z ∼ 64). Thus, even if this component is not required in June 2022 and December 2023, we keep it in the best-fit model so that uncertainties on all parameters are marginalized over this potential spectral component. Finally, we note that we model this residuals as Comptonization given its ubiquity in accreting systems, although its nature in QPE sources is still up to debate. As further discussed in Sect. 2, the time evolution of eRO-QPE2 is remarkable because there is no significant evolution: the accretion disk and eruption component remain remarkably stable in flux and temperature over the course of the ∼3.3 yr baseline (Table A.1 and Fig. 3).

All Tables

All Figures

Fig. 1. XMM-Newton 0.2–2.0 keV light curves of the four epochs, scaled by the arrival time of the first eruption in each observation. Time intervals contaminated by the flaring background are shown in grey. In the top panel, we additionally show the rise, peak, and decay phases selected for August 2020, as an example. In the second panel from the top, the first eruption is only partially resolved in EPIC-pn data and renormalized MOS1 data are shown with empty symbols. As the x-axis range is shared, we note the shortening of recurrence times over the ∼3.3 yr baseline, with vertical dashed lines showing the average arrival time interval of the previous epoch to guide the eye.

In the text

	Fig. 2. Evolution of the recurrence time between eruptions over the ∼3.3 yr baseline. The start time t0 represents the start of the August 2020 observation. Data points with an orange contour contain an additional systematic uncertainty, as is described in the text. The mean value in each epoch, with an associated standard error of the mean, is shown with a dashed line and shaded contours.
In the text

Fig. 3. Evolution of the emission properties. Left: Evolution of the 0.2–8.0 keV luminosity of individual spectral components (as is labeled in the top panel). Histograms represent the fitted flux chains, converted to luminosity and corrected for absorption. Vertical lines and shaded area (following the same coding) highlight the median and 16th–84th interquintile range. There is no significant evolution either of the accretion disk (dashed) or the eruptions (dash-dotted), while the additional spectral component required for the quiescence spectrum in the first epoch becomes less significant and ultimately undetected at later epochs. Right: Same as the left panels, but showing posterior chains of the fitted temperature values. There is no significant evolution of the accretion disk temperature (solid) over the ∼3.3 yr baseline. For the eruption temperature (dash-dotted), despite the median increasing over time and the first and last epochs not being compatible at 1σ, all values are consistent within 2σ.

In the text

Fig. 4. Spectral evolution of QPEs. Darker to lighter colors represent the evolution from the start to end of the eruptions, color-coded as in the legend as a function of the epoch. The top right (medium right) panel shows the QPE bolometric luminosity (temperature) evolution with respect to the eruptions phase. The top left panel shows their coevolution. The bottom right panel shows the evolution of the emitting radius, assuming a thermal spectrum and spherical geometry.

In the text

Fig. 5. Observed recurrence time evolution, overplotting some predicted evolution tracks due to gas drag and slowly evolving disk precession. Different lines are shown for different sets of parameters (as is shown in the legend), such as the spin (a•), the proportionality constant between disk and star nodal precessions (αLT), and the parameter η ≡ ΣR⋆2/m⋆. Variations in the slope of the calculated curves correspond to variations in the relative star-disk inclination, such that the plateaus (low Ṗ) correspond to when the star and disk are near alignment.

In the text

	Fig. A.1. Median light curve profile (with the associated 16th and 84th percentile contours) for the August 2020 (blue) and December 2023 (red, dashed) epochs, folded at the eruption peaks.
In the text

Fig. A.2. XMM-Newton spectra for different phases, namely quiescence (circles), rise (up triangles), peak (stars), and decay (down triangles), extracted following the same colors and symbols as in the top panel of Fig. 1. The median models from the best-fit posterior are shown with a dashed (disk component, here diskbb), dotted (additional component, here nthComp) and dash-dotted (eruption component, here zbbody) line.

In the text