Quasi-periodic X-ray eruptions years after a nearby tidal disruption event

Observations and data analysis

X-ray data

Chandra

We downloaded processed Chandra images and event files and associated calibration data from the Chandra archive. We carried out analysis using CIAO version 4.16 (ref. ³⁶) and CALDB version 4.11.0. We checked for pileup using the pileup_map task, finding a pileup fraction of about 1% only for the central 4 pixels of the middle exposure. Therefore, pileup has negligible impact on our analysis. Count rates were extracted using the srcflux task. We used a 2-arcsec (4-pixel) circular radius and the default point-spread function (PSF) model. The background was estimated using an annular region with inner and outer radii of 15 and 60 arcsec, respectively, centred on AT2019qiz. This excludes other point sources, including the SE source (see below). The CIAO srcflux task includes the Bayesian Gregory–Loredo algorithm³⁷ to determine the optimal number of bins for investigating a time-varying (or, more formally, periodic) signal. The algorithm provides an odds ratio for variability (2.5 for AT2019qiz) and a light curve with the number of bins that maximizes this odds ratio. None of the other five sources in Extended Data Fig. 1 show an odds ratio >1.

We extract the spectrum in both eruption and quiescence (see below) using the specextract task. The spectrum of the eruption is soft and can be reasonably fit with a blackbody of about 100 eV. We perform a more detailed spectral analysis of AT2019qiz using the later eruptions and quiescent-phase data from instruments with greater sensitivity to softer (0.3–0.7 keV) X-rays (see sections ‘Swift/XRT and the quiescent spectrum of AT2019qiz’ and ‘NICER’).

The nature of the SE X-ray source

The Chandra images show a nearby source approximately 7 arcsec to the SE (labelled ‘SE source’ in Fig. 1). It overlaps with the PSF of AT2019qiz in all instruments other than Chandra. We extracted individual X-ray (0.5–7.0-keV) spectra from all three Chandra obsIDs to characterize the SE source. We perform spectral analysis with the Bayesian X-ray Analysis (BXA) software version 4.0.7 (ref. ³⁸), which connects the nested sampling algorithm UltraNest³⁹ with the fitting environment XSPEC version 12.13.0c (ref. ⁴⁰), in its Python version PyXspec. To improve the quality of the spectrum, we jointly fit all three Chandra obsIDs. The source can be fit with a simple power-law model with foreground absorption (tbabs × cflux(pow)) and is consistent with being constant over all three obsIDs. The neutral column density was fixed at the Milky Way value of 6.6 × 10²⁰ cm^−2. The 0.5–3.0-keV flux in the model is \({2.1}_{-0.9}^{+1.6}\times 1{0}^{-14}\,{\rm{erg}}\,{{\rm{s}}}^{-1}\,{{\rm{cm}}}^{-2}\) (90% posterior) and the photon index of the power law is Γ = 1.8 ± 0.5 (90% posterior). The fit is shown in Extended Data Fig. 4a.

Swift/XRT and the quiescent spectrum of AT2019qiz

We obtained Target of Opportunity time to follow up AT2019qiz with Swift/XRT. Eleven observations were obtained from 12 March 2024 to 14 March 2024, with a typical exposure time of about 1,200 s per visit and cadence of 4.5 h. We clearly detect one eruption in the new data (Fig. 1). We also reanalysed all previous XRT data for this source obtained under previous programmes, using the online tools available through the UK Swift Science Data Centre^41,42.

Owing to the better sensitivity at soft energies compared with Chandra, we are able to model the underlying disk spectrum using the XRT observations during the quiescent phase. For this, we use a colour-corrected thermal disk model (tdediscspec)⁴³, to be consistent with the full spectral energy distribution (SED) fit (see section ‘Disk modelling’). Given the larger PSF of XRT, we simultaneously model the AT2019qiz and the SE source contributions to the total spectrum. We use the model tbabs × (zashift(tdediscspec) + cflux(pow)), in which zashift(tdediscspec) is the contribution from AT2019qiz and cflux(pow) is the contribution from the SE source. The fit does not require a redshifted absorption component. We use PyXspec and BXA. For the disk parameters (that is, AT2019qiz), we assume flat priors; however, for the SE source, we use the posteriors from fitting its spatially resolved Chandra spectrum (see section ‘The nature of the SE X-ray source’) as the priors. Extended Data Fig. 4b shows their individual contributions to the observed spectrum, confirming that AT2019qiz dominates at energies below ≃ 1.0 keV. The posteriors of the fit indicate a peak disk temperature kT_p = 67 ± 10 eV (90% posterior), in agreement with the bulk TDE population⁴⁴.

Archival data from Swift/XRT

The X-ray spectrum of AT2019qiz observed by Swift/XRT in 2019–2020 was reported to be hard^16,21, suggesting a possible contribution from the SE source. To test this, we fit the combined spectrum (MJD 58714 to 59000) with the same power law plus disk model. We again use our power-law-fit posteriors for the SE source from Chandra as a prior in BXA and this time fix the temperature of the disk component while letting its flux vary freely. The early-time XRT spectrum is entirely consistent with the SE source, with no statistically significant contribution from the disk component (Extended Data Fig. 4c). This results in a 3σ upper limit on the flux (0.3–1.0 keV) from AT2019qiz at early times of ≤1.4 × 10^−14 erg s^−1 cm^−2, or a luminosity ≤7.2 × 10³⁹ erg s^−1.

By contrast, AT2019qiz is brighter and detected at high significance in data from 2022 onwards, with a spectrum dominated by the thermal component²¹. The luminosity of AT2019qiz measured during all quiescent phases with XRT and Chandra is roughly 10⁴¹â€‰erg s^−1, more than an order of magnitude fainter than the eruptions. Extended Data Fig. 6 shows the observation from 2022 in bins of 5 ks. The final bin shows an increase in flux, but the temporal baseline is too short to confirm or rule out that this represents the onset of a QPE (see also Fig. 4). The spectral fit from ref. ²¹ is consistent with a blackbody with kT_BB = 130 ± 10 eV, dominated by the final bin. We use the blackbody spectrum to calculate the luminosity in the final bin and exclude this bin from the disk model fit in Fig. 4a. We stack the remaining counts in a single bin and compute the quiescent luminosity using the fit from Extended Data Fig. 4.

NICER

NICER^45,46 observed AT2019qiz in two distinct campaigns, first at early times (around optical peak) from 25 September 2019 to 5 November 2019 and another at late times (about 1,600 days after optical peak) from 29 February 2024 to 9 March 2024.

The cleaned events lists were extracted using the standard NICER Data Analysis Software (HEASoft 6.33.2) tasks nicerl2 using the following filters: nicersaafilt=YES, saafilt=NO, trackfilt=YES, ang_dist=0.015, st_valid=YES, cor_range=“*-*”, min_fpm=38, underonly_range=0-80, overonly_range=“0.0-1.0”, overonly_expr=“1.52*COR SAX**(-0.633)”, elv=30 and br_earth=40. The whole dataset was acquired during orbit night-time and hence the daytime optical light leak (https://heasarc.gsfc.nasa.gov/docs/nicer/data_analysis/nicer_analysis_tips.html#lightleakincrease) does not apply to our data analysis. The latest NICER calibration release xti20240206 (6 February 2024) was used. Light curves in the 0.3–1.0-keV range were extracted using the nicerl3-lc task with a time bin size of 100 s and the SCORPEON background model.

The data obtained in the first campaign show no evidence for QPEs. Although the cadence is lower than that of the late-time data, it should be sufficient to detect QPEs occurring with the same frequency and duration as at late times, with a probability of detecting no QPEs of about 0.02 (using binomial statistics with a 20% duty cycle). We can therefore probably rule out QPEs within the first approximately two months after TDE fallback started (estimated to have occurred around 11 September 2019 (ref. ¹⁶)). However, we note that we would not expect QPEs during this phase in any model, as AT2019qiz was found to have an extended debris atmosphere¹⁶, which remained optically thick to X-rays until much later²¹.

During the second observing campaign, we clearly detect QPEs. The field of view of NICER is shown in Extended Data Fig. 1, overlaid on the Chandra image. All of the sources detected by Chandra have intensities (at energies lower than 1 keV) more than a factor of 10 below the measured peak of the QPE. Any contributions from these sources to the NICER spectra are further diminished by their offset angles from the centre of the field. We conclude that the NICER counts during eruptions are completely dominated by AT2019qiz. The six consecutive eruptions detected by NICER were modelled using a skewed Gaussian fit to each peak (Extended Data Fig. 2). We measure rest-frame delay times of 39.3 ± 0.3, 56.3 ± 0.3, 42.1 ± 0.3, 51.2 ± 0.2 and 53.5 ± 0.2 h between successive eruptions.

Given the high count rate and good coverage, we extracted time-resolved X-ray spectra from the second NICER eruption (Fig. 1) in the 0.3–0.9-keV band. We created Good Time Intervals (GTIs) with nimaketime for four intervals representing the rise, peak and decay (two phases) of the eruption. We extracted these spectra using the nicerl3-spec task and produced SCORPEON background spectra in ‘file mode’ (bkgmodeltype=scorpeon bkgformat=file) for each of the four GTIs. We simultaneously fit the four spectra using PyXspec and BXA, assuming the model tbabs × zashift(bbody). We fixed the redshift to z = 0.0151 and included foreground absorption, with a neutral hydrogen column density fixed to n_H = 6.6 × 10²⁰ cm^−2 (ref. ⁴⁷). We initially included a redshifted absorber, but the model preferred zero contribution from this component, so we excluded it for simplicity. The full posteriors of the parameters are shown in Extended Data Fig. 3.

AstroSat/SXT

We observed AT2019qiz with AstroSat⁴⁸ for four days starting on 12 March 2024 UT using the Soft X-ray Telescope (SXT)⁴⁹ and the Ultra-Violet Imaging Telescope (UVIT)^50,51. We used the level2 SXT data processed at the Payload Operation Centres using sxtpipeline v1.5. We merged the orbit-wise level2 data using SXTMerger.jl. We extracted the source in 200-s bins using a circular region of 12 arcmin. The broad PSF of the SXT does not leave any source-free regions for simultaneous background measurement. However, the background is low (0.025 ± 0.002 counts s^−1) and steady. As the quiescent flux measured by Chandra is below the SXT detection limit, we take this count rate as our background estimate and subtract it from the light curve. SXT detected one eruption (MJD 60383.548).

Optical/UV observations

HST

We observed AT2019qiz using HST on 21 December 2023 UT (MJD 60299.55), obtaining one orbit with the Wide-Field Camera 3 (WFC3) UVIS channel in the F225W band. We downloaded the reduced, drizzled and charge-transfer-corrected image from the HST archive. We clearly detect a UV source coincident with the nucleus of the host galaxy. We verify that this source is consistent with a point source both by comparing the profile with other point sources in the image using the RadialProfile task in photutils and by confirming that the fraction of counts within apertures of 3 and 10 pixels are consistent with published encircled energy fractions in the UVIS documentation.

We perform aperture photometry using a 10-pixel (0.396-arcsec) circular aperture, measuring the galaxy background per square arcsecond using a circular annulus between 20 and 40 pixels and subtracting this from the source photometry. Although we cannot measure the galaxy light at the precise position of AT2019qiz, having no UV images free from TDE light, the estimated background within our aperture is <2% of the transient flux, so our results are not sensitive to this approximation. We correct to an infinite aperture using the encircled energy fraction of 85.8% recommended for F225W. The zero point is derived from the image header, including a chip-dependent flux correction. We measure a final magnitude of 20.63 ± 0.03 (AB).

Although the angular scale of about 25 pc is not small enough to rule out a nuclear star cluster (NSC), the UV source is an order of magnitude brighter than known NSCs⁵². Moreover, NSCs are generally red⁵³ and many magnitudes fainter than their host galaxies in bluer bands. The magnitude of the source we detect is comparable with the total UV magnitude of the galaxy¹⁶. An unresolved nuclear source was also detected in the QPE source GSN 069 (ref. ⁵⁴).

Ground-based photometry

Numerous observations of this galaxy have been obtained by all-sky optical surveys both before and after the TDE. The optical emission was independently detected by ZTF^55,56, the Asteroid Terrestrial-impact Last Alert System (ATLAS)⁵⁷, Pan-STARRS⁵⁸ and the Gaia satellite⁵⁹.

Pan-STARRS reaches a typical limiting magnitude of about 22 in the broad w filter (effective wavelength of 6,286 Å) in each 45-s exposure. All observations are processed and photometrically calibrated with the PS image-processing pipeline^60,61,62. We downloaded and manually vetted all w-band observations of AT2019qiz since September 2019 and, in most cases, confirm a clean subtraction of the host galaxy light. We also retrieved ZTF forced photometry⁶³ in the r band (with a similar effective wavelength of 6,417 Å). Owing to the shallower limiting magnitude of about 20.5, we stack the fluxes in 7-day bins. Both surveys clearly detect a continuing plateau, persisting for >1,000 days with a luminosity νL_ν ≈ 7 × 10⁴⁰ erg s^−1. All Pan-STARRS and ZTF photometry was measured after subtraction of pre-TDE reference images using dedicated pipelines and hence include only light from AT2019qiz.

Although the optical light curves show scatter consistent with noise, they do not seem to exhibit the intense flaring behaviour seen in the X-rays. An order-of-magnitude flare in the optical would easily be detected even in the unbinned ZTF photometry. Assuming a duty cycle of 20%, and conservatively restricting to data since January 2022 (when we first see signs of day-timescale X-ray variability with XRT), the probability of never detecting an eruption simply because of gaps in cadence is ≲10^−13.

To test for optical variability on shorter timescales, we conducted targeted observations with the 1.8-m Pan-STARRS2 telescope in Hawaii on 11 February 2024, with the IO:O instrument on the 2.0-m Liverpool Telescope⁶⁴ in La Palma on 15 February 2024 and with ULTRACAM⁶⁵ on the 3.5-m New Technology Telescope at the European Southern Observatory (La Silla) in Chile on 10 February 2024. Pan-STARRS images were obtained in the w band (50 × 200-s exposures) and Liverpool Telescope in the r band (32 × 120 s), whereas ULTRACAM observed simultaneously in the u_s, g_s and r_s bands⁶⁶ (384 × 20 s, with only 24 ms between exposures). All images were reduced through standard facility pipelines. For Pan-STARRS, this included subtraction of a pre-TDE reference image and forced photometry at the position of AT2019qiz. In the case of Liverpool Telescope and ULTRACAM, we performed photometry using psf⁶⁷, an open-source Python wrapper for photutils and other image-analysis routines. We excluded 17 ULTRACAM images affected by poor seeing. We attempted manual subtraction of the Pan-STARRS reference images using psf; however, we found that the extra noise introduced by the subtraction was larger than any detectable variability. As shown in Extended Data Fig. 5, there is no strong evidence for variability on timescales on the order of hours.

Swift/UVOT

UV observations were taken with Swift/UVOT in the uvm2 filter contemporaneously with the XRT observations. We used the uvotsource package to measure the UV photometry, using an aperture of 12″. We subtracted the host galaxy contribution by fitting archival photometry data with stellar population synthesis models using Prospector⁶⁸. This standard procedure has been used to analyse previous UVOT observations of TDEs⁵⁶. We apply Galactic extinction correction to all bands using a E(B–V) value of 0.094 (ref. ⁶⁹).

The UVOT photometry is shown in Extended Data Fig. 5. Although lacking the resolution of HST to separate the central point source from the host light, the mean measured magnitude of 20.1 is about 0.5 mag brighter than the host level estimated by SED modelling¹⁶. The individual measurements exhibit root-mean-square variation of 0.27 mag (Extended Data Fig. 5), possibly indicating variability that would further exclude a nuclear star cluster. The timing of the XRT QPE is marked, coinciding with a possible (but not statistically significant) dip in UV flux as seen in the QPE candidate XMMSL1 J0249 (ref. ¹³).

AstroSat/UVIT

We observed AT2019qiz with UVIT using the broad filter CaF2 (F148W)⁵⁰. We processed the level1 data with the CCDLAB pipeline⁷⁰ and generated orbit-wise images, detecting a bright nuclear source. We performed aperture photometry using the UVITTools.jl package and the latest calibration⁵¹, in a circular region of 20 pixels (8.2 arcsec). We also extracted background counts from a source-free area of the image. The background-corrected count rate in the merged image corresponds to a flux density f_λ = 3.16 ± 0.97 × 10^−16 erg cm^−2 s^−1 Å^−1 or magnitude m = 20.49 ± 0.03 (AB). We found no statistically significant far-ultraviolet variability between the orbit-wise images. We do not attempt to remove host galaxy flux for the UVIT data, as the field has not been covered by previous far-ultraviolet surveys. SED modelling would require a large extrapolation. Regardless, we expect that the galaxy flux should be negligible at these wavelengths²⁰.

Analysis

Assessing variability

We perform two checks that the X-ray variability corresponds to QPEs rather than random variation. First we compare with physically motivated models of stochastic variability. Reference ⁷¹ demonstrated a mechanism to produce order-of-magnitude X-ray variability through Wien-tail amplification of accretion-disk perturbations. Their Fig. 3 shows the X-ray light curve of a model with a SMBH mass of 2 × 10⁶ M_⊙, consistent with AT2019qiz. The light curves are of a visibly different character to our data, with random variability rather than flares of consistent duration and no obvious ‘quiescent’ level. We ran further simulations using their model and never found a light-curve segment resembling AT2019qiz.

We also take a model-agnostic approach and assume the null hypothesis that the times of the X-ray peaks are random. Drawing a list of 10⁵ delay times from a flat probability distribution between 0 and 60 h and examining every consecutive sequence of eight, we ‘measure’ the standard deviation in delay times to be ≤15% of the mean in only ≲0.1% of trials. This is not sensitive to where we place the upper and lower bounds of the distribution. Therefore, we can exclude random peak times at >3σ confidence.

QPE duration–recurrence time correlation

The data in Fig. 3a show an apparent correlation between the mean duration and mean recurrence time of QPEs from a given source⁵. An equivalent statement is that QPEs seem to show a constant duty cycle across the population, with previous work indicating a duty cycle of 0.24 ± 0.13 (ref. ⁵). We reanalyse this correlation including AT2019qiz by performing Bayesian regression with a linear model T_duration = αT_recurrence + β. We find \(\alpha =0.2{2}_{-0.04}^{+0.11}\) (95% credible range), consistent with previous findings⁵. Comparing this model with the null hypothesis (α = 0), we find a change in the Bayesian Information Criterion ΔBIC ≈ 50, indicating a strong preference for a positive linear correlation over the null hypothesis of no correlation.

Disk modelling

We use the time-dependent relativistic thin disk model developed in refs. ^19,25. This computes the spectrum of an evolving accretion flow, produced at early times by the circularization of some fraction of the TDE stellar debris. To generate light curves, we follow the procedure of ref. ¹⁹ (their Fig. 2). The important input parameters are the mass and spin of the SMBH, the initial disk mass, the disk–observer inclination angle and the turbulent evolutionary timescale. Also, there are nuisance parameters relating to the initial surface density profile of the disk, which is generally unknown and has minimal effect on the late-time behaviour. As this initial condition is so poorly constrained, we simply consider an initial ring of material (as in ref. ²⁵).

For each set of parameters {Θ}, we compute the total (log-)likelihood

in which O_i,j, M_i,j and E_i,j are the observed flux, model flux and flux uncertainty of the jth data point in the ith band, respectively. For the X-ray data, we compute the integrated 0.3–1.0-keV flux using the best-fit models to the quiescent Swift/XRT and Chandra data, whereas for optical/UV bands, we compute the flux at the effective frequency of the band. We correct all data for foreground extinction/absorption^47,69.

The early optical and UV observations do not examine direct emission from the accretion flow, because of either reprocessing⁷² or shock emission from streams⁷³. We add an early-time component to model out this decay¹⁹, with functional form

in which B(ν, T) is the Planck function and ν₀ = 6 × 10¹⁴ Hz is a reference frequency. We fit the amplitude L₀, temperature T and decay timescale τ_dec, as well as the disk parameters. We only include data taken after the peak of the optical light curves.

The fit was performed using Markov chain Monte Carlo techniques, using the emcee formalism⁷⁴. To speed up computations, analytic solutions of the relativistic disk equations⁷⁵ were used. The model satisfactorily reproduces all data. The model X-ray light curve shows a slow rise; however, this is completely unconstrained by data and is therefore very sensitive to the uncertain initial conditions of the simulation. After a few hundred days (by the time of the earliest X-ray data in Fig. 4), the disk has spread to large radii and is no longer sensitive to initial conditions. We present the posterior distributions of the physically relevant free parameters in Extended Data Fig. 8. The best-fitting SMBH mass is consistent with all other observational constraints.

We note that a dimensionless SMBH spin parameter a_• > 0 is favoured by the model (although see caveats below), with a peak in the posterior around a_• ≈ 0.4–0.5. This constraint originates from the relative amplitudes of the optical/UV and X-ray luminosities, as highlighted in Extended Data Fig. 9. As the optical and UV light curves are well separated in frequency, the properties of the disk at scales r ≳ 20r_g are tightly constrained. The amplitude of the X-ray luminosity is controlled by the temperature of the inner disk, close to the innermost stable circular orbit. For a given large-scale structure, this radius is determined by a_•.

Our disk model parameterizes the colour correction factor f_col in terms of the local disk temperature⁷⁶, but our posteriors do not marginalize over its unknown uncertainty. Recognizing that modest uncertainties in f_col lead to substantial uncertainties in spin (for non-maximal black hole spins)⁷⁷, we do not claim a spin measurement here but simply note that a modest spin is consistent with our data. The spin estimates in this model also assume a planar disk that is aligned with the SMBH spin, which is not true in the case of a precessing disk (see next section).

Although the disk temperature profile (and therefore the location of the disk’s outer edge) is tightly constrained from the multiband late-time observations, it is well known that disk temperature constraints only scrutinize the product \({W}_{\phi }^{r}\varSigma \), in which \({W}_{\phi }^{r}\) is the turbulent stress and Σ is the surface mass density. As the functional form of the turbulent stress cannot be derived from first principles, and must be specified by hand, there is some uncertainty in the mid-disk density slope. Our choice of \({W}_{\phi }^{r}\) parameterization is optimized for computational speed⁷⁵ and is given by \({W}_{\phi }^{r}=w={\rm{constant}}\). Rather than fit for w, we fit for the evolutionary timescale of the disk (which has a more obvious physical interpretation), given by \({t}_{{\rm{evol}}}\equiv 2\sqrt{G{M}_{\bullet }{r}_{0}^{3}}/9w\). We emphasize that this uncertainty has no effect on our constraints on the size of the disk.

With this choice of parameterization for the turbulent stress, the disk density profile (Fig. 4) can be approximated as Σ ∝ r^−ζ, with ζ = 1/2, for r = (2–600)GM_•/c². The density slope is not very sensitive to modelling assumptions, with the (potentially) more physical radiation-pressure-dominated α-disk model having ζ = 3/4.

Precession timescales

If the SMBH is rotating, any orbit or disk that is misaligned with the spin axis will undergo Lense–Thirring precession. This is a possible cause of timing variations in QPEs³⁰. Changes in QPE timing in AT2019qiz are seen over the course of ≲8 observed cycles, which would require that the precession timescale T_prec is approximately several T_QPE, in which T_QPE ≈ 48.4 h is the QPE recurrence time.

The precession timescale can be calculated following²⁹:

in which r_in and r_out are the inner and outer radii of the disk or orbit, respectively, in Schwarzschild units (see also ref. ⁷⁸). We assume log(M_•/M_⊙) = 6.3 and investigate the plausible precession period for different values of a_•.

The nodal precession timescale for an orbiting body can be estimated by calculating T_prec at the orbital radius (setting R_in ≈ R_out ≈ R_orb). For a_• = 0.1–0.9, this gives T_prec,orbit ≈ (10³â€“10⁴) × T_QPE, independent of ζ. Therefore, in the EMRI model, nodal precession is too slow to account for changes in QPE timing over several orbits.

The precession timescale of the disk can be calculated by assuming that it behaves as a rigid body with r_in = 2GM_•/c², r_out = 600GM_•/c² and a density slope ζ = 1/2 from our disk model. We use the above equation to find T_prec,disk ≈ (70–200) × T_QPE (for the same range of spins). With a steeper density profile having ζ = 1, this would reduce to T_prec,disk ≈ (8–70) × T_QPE (because more mass closer to the SMBH enables stronger precession). Therefore, precession can explain detectable changes in QPE timing over the course of several orbits only in the case of a rapidly spinning SMBH (a_• ≳ 0.5–0.7) and a steep disk density profile.

With these constraints, attributing the timing residuals primarily to disk precession becomes challenging. The larger the SMBH spin magnitude, the faster an initially inclined disk will come into alignment with the black hole spin axis, damping precession on a timescale ≲100 days for a_• > 0.6 and M_• ≈ 10⁶ M_⊙ (ref. ⁷⁹). To maintain precession for more than 1,000 days requires a spin a_• ≲ 0.2, in which case the precession is not fast enough to fully explain the timing variations in our data.

We also note that the disk inner radius used in our precession calculation was derived from a planar disk model. In a tilted disk around a spinning SMBH, the radius of the innermost stable circular orbit will differ from the equatorial case. Understanding the effect of disk precession in AT2019qiz will probably require both continued monitoring to better understand the QPE timing structure and a self-consistent model of an evolving and precessing disk that can explain both the multiwavelength light curve and the timing residuals.

Constraints on QPE models

Many models have been proposed to explain QPEs. Disk tearing owing to Lense–Thirring precession has been suggested⁸⁰. This effect has plausibly been detected in the TDE AT2020ocn (ref. ⁸¹). However, its X-ray light curve did not resemble that of AT2019qiz or those of other QPEs. As discussed above, it is also unclear whether strong precession will persist until such late times. The X-ray variability in AT2020ocn occurred only in the first months following the TDE.

Gravitational lensing of an accretion disk by a second SMBH in a tight binary could cause periodic X-ray peaks for the right inclination⁸². However, in the case of AT2019qiz, no signs of gravitational self-lensing were detected during the initial TDE. In this model, a QPE magnification by a factor ≳10 requires an extremely edge-on view of the disk, which leads to a shorter duration of the QPE flares. This was already problematic for previous QPEs⁸² and is more so for the longer-duration flares in AT2019qiz. Moreover, finding a TDE around a close SMBH binary within a very narrow range of viewing angles (≳89.5°) is very unlikely in the small sample of known TDEs, so a strong TDE–QPE connection is not expected in this model.

Limit-cycle instabilities are an appealing way to explain recurrent variability^7,83. The recurrence timescale for disk-pressure instabilities depends on whether the disk is dominated by radiation pressure or magnetic fields⁸, as well as the accretion rate. Our disk model, which is well constrained by the multiwavelength data, gives an Eddington ratio \(\dot{M}/{M}_{{\rm{Edd}}}\approx L/{L}_{{\rm{Edd}}}=0.04\pm 0.01\). Reference ²⁶ gives formulae to interpolate the recurrence time for radiation-pressure instabilities, for a given amplitude relative to quiescence. We assume a peak-to-quiescence luminosity ratio of 60, although our analysis is not sensitive to this. Using the prescription for either the intermediate-mass black holes (their equation 33) or SMBHs (their equation 34), we find a recurrence time of about 5,000 days.

In the magnetic case, we use equation 14 from ref. ⁸. Matching the observed recurrence time requires a dimensionless magnetic-pressure scaling parameter p₀ ≈ 10. However, at this Eddington ratio, the disk should be stable⁸ if p₀ ≳ 1. This leaves no self-consistent solution in which magnetic-pressure instabilities cause the QPEs in AT2019qiz. The possibility of a long–short cycle in recurrence time, and the asymmetric profile of the eruptions³, also disfavour pressure instabilities. We also note that, in disk-instability models, the recurrence time of the instability correlates with SMBH mass. For the known QPEs, there is no apparent correlation in recurrence time with mass (Fig. 3).

The final class of models to explain QPEs involves an orbiting body (EMRI) either transferring mass to an accretion disk or colliding with it repeatedly^{9,10,11,27,28,30,35,84,85,86}. Note that this is very unlikely to be the same star that was disrupted during the TDE: if a bound remnant survived the disruption, it is expected to be on a highly eccentric orbit with a much longer period than the QPEs¹¹. The fundamental requirement for star–disk collisions to explain QPEs is that the disk is wider than the orbit of the EMRI. The size of the disk in AT2019qiz is well constrained by our analysis and the posteriors of our fit fully satisfy this requirement, at least in the case of a circular disk.

For an orbit with the QPE period to avoid intersecting the disk would require a disk ellipticity e > 0.7 (assuming that the semimajor axis of the disk is fixed) and an appropriately chosen orbital inclination. Although some TDE spectra support a highly elliptical disk in the tens of days after disruption⁸⁷, most can be explained with an approximately circular disk^88,89,90. Simulations of TDE accretion disks show a high ellipticity in the days after disruption⁹¹, but shocks are expected to circularize the disk over the course of a few debris orbital periods⁹² (days to weeks), whereas we observe QPEs on a timescale of years after AT2019qiz. An initially highly eccentric disk becomes only mildly elliptical (e ≈ 0.6) on a timescale of several days (refs. ^{93,94,95,96,97,98,99,100}). Once notable fallback has stopped (before the plateau phase), no more eccentricity will be excited in the disk, whereas turbulence will act to further circularize it, so we expect that the disk in AT2019qiz will be circular to a good approximation.

The case of an EMRI interacting with a TDE disk was specifically predicted by refs. ^11,30. The formation rate of EMRIs by the Hills mechanism is about 10^−5 year^−1 galaxy^−1, about one-tenth of the TDE rate. Because the time for inspiral resulting from gravitational-wave emission (approximately 10⁶ years) is longer than the time between TDEs (approximately 10⁴ years), theory predicts that ≳1 in 10 TDEs could host an EMRI capable of producing QPEs^11,35. This is consistent with recent observational constraints on the QPE rate²⁴.