© The Authors 2023

1. Introduction

GSN 069 was first detected in July 2010 during an XMM-Newton slew (Saxton et al. 2011) at a flux level more than a factor 240 above previous upper limits from ROSAT observations performed 16 years earlier. Subsequent observations with the Neil Gehrels Swift observatory (hereafter Swift) show a relatively constant X-ray flux for the first ∼1 yr (Miniutti et al. 2013), and further observations with Swift and XMM-Newton reveal a smooth flux decay during the following ∼7–8 yr. This X-ray long-term evolution is best interpreted as the result of a tidal disruption event (TDE) whose long-lived nature and UV spectral properties suggest the disruption of an evolved star (Shu et al. 2018; Sheng et al. 2021). During XMM-Newton and Chandra observations between December 2018 and February 2019, the X-ray light curve of GSN 069 exhibits high-amplitude, short-lived X-ray flares recurring every ∼9 h (Miniutti et al. 2019). These quasi-periodic eruptions (QPEs) produce an increase in the X-ray count rate by up to two orders of magnitude in the hardest energy bands. They are characterised by thermal-like X-ray spectra with a typical temperature evolving from ≃50 eV up to ∼100–120 eV, superimposed to an otherwise stable blackbody-like quiescent level, most likely due to the emission from the accretion disc resulting from the initial TDE.

Following the discovery of QPEs in GSN 069, this peculiar phenomenon has been identified in a select number of additional sources: RX J1301.9+2747 (Giustini et al. 2020), eRO-QPE1 and eRO-QPE2 (Arcodia et al. 2021), and most likely XMMSL1 J024916.6-041244 (Chakraborty et al. 2021). So far, QPEs have been consistent with being associated with the nuclei of galaxies harboring relatively small mass black holes (10⁵ − few × 10⁶ M_⊙) as derived from a variety of different methods (Miniutti et al. 2019; Arcodia et al. 2021; Wevers et al. 2022). None of the sources show the broad optical emission lines that are associated with active nuclei; however, the presence of a narrow line region (most evident in GSN 069), combined with the lack of intrinsic absorption capable of obscuring the broad line region, suggests some level of past (or present, but weak) nuclear activity (Miniutti et al. 2013; Wevers et al. 2022).

As discussed by Miniutti et al. (2019) the X-ray and optical spectra of the Seyfert 2 galaxy 2XMM J123103.2+110648 (Terashima et al. 2012) closely match those of GSN 069. Its X-ray variability is also peculiar with a ∼3.8 h quasi-periodic oscillation detected in different XMM-Newton exposures (Terashima et al. 2012; Lin et al. 2013), and with overall long- and short-timescale variability properties that suggest a possible connection with both TDEs and QPEs (Lin et al. 2017; Webbe & Young 2023). Recurrent X-ray flares sharing some common characteristics with QPEs have also been recently reported in the ultra-luminous X-ray source XMMU J122939.7+075333 located in the globular cluster RZ 2109 in the Virgo galaxy NGC 4472 (Tiengo et al. 2022). Other sources that may perhaps be associated with QPEs or with the associated phenomenology are ESO 249-39 HLX-1 (Farrell et al. 2009; Godet et al. 2014), ASASSN-14ko (Payne et al. 2021; Liu et al. 2022), and eRASSt J045650.3-203750 (Liu et al. 2023).

Several theoretical models have been considered to explain the QPE phenomenon. Proposed scenarios include (i) modified disc instability models (Sniegowska et al. 2020; Raj & Nixon 2021; Pan et al. 2022; Kaur et al. 2022); (ii) gravitational lensing in a nearly equal-mass supermassive black hole (SMBH) binary (Ingram et al. 2021; iii) mass transfer from one or more orbiting bodies in a variety of different configurations as discussed by King (2020, 2022), Chen et al. (2022), Metzger et al. (2022), Zhao et al. (2022), Wang et al. (2022), Krolik & Linial (2022), and Lu & Quataert (2022); or (iv) collisions between an orbiting secondary body and the accretion disc that is formed following the initial TDE (Suková et al. 2021; Xian et al. 2021). Most proposed models focus on QPEs ignoring the long-term evolution of the quiescent level emission which is particularly well monitored in GSN 069.

In this work, we present results obtained from a series of X-ray observations of GSN 069 (11 by XMM-Newton, one by Chandra, and 34 by Swift) and we discuss the short- and long-timescale properties of both QPEs and continuum (quiescent) emission over the past 12 yr. Our results represent a significant step in deriving observational constraints that can inform further theoretical and numerical work on the QPE formation process and long-term evolution in GSN 069, as well as in the other QPE-emitting galactic nuclei identified so far.

2. Observations and data analysis

GSN 069 was observed in the X-rays on several occasions with the XMM-Newton, Swift, and Chandra X-ray observatories. Most of our work is based on the 11 pointed XMM-Newton observations that are reported in Table 1. Data also comprise one 73 ks Chandra observation performed on 14 February 2019; we only used Chandra data for X-ray variability analysis (light curve). Due to the severe contamination affecting the low-energy response of the ACIS detector, we considered Chandra data above 0.4 keV only, and we ignored them in our spectral analysis as the X-ray emission in GSN 069 is super-soft and cannot be reliably constrained by Chandra. XMM-Newton and Chandra source and background products were extracted from circular regions on the same detector chip using the latest versions of the SAS (XMM-Newton) and CIAO (Chandra) dedicated software. X-ray light curves were background subtracted, as well as corrected for various effects (bad pixels, quantum efficiency, vignetting, dead time) using the SAS epiclccorr and CIAO dmextract tasks. For XMM-Newton, we used data from the EPIC-pn camera only for simplicity. We also analysed 34 observations of GSN 069 made with the Swift X-ray telescope, XRT (Burrows et al. 2005) in photon counting mode. The Swift-XRT observations were analysed following the procedure outlined in Evans et al. (2009), which uses fully calibrated data and corrects for effects such as pile-up and the bad columns on the CCD, to obtain count-rates on an observation-by-observation basis. Occasionally, when the source counts in two neighbouring observations were insufficient to obtain detection of the source, the observations were combined. Finally, although only a minuscule effect, the photon arrival times from all observations were barycentre-corrected in the DE405-ICRS reference system.

3. The transient nature and properties of QPEs in GSN 069

Background-subtracted light curves from all observations with QPEs are shown in Fig. 1 in a common 0.4–1 keV band. The Chandra light curve has been re-scaled to match the XMM-Newton EPIC-pn effective area to better show the QPE intensity evolution regardless of the detector in use.¹ QPEs are consistently detected from December 2018 (XMM3) to January 2020 (XMM6), although they appear much weaker and more irregularly spaced in the latter observation. No QPEs are detected during XMM2, a ∼92 ks XMM-Newton observation in December 2014, that is QPEs first appeared in GSN 069 sometime between December 2014 and December 2018. No clear QPEs are detected in a ∼125 ks observation performed in May 2020 (XMM7) and in none of the subsequent observations. We hence conclude that QPEs in GSN 069 are a transient phenomenon with an observed life-time ≳1.05 yr. Weak QPEs might have been present at other epochs, but undetected because of the reduced contrast against a higher X-ray flux quiescent level, which may increase the actual QPE life-time.

Fig. 1. XMM-Newton (EPIC-pn) and Chandra (ACIS-S) background-subtracted light curves from all observations with QPEs in a common 0.4–1 keV band. The Chandra light curve has been re-scaled to the XMM-Newton pn effective area to better display the QPE amplitude evolution and to ease comparison between data from different detectors. We use time bins of 200 s and 500 s duration for XMM-Newton and Chandra data respectively. All panels have a common y-axis range.

We fitted each 0.4–1 keV light curve with the simplest possible model comprising an observation-dependent constant C, representing the average quiescent level during each exposure, and a series of Gaussian functions with normalisation N and width σ describing QPEs. We define the recurrence time T_rec as the time interval between the peak of two consecutive QPEs. The best-fitting results for all observations are reported in Table A.1. We point out that results are only valid in the 0.4–1 keV band, all quantities being energy-dependent as discussed by Miniutti et al. (2019). In Fig. 2, we show the evolution of the QPE intensity, as measured from the best-fitting normalisation of the Gaussian functions (see Table A.1), averaged over each observation. The QPE intensity decays monotonically by 0.22 cts s⁻¹ per 100 d in the 0.4–1 keV band. We also show the continuum level of the first observation with no QPEs (XMM7). If the decaying trend continued after XMM6, our ability to detect QPEs was seriously compromised already ≃20 d before the XMM7 observation.

Fig. 2. Time evolution of the QPE intensity since first detection. The QPE intensity is measured as the average of the best-fitting Gaussian normalisation in each observation (i.e. we ignore here the difference between strong and weak QPEs). The colour scheme is the same as in Fig. 1. The dotted line represents the best-fitting linear relation (corresponding to a decay of about 0.22 cts s−1 per 100 d). We also show (cross) the continuum level during the first observation with no detected QPEs (XMM7).

As is visually clear from Fig. 1 and quantitatively shown in Table A.1, the QPE intensity oscillates in long enough observations allowing for alternating strong and weak QPE types to be defined. The same is true for the recurrence time between QPEs, with strong (weak) QPEs being systematically followed by longer (shorter) recurrence times to the next QPE (with the exception of the more irregular XMM6 observation). The QPE intensity and recurrence time evolution during all observations comprising QPEs is shown in Fig. 3. The dashed lines are sine functions with period fixed at the observation-dependent average separation between QPEs of the same type (strong or weak), that is at about twice the average recurrence time between consecutive QPEs. As the difference between strong and weak QPE intensities and between long and short T_rec decreases in all observations, we also impose an exponentially decaying amplitude of the sine functions. This signals that the intensity and recurrence time are likely modulated on a longer timescale and in a similar way. In fact, all exponential folding times are consistent with each other (although with very large uncertainties) and are of the order of 1–2 d. The origin of such a longer-term modulation is being analysed, and results will be presented in a forthcoming publication (Miniutti et al., in prep.). The dashed lines in Fig. 3 are not to be considered as actual fits, and they are only included to guide the eye.

Fig. 3. Alternating QPE intensities and recurrence times. The QPE intensity (N) for all observations with QPEs is shown in the upper panels. The recurrence time between consecutive QPEs is shown in the lower panels (each data point is placed at a time corresponding to half the separation between consecutive QPEs). The y-axis is always in the same range, except for the irregular XMM6 observation. Dashed lines in the XMM4, Chandra, and XMM5 panels are sine functions with period equal to twice the averaged, observation-dependent recurrence time, and with an exponentially decaying amplitude that most likely indicates a common long-term modulation. The lines are intended to guide the eye rather than to provide fits to the data.

The behaviour is different during XMM6: while intensities still oscillate, they do so in a much less regular manner to the extent that it is difficult to unambiguously define the QPE type (especially for the first and last QPE). Most importantly, the recurrence time during XMM6 does not oscillate as in all other cases, but increases monotonically. As is reported in Table A.1, the typical QPE duration, defined here arbitrarily as twice the best-fitting Gaussian FWHM, is ≃3780 s (≃1.05 h), and there is no clear difference between strong and weak QPEs duration, nor evident long-term evolution.

In Fig. 4, we plot the ratio between the intensity of consecutive QPEs as a function of the recurrence time between them. The recurrence time is shown as fraction of the average separation between QPEs of the same strong or weak type (or, equivalently, of the sum of consecutive long and short recurrence times), T_sum ≃ 64 ks. The two quantities are well correlated, although significant scatter is present. The upper-left and lower-right quadrants in Fig. 4 are not populated. This means that (at least for the events that we have observed) long recurrence times always follow strong QPEs. Remarkably, the ratio between the intensity of consecutive QPEs is consistent with being the same (N_i /N_i + 1 = 1) when long and short recurrence times have equal duration ($T_{\mathrm{rec}}^i/T_{\mathrm{sum}}=0.5$), as is shown by the intersection of the vertical and horizontal lines with the best-fitting linear relation. Figure 4 suggests an oscillatory behaviour of consecutive QPE intensity and separation around a mean in which intensities and separations of consecutive QPEs are equal. This configuration is never observed exactly, and is almost reached for the 3rd-4th QPEs during XMM4 when N_i /N_i + 1 = 1.18 ± 0.07 and $T_{\mathrm{rec}}^i/T_{\mathrm{sum}}=0.506\pm 0.002$. The QPE X-ray spectral evolution is presented in Sect. 7 together with estimates of the total radiated energy from QPEs in GSN 069.

Fig. 4. Consecutive QPEs intensity ratio and recurrence time correlation. The ratio between the intensity of consecutive QPEs is shown as a function of the recurrence time between them (normalised to the sum of consecutive long and short intervals Tsum). The horizontal line (Ni /Ni + 1 = 1) separates events where strong QPEs precede weak ones (upper half of the figure), from those where weak QPEs precede strong ones (lower half). The $T_{\mathrm{rec}}^i/T_{\mathrm{sum}}=0.5$ vertical line separates short from long recurrence times. The (red) dotted line show the best-fitting linear relation between the two quantities.

4. Quasi-periodic oscillation of the quiescent level (or secondary QPEs)

The statistical quality of the fits to the individual light curves with the simple constant plus Gaussian model is fair, but not excellent (see Table A.1). This is not due to inaccurate modelling of the QPE profiles, but rather to residual quiescent level variability. Visual inspection of the best-fitting residuals reveals that excess emission is systematically present ∼10 ks after most QPEs. As discussed in Appendix B and Fig. B.1, and excluding the irregular XMM6 observation, the excess is quite clear after nine of the ten QPEs observed with XMM-Newton for which sufficient data following a QPE exist (the Chandra data do not have sufficient quality to reveal variability of the quiescent level, as this is barely detected). In fact, the quiescent level appears to vary systematically with a characteristic timescale similar to the average observation-dependent recurrence time between consecutive QPEs, except during the XMM6 observation where excess emission is tentatively seen following QPEs, but with distinct properties (e.g. duration) with respect to all others (see Fig. B.1).

In order to search for systematic trends of the quiescent level variability with better signal-to-noise, we folded the light curves from the XMM3, XMM4, and XMM5 observations at the average observation-dependent recurrence time. We ignored the XMM6 observation whose increasing recurrence times prevented us from defining a clear folding timescale (see Fig. 3). As Chandra data were ignored, there was no need to restrict the energy band above 0.4 keV, and we used the full 0.2–1 keV band in this analysis.

The 0.2–1 keV folded light curves are shown in the upper panels of Fig. 5 together with their best-fitting baseline (constant plus Gaussian functions) model. The middle panels show the resulting residual light curves where a sinusoidal trend is evident, strongly suggesting the presence of a quiescent level quasi-periodic oscillation (QPO) with period equal to the average observation-dependent recurrence time. Re-fitting the original folded light curves with the addition of a sine function (with period fixed at the folding timescale) improves very significantly the fits in all cases, and best-fitting residuals are shown in the lower panels. The best-fitting statistics for the model comprising the QPO provides and improvement of Δχ² = −134 (XMM3), Δχ² = −302 (XMM4), and Δχ² = −133 (XMM5) for 2 degrees of freedom. The remaining residuals (see lower panel of Fig. 5) can be attributed to the fact that QPE profiles are not exactly Gaussian, especially once folded at the average recurrence time since the periodicity is not perfect and the folding process introduces small distortions. We conclude that a QPO is highly significant in the folded light curves. Although deriving periodicity from arbitrarily folded light curves is not always robust, the extremely similar shape of the residual light curves in all observations makes it highly unlikely that the derived sinusoidal trend is spurious. The vertical lines in Fig. 5 mark the phase of the QPE (dotted) and of the peak of the QPO (dashed).

Fig. 5. Quasi-periodic oscillation of the quiescent level. In the upper panels, we show the 0.2–1 keV light curves from the XMM3, XMM4, and XMM5 folded at the average observation-dependent recurrence time, and the corresponding best-fitting baseline model (constant plus Gaussian functions, dotted). We show the 0–2 phase interval for visual clarity. In the middle panels, we show the resulting residual light curves together with a sine function with period fixed to 1 (in phase space) to guide the eye. The vertical lines represent the phase of the primary QPE (dotted) and of the peak of the sine function (dashed). In the lower panel, we show the residual for a fit of the original folded light curve that comprises a sine function with period fixed to 1 (i.e. to the folding timescale) in addition to the baseline model.

In Fig. 6, we show the time-evolution of the QPE intensity (here the average between strong and weak ones, as we folded on the recurrence time between consecutive QPEs) as well as of the ratio between the QPE and QPO best-fitting normalisation. The intensity of both QPEs and QPOs decays with time, but their ratio is consistent with being constant (and equal to ∼37.2). The peak of the QPO emission lags the preceding QPE by 8.3 ± 0.3 ks, 10.6 ± 0.3 ks, and 10.5 ± 0.3 ks in XMM3, XMM4, and XMM5 respectively. Although the number of data points is limited, it is tempting to associate shorter time delays to observations with shorter recurrence times. In fact, the delay and observation-dependent recurrence (folding) time appear to obey a 1:1 correlation, as shown in Fig. 7.

Fig. 6. Evolution of the QPE intensity (upper panel) and of the ratio between the QPE and QPO intensities (lower panel). The latter show that the QPE and QPO decay together preserving a ratio of 37.2 ± 1.3 (horizontal lines).

Fig. 7. QPE-QPO time delay. We show the time delay of the QPO peak with respect to the preceding QPE as a function of the average observation-dependent recurrence time. The dotted line is a linear function with slope fixed to 1, the best-fitting one being 0.95 ± 0.15.

The detection of a quiescent level QPO with characteristic timescale of ∼9 h, equal to the QPE recurrence time, suggests that the QPO production mechanism is linked to every single QPE regardless of its intrinsic type (strong or weak), and likely indicates that every QPE induces either an oscillation of the quiescent (disc) emission, or the production of a weaker and broader (longer-lasting) secondary QPE. A common production mechanism (or at least a link between the QPE and the QPO) is also suggested by their common time evolution which preserves their intensity ratio as well as by the 1:1 correlation between the QPO time delay (with respect to the preceding QPE) and the observation-dependent recurrence time.

A QPO of the quiescent level was also reported for another QPE-source, RX J1301.9+2747 (Giustini et al. 2020) by Song et al. (2020). In RX J1301.9+2747, the QPO is detected in two observations ∼18, 5 yr apart (both of which presenting X-ray QPEs), with a stable timescale of ∼1 500 s. The QPO timescale in GSN 069 (∼30–32 ks) and RX J1301.9+2747 (∼1.5 ks) are remarkably different despite similar QPE recurrence times in the two sources (within a factor of 2) which may indicate a different production mechanism. In eRO-QPE1 no QPO is detected, although QPEs there are significantly more complex than in GSN 069 and it may be difficult to separate QPOs from multi-component QPEs (Arcodia et al. 2022). In eRO-QPE2, the quiescent level does not have high enough signal-to-noise to enable us to investigate the presence of any QPO in detail (Arcodia et al., in prep.).

5. Long-term evolution of the quiescent level emission

Following its first X-ray detection during an XMM-Newton slew (14 July 2010), GSN 069 has been consistently detected in the soft X-rays. Its X-ray flux evolution during the first few yr shows a smooth decay by a factor of a few, consistent with a long-lived TDE, as first noted by Shu et al. (2018). UV spectra from Hubble Space Telescope observations obtained in 2014 and 2018 show anomalous N-rich abundance (Sheng et al. 2021), also suggesting an association with TDEs (Kochanek 2016).

5.1. XMM-Newton spectral analysis

We extracted EPIC-pn spectra from the available 11 XMM-Newton observations excluding time intervals comprising QPEs, as our goal here is to present the properties and evolution of the quiescent emission, most likely due to the accretion flow formed by the initial TDE. We assumed a model comprising thermal disc emission using the DISKBB model in XSPEC (Mitsuda et al. 1984) and a power law component, both at the redshift of GSN 069 (z = 0.0181). We also included neutral absorption with N_H = 2.3 × 10²⁰ cm⁻² fixed at the Galactic value (HI4PI Collaboration 2016).

Most spectra exhibit an absorption structure around 0.7 keV, first noted by Miniutti et al. (2013) in the XMM1 spectrum. Indeed, the feature is clearer in the XMM1 highest-flux observation, but appears possibly present in all others as well. In order to first establish whether the quiescent level in observations with QPEs is different from the others in terms of absorption, we combined X-ray spectra from different observations to increase the signal-to-noise and we produced: (i) a merged spectrum of the quiescent level from all observations comprising QPEs (XMM3 to XMM6); (ii) a merged spectrum from all other observations except XMM1 (where the feature is already clear). The two merged spectra include observations with similar flux levels, so that none of the two spectra is dominated by one particular observation, but represents a sensible average in both cases. Both merged spectra show clear signs of an absorption feature around 0.7 keV, consistent with being the same.

As the two merged spectra (from observations with and without QPEs) are consistent with each other in terms of absorption, we then produced a single merged spectrum by combining the XMM2 to XMM11 spectra, and we compared it to the XMM1 spectrum (the one with the clearest absorption feature). In order to account for the extra-absorption (most likely from ionised gas), we included a warm absorber model using a custom-built XSTAR grid (Kallman & Bautista 2001) that assumes a black body with kT = 50 eV as irradiating spectral energy distribution (SED), as appropriate for GSN 069 (Miniutti et al. 2019). We assumed a turbulent velocity of 100 km s⁻¹. Residuals for the merged XMM2 to XMM11 spectrum with and without the warm absorber model are presented in Fig. 8. The model with no warm absorber (upper panel) does require a higher neutral column density (≃5.3 × 10²⁰ cm⁻²) than the Galactic one (≃2.3 × 10²⁰ cm⁻²), while the model with warm absorption does not. The warm absorber best-fitting parameters for the merged spectrum are N_H = (6.5 ± 0.1)×10²¹ cm⁻² and log ξ = 0.49 ± 0.01, consistent with those derived from XMM1, N${}_{\mathrm{H}}^{(\mathrm{XMM}1)}=(7.0\pm 2.0)\times {10}^{21}$ cm⁻² and log ξ^(XMM1) = 0.52 ± 0.03, although the latter have larger uncertainties due to the much lower signal-to-noise. We conclude that a warm absorber is always present in GSN 069 and that its properties (column density and ionisation) are consistent with being the same at all epochs.

Fig. 8. Residuals (in terms of σ) when the merged XMM2 to XMM11 spectrum is fitted without (upper panel) and with (lower panel) a warm absorber component. We only show the most interesting energy range close to the most prominent absorption feature at ≃0.7 keV for better visual clarity. Data have been slightly re-binned for visualisation purposes only.

We then fitted jointly the 11 individual spectra (XMM1 to XMM11) including the warm absorber in the spectral model. After a few initial tests, and as suggested by the previous analysis, we found that the quality of the fits did not worsen significantly by forcing the warm absorber parameters to be the same at all epochs. The final fit is acceptable with χ² = 922 for 814 degrees of freedom. The (common) warm absorber column density is N_H = (6.6 ± 0.1)×10²¹ cm⁻² with log ξ = 0.49 ± 0.01. Allowing for an outflow/inflow velocity of the absorber does not improve the fit, as expected because of the moderate resolution of the EPIC-pn detector. The hard power law component is always detected, but it is so weak that its contribution below 1 keV is negligible (Miniutti et al. 2019). Excluding the power law component from the spectral model slightly worsen the fitting statistics but does not change significantly any of the relevant parameters.

The resulting observed 0.3–1 keV X-ray flux, intrinsic X-ray luminosity, and (rest-frame) disc temperature evolution are shown in Fig. 9. As is clear from the upper two panels, the X-ray emission of GSN 069 exhibits an initial decay during the first ∼9 yr, followed by a highly significant X-ray re-brightening that appears to start just after XMM5, that is the last observation where regular QPEs are detected. Here we focus first on the overall spectral evolution, deferring the discussion of this interesting behaviour to Sect. 6. Besides the XMM-Newton slew data point that is affected by very large uncertainties, the highest 0.3–1 keV X-ray luminosity, L_X ≃ 8.7 × 10⁴² erg s⁻¹, is observed during XMM1 (∼140 d after the initial slew detection), while the lowest is reached during XMM5 (∼3 240 d after the slew detection) and is L_X ≃ 1.5 × 10⁴² erg s⁻¹, namely a factor of ≃5.8 lower. We point out that the bolometric luminosity variation is likely to be significantly smaller because part of the variability amplitude in the L_X (and F_X) light curve is due to the shift of the overall SED to lower energies as the temperature decays (see lower panel in Fig. 9).

Fig. 9. X-ray flux, luminosity, and temperature evolution. From upper to lower panel, we show the observed 0.3–1 keV X-ray flux, the 0.3–1 keV intrinsic X-ray luminosity, and the (rest-frame) disc temperature evolution, as derived from our best-fitting DISKBB models to the 11 XMM-Newton pointed observations. The first data point in the upper and middle panels refers to the slew observation whose observed count rate has been converted into X-ray flux and luminosity by using the best-fitting model from the nearest XMM-Newton pointed observation (XMM1). In all panels, the shaded areas highlights the period during which regular (red, XMM3 to XMM5) and irregular (blue, XMM6) QPEs are detected.

The relation between intrinsic disc X-ray luminosity and rest-frame temperature is shown in Fig. 10. We show the relation for the observed 0.3 − 1 keV L_X (grey squares) as well as that for the 0.2 − 2 keV band (black circles), obtained by extrapolating the best-fitting model. The dotted lines are best-fitting relations of the form L_X ∝ T^q. The 0.3 − 1 keV L_X − T relation is steeper (q ≃ 5.6) than the 0.2 − 2 keV one (q ≃ 4.1). The different slope is expected because the narrower 0.3 − 1 keV band is more severely affected by the shift of the SED to lower energies as the temperature decays. What is surprising, however, is that the 0.2 − 2 keV luminosity is consistent with the L ∝ T⁴ relation that, for constant-area blackbody emission, is naturally expected to be satisfied by the bolometric, rather than narrow-band, luminosity. As, for a standard disc, the X-ray bolometric correction (K_{0.2−2 keV} = L_bol/L_{0.2−2 keV}) is a function of temperature² and varies in the range ≃[4.9 − 11.6] for the observed range of kT ≃ [43 − 60] eV, the L ∝ T^q relation cannot be satisfied by L_bol with q = 4, and the slope is necessarily shallower. Forcing the L_bol ∝ T⁴ relation instead,³ produces a lower quality fit with χ² = 1002 for 824 degrees of freedom (to be compared with χ² = 922 for 814 degrees of freedom for the case in which the DISKBB normalisation is free to vary).

Fig. 10. Intrinsic 0.2 − 2 keV (black circles) and 0.3 − 1 keV (grey squares) X-ray luminosity as a function of rest-frame temperature for all 11 XMM-Newton observations. The dotted lines are best-fitting relations assuming LX ∝ Tq.

The L_{0.2−2 keV} ∝ T⁴ relation suggests that L_{0.2−2 keV} is a good proxy of the bolometric luminosity, that is that the bolometric correction is small and approximately constant in the observed range of temperatures. This cannot be achieved with a standard disc, but it can approximately be obtained with a compact accretion flow with small outer radius (R_out) significantly suppressing any disc optical/UV emission, that is a disc emitting most of its energy in the soft X-rays. We point out that, while a small R_out is implausible in active galactic nuclei (AGNs), the accretion flow in TDEs is at least initially expected to be compact, as the stellar debris circularise at radii comparable (within a factor of a few) to the star’s pericenter distance that can be as small as few-tens gravitational radii (R_g = GM_BH/c²).

5.2. UV emission

In principle, the disc contribution in the UV can be derived from data obtained with the Optical Monitor (OM) on board XMM-Newton. Most UV fluxes from the OM, obtained simultaneously with the X-ray data, do not show any variability despite the large X-ray variation, which suggests that the UV emission in GSN 069 is quite severely contaminated by stellar light or not directly associated with disc emission, but rather to another process (for example, UV reprocessing). The UV flux light curves for the three available OM filters are shown in Fig. 11. The UVM2 filter (at 2 310 Å) is the only one where some variability is observed with higher fluxes during the first two observations (XMM1 and XMM2) than during all others (XMM3 to XMM11). Fitting a constant to the light curves results in ${\chi }_{\nu }^2 \simeq 0.2$ (0.4) for UVW1 (UVW2), and ${\chi }_{\nu }^2\simeq 2.2$ for UVM2 (due to the first two data points only, XMM1 and XMM2). We must however point out that the XMM7 to XMM11 observations (last five data points in Fig. 11) have very similar X-ray flux to XMM2, but the UVM2 flux does not recover, casting some doubts on the reliability of the UVM2 fluxes (and therefore putative variability) in XMM1 and XMM2. Moreover, UV fluxes from the bluer UVW2 filter (2120 Å) and the redder UVW1 one (2910 Å) do not show any significant variability, which is surprising if the UV are dominated by disc emission.

Fig. 11. UV light curves from OM data obtained during the XMM-Newton observations. The UVW2 (λeff ≃ 2 120 Å), and UVW1 (≃2 910 Å) fluxes are re-scaled by factors of 1.1 and 0.8 respectively for visual clarity, while the UVM2 (≃2 310 Å) remains as is.

Another estimate of the UV flux can be obtained by making use of two Hubble Space Telescope spectroscopic observations that have been performed quasi-simultaneously with XMM2 and XMM3. The 2014 (XMM2) flux is slightly higher than that measured in 2018 (XMM3) by ≲7%, while the X-ray flux decayed by ∼60%. Milky Way’s absorption appears to be sufficient to account for the slight reddening in the Hubble Space Telescope spectra, and the de-reddened UV flux in the relatively features-free range around 1400 Å is ≃1.4 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ which, after subtracting appropriate stellar continuum, reduces to ≃1 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ (Sheng et al. 2021). Our best-fitting DISKBB model over-predicts the intrinsic 1400 Å flux by a factor of 4 − 5 during XMM2 and XMM3, which may be taken as an indication of a compact accretion flow with suppressed UV emission (as suggested by the L_{0.2−2 keV} ∝ T⁴ relation), or of the overall inadequacy of the model to reproduce the broadband SED.

A conservative analysis should thus look into solutions in which the UV disc emission is constrained between two extremal cases, namely: (i) the UV flux is entirely due to disc emission and the disc model only needs not to over-predict it; (ii) the UV flux is entirely due to a different phenomenon (e.g. reprocessing) and the disc model is allowed to have negligible contribution in the UV. These extremal cases will set the systematic uncertainty on the disc bolometric luminosity derived below. In our analysis, we took as UV flux that measured by Hubble Space Telescope after correcting for stellar contamination and Milky Way reddening, that is ≃1 × 10⁻¹⁵ erg s⁻¹ cm⁻² Å⁻¹ at 1 400 Å, and compared it with the low-energy extrapolation of unabsorbed X-ray accretion disc models at the same epoch (XMM2 and XMM3).

5.3. Bolometric luminosity evolution

We considered an accretion disc model that is constructed with the goal of reproducing the broadband SED of AGNs, namely the OPTXAGNF model in XSPEC (Done et al. 2012). As we are only interested in accretion disc emission, we switched off all Comptonisation zones that are originally introduced to account for the high-energy power law and soft excess X-ray components in AGNs. The treatment of the disc emission is more physical than in the DISKBB model used so far. In particular, each annulus on the disc thermalises to the local blackbody temperature, but a temperature-dependent colour-correction factor is introduced to account for electron scattering within the disc. Such correction factor basically shifts the emission towards higher energies and produces higher X-ray-to-UV flux ratios with respect to models where the correction is ignored, which may potentially account for the depressed UV flux discussed above. The model naturally enforces the L_bol ∝ T⁴ relation. When all Comptonisation components are switched off, the model only depends on black hole mass M_BH and spin a, Eddington ratio, and outer disc radius R_out. The black hole spin sets the inner disc radius, assumed to be coincident with the innermost stable circular orbit at radius R_isco. In our analysis, we assumed a = 0, because the corresponding R_isco = 6 R_g is intermediate between those associated with accretion around a maximally spinning black hole for a prograde or retrograde accretion flow (R_isco ≃ 1.24 R_g and R_isco ≃ 9 R_g respectively). We initially fixed R_out = 10³ R_g.

We fitted all XMM-Newton spectra jointly, keeping the warm absorber and power law components as before (with warm absorber parameters forced to be the same in all spectra, and power law photon index fixed at Γ = 1.9). The fit is relatively good with χ² = 999 for 823 degrees of freedom and has similar quality to that obtained with DISKBB forcing the normalisation to always be the same (i.e. enforcing the L_bol ∝ T⁴ relation), which resulted in χ² = 1002 for 824 degrees of freedom. The unabsorbed best-fitting model, extrapolated to the UV, still over-predicts the 1 400 Å flux by a large factor. We then decreased the outer disc radius R_out until a good match with the measured Hubble Space Telescope flux was found, and we derive R_out ≃ 140 R_g, that is the maximum allowed disc outer radius that is needed not to over-predict the UV flux. As for the minimum one, corresponding to negligible UV contribution, it cannot be constrained by using the UV flux level. Hence, we re-fitted the X-ray data by allowing R_out to vary, forcing it to be the same for all observations (i.e. at all epochs). The statistical quality of the fit improves significantly, by Δχ² = −30 for 1 degree of freedom, and we obtain R_out = 15 ± 5 R_g, that is the X-ray data alone tend to prefer a compact accretion disc with small outer radius.

The L_bol evolution from the 11 XMM-Newton pointed observations is shown in Fig. 12 where we also indicate the XMM-Newton observation number corresponding to each data point for reference. We use two different y-axis scales: the left axis is L_bol as measured by the compact disc model (R_out ≃ 15 R_g), while the right axis refers to the more extended one (R_out ≃ 140 R_g). The bolometric luminosity ranges from L_bol = 0.5 − 1.3 × 10⁴⁴ erg s⁻¹ (XMM1) to 0.2 − 0.5 × 10⁴⁴ erg s⁻¹ (XMM5) where the lower luminosity is associated with the compact disc.

Fig. 12. Bolometric luminosity evolution as inferred from OPTXAGNF best-fitting models to the 11 XMM-Newton observations. The left y-axis corresponds to the compact disc (Rout ≃ 15 Rg), and the right y-axis to the more extended one (Rout ≃ 140 Rg). The conversion between the two is a factor of ≃2.5. For each data point, we also show the corresponding XMM-Newton observation number (XMM1 to XMM11) as reference.

As for the remaining best-fitting parameters, the warm absorber and power law are unchanged with respect to previous fits using the DISKBB model, while we obtain a black hole mass of M_BH = 2 − 4 × 10⁶ M_⊙ depending on the adopted model, the lower mass being associated with the larger R_out. The derived M_BH is significantly higher than that estimated by Miniutti et al. (2019) using DISKBB fits (≃4 × 10⁵ M_⊙), and this is due to the different disc model used here that includes the colour-correction factor, shifting the SED to higher energies and therefore requiring higher black hole mass. Indeed, if we remove the colour-correction, adopting the similar but not colour-corrected AGNSED model (Kubota & Done 2018), we recover a best-fitting black hole mass of (4 ± 1)×10⁵ M_⊙. Hence, while our analysis allows us to estimate the disc bolometric luminosity of GSN 069, at least within the systematic uncertainty associated with the different possible R_out, the black hole mass (and hence Eddington ratio evolution) remains highly model-dependent and thus largely unconstrained.

As a final note on black hole mass, the stellar velocity dispersion in GSN 069 is σ_⋆ = (63 ± 4) km s⁻¹ which can be used to infer a black hole mass of M_BH = 3 − 30 × 10⁵ M_⊙ (Wevers et al. 2022), consistent with the range of black hole masses that we infer from the X-ray spectral analysis using models with and without colour-correction. When needed, hereafter we consider the full range of plausible black hole masses as obtained from our spectral analysis using models with (OPTXAGNF model) and without (AGNSED model) colour-correction, namely M_BH = 0.3 − 4 × 10⁶ M_⊙. The estimated range assumes a non-rotating Schwarzschild black hole. Significant black hole spin associated with either prograde or retrograde accretion induces an even broader range of allowed black hole mass. For reference, assuming (quite arbitrarily) that the historical peak luminosity of GSN 069 corresponds to its Eddington luminosity, the black hole mass range shrinks to 0.4 − 1 × 10⁶ M_⊙ depending on the assumed outer disc radius.

As a caveat, we note that the X-ray luminosity in Figs. 9 and 10 are different from those reported in Miniutti et al. (2019) for the common observations. This is due to the different absorption model adopted here. The warm absorber induces a significantly higher X-ray luminosity than the simpler neutral absorber used by Miniutti et al. (2019) as it is associated with a higher column density and only moderate ionisation. While more physical, our specific absorption model uses an XSTAR grid obtained by assuming a standard disc SED (with significant optical/UV emission) as well as Solar abundances, which may not be accurate if the accretion flow if confined to small radii or if the absorbing matter is associated with the unbound debris from the TDE which may well have different abundances, depending on the properties of the disrupted star. While the relationship between luminosity and temperature shown in Fig. 10 is preserved for whatever specific absorption model (provided that its parameters do not evolve with time, as the data suggest), the absolute value of the reported X-ray luminosities, as well as the estimate on radiated energy and accreted mass given below, should then be taken with caution.

6. GSN 069 as a repeating TDE candidate

In order to study the long-term evolution of GSN 069, we considered, besides the 11 XMM-Newton pointed observations discussed above, an additional set of data provided by long-term X-ray monitoring observations of GSN 069 with the Swift observatory. As Swift observations (as well as the initial XMM-Newton slew detection) are shallow, only few-tens of X-ray counts per pointing are typically collected, so that spectral information is limited or absent. We then used best-fitting models from appropriate (i.e. close in time) deep XMM-Newton exposures to convert count rates into luminosities. Figure 13 shows the L_bol evolution over the whole ∼12 yr period probed by our X-ray campaign. The upper panel shows the original light curve from all XMM-Newton and Swift observations. The left (right) y-axis refers to L_bol as measured with the compact (extended) disc model. Hereafter, when discussing the L_bol evolution, we refer to measurements obtained with the compact disc. However, in our estimates of physical quantities (e.g. radiated energy and accreted mass), we always consider the full range of possible L_bol. As shown in the upper panel of Fig. 13, significant short-term variability is present on timescales of the order of few-tens of days. In order to study the long-term evolution reducing this source of noise as well as the uncertainties of the Swift data points, we slightly re-binned the original L_bol light curve by combining observations that were performed close in time (we combined Swift and XMM-Newton observations separately, while the XMM-Newton slew detection remained as is). The slightly re-binned L_bol light curve is shown in the lower panel of Fig. 13.

Fig. 13. Long-term evolution of the (quiescent level) bolometric luminosity. The y-axis scales have been defined in the caption of Fig. 12. In the upper panel we show the original light curve from all XMM-Newton and Swift observations. The lower panel shows a re-binned version of the same light curve (see text for details). In both panels, a vertical dashed line separates two different phases of the evolution: Phase1 comprises the initial X-ray peak and subsequent decay during the first ≃9 yr, while Phase2 is associated with the late-time X-ray re-brightening.

The L_bol light curve in Fig. 13 (lower panel) is characterised by an initial rise peaking at L_bol ≃ 6 − 7 × 10⁴³ erg s⁻¹, and subsequent decay lasting ≃9 yr. A significant X-ray re-brightening is seen after that, peaking at ≃3 360 d with L_bol ≃ 4 − 5 × 10⁴³ erg s⁻¹. The re-brightening peak is followed by a fast luminosity drop by a factor of ≃1.5 in less than ∼100 d. The final evolution is characterised by a further rise and subsequent decay. We define Phase1 as the initial 9 yr-long decay, and Phase2 as the subsequent X-ray re-brightening.

6.1. Phase1: The initial decay

The Swift and XMM-NewtonL_bol light curve during Phase1 is shown in Fig. 14. The evolution is highly suggestive of the typical rise and decay light curve of TDEs. We adopted a simple parametrisation of the L_bol = L_bol(t) light curve in which the rise is described by a Gaussian function, while the decay has a power law form (van Velzen et al. 2021):

$$\begin{array}{cc}\hfill L_{\mathrm{bol}}(t)& =L_{\mathrm{peak}}\times \{\begin{array}{cc}{e}^{-{(t-t_{\mathrm{peak}})}^2/2{\sigma }^2}\hfill & t\le t_{\mathrm{peak}}\hfill \\ {[(t-t_{\mathrm{peak}}+t_0)/t_0]}^n\hfill & t>t_{\mathrm{peak}}\hfill \end{array}\hfill \end{array}$$(1)

Fig. 14. Lbol light curve during Phase1. The y-axis scales have been defined in the caption of Fig. 12. We show a simple model comprising a Gaussian rise and a power law decay with t−9/4.

where t_peak is the time when the peak luminosity L_peak is reached, σ is the Gaussian rise width, t₀ is the normalisation of the power law decay, and n is its power law index.

The decay power law index cannot be constrained well by the data. In particular, both canonical values of n = −5/3 and n = −9/4, describing the decay for full or partial TDEs respectively based on fallback rate calculations (Rees 1988; Phinney et al. 1989; Coughlin & Nixon 2019; Miles et al. 2020; Nixon et al. 2021), provide a good representation of the decay and cannot be distinguished on a statistical basis. We point out, however, that even an exponential decay with e-folding timescale τ ≃ 7.6 yr can describe the data well with similar statistical quality. For simplicity, and for reasons that will become clear in the subsequent sections, we assume n = −9/4, as expected from the fallback rate in the case of a partial TDE (pTDE). As is clear from Fig. 14, our description is not unique. In particular, the rise and decay phases may be connected by a plateau that we ignore here for simplicity. As for the rise phase, the Gaussian width is σ ≃ 500 − 700 d, peaking at t_peak ≃ 500 − 800 d with L_bol ≃ 6 − 7 × 10⁴³ erg s⁻¹ cm⁻². The best-fitting model is shown as a (blue) solid line in Fig. 14. We note that the first two data points (from the XMM-Newton slew and first two Swift observations combined) lie above the model. We shall discuss this small discrepancy further in Sect. 8.4.

By integrating L_bol(t), the total radiated energy during Phase1 is E^(obs) = (2.3 ± 1.0)×10⁵² erg, where the uncertainty (here and hereafter) reflects the systematic uncertainty on L_bol associated to the compact or extended disc. By assuming a radiative efficiency η = 0.1, this corresponds to an accreted mass of $M_{\mathrm{accr}}^{(\mathrm{obs})}=0.13\pm 0.06M_{\odot }$. Extending to t = ∞, one has E^(∞) = (4.1 ± 1.8)×10⁵² erg and $M_{\mathrm{accr}}^{(\infty )}=0.23\pm 0.10M_{\odot }$. As bound and unbound debris in TDEs likely represent ∼50% of the disrupted mass, the TDE in Phase1 likely disrupted of the order of 0.46 ± 0.20 M_⊙.

6.2. Phase2: The X-ray re-brightening

The L_bol light curve during Phase2 is shown in Fig. 15. The first three data points in Fig. 15 are the last three of the Phase1 decay in Fig. 14, and we also re-plot the best-fitting decay law (n = −9/4) until Phase2 begins. Besides an apparent high-amplitude X-ray flare at the beginning of Phase2 (red dotted line in Fig. 15), the evolution is very similar to that during Phase1. As a first step, we treated Phase2 as independent from Phase1. This is likely inaccurate, as part of the observed luminosity is still associated with residual accretion of debris from Phase1. We adopted a simple parametrisation of Phase2 by describing the initial X-ray flare with a Gaussian function, and the subsequent evolution with the same model that successfully describes Phase1, namely a Gaussian rise and subsequent power law decay with index n. We first discuss the TDE-like evolution (Gaussian rise and power law decay), ignoring for a moment the X-ray flare precursor.

Fig. 15. Lbol light curve during Phase2. The y-axis scales have been defined in the caption of Fig. 12. The first three data points are from the final decay of Phase1 and we also show the best-fitting decay model until Phase2 begins. A simple model comprising a Gaussian rise and a power law decay with t−9/4 for Phase2 is shown as a solid blue line, together with a Gaussian X-ray flare precursor (red dotted line).

As was the case for Phase1, power law indices of n = −5/3 and n = −9/4 are both acceptable for the decay, and produce the same statistical quality with reduced χ² ≃ 1. We adopt n = −9/4 as during Phase1. An exponential is also an excellent fir to the decay of Phase2 and results in an e-folding timescale of ≃4.8 yr, slightly shorter than that characterising Phase1 (≃7.6 yr). The Gaussian rise has width σ ≃ 400 d, slightly shorter than the rise of Phase1. The peak luminosity is ≃4 − 5 × 10⁴³ erg s⁻¹ (the luminosity for the extended disc being ≃2.5 times higher) and is reached at t_peak ≃ 3750 − 3900 d. By integrating over the observed time-interval, and accounting for the systematic uncertainty in L_bol, the total radiated energy is E^(obs) = (4.5 ± 0.9)×10⁵¹ erg, corresponding to $M_{\mathrm{accr}}^{(\mathrm{obs})}=0.025\pm 0.005M_{\odot }$. Extending to t = ∞, E^(∞) = (2.3 ± 1.0)×10⁵² erg with $M_{\mathrm{accr}}^{(\infty )}=0.13\pm 0.06M_{\odot }$.

However, part of the emitted radiation in Phase2 is likely associated with the accretion of residual mass disrupted in Phase1 and not yet accreted. In order to account for it (at least qualitatively), we re-fitted Phase2 by including, as baseline level, the extrapolation of the best-fitting Phase1 model. This significantly reduces the luminosity associated solely to the Phase2 event, as well as the associated timescales. The decay e-folding time shortens from ∼4.8 yr to ≃3.4 yr, and the Gaussian rise width from ≃400 d to ≃270 d. The total radiated energy lowers to E^(obs) = (2.1 ± 0.4)×10⁵¹ erg, and E^(∞) = (6.5 ± 1.3)×10⁵¹ erg, corresponding to $M_{\mathrm{accr}}^{(\mathrm{obs})}=0.0125\pm 0.002M_{\odot }$ and $M_{\mathrm{accr}}^{(\infty )}=0.036\pm 0.007M_{\odot }$. The disrupted mass during Phase2 is therefore highly uncertain and strongly depends on whether we consider or not that part of the X-ray luminosity comes from the accretion of debris from the previous Phase1. We conclude that sensible lower and upper limits are set by the two extremal situations we have considered, so that we conservatively estimate a disrupted mass in the range of 0.06 − 0.38 M_⊙ during Phase2 (assuming M_disr ≃ 2 M_accr). We point out, however, that the disrupted mass is likely closer to the lower rather than upper limit because the luminosity from the accretion of the remaining debris from Phase1 cannot be completely ignored.

As is shown in Fig. 15, an X-ray precursor flare is seen before the peak of Phase2. When described with a Gaussian function, the precursor flare peaks at t ≃ 3 380 d with peak luminosity in the range of of 3 − 5 × 10⁴³ erg s⁻¹ cm⁻² (depending on whether we include the Phase1 model extrapolation in the fit or not). The Gaussian width is 50 − 90 d, and the radiated energy is E = 0.3 − 1.8 × 10⁵¹ erg once all uncertainties are included. Under the (strong) assumption that the precursor flare is associated with a radiative efficiency η ≃ 0.1, one has M_accr = 0.002 − 0.01 M_⊙. In principle, as the precursor flare is only defined by one single data point in the re-binned light curve of Fig. 15, one should consider its detection as merely tentative. However, in Fig. 16 we re-plot the same Phase2 re-binned light curve shown in Fig. 15 but with the original Swift data points around the precursor flare superimposed. The precursor flare appears to be confirmed by at least another couple of data points, and we therefore consider it as a real feature of the Phase2 light curve.

Fig. 16. Re-binned light curve during Phase 2 (black circles), together with its best-fitting model (see text for details). The green squares are the original Swift data points (before re-binning) around the precursor flare.

7. QPE spectral evolution

Having studied in detail the spectral properties of the quiescent emission of GSN 069, we are now in a position to describe the QPE spectral evolution. We considered as case study QPEs during the XMM4 observation. The same conclusions discussed here were reached when the XMM3 and XMM5 observations were studied, the only difference being that the QPE peak is slightly more (less) luminous in XMM3 (XMM5) with respect to XMM4 (see e.g. Fig. 2). As for the weak, irregular QPEs during the XMM6 observation, the low signal-to-noise prevented us from reaching a sufficiently high spectral quality to perform a detailed spectral analysis, and we can only report that the typical peak luminosity if a factor of ∼2 lower than during XMM4.

To increase the signal-to-noise during the XMM4 observation, we accumulated spectra during the rise, peak, and decay of QPEs without distinguishing between strong and weak ones⁴. We defined five different spectra during the QPE evolution, namely: the quiescent spectrum (Q), two spectra defined during the QPE rise (hereafter R1 and R2), the peak spectrum (P), and two spectra during the QPE decay (D1 and D2). The phase-resolved spectra were accumulated during the intervals shown in Fig. 17, as in Miniutti et al. (2019). We then subtracted the Q spectrum from all others to study the evolution of the intrinsic QPE spectrum. The resulting difference spectra were fitted jointly, and are all well described by a blackbody model (at the redshift of GSN 069) with only Galactic absorption. The blackbody model we used is not unique, and the spectral evolution could also be discussed in terms of other (thermal-like) spectral models. However, we used the simplest possible spectral model as it allows us to highlight the main properties of the QPE spectral evolution.

Fig. 17. XMM4 folded light curve (at the average recurrence time between QPEs). The light curve is normalised to the peak amplitude and shifted so that the peak is at phase 1. We show with different colours the phase intervals used to accumulate the spectra during the rise (R1 and R2), peak (P), and decay (D1 and D2). The quiescent spectrum (Q) is accumulated throughout the exposure (excluding time-intervals when QPE are present) and only a small fraction of the corresponding phase-interval is shown here. A phase interval Δp = 0.05 corresponds to a time interval Δt ≃ 1.6 ks.

A joint fit to the five QPE difference spectra results in χ² = 342 for 305 degrees of freedom. No signatures for a warm absorber (nor for any excess absorption with respect to the Galactic column density) is seen during the QPE evolution. As already discussed in Miniutti et al. (2019) QPEs exhibit a rather smooth temperature and luminosity evolution peaking at L_bol ≃ 3.3 × 10⁴² erg s⁻¹. The maximum temperature (kT = 111 ± 4 eV) is reached slightly earlier than the peak (kT = 103 ± 2 eV). In the upper panel of Fig. 18 we show the resulting relation between the QPE bolometric luminosity and the derived blackbody (rest-frame) temperature. We also show, as dotted lines, two L_bol ∝ T⁴ relations normalised to the extremal data points in the evolution, R1 and D2. While the initial rise (R1 to R2) lies on the L_bol ∝ T⁴ relation, the overall evolution is not consistent with constant-area blackbody emission but shows instead hysteresis. The hysteresis behaviour in the L_bol − T plane is fully consistent with that observed in the other QPE source eRO-QPE 1 (Arcodia et al. 2022) despite its very different timing properties, as well as in RX J1301.9+2747 (Giustini et al. in prep.), suggesting it is a common property of QPEs.

Fig. 18. QPE spectral evolution during XMM4. In the upper panel we show Lbol as a function of blackbody rest-frame temperature. The colour-code is that defined in Fig. 17. The evolution starts at R1 (lowest Lbol at ≃83 eV, red data point) and proceeds counter-clockwise. The two dotted lines are Lbol ∝ T4 relations normalised to the extremal of the QPE evolution (R1 and D2). The lower panel shows the estimated radius of the blackbody emitting region as a function of QPE phase. The re-scaled QPE profile is also shown for reference (light grey).

Under the assumption of blackbody emission, the area A of the emitting region can be estimated for each data point (e.g. at each phase during the QPE evolution) as A = L_bol σ⁻¹ T⁻⁴, where σ is the Stefan-Boltzmann constant. In the lower panel of Fig. 18, we show the estimated blackbody radius (assuming A = 2πR²) as a function of QPE phase. The derived blackbody radius is extremely small at all phases, as noted by King (2020) and Krolik & Linial (2022). The range of observed radii (4 − 12  ×  10¹⁰ cm) corresponds to ≃0.04 − 0.13 R_g (M_BH ≃ 4 × 10⁶ M_⊙) and ≃1 − 3 R_g (M_BH ≃ 3 × 10⁵ M_⊙) for the extremal M_BH values we derive from the spectral analysis in Sect. 5.1. We point out, however, that the derived blackbody size could be under-estimated by as much as ∼ one order of magnitude as we use a pure blackbody spectral model ignoring scattering, that is ignoring any X-ray colour-correction.

As was noted, the QPE spectra do not require any absorption in excess of the Galactic one, that is the warm absorber that is clearly present during the quiescent period disappears during QPEs. This could in principle be attributed to a sudden increase in ionisation during QPEs, or to the lower signal-to-noise of the phase-resolved QPE spectra with respect to the quiescent level one. In the former case, the bolometric luminosity associated with QPEs is (even at peak) only a tiny fraction (≤15%) of that of the quiescent level. However, the QPE luminosity fraction is drastically higher above ∼0.7 keV, which is the relevant band in terms of ionising luminosity as far as absorption features around that energy are concerned (see Fig. 8), with L_QPE ≃ 200 × L_quiescence at QPE peak. Such a sudden increase in ionising luminosity could then be responsible for the disappearance of the warm absorber during QPEs.

The total energy radiated during one QPE is of the order of E_QPE ≃ 6.7 × 10⁴⁵ erg. We do not have enough observations with QPEs to constrain the peak luminosity evolution on timescales longer than ∼1 yr. Hence, we assume that E_QPE is representative of the average QPE radiated energy at all times. Although QPEs cannot be interpreted as accretion-related without ambiguity (see discussion in Sect. 8.2), we estimate, for future reference, the corresponding accreted mass assuming an efficiency η = 0.1 which results into $M_{\mathrm{accr}}^{(\mathrm{QPE})}\simeq 3.7\times {10}^{-8}{M}_{\odot }$ per QPE or, considering ∼971 QPEs per year, $M_{\mathrm{accr}}^{(\mathrm{QPE})}\simeq 3.6\times {10}^{-5}{M}_{\odot }$ yr⁻¹. If QPEs were always present but undetected against higher quiescent level emission, the total mass accreted from QPEs in ∼12 yr would be ≤4.3 × 10⁻⁴ M_⊙, a very conservative upper limit.

8. Discussion

Before discussing some of the physical consequences of our findings, we provide a summary of all results reported here and in previous studies of GSN 069 to inform the discussion. The body of results presented here can be used in future studies to constrain the different theoretical models that have been (and likely will be) proposed for QPEs.

8.1. Summary of the main results

• QPEs in GSN 069 are abrupt flares occurring approximately every ∼9 h and lasting ∼1 h during which the X-ray count rate increases by up to two orders of magnitude with respect to the quiescent level, depending on the considered energy band (Miniutti et al. 2019). They were first discovered in December 2018 and last detected in January 2020. QPEs are then a transient phenomenon in GSN 069 with an observed life-time of ∼1.05 yr. However, considering that QPEs may have been missed in observations with high (quiescent) flux level, as well accounting for periods that were not covered by long enough exposures, the actual QPE life-time could be significantly longer.

• With the exception of the last observation where they are detected (XMM6), QPEs are characterised by alternating recurrence times and intensities. The contrast between the intensity of consecutive QPEs is correlated with the recurrence time between them, as shown in Fig. 4. Strong QPEs are always followed by long recurrence times (and vice versa). Remarkably, the correlation implies that consecutive QPEs have the same intensity when long and short recurrence times have equal duration. During the XMM6 observation, the alternating intensity is (roughly) preserved, but recurrence times increase monotonically suggesting that the system was experiencing some major change at the epoch of the XMM6 observation (see Fig. 3).

• QPEs measured in harder energy bands are stronger (with respect to the quiescent level), peak earlier, and have shorter duration than when measured at softer energies (Miniutti et al. 2019). When modelled with blackbody emission, the QPE spectral evolution shows hysteresis in the L_bol − kT plane (Fig. 18) and the maximum temperature, of the order of ≃110 eV, is reached slightly earlier than the luminosity peak. The blackbody emitting area is extremely small with radius evolving from 1 to 3 R_g at most if M_BH ∼ few × 10⁵ M_⊙ (although the size may be larger if scattering is important). The total energy radiated away is ≃6.7 × 10⁴⁵ erg per QPE during XMM4, and slightly higher (lower) during XMM3 (XMM5).

• During observations with QPEs, the quiescent level emission varies with a characteristic timescale equal to the average observation-dependent recurrence time, and folded light curves exhibit a characteristic quasi-periodic oscillation, or QPO (see Fig. 5). The fact that the quiescent level QPO period is equal to the recurrence time between consecutive QPEs suggests that the QPO is associated with every single QPE, independently on its strong or weak type. The peak of the quiescent level QPO lags the preceding QPE by 8–10 ks, and the time delay is well correlated with the observation-dependent average recurrence time. The QPE-QPO intensity ratio in the folded light curves is constant over the XMM3 to XMM5 observations spanning ≃160 d (and is of the order of ∼37, see Fig. 6).

• The quiescent level X-ray spectrum is consistent with accretion disc emission satisfying L_X ∝ T⁴ in the 0.2–2 keV band which suggests that L_X is a good proxy of the bolometric luminosity. This most likely implies a compact accretion flow with suppressed optical/UV emission. Indeed, by combining the X-ray data with two Hubble Space Telescope observations performed ∼4 years apart and quasi-simultaneously with two XMM-Newton exposures, we infer an outer disc radius R_out ≤ 140 R_g, with a statistical preference for a very compact disc with R_out ≃ 15 R_g. If the disc is really only ∼15 R_g in size, the UV emission must be due to a process different from intrinsic disc emission such as, for example, reprocessing. The black hole mass can only be constrained in the range M_BH = 0.3 − 4 × 10⁶ M_⊙ (for a non-spinning Schwarzschild black hole) and mostly depends on the assumed colour-correction for the disc emission.

• A warm absorber (mostly detected through a broad absorption feature around 0.7 keV) is present in the quiescent spectra throughout the 12 yr evolution, and it is consistent with being constant in both column density and ionisation. The warm absorber disappears during QPEs, possibly because of the ∼200 times higher X-ray luminosity above ∼0.7 keV at QPE peak, likely inducing a much higher ionisation and thus reduced absorption features equivalent width.

• The 12 yr-long evolution of the quiescent emission is consistent with repeating TDEs in GSN 069, with two events ∼9 yr apart peaking at L_bol ≃ 6 − 7 × 10⁴³ erg s⁻¹ and ≃4 − 5 × 10⁴³ erg s⁻¹ respectively when assuming a compact disc with R_out ∼ 15 R_g (the luminosity for the more extended case being ∼2.5 times higher). Assuming that ∼50% of the disrupted mass is bound to the SMBH and subsequently accreted, we estimate a disrupted mass of 0.46 ± 0.20 M_⊙ for the first TDE. The disrupted mass during the second event is more uncertain and it depends on how we consider the effect of the debris from the first event that still have to be accreted. We conservatively estimate a disrupted mass in the range of 0.06 − 0.38 M_⊙ for the second TDE.

• We detect a precursor X-ray flare before the second TDE, peaking at L_bol ≃ 3 − 5 × 10⁴³ erg s⁻¹. The precursor flare lasts ≃150 − 200 d in total, and precedes the peak of the second TDE by 400 − 500 d. A possible precursor flare is present also before the first TDE, but its detection is only tentative. If interpreted as an accretion flare, the radiated energy of the precursor flare in TDE 2 corresponds to an accreted mass of ≃0.006 ± 0.004 M_⊙.

8.2. Constraints on the QPE emission mechanism

As is discussed by Arcodia et al. (2022), models based on pure gravitational lensing (Ingram et al. 2021) are disfavoured due to the QPE energy-dependence (lensing being achromatic) and the difficulties in explaining simultaneously both the duration and amplitude of QPEs. Classical disc instability models are difficult to reconcile with the fast QPE timescales even for a black hole mass at the lower end of the allowed range (few × 10⁵ M_⊙) unless the disc properties are significantly different than standard. Moreover, the alternating pattern of QPE intensities and recurrence times is not naturally predicted, and may need to be introduced by an ad-hoc feedback mechanism. Instability-based models generally produce light curves with slow rise and fast decay, while QPEs in GSN 069 appear symmetric and, if anything, with a slower exponential decay as shown in Fig. 2 of Miniutti et al. (2019). However, modified models that significantly reduce the size of the unstable region and include disc-corona interaction or the presence of magnetic fields are worth exploring in the future (Sniegowska et al. 2020; Pan et al. 2022; Kaur et al. 2022) and, in fact, the QPE spectral evolution is consistent with expectations from instability models (Pan et al. 2022).

Another type of instability has been suggested by Raj & Nixon (2021) as a potential source of AGN variability, including QPEs. The idea is that accretion flows that are misaligned to the black hole spin can become warped due to Lense-Thirring precession, to the point that strongly warped discs may be unstable and break up into discrete rings. The rings then precess basically independently from each other, and their misalignment can give rise to shocks which may result in radiative losses and also in an enhancement of mass accretion rate on timescales that can be shorter than the viscous one. As shown by Raj & Nixon (2021), if the warp is unstable in the innermost regions, quasi periodic shocks as well as mass accretion rate variations are expected, which may give rise to QPEs (and, possibly, to QPOs). In order for the predicted timescales of the instability to match the observed QPE recurrence time in GSN 069, relatively large values of the black hole mass are favoured (∼10⁶ M_⊙), which is consistent with our findings (see Sect. 5.1). As the properties of the instability depend on the disc structure (e.g. viscosity α and scale-height H/R), the proposed mechanism might explain why QPEs are not always present in GSN 069 by assuming, for instance, that the accretion flow is stabilised against disc tearing at high mass accretion rates.

Most other QPE models invoke the interaction between a stellar-like body and the SMBH (or SMBH plus accretion flow) identifying the QPE quasi-periodicity with the orbital period of the stellar companion. The proposed models can be roughly separated into two classes in which QPEs are either due to collisions between the stellar companion and the accretion flow in a nearly circular orbit (Xian et al. 2021), or to episodic mass transfer events at each pericenter passage in a mildly or highly eccentric orbit (King 2020, 2022; Chen et al. 2022; Wang et al. 2022; Zhao et al. 2022; Krolik & Linial 2022; Lu & Quataert 2022). Below we discuss, with no ambition of completeness, some of the most relevant aspects of the two classes of proposed models in comparison with the rich phenomenology of GSN 069.

As shown in Fig. 18, the QPE bolometric luminosity does not follow the L ∝ T⁴ relation except, possibly, during the initial rise. Rather, hysteresis is clearly present in the L_bol − T plane, a likely general QPE property (Arcodia et al. 2022). Moreover, assuming blackbody emission, the resulting emitting area is very small, although we must point out that the inferred sizes may be underestimated by a factor of a few as the spectral model (a simple blackbody one) does not include scattering. At face value, this suggests that QPEs originate from a compact region that heats up rapidly, inflates, and cools down via X-ray emission while expanding. Qualitatively, QPEs appear then to be more consistent with shocks than with accretion onto the SMBH.

One class of QPE models where the X-ray emission is powered by shocks is that of collisions between an orbiting stellar body and the accretion flow occurring twice per orbit. Xian et al. (2021) show that such a model can explain the irregularities of the QPE arrival times and the alternating pattern of their recurrence times via the combination of small eccentricity and relativistic precession. Their orbital model successfully reproduces the QPE arrival times for an orbit with semi-major axis a ≃ 350 − 400 R_g and small eccentricity (e ≃ 0.06) crossing the accretion disc around a 3 × 10⁵ M_⊙ SMBH twice per orbit. The derived SMBH mass is consistent with the very low end of the allowed mass range in Sect. 5.1. Xian et al. (2021) propose that the orbiter is the remnant of the pTDE of a red giant star, in line with our analysis (see Sect. 8.4).

Further work is needed to account for the time delay between TDE 1 and the appearance of QPEs 4.5 − 8.5 yr later. Moreover, that shocks produced by impacts at such large radii can produce relatively high soft X-rays luminosity with negligible associated UV emission (Miniutti et al. 2019) needs to be studied in detail. Xian et al. (2021) do not discuss the alternating QPE intensity pattern, nor its relation with the alternating recurrence times. If, as would appear natural, strong (weak) QPEs are due to impact towards (away from) the observer, precession inevitably means that, depending on precession phase, strong QPEs would be followed by either long or short recurrence times as shown in the right panel of Fig. 4 of Xian et al. (2021). This behaviour is not consistent with QPE data, as strong QPEs are always followed by long recurrence times, at least for the events that we have observed (see Fig. 4).

A different class of models is based on the idea that QPEs are instead due to episodic mass transfer from a donor star at each pericenter passage through Roche lobe overflow, which can be further differentiated in models assuming QPE emission from shocks or from accretion onto the SMBH. Circularisation shocks are proposed as a potential QPE emission mechanism by Krolik & Linial (2022) and Lu & Quataert (2022). While Krolik & Linial (2022) ignore the presence of the accretion flow that gives rise to the quiescent level emission and invoke stream-stream interactions as a source of QPEs, Lu & Quataert (2022) consider the stream-flow interaction where they assume that the flow is due to previous mass transfer events. Stream self-interactions, as proposed by Krolik & Linial (2022), may also play an important role in this context, but it seems difficult that the incoming streams would always miss the existing accretion flow and interact first between themselves. As shown by Krolik & Linial (2022), shocks involving velocities of the order of the orbiter’s velocity at pericenter produce emission with a typical temperature that is consistent with that observed in QPEs. For both scenarios (stream-flow and stream-stream interaction) the energy budget also appears to be consistent with observations. Both proposed models invoke episodic mass transfer of a Sun-like star with M_* ≃ 0.5 − 1 M_⊙ in a mildly eccentric orbit. As discussed in Sect. 8.4, however, the properties of TDE 1 suggest that the remnant is unlikely to be a Sun-like star, and that a compact core, potentially surrounded by a remaining envelope, is more likely. Hence, if a (low-mass) Sun-like star is responsible for QPEs, it is likely unrelated to TDE 1 and TDE 2 and the system may be a triple.

QPEs due to episodic mass transfer from a compact object (e.g. a white dwarf or a compact stellar core) in a highly eccentric orbit, possibly captured through the Hills mechanism (Cufari et al. 2022), are proposed instead by King (2020, 2022), Chen et al. (2022), and Wang et al. (2022), while Zhao et al. (2022) consider mass transfer from the He envelope of a H-deficient post-AGB star. In all of these cases cases, the QPE X-ray emission is assumed to be powered directly by accretion onto the SMBH of the tidally stripped mass, leading to self-consistent solutions for the stellar structure and orbit.

The hysteresis in Fig. 18 is not necessarily inconsistent with an accretion origin for QPEs considering their highly dynamical nature. If a mass accretion rate perturbation propagates inwards, soft X-rays respond first on timescales associated with their typical emitting radii, while harder X-rays come later but evolve on faster timescales. Hence, at any given time, we observe the superposition of X-ray emission from different annuli on the disc, with hotter region having already responded, and colder ones still rising, which does not produce an obvious, global L_bol ∝ T⁴. The accretion-propagation scenario is in fact consistent with the energy dependence of QPEs, with higher energy X-rays peaking earlier and evolving faster than lower energy ones during QPEs as shown in Fig. 2 of Miniutti et al. (2019). On the other hand, the accretion timescale must be extremely fast, which is only possible if QPEs are induced by accretion from a narrow ring. This is consistent with a compact orbiter with very high eccentricity (and hence small pericenter distance) and, possibly, with a compact accretion flow whose extension is limited by the repeated passages of the orbiting star (King 2020, 2022). A compact donor on a highly eccentric orbit requires the black hole mass to be of few × 10⁵ M_⊙ at most for the object not plunge directly into the black hole. As a side note, we point out that X-ray flares with soft X-ray emission lasting longer than harder X-ray one (as observed in QPEs) may originate from magnetic energy release in the plunging region, that is at radii smaller than R_isco, as discussed by (Machida & Matsumoto 2003) making use of three-dimensional magnetohydrodynamical simulations, see in particular their Fig. 17.

Within the context of mass transfer scenarios, the alternating strong-weak intensities and long-short recurrence times have been addressed by King (2022) who propose that the very fast mass transfer timescale at pericenter induce an oscillatory behaviour around the evolutionary mean (basically set by gravitational radiation). We point out that such an idea is consistent with the properties shown in Fig. 4 as intensities and recurrence times appear to indeed oscillate around the equal-intensity, equal-recurrence time point that is naturally associated with the evolutionary mean. Another plausible mechanism, where changes in intensities and recurrence times are associated to the X-ray irradiation of the stellar surface at QPEs, has been instead proposed by Krolik & Linial (2022).

Finally, we note that disc tearing resulting from instabilities in a warped accretion flow may also naturally give rise to shocks between misaligned neighbouring rings in the innermost regions (Raj & Nixon 2021). Moreover, shocks are expected to drive mass inflow on fast timescales associated with narrow rings which, as discussed above, may be consistent with the X-ray spectral evolution of QPEs. As mentioned, the variability timescales can match the observed QPE period if the black hole mass is of the order of ∼10⁶ M_⊙, consistent with the black hole mass range in GSN 069. The dependency of the instability on the accretion flow properties (as well as on black hole spin) may produce enough diversity to account for the different characteristics of the observed QPE sources and could explain why QPEs are not always present during the evolution of GSN 069.

8.3. The quiescent level QPO

Proposed models should also account for the quiescent level QPO discussed in Sect. 4. Disturbances in the accretion disc are naturally produced in both the mass transfer and impacts scenarios with a characteristic timescale of ∼9 h. In this context, the QPE-QPO time delay we measure is not necessarily relevant, as the QPO observed at a given epoch may be due to events having occurred long time before. On the other hand, the stable QPE-QPO intensity ratio (Fig. 6), and the 1:1 correlation (Fig. 7) suggest that each QPO is indeed associated with the immediately preceding QPE.

Here we discuss one of the possible mechanisms that can produce the observed QPO: we assume (quasi-periodic) shocks on the accretion flow at radius R_sh producing QPE emission locally, and also driving a perturbation that propagates in the flow giving rise to enhanced quasi-periodic X-ray emission (the QPO). That outlined above is a plausible but not unique QPO origin, so that the conclusions reached below are model-dependent. Shocks at R_sh may be produced either by impacts on the accretion flow (twice per orbit) or by the interaction between incoming streams and the flow in the mass transfer scenario (with mass transfer events occurring once per orbit). In both cases, shocks are expected to roughly occur at radii of the order of R_sh ≃ R_p = (1 − e) a on the flow, where a is the orbit semi-major axis, and e the orbital eccentricity.

We shall first discuss mass transfer scenarios in which the QPE and QPO periods are equal to the orbital one, namely P ≃ 32 ks (see Table A.1). We note that the timescale for the propagation of the perturbation from R_sh inwards, t_prop, must be shorter than the period P of the perturbation in order not to smear out the QPO. The observed QPE-QPO time delay of 8 − 11 ks could perhaps be taken as a plausible estimate of t_prop, although the delay may also be affected by the time it takes to produce the QPE from the shock. The propagation of perturbations in the accretion flow is associated with the viscous timescale t_visc. Faster timescales are reached by propagating cooling or heating fronts in the disc (Ross et al. 2018), with t_front ≃ (H/R) t_visc. Although such a fast timescale is unlikely to be associated with the propagation of mass accretion rate fluctuations, we adopt conservatively t_front ≤ t_prop ≤ t_visc, baring in mind that we are only considering order-of-magnitude estimates here.

All accretion timescales are proportional to the dynamical one t_dyn = (G M_BH)^−1/2R^3/2 at radius R. When evaluated at R = R_sh ≃ (1 − e) a, t_dyn does not depend anymore on M_BH, so that, for a given period P, t_visc and t_front are a function of the orbital eccentricity e only. In order to estimate t_visc and t_front at R_sh, we set the viscosity parameter to α = 0.1 and the flow scale-height to H/R = 0.1, as seems reasonable for accretion flows produced by TDEs. The resulting t_visc and t_front are shown in Fig. 19 as a function of orbital eccentricity, where the upper (green) horizontal shaded area is the perturbation period (P ≃ 32 ks) while the lower (magenta) one is the shortest observed QPE-QPO time delay (8.3 ± 0.3 ks in XMM3) shown for completeness. The less stringent constraint on orbital eccentricity is obtained by requiring that the propagation timescale is shorter than the perturbation (in this case orbital) period P in order not to smear out the QPO (i.e. t_prop ≲ P). Even using the faster of the two timescales (t_front), the orbital eccentricity must be e ≥ 0.84. If, perhaps more reasonably, we adopt t_prop ∼ t_visc, one has e ≥ 0.965. A geometrically thinner accretion flow (H/R ≤ 0.1) requires even higher e. If the QPE-QPO time delay is assumed to represents the actual t_prop, higher eccentricities are required, see Fig. 19.

Fig. 19. Viscous and front propagation timescales tvisc and tfront evaluated at Rsh as a function of orbital eccentricity. The upper (green) horizontal shaded area represents the period of the shock-induced perturbation (equal to that of the QPE or QPO, i.e. to the orbital period in the assumed scenario), while the lower (magenta) one is the shortest observed time delay of the QPO with respect to the immediately preceding QPE (XMM3).

Hence, under the assumption that the QPO is due to a propagating perturbation caused by each QPE and that QPEs are associated with an orbiting body, the QPE-inducing body must describe a highly eccentric orbit. In mass transfer scenarios, the orbiting object has to fill its Roche lobe at pericenter. As discussed by King (2020, 2022), Chen et al. (2022), and Wang et al. (2022) such a condition, when combined with a high eccentricity (i.e. small R_p), implies a compact white-dwarf like object. In particular, assuming a typical stellar core mass M_* ≃ 0.3 M_⊙, and the less stringent constraint on eccentricity (e ≃ 0.84), one has R_* ≃ 0.1 R_⊙, allowing for the presence of a significant envelope bound to the core which may be then tidally stripped at pericenter passage, giving rise to QPEs. Higher eccentricities result in smaller R_*. Such a compact object in a highly eccentric orbit has pericenter distance of only few−few tens of R_g (Maguire et al. 2020; Cufari et al. 2022) which, as mentioned, favour a low black hole mass in order for the object not to be swallowed by the black hole. In any case, the orbit is highly relativistic, and the flash of X-ray emission at each pericenter passage (QPE) represents a probe sending back to us a signal (that can be accurately timed) from a highly distorted region of spacetime. QPEs may thus enable us to constrain the strong field regime of General Relativity at distances from the SMBH event horizon that appear difficult to probe by other means.

On the other hand, if QPEs are produced by impacts, the same arguments apply with the exception that the orbital period is twice as large, increasing all accretion timescales. Within this scenario, the periodic perturbation (each impact) has still the same period of ≃32 ks (as the model assumes two impacts per orbit), so that the eccentricity needs to be even higher than in the mass transfer scenario. However, one of the most attractive aspects of the impacts model, is that it can explain the alternating recurrence times of QPEs by invoking a small orbital eccentricity (Xian et al. 2021), which is therefore inconsistent with the proposed QPO interpretation. Moreover, and independently from the nature of the QPO, we point out that, for small eccentricities, the pericenter distance (which is also roughly the radial distance where impacts occur) is consistent with the maximum allowed disc radius that is needed not to over-predict the UV emission (R_out ≃ 140 R_g) only for M_BH ≳ 1 × 10⁶ M_⊙, while the orbital solution of Xian et al. (2021) implies a best-fitting black hole mass of 3 × 10⁵ M_⊙. One further potential weakness of this model is that, while it can likely explain the regular QPE patterns observed in GSN 069 and eRO-QPE2, it appears difficult that it alone can provide the diversity that is needed to account for the less regular QPEs in RX J301.9+2747 and eRO-QPE1 (Giustini et al. 2020; Arcodia et al. 2022). Mass transfer scenarios may be more naturally associated with different stellar and orbital properties leading to a broader range of observable properties, see for example Metzger et al. (2022).

As mentioned, the constraints derived above strongly depend on the adopted interpretation of the QPO. As such, our conclusions should be treated as model-dependent and therefore do not represent a robust, conclusive result. Other possible QPO origins should be explored and contrasted with the data. One possibility is that QPEs and QPOs are produced by a combination of the QPE models proposed so far. For instance, if the orbit is eccentric and the accretion flow does not extend much beyond R_p, star-disc collisions may occur only once per orbit, depending on the relative inclination, and could give rise to one of the two phenomena (QPE or QPO), with the other being associated with mass transfer at pericenter. In this case, the constraint on e is relaxed, and mild eccentricities are likely possible as well. Within this scenario, and depending on eccentricity and relative orientation of the orbital and accretion planes, up to three events per orbit might be expected (two impacts and one mass transfer event) which may help explaining the complex timing behaviour of some QPE sources, most notably RX JJ301.9+2747 (Giustini et al. 2020). In mass transfer scenarios, we also note that low eccentricities would likely increase the fraction of the orbit where the star over-fills its Roche lobe, perhaps accounting for the much longer-lived QPEs observed in eRO-QPE 1 (Arcodia et al. 2022) and, potentially, also for the QPO-like variability in 2XMM J123103.2+110648 (Webbe & Young 2023). Tidal interactions close to pericenter may also induce a time-dependent warp of the existing accretion flow which then precesses and possibly accounts for the observed QPO. In this context, disc tearing may be relevant (Raj & Nixon 2021), and solutions in which QPEs are produced, for example, by shocks between neighbouring rings should be explored further. The QPO may be due to the modulation of the X-ray emission from a dominant ring at its own precession timescale, or to the propagation of mass accretion rate fluctuations in the inner regions caused by ring-ring interactions. We encourage exploring these and other possibilities theoretically and numerically to constrain the different proposed models, deepening our understanding of the QPE phenomenon in GSN 069 and the other QPE sources.

8.4. The TDE-QPE connection

The overall evolution of GSN 069 over the ≃12 yr period probed by our campaign is shown in Fig. 20 together with the best-fitting model discussed in the Sects. 6.1 and 6.2. As shown there, Phase1 and Phase2 are well described by the typical rise-and-decay light curve of TDEs, and we therefore re-name Phase1 and Phase2 as TDE 1 and TDE 2 respectively. We also include in Fig. 20 a possible precursor X-ray flare before TDE 1: this is not a fit, but rather a re-scaled version of the flare seen at the beginning of TDE 2 that accounts well for the first two data points that, as noted in Sect. 6.1, lie above the best-fitting model for TDE 1.

Fig. 20. Overall Lbol evolution over the ∼12 yr probed by our Swift and XMM-Newton campaign, together with the best-fitting models during Phase1 and Phase2, re-named here as TDE 1 and TDE 2 respectively. The y-axis scales are the same as those defined in the caption of Fig. 12.

The evolution in Fig. 20 strongly suggests repeating TDEs in GSN 069 whose peaks (we ignore here the precursors) are separated by ≃3340 d (≃9.1 yr). Similar recurrent optical or X-ray flares that can be interpreted as repeating TDEs have been recently reported in ASASSN-14ko (Payne et al. 2021), eRASSt J045650.3-203750 (Liu et al. 2023), and AT2018fyk (Wevers et al. 2023), with source-dependent timescales ranging from ≃100 − 200 d to ≃3 yr. One striking property of both TDEs in GSN 069 is that their X-ray emission evolves on much longer timescales than typical TDEs. The X-ray flux decays by a factor of ∼6.8 in ∼9 yr, and the corresponding bolometric luminosity drops by an even smaller factor of ∼3 (see Figs. 9 and 13). This very slow decay is rarely seen in TDEs (Lin et al. 2022). When the TDE 1 (TDE 2) decay is described with an exponential, its e-folding timescale is τ ≃ 7.6 yr (≃3.4 yr), more than one order of magnitude longer than typical for both X-ray and optical TDEs (Komossa 2015; van Velzen et al. 2021). The width around the peak is also, in both cases, much longer than in typical TDEs (about ≃2.5 − 4 yr and ≃1.5 − 3 yr for TDE 1 and TDE 2 respectively). The years-long timescales characterising TDE 1 and TDE 2 can be compared with the so-called fallback time t_fb that sets the typical timescale of the fallback rate evolution and that can be rewritten as (Gezari 2021):

$$\begin{array}{c}\hfill t_{\mathrm{fb}}\simeq 41{\left(\frac{R_{\ast }}{R_{\odot }}\right)}^{3/2}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1}{\left(\frac{M_{\mathrm{BH}}}{{10}^6{M}_{\odot }}\right)}^{1/2}\mathrm{d},\end{array}$$(2)

which, for the disruption of a Sun-like star, is ≃20 − 80 d in the allowed range of black hole masses in GSN 069 (0.3 − 4 × 10⁶ M_⊙). Although the timescales for partial disruptions are likely a factor of a few longer than those characterising full TDEs (Guillochon & Ramirez-Ruiz 2013; Nixon et al. 2021), the disruption of a Sun-like star does not appear to be consistent with the observed timescales. On the other hand, the disruption of an evolved star comprising a compact core surrounded by a tenuous envelope may be viable. Indeed, assuming a typical mass M_* ∼ 1 M_⊙, the disruption of a a star with R_* ≃ 5 − 12 R_⊙ (depending on the assumed M_BH), gives t_fb ≃ 1 000 d, roughly consistent with the observed timescales. Hence, the years-long timescales associated with TDE 1 indicate the likely disruption of an evolved star extending out to at least a few times the Solar radius (MacLeod et al. 2012).

The presence of TDE 2 suggests that TDE 1 was only partial and that a remnant survived TDE 1. As shown by MacLeod et al. (2012), the mass loss for the disruption of the envelope of an evolved star ranges from ∼0.1 M_⊙ to ∼0.7 M_⊙ depending on stellar structure and penetration factor β = R_t/R_p, where R_p is the pericenter distance of the orbit, and R_t is the tidal radius

$$\begin{array}{c}\hfill R_{\mathrm{t}}\simeq 47{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1/3}\left(\frac{R_{\ast }}{R_{\odot }}\right){\left(\frac{M_{\mathrm{BH}}}{{10}^6{M}_{\odot }}\right)}^{-2/3}{R}_g.\end{array}$$(3)

The disrupted mass we derive at TDE 1 (0.46 ± 0.20 M_⊙) is thus fully consistent with the partial disruption of the envelope of an evolved star.

As discussed in Sect. 6.2, TDE 2 is also characterised by years-long timescales, only slightly shorter than those associated with TDE 1. At face value, this would suggest the disruption of an extended object also at TDE 2, consistent with partial/full disruption of the remnant of TDE 1, if the remnant still comprises a significant envelope. TDE 1 and TDE 2 might then be due to the repeated tidal stripping of the envelope of an evolved giant star at each pericenter passage (MacLeod et al. 2013). If so, the star is on a ∼9 yr, likely eccentric orbit around the central SMBH.

If the repeating TDEs in GSN 069 are due to tidal stripping of the envelope of a star on a ∼9 yr orbit, the remnant that survived TDE 1 cannot be responsible for QPEs because of the different required orbital periods (∼9 yr versus 9 − 18 h, depending on the adopted QPE model). The system would then be a triple comprising the SMBH, a giant star on a ∼9 yr orbit whose envelope is tidally stripped at each pericenter passage, and a further stellar companion on a much shorter-period (9 − 18 h) orbit producing QPEs. One further possibility is that QPEs are not associated with an orbiting companion and are due, for instance, to shocks between misaligned precessing rings in the accretion flow, formed as a consequence of disc tearing (Raj & Nixon 2021). In such a case, the system only comprises the SMBH and an orbiting evolved star on a ∼9 yr orbit. The scenario in which the two repeating TDEs are associated with the partial TDE of the envelope of an evolved star on a ∼9 yr orbit is testable, as it predicts a third TDE ∼9 yr after TDE 2. Monitoring GSN 069 in the soft X-rays during the next few years is therefore crucial to confirm or falsify this possibility.

However, the partial nature of TDE 1, that is the presence of a surviving remnant, suggests that a different scenario in which TDE 1, QPEs, and TDE 2 are all associated with a common progenitor should be explored. An evolved giant star could be partially disrupted at TDE 1, leaving a remnant able to survive thousands of pericenter passages and to give rise to QPEs before being driven to its own partial/full disruption at TDE 2 by orbital evolution. If TDE 1 is the pTDE of an evolved star, as argued above, this is in fact possible. The condition to have the envelope (or part of it) disrupted at TDE 1 leaving a surviving remnant is that the penetration factors of the progenitor star (β_*) and of the remnant (β_rem) satisfy β_* ≥ 1 and β_rem < 1, although this is only approximately true as the critical β is slightly different from unity in the case of pTDEs and, more generally, β also depends on the stellar structure (Guillochon & Ramirez-Ruiz 2013; Mainetti et al. 2017). Assuming that the orbit of the remnant is unperturbed (i.e. the progenitor star and the remnant have the same pericenter distance), one has

$$\begin{array}{c}\hfill {\beta }_{\mathrm{rem}}/{\beta }_{\ast }=\left(\frac{R_{\mathrm{rem}}}{R_{\ast }}\right){\left(\frac{M_{\ast }}{M_{\mathrm{rem}}}\right)}^{1/3}.\end{array}$$(4)

In the limiting case of a bare compact core remnant (with say M_rem ≃ 0.3 M_⊙ and R_rem ≃ 0.014 R_⊙) and a relatively common red giant progenitor star (M_* ≃ 1 M_⊙ and R_* ≃ 12 R_⊙), one has β_rem/β_* ≃ 1/574, so that there is a significant range of parameters that can satisfy the condition for which a remnant survives the next pericenter passage even if part of the envelope is still bound to the core. This also means that the disruption of the envelope alone (or part of it) can occur even for β_* ≫ 1, namely long before the star reaches R_p. The number of orbits that the remnant can complete before being driven to disruption by further orbital (and stellar size) evolution depends on the initial β_rem as well as on the details of its subsequent evolution (gravitational wave emission, mass transfer, tidal heating, etc).

However, such an interpretation implies that the initial evolved star is on a ∼9 h (or ∼18 h) around the black hole. As noted by Cufari et al. (2022), it is difficult (basically impossible) to place a large star on an orbit with such a short period even invoking the break-up of an initial binary system on a very tight orbit (Hills 1988). One possibility is that the evolved star was captured long ago in a significantly wider (longer-period) orbit that hardened due to orbital evolution leading the star to its tidal radius on an orbit with a ∼9 − 18 h period long after being captured by the SMBH.

On the other hand, a significantly more compact object, such as a white dwarf (WD) may be placed on a very short period, highly eccentric orbit via the Hills mechanism (Cufari et al. 2022), although Lu & Quataert (2022) argue that the rate of delivering WDs to the required orbit is too small to explain the observed QPEs. Within this context, Wang et al. (2022) have proposed that TDE 1 may represent the disruption of a WD envelope significantly inflated due tidal heating that can induce runaway fusion in the envelope (Fuller & Lai 2012). The inflated envelope may produce TDE 1 at larger distances (and hence with longer timescales) than expected from the partial disruption of an unperturbed WD. While attractive because it potentially provides a link between TDE 1, QPEs, and TDE 2, the model faces at least two challenges: firstly, the disrupted mass that we estimate at TDE 1 (0.46 ± 0.20 M_⊙) is much too large for the disruption of a typical WD envelope, and, secondly, the surviving remnant would likely be just the WD itself, whose further disruption at TDE 2 should in principle proceed on much shorter timescales than observed (Maguire et al. 2020). The latter inconsistency may perhaps be cured by noting that TDE 2 does not occur in a clean environment because of the presence of the accretion flow from TDE 1. Hence, the TDE 2 timescales may not reflect accurately the structure of the disrupted body, but could rather still be dominated by the fallback rate from TDE 1 which is likely dominant in terms of bound mass, as discussed in Sects. 6.1 and 6.2.

As discussed above, weak and irregular QPEs are observed during XMM6. As shown in Fig. 21, XMM6 is performed after TDE 2, basically at the end of the precursor X-ray flare. The irregular and weak nature of QPEs may be understood in terms of the strong perturbation of the system (including the accretion flow) caused by TDE 2. For example, if QPEs are produced by shocks on the accretion flow (see e.g. Sect. 8.2), the disruption (or, at least, disturbance) and subsequent re-formation of the accretion flow due to TDE 2 is bound to perturb the regularity of QPEs. Moreover, if QPEs are due to mass transfer at pericenter, they are expected to be initially weak and then grow in amplitude with time, as discussed by Zalamea et al. (2010) in the case of a white dwarf donor. The increase of the continuum level during the rise of TDE 2 might have reduced the contrast of weak, XMM6-like QPEs with respect to the quiescent continuum level, possibly preventing us from detecting QPEs after XMM6.⁵ However, as the continuum is now decaying again (and the QPE amplitude is possibly growing), high-contrast QPEs may soon re-appear in GSN 069.

Fig. 21. Portion of the Lbol light curve comprising the final part of TDE 1 and the overall TDE 2 evolution. Black circles (green squares) denote XMM-Newton (Swift) observations. Only XMM-Newton observations (whose observation number is also indicated) are long enough to detect QPEs. The shaded area highlight time-intervals during which high-amplitude, regular QPE are detected (XMM3 to XMM5, red) and those in which weak QPEs with irregular recurrence times are observed (XMM6, blue). No QPEs are detected after XMM6, although the higher level of the quiescent emission might have decreased significantly their contrast in subsequent observations.

In mass transfer scenarios for QPEs, the detection of of a new QPE phase in the future (of which QPEs in XMM6 may be considered the beginning) would then imply that an orbiter is still present. The QPE properties of the new (putative) QPE phase could be used to infer whether the orbiter producing them is related to that inducing the previous set of QPEs. If TDE 2 partially disrupted the object that is responsible for previous QPEs, mass exchange at TDE 2 should result into different QPE recurrence times (associated with the orbital period), while if TDE 2 is due to tidal stripping of the envelope of a giant star on a ∼9 yr orbit and QPEs are associated with an independent orbiting companion, the new QPE phase should have retained the same recurrence times as before. Note, however, that the latter scenario implies a triple system whose dynamical stability should be further studied.

8.5. The precursor X-ray flare

The X-ray flare detected at the beginning of TDE 2 is a novel addition to the rich phenomenology of X-ray TDEs evolution. Its peak luminosity is similar to the subsequent TDE 2 peak (most likely associated with the accretion phase), but its evolution timescales is much shorter (50 − 90 d as opposed to ≃1 − 2 yr). As mentioned, a precursor X-ray flare preceding TDE 1 is also tentatively observed as shown in Fig. 20.

Under the assumption that the repeating TDEs are due to the tidal stripping of a giant star’s envelope at pericenter, the two TDE light curve should be similar, as observed. In this case, it is tempting to associate the precursor flare (or flares) to the circularisation phase in which the TDE tidal streams self-interact and circularise during the formation of the accretion disc. As shown by Hayasaki et al. (2016), if radiative cooling is efficient, the debris circularise in a geometrically thin ring-like accretion flow and the dissipated energy during circularisation can give rise to observable precursor flares (with respect to the subsequent accretion phase) that can be luminous enough to be detectable, especially for parabolic TDEs (and likely highly eccentric ones), with an upper limit comparable to the Eddington luminosity. A similar conclusion is also reached by Bogdanović et al. (2014) in the case of the pTDE of a red giant envelope (which is likely relevant for GSN 069), showing that the luminosity due to circularisation can exceed the fallback luminosity at early times. Precursor X-ray flares are also found in kinematic simulations by Rossi et al. (2021) if the circularisation timescale is significantly shorter than the characteristic fallback timescale at which the accretion phase peaks.

On the other hand, assuming instead that TDEs and QPEs are all related to the same progenitor, we note that the precursor flare relatively short timescales are roughly consistent with the partial or full disruption of a compact stellar core plus small remaining envelope (Gezari et al. 2012; Strubbe & Murray 2015; Law-Smith et al. 2017). In fact, timescales ≲100 d, as observed for the precursor flare, are also compatible with typical optical TDEs in general (van Velzen et al. 2021). The subsequent (accretion-dominated) evolution of TDE 2 might evolve on timescales still dictated by the fallback rate from TDE 1 (dominant with respect to TDE 2 in terms of disrupted mass), which could account for the years-long TDE 2 overall timescales once the fact that TDE 2 does not occur in a clean environment is considered. It is also worth to point out that the disruption of a WD-like remnant may be associated with tidal compression and subsequent ignition, leading to a thermonuclear explosion accompanied by a relatively short-lived thermal-like soft X-ray flare (Rosswog et al. 2009).

As mentioned, the two scenarios outlined above may be tested with future observations. A third TDE ∼9 yr after TDE 2 would confirm the tidal stripping of the giant star’s envelope at pericenter, and the properties of a putative new QPE phase may be used to infer the origin of the stellar-like companion that is responsible for QPEs, as discussed in the previous Sect. Other theoretical origins of the precursor flare are possible, and they should be explored in future theoretical and numerical works adopting parameters that can account for the overall evolution (and QPEs) in GSN 069. As discussed here, understanding the origin of the precursor flare(s) may help elucidate the relationship between QPEs and repeating TDEs in GSN 069. We stress, once more, that the presence of a pre-existing accretion flow due to the accretion of debris from TDE 1 should be properly taken into account when describing the overall evolution and, in particular, the properties of TDE 2.

9. Conclusions

We study the long-term evolution of GSN 069 and the properties of its QPEs using 11 pointed XMM-Newton, one Chandra, and 34 Swift X-ray observations obtained in the 12 yr since its first X-ray detection in July 2010 during an XMM-Newton slew. QPEs in GSN 069 are a transient phenomenon observed between December 2018 and January 2020. However, considering that QPEs might have been missed in observations with high quiescent flux level, and accounting for observational gaps, the QPE life-time could be considerably longer.

We report repeating TDEs in GSN 069 with two events separated by about 9 yr. The presence of TDE 2 suggests that the first event (TDE 1) was a partial TDE. We show that TDE 1 is most likely produced by the partial disruption of an evolved star, leaving a compact stellar core possibly surrounded by a remaining envelope. We detect a precursor X-ray flare 400 − 500 d before the TDE 2 peak. The precursor might be related to the circularisation phase during the formation of the accretion flow. A similar precursor flare may be present ≃500 d before TDE 1, but its detection is only tentative.

The X-ray spectral evolution during QPEs suggests an origin in shocks, consistent with a compact region that heats up rapidly, and then cools down by emitting X-rays while expanding by a factor of ∼3 in radius in ≲1 h. Shocks may be produced by impacts between the accretion flow and an orbiting stellar companion, the interaction between the incoming stream(s) and the flow in mass transfer scenarios, or the interaction between neighbouring misaligned rings in a warped, unstable disc. A purely accretion-based origin of QPEs cannot be discarded, and would imply a very compact, ring-like emitting region for QPEs. In any case, we stress that the presence of a pre-existing accretion flow from TDE 1 should be taken into account in future realisations of QPE models.

During observations with QPEs, the quiescent level emission exhibits a QPO with a period equal to the separation between consecutive QPEs. Assuming that the QPO is produced by mass accretion rate fluctuations caused by shocks and that QPEs are associated with an orbiting body, the eccentricity of the QPE-inducing orbiter must be e ≳ 0.84, likely ruling out star-disc collisions scenarios associated with a low-eccentricity orbit. However, this conclusion strongly depends on the assumed QPO origin and cannot be considered as conclusive yet. We encourage further theoretical work that, informed also by our observational results, should ideally be able to account for the QPEs and QPO properties in a self-consistent manner.

As for the two TDEs that we observe ∼9 yr apart, and their connection with QPEs, we propose two possible scenarios. In the first, all events (TDE 1, QPEs, and TDE 2) are associated with the same progenitor (an evolved star), that is QPEs are due to the remnant of TDE 1 that survives thousands of pericenter passages before being partially or fully disrupted at TDE 2. In this case, as it is difficult to place an evolved star on a very short-period orbit, the progenitor could have been captured long before TDE 1, and been driven to its tidal radius by subsequent orbital evolution. The long timescales associated to TDE 2 (only marginally shorter than the TDE 1 timescales) are somewhat difficult to account for if TDE 2 disrupted a significantly more compact object than TDE 1, as expected within this scenario. However, the disrupted mass at TDE 2 likely interacts with the existing accretion flow (from TDE 1). Since the disrupted mass at TDE 2 is likely comparable to that of the region of the disc with which it interacts, the disc can be significantly disturbed, or even disrupted. The disc has then to re-stabilise, or re-form, by accreting not only the debris from TDE 2 after circularisation, but also the remaining bound mass from TDE 1. It is therefore possible that the TDE 2 timescales are still dominated by the TDE 1 fallback rate, as the disrupted mass at TDE 1 is likely dominant.

A second scenario assumes that the two TDEs represent repeating partial disruptions of the envelope of an evolved star at each pericenter passage so that the star is on a ∼9 yr orbit, thus removing the difficulty to place a large star on a very short-period orbit and also accounting naturally for the similar shape and timescales of TDE 1 and TDE 2. If so, QPEs cannot be produced by the remnant surviving TDE 1. A QPE-inducing orbiter could have been captured, possibly through the Hills mechanism, at a different epoch on its much shorter-period orbit (P ≃ 9 − 18 h, depending on the adopted QPE model). The stability of the putative triple system needs to be studied. On the other hand, QPEs may be induced not by an orbiter, but rather by shocks associated with disc tearing in a warped, unstable disc, thus removing the need for a triple system in GSN 069.

We note that, in all cases, as the tidal radius at TDE 2 is reached smoothly from larger radii, TDE 2 is most likely another partial rather than full TDE. If so, a third event is expected in the future. Would a third TDE occur ∼9 yr after TDE 2, the scenario in which the two repeating TDEs in GSN 069 represent the partial disruption of the envelope of an evolved star on a ∼9 yr orbit will be proved right. Future observations may also detect a new QPE phase whose beginning could be identified with the last observation in which QPEs are seen (XMM6, just after TDE 2). In mass transfer scenarios, the detection of clear new QPE phase would signal than an orbiter is still present after TDE 2. The properties of the putative new QPE phase (most notably the QPE recurrence time, i.e. orbital period) will significantly help to constrain the origin of the orbiter, namely whether it is a new remnant that survived TDE 2, or the same object that was producing previous QPEs. Changes in the QPE recurrence time may also help constraining the alternative disc tearing scenario, as the instability strength, location, and timescales are expected to vary with disc properties.

Future X-ray observations of GSN 069 on both long and short timescales will most likely allows us to pinpoint the origin of repeating TDEs, of QPEs, and of the possible TDE-QPE connection in this galactic nucleus, with consequences for QPE models in the other galaxies where QPEs have been identified so far. Our work provides a comprehensive summary of the rich phenomenology of GSN 069 that can be used to inform further theoretical and numerical studies on the origin of QPEs. QPE models are evolving fast, and will ideally be able to not only reproduce the most relevant properties of one specific source, but also to account for the diversity of QPE phenomenology that is starting to be revealed as new QPE sources are discovered and new data are accumulated.

Acknowledgments

We thank the referee for a constructive review of our work that contributed to improve our paper. GM warmly thanks Andrew King for valuable discussions. GM, MG, and RA also thank A. Sesana, E. Bortolas, A. Franchini, M. Bonetti, and A. Lupi for a pleasant and useful QPE brain-storming session in Milan. Based on observations obtained with XMM-Newton an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester, and we thank the Swift science team for organising and scheduling a series of unanticipated TOO observations of GSN 069. The scientific results reported in this article are also based on observations made by the Chandra X-ray Observatory. MG is supported by the “Programa de Atracción de Talento” of the Comunidad de Madrid, grant number 2018-T1/TIC-11733. SB acknowledges financial support from the Italian Space Agency under grant ASI-INAF 2017-14-H.O.

References

Appendix A: Light curves model

Appendix B: Quiescent level variability

In the upper panels of Fig. B.1, we show light curves from all XMM-Newton observations in which QPEs are observed. We ignored the Chandra data as the quiescent level was barely detected, making it difficult to study its variability. As all data are from the XMM-Newton EPIC-pn detector, there was no need to restrict the analysis to energies above 0.4 keV, and we used the full 0.2-1 keV band. The light curves were modelled with a constant and a series of Gaussian functions describing the quiescence and QPEs respectively. The lower panels of Fig. B.1 show the resulting residual light curves that have been re-binned for visual clarity to highlight the relatively long-term variations (tens of ks timescales) of the quiescent level emission. In all cases but the XMM6 observation, the quiescent level exhibits broad excess emission following most QPEs, regardless of their strong or weak type. Excluding XMM6, an excess can be clearly seen after nine out of the ten QPEs that are followed by a significant exposure. The only QPE that is not followed by any excess if the first during XMM5, which is also the weakest QPE detected in the XMM3 to XMM5 observations (see Table A.1). The dashed and solid lines represent models describing the residual light curves either with a constant (dashed lines) or with an additional periodic function (solid lines) with period fixed at the average observation-dependent recurrence time between consecutive QPEs. The sinusoidal models represent a very significant improvement with respect to the constant one in all cases. In particular, re-fitting the original light curves with the addition of a sine function improves the statistical result by Δχ² = −114,   − 261,   − 63 for the XMM3, XMM4, and XMM5 observations for −2 degrees of freedom, see Table A.1 for the statistical quality of the baseline model fits when the sine function is not included.

Fig. B.1. Quiescent level variability. In the upper panels, we show the 0.2-1 keV light curves from all XMM-Newton observations with QPEs, together with their best-fitting models comprising a constant and a series of Gaussian functions to describe the quiescent level and QPEs respectively. The lower panels show the corresponding residual light curves re-binned by a factor of 8, together with two models, a constant (dashed lines) and a sine function (solid lines) with a period equal to the average observation-dependent recurrence time.

All Tables

All Figures

Fig. 1. XMM-Newton (EPIC-pn) and Chandra (ACIS-S) background-subtracted light curves from all observations with QPEs in a common 0.4–1 keV band. The Chandra light curve has been re-scaled to the XMM-Newton pn effective area to better display the QPE amplitude evolution and to ease comparison between data from different detectors. We use time bins of 200 s and 500 s duration for XMM-Newton and Chandra data respectively. All panels have a common y-axis range.

In the text

Fig. 2. Time evolution of the QPE intensity since first detection. The QPE intensity is measured as the average of the best-fitting Gaussian normalisation in each observation (i.e. we ignore here the difference between strong and weak QPEs). The colour scheme is the same as in Fig. 1. The dotted line represents the best-fitting linear relation (corresponding to a decay of about 0.22 cts s−1 per 100 d). We also show (cross) the continuum level during the first observation with no detected QPEs (XMM7).

In the text

Fig. 3. Alternating QPE intensities and recurrence times. The QPE intensity (N) for all observations with QPEs is shown in the upper panels. The recurrence time between consecutive QPEs is shown in the lower panels (each data point is placed at a time corresponding to half the separation between consecutive QPEs). The y-axis is always in the same range, except for the irregular XMM6 observation. Dashed lines in the XMM4, Chandra, and XMM5 panels are sine functions with period equal to twice the averaged, observation-dependent recurrence time, and with an exponentially decaying amplitude that most likely indicates a common long-term modulation. The lines are intended to guide the eye rather than to provide fits to the data.

In the text

Fig. 4. Consecutive QPEs intensity ratio and recurrence time correlation. The ratio between the intensity of consecutive QPEs is shown as a function of the recurrence time between them (normalised to the sum of consecutive long and short intervals Tsum). The horizontal line (Ni /Ni + 1 = 1) separates events where strong QPEs precede weak ones (upper half of the figure), from those where weak QPEs precede strong ones (lower half). The $T_{\mathrm{rec}}^i/T_{\mathrm{sum}}=0.5$ vertical line separates short from long recurrence times. The (red) dotted line show the best-fitting linear relation between the two quantities.

In the text

Fig. 5. Quasi-periodic oscillation of the quiescent level. In the upper panels, we show the 0.2–1 keV light curves from the XMM3, XMM4, and XMM5 folded at the average observation-dependent recurrence time, and the corresponding best-fitting baseline model (constant plus Gaussian functions, dotted). We show the 0–2 phase interval for visual clarity. In the middle panels, we show the resulting residual light curves together with a sine function with period fixed to 1 (in phase space) to guide the eye. The vertical lines represent the phase of the primary QPE (dotted) and of the peak of the sine function (dashed). In the lower panel, we show the residual for a fit of the original folded light curve that comprises a sine function with period fixed to 1 (i.e. to the folding timescale) in addition to the baseline model.

In the text

	Fig. 6. Evolution of the QPE intensity (upper panel) and of the ratio between the QPE and QPO intensities (lower panel). The latter show that the QPE and QPO decay together preserving a ratio of 37.2 ± 1.3 (horizontal lines).
In the text

	Fig. 7. QPE-QPO time delay. We show the time delay of the QPO peak with respect to the preceding QPE as a function of the average observation-dependent recurrence time. The dotted line is a linear function with slope fixed to 1, the best-fitting one being 0.95 ± 0.15.
In the text

	Fig. 8. Residuals (in terms of σ) when the merged XMM2 to XMM11 spectrum is fitted without (upper panel) and with (lower panel) a warm absorber component. We only show the most interesting energy range close to the most prominent absorption feature at ≃0.7 keV for better visual clarity. Data have been slightly re-binned for visualisation purposes only.
In the text

Fig. 9. X-ray flux, luminosity, and temperature evolution. From upper to lower panel, we show the observed 0.3–1 keV X-ray flux, the 0.3–1 keV intrinsic X-ray luminosity, and the (rest-frame) disc temperature evolution, as derived from our best-fitting DISKBB models to the 11 XMM-Newton pointed observations. The first data point in the upper and middle panels refers to the slew observation whose observed count rate has been converted into X-ray flux and luminosity by using the best-fitting model from the nearest XMM-Newton pointed observation (XMM1). In all panels, the shaded areas highlights the period during which regular (red, XMM3 to XMM5) and irregular (blue, XMM6) QPEs are detected.

In the text

	Fig. 10. Intrinsic 0.2 − 2 keV (black circles) and 0.3 − 1 keV (grey squares) X-ray luminosity as a function of rest-frame temperature for all 11 XMM-Newton observations. The dotted lines are best-fitting relations assuming LX ∝ Tq.
In the text

	Fig. 11. UV light curves from OM data obtained during the XMM-Newton observations. The UVW2 (λeff ≃ 2 120 Å), and UVW1 (≃2 910 Å) fluxes are re-scaled by factors of 1.1 and 0.8 respectively for visual clarity, while the UVM2 (≃2 310 Å) remains as is.
In the text

Fig. 12. Bolometric luminosity evolution as inferred from OPTXAGNF best-fitting models to the 11 XMM-Newton observations. The left y-axis corresponds to the compact disc (Rout ≃ 15 Rg), and the right y-axis to the more extended one (Rout ≃ 140 Rg). The conversion between the two is a factor of ≃2.5. For each data point, we also show the corresponding XMM-Newton observation number (XMM1 to XMM11) as reference.

In the text

Fig. 13. Long-term evolution of the (quiescent level) bolometric luminosity. The y-axis scales have been defined in the caption of Fig. 12. In the upper panel we show the original light curve from all XMM-Newton and Swift observations. The lower panel shows a re-binned version of the same light curve (see text for details). In both panels, a vertical dashed line separates two different phases of the evolution: Phase1 comprises the initial X-ray peak and subsequent decay during the first ≃9 yr, while Phase2 is associated with the late-time X-ray re-brightening.

In the text

	Fig. 14. Lbol light curve during Phase1. The y-axis scales have been defined in the caption of Fig. 12. We show a simple model comprising a Gaussian rise and a power law decay with t−9/4.
In the text

Fig. 15. Lbol light curve during Phase2. The y-axis scales have been defined in the caption of Fig. 12. The first three data points are from the final decay of Phase1 and we also show the best-fitting decay model until Phase2 begins. A simple model comprising a Gaussian rise and a power law decay with t−9/4 for Phase2 is shown as a solid blue line, together with a Gaussian X-ray flare precursor (red dotted line).

In the text

	Fig. 16. Re-binned light curve during Phase 2 (black circles), together with its best-fitting model (see text for details). The green squares are the original Swift data points (before re-binning) around the precursor flare.
In the text

Fig. 17. XMM4 folded light curve (at the average recurrence time between QPEs). The light curve is normalised to the peak amplitude and shifted so that the peak is at phase 1. We show with different colours the phase intervals used to accumulate the spectra during the rise (R1 and R2), peak (P), and decay (D1 and D2). The quiescent spectrum (Q) is accumulated throughout the exposure (excluding time-intervals when QPE are present) and only a small fraction of the corresponding phase-interval is shown here. A phase interval Δp = 0.05 corresponds to a time interval Δt ≃ 1.6 ks.

In the text

Fig. 18. QPE spectral evolution during XMM4. In the upper panel we show Lbol as a function of blackbody rest-frame temperature. The colour-code is that defined in Fig. 17. The evolution starts at R1 (lowest Lbol at ≃83 eV, red data point) and proceeds counter-clockwise. The two dotted lines are Lbol ∝ T4 relations normalised to the extremal of the QPE evolution (R1 and D2). The lower panel shows the estimated radius of the blackbody emitting region as a function of QPE phase. The re-scaled QPE profile is also shown for reference (light grey).

In the text

Fig. 19. Viscous and front propagation timescales tvisc and tfront evaluated at Rsh as a function of orbital eccentricity. The upper (green) horizontal shaded area represents the period of the shock-induced perturbation (equal to that of the QPE or QPO, i.e. to the orbital period in the assumed scenario), while the lower (magenta) one is the shortest observed time delay of the QPO with respect to the immediately preceding QPE (XMM3).

In the text

	Fig. 20. Overall Lbol evolution over the ∼12 yr probed by our Swift and XMM-Newton campaign, together with the best-fitting models during Phase1 and Phase2, re-named here as TDE 1 and TDE 2 respectively. The y-axis scales are the same as those defined in the caption of Fig. 12.
In the text

Fig. 21. Portion of the Lbol light curve comprising the final part of TDE 1 and the overall TDE 2 evolution. Black circles (green squares) denote XMM-Newton (Swift) observations. Only XMM-Newton observations (whose observation number is also indicated) are long enough to detect QPEs. The shaded area highlight time-intervals during which high-amplitude, regular QPE are detected (XMM3 to XMM5, red) and those in which weak QPEs with irregular recurrence times are observed (XMM6, blue). No QPEs are detected after XMM6, although the higher level of the quiescent emission might have decreased significantly their contrast in subsequent observations.

In the text

Fig. B.1. Quiescent level variability. In the upper panels, we show the 0.2-1 keV light curves from all XMM-Newton observations with QPEs, together with their best-fitting models comprising a constant and a series of Gaussian functions to describe the quiescent level and QPEs respectively. The lower panels show the corresponding residual light curves re-binned by a factor of 8, together with two models, a constant (dashed lines) and a sine function (solid lines) with a period equal to the average observation-dependent recurrence time.

In the text