Recurring tidal disruption events a decade apart in IRAS F01004-2237

The ATLAS has been observing F01004-2237 in the $c$ and $o$ bands since August 2015. We obtained the photometries on the subtracted images from the forced photometry server¹¹1https://fallingstar-data.com/forcedphot/ (Shingles et al. 2021). We only selected data taken under good weather conditions with the criterion that the limit magnitude is fainter than 18 for both bands. We then binned the LC data nightly²²2The survey strategy of ATLAS results in 2–8 consecutive exposures per night of observation.. The flux in the binned LC was calculated using the weighted average value of the single-exposure fluxes in the time bin, where the weight is the reciprocal of the square of the flux errors, and its error was calculated according to the law of propagation of uncertainties. We finally obtained reference fluxes from the forced photometry server, 249 and 329 $\mu$Jy for $c$ and $o$ bands, respectively, and added them to the photometries.

The ASASSN has been observing F01004-2237 in the $g$ band since September 2017³³3Although ASASSN observed F01004-2237 in the $V$ band between 2013 and 2018, we did not use the data because the V-band observation covered neither of the two flares and the data quality is poor.. We obtained the original light curve data from the ASASSN Sky Patrol⁴⁴4https://asas-sn.osu.edu/ (Kochanek et al. 2017). We used the photometries on the subtracted images with reference flux added and selected good data with a criterion that the limit magnitude is fainter than 18. We binned the LC data nightly as we did for the ATLAS data.

We collected the $G$-band LC of F01004-2237 from the Gaia Photometric Science Alerts website⁵⁵5http://gsaweb.ast.cam.ac.uk/alerts/alert/Gaia22ahd/. We did not analyze this LC data in detail and only used it for display purposes because of the sparse sampling and the lack of error data. To study the 2010 flare, we also collected the CRTS $V$-band LC of F01004-2237 from July 2005 to January 2014 from the data release 2⁶⁶6http://nesssi.cacr.caltech.edu/DataRelease/. CRTS observed the source 4–10 nights per year by taking 2–4 exposures per night, and we also binned the LC data nightly.

We made follow-up observations in $g$, $r$ and $i$ bands from December 2021 to January 2023 using the Las Cumbres Observatory Global Telescope Network (LCOGT, Brown et al. 2013). Using the images taken before the flare by the SkyMapper Southern Survey (Wolf et al. 2018) as the reference images, we made image subtraction with the software Saccadic Fast Fourier Transform (SFFT, Hu et al. 2022). We then performed PSF photometries on the difference images using the Photutils package of Astropy. We finally made flux calibration using nearby stars with the source catalog⁷⁷7https://catalogs.mast.stsci.edu/panstarrs/ of the Panoramic Survey Telescope & Rapid Response System (PanSTARRS, Kaiser et al. 2002). We show all the original light curve data in Fig. 1.

After we noticed the 2021 flare, we triggered target of opportunity observations of Swift telescope. We present the information of all the observations in Table A1 in the appendix.

We obtained the XRT count rate in 0.3–10 keV of F01004-2237 from the UK Swift Science Data Centre⁸⁸8https://www.swift.ac.uk/user_objects/ (UKSSDC). Their server calculated the count rate following (Evans et al. 2007, 2009) as:

where $N_{\rm tot}$ is the total photon count in the source region, $N_{\rm bkg}$ is the expected count of background photons in the source region, $t_{\rm exp}$ is the exposure time, and $f_{c}$ is the correction factor for pile-up, dead zones on the CCD, and source photons falling outside of the source extraction region. We list the count rates with 1$\sigma$ errors (or 3$\sigma$ upper limits) and data for calculating them in Table A1.

We also extracted the photometry of F01004-2237 on the UVOT images using uvotsource in High Energy Astrophysics Software (HEASoft). We used an aperture with a radius of 5$\arcsec$, and calculated background in an annular region with inner and outer radii of 35$\arcsec$ and 40$\arcsec$, respectively. The results are listed in Table A2.

We triggered a target of opportunity observation of XMM-Newton telescope (ID 0911791101) on June 8, 2022. In addition, there is an archival XMM-Newton observation (ID 0605880101) on December 8, 2009 (see details in Nardini & Risaliti 2011). This observation was taken before the 2010 flare occurred (later than Jan 2010; see section 3.2) and can be used to measure pre-flare fluxes from the host galaxy. So we analyze both the two data. We process the XMM-Newton data using the Science Analysis System (SAS) version 20.0.0. We included event patterns 0–4 for the PN camara and 0–12 for the MOS cameras. The time intervals of high-flaring background contamination were identified and excluded by examining the light curves in the 10–12 keV energy range. The total cleaned exposure times for the PN and MOS cameras are 28 and 33 ks, respectively, in 2009, and are 12 and 20 ks in 2022. We extracted the source+background spectra from a circular region centred on the target with a 32$\arcsec$ radius, and the background spectra from an adjacent source-free circular region with a larger radius. We combined the spectra taken from PN and two MOS cameras using the task epicspeccombine, during which the background spectra and the response matrices were also combined. The combined EPIC spectra were then grouped so that there are at least 20 net counts in each energy bin to ensure a $\chi^{2}$ statistic. We also extracted the photometry of F01004-2237 on the OM images using the task omichain, and the results are listed in Table A2.

We made follow-up optical spectroscopic observations of F01004-2237 using Lijiang 2.4-meter telescope (LJT) on December 20, 2021 and January 3, 2022, Southern Africa Large Telescope (SALT) on December 30, 2021, Magellan telescope on January 2, 2022, Gemini South telescope on July 3, 2022 (Proposal ID GS-2022A-DD-104), and Palomar 200-inch Hale telescope (P200) on September 3, 2022.

For the SALT observation, we used the Robert Stobie Spectrograph (RSS, Burgh et al. 2003) with the PG0900 grating, resulting in a mean resolving power of $R\sim 800-1200$. The exposure time is 1000 seconds. To reduce the spectra, we used the PyRAF-based PySALT package⁹⁹9https://astronomers.salt.ac.za/software/ (Crawford et al. 2010), which includes corrections for gain and cross-talk and bias subtraction. We extracted the science spectrum using standard IRAF tasks, including wavelength calibration, background subtraction, and 1D spectra extraction. We obtained flux calibration using observations of standard stars during twilight. For the Magellan observation, we used the Low Dispersion Survey Spectrograph (LDSS-3). The total exposure time is 900 seconds. We reduced the spectra using IRAF following the standard routine. For the Gemini South observation, we used the GEMINI Multiple Object Spectrographs (GMOS, Davies et al. 1997) with an R400 grating, centring at 655 nm, resulting in a continuous wavelength coverage of 4860–9200 Å. Three exposures with 600 seconds exposure each were taken. Following the standard routine, we reduced the data using Pyraf, during which the flux calibration was made using a standard star spectrum taken on the same night. For the P200 observation, we used the Double Spectrograph (DBSP) with a D55 dichroic, a 600/4000 grating for the blue side, and a 316/7500 grating for the red side, resulting in a continuous wavelength coverage of 3200–10200 Å. The total exposure time is 700 seconds. We reduced the data using Pyraf following the standard routine.

To model the starlight and identify narrow emission lines (NELs), which are crucial for analyzing the flare’s spectra, we obtained the pre-flare VLT/XShooter spectrum observed in August 2018 (see details in Tadhunter et al. 2021). We downloaded the two-dimensional spectra from the two exposures from ESO database¹⁰¹⁰10http://archive.eso.org/wdb/wdb/adp/phase3_spectral/form , which had been processed by the XShooter pipeline, and then combined them to remove cosmic rays, and finally extracted one-dimensional spectra with an aperture of 1.8”. We recalibrated the fluxes of all the spectra using quasi-simultaneous photometric data. All the spectra used in this work are shown in Fig. 2.

We triggered a follow-up UV spectroscopic observation in July 2022 with Space Telescope Imaging Spectrograph (STIS) mounted on Hubble Space Telescope (HST, proposal ID 16943). The observations using G140L and G230L gratings were made on July 30 and July 14 with exposure times of 3766 and 5947 seconds, respectively. A slit with an aperture of 52$\times$0.2 arcsec² was used and positioned through the galaxy centre. We downloaded the data from the Mikulski Archive for Space Telescopes¹¹¹¹11https://archive.stsci.edu (MAST). We obtained the two-dimensional spectra that the STIS pipeline had processed. We then combined the exposures to remove cosmic rays and extracted one-dimensional spectra with an aperture of 0.2$\arcsec$. As a result, the spectrum is from a region in the galaxy centre with an area of 0.2$\times$0.2 arcsec², taking into account the slit width.

To obtain the host galaxy’s contribution to UV emission, we collected an archival HST/STIS spectrum (Proposal ID 8190) observed on February 9, 2000 (see details in Farrah et al. 2005). Observations with G140L, G430L and G750L gratings were taken. A 52$\times$0.2 slit was also used and positioned through the galaxy centre. With a similar procedure, we obtained a pre-flare spectrum, which is also from a 0.2$\times$0.2 arcsec² region in the galaxy centre. We made subtraction with the two G140L spectra to obtain a difference spectrum covering a wavelength range of 1165 to 1710 Å. All the STIS spectra used in this work and the difference spectrum are shown in Fig. 3. Because the aperture sizes of the two spectra are the same and the flux calibration of STIS is stable, we adopted the difference spectrum as emission from the 2021 flare.

We checked the position of the 2021 flare using Gaia data with the best positional accuracy. According to the Gaia alert, when the 2021 flare was detected, the coordinate of the target was RA = 01h 02m 49.990s and Dec = -22d 21m 57.24s (J2000). The coordinate before the flare occurred was RA = 01h 02m 49.991s and Dec = -22d 21m 57.27s, which was obtained from the source catalogue of Gaia data release 2 (DR2, Gaia Collaboration et al. 2018). The offset between the two positions is 33 mas. Hodgkin et al. (2021) found that the offset between Gaia alerts and their Gaia DR2 counterparts obeys a Rayleigh function with an average value of 55 mas. Thus, we found no evidence that the flare’s position deviates from the centre of the galaxy. Using the Rayleigh function above, we can limit the offset of the two positions to be within 160 mas (99.73% probability).

The constraint to the position of the Gaia alert is not equivalent to the constraint to the flare’s position. This is because the Gaia position is the average position based on the brightness distribution. When an off-centre flare occurs, the average position should be between the flare’s position and the centre of the galaxy. The alerting magnitude of Gaia22ahd is 0.8 magnitudes brighter than historical values, which means that the flare is 1.1 times brighter than the host galaxy. A simulation shows that the offset of a hypothetical off-centre flare should be 1.9 times the offset of the average position. Thus, we can limit the flare’s offset relative to the galaxy centre to be within 300 mas.

We also measured the flare’s position using the LCOGT $g$-band data. We measured the positions on the reference-subtracted images with high signal-to-noise ratios. We then calculated an averaged position of RA = 01h 02m 49.987s and Dec = -22d 21m 57.21s, which has an 82 mas offset from the Gaia position. We estimated a positional error of $\sim$80 mas from the scatter of the positions. Thus, the analysis of LCOGT images yields no evidence of offset between the flare’s position and the galaxy centre and limits the offset to be within 240 mas.

The results using Gaia and LCOGT data are consistent. Taking them together, we can limit the offset to be within $\sim$0.2$\arcsec$, corresponding to a physical size of 440 pc.

In Fig. 1, we display the long-term variability of F01004-2237 with CRTS/$V$, ATLAS/$o$, ASASSN/$g$ and Gaia/$G$ LCs, which clearly show the two flares. The CRTS/$V$ LC shows that the galaxy’s emission remained quiescent until January 23, 2010 (MJD$=$55219.5) with $V=17.54$, and rose by 0.5 magnitude to $V=17.04$ on September 17, 2010 (MJD$=$55456.4), and then slowly declined over the next three months, and then dropped to $V\sim 17.3$ in 2011, and finally fading gradually. With this LC, we inferred that the flare peaked between January and September 2010 (MJD$=55338\pm 118$). We fit the LC with a $t^{-5/3}$ curve expressed as:

where $t_{\rm peak}$ is the peak time and is set to be in the range of MJD$=55338\pm 118$, and $f_{\rm peak}$ is the peak flux. As can be seen in Fig. 4, the curve well fit the data. Using this curve, we estimated a peak $V$-band magnitude of $17.4\pm 0.5$. The Gaia/$G$ LC shows a long-term fading trend from 2014 to 2018, so the 2010 flare might have lasted seven years. The galaxy’s emission went quiescent again from 2019 to August 2021, during which no variation was detected in any band. Using data in this period, we derived the quiescent levels before the 2021 flare that $o=17.59$, $c=17.87$, and $g=17.81$.

In Fig.5, we present the LCs of the 2021 flare after subtracting quiescent fluxes. For the six Swift/UVOT bands, we adopted the fluxes observed in May 2024 as the quiescent fluxes, because they are close to those observed by XMM-Newton/OM in 2009 (the difference could be explained as statistical fluctuation, difference in filters’ response, and variation of possible weak underlying AGN). The first $>3\sigma$ detection was on September 4, 2021 (MJD=59461.6, the discovery date hereafter) in the ATLAS/$c$ band with $c=20.22\pm 0.16$. Then, the flare rapidly rose, reaching a peak about 50 days after discovery. We fit the $o$-band LC near the peak with polynomial and obtained a $t_{\rm peak}$ of MJD$=59519.8\pm 0.7$, and a rising time $t_{\rm rise}$ of $52$ days in the rest frame. At $t_{\rm peak}$, the flare has $o=17.64\pm 0.02$, $c=17.42\pm 0.02$, and $g=17.36\pm 0.07$. Here, the $c$-band and $g$-band magnitudes were calculated by shifting the best-fitting $o$-band LC to match the data points in these bands. These apparent magnitudes correspond to $M_{g}=-21.16$ and $M_{V}=-21.05$ after K correction. The flare declined slowly after the peak in all bands, although some fluctuations could be seen. After $+$500 days, the flare was only barely detected in the $c$ band and was swamped by host galaxy emission at other bands. After $+$650 days, the flare dropped to an undetectable level with $c>21$.

The LC that rises rapidly in 50–60 days, peaks with absolute magnitudes of $\sim-21$, and fades in years is reminiscent of superluminous supernovae (SLSNe) and TDEs. To distinguish between the two possible origins, we measured the variation of the photospherical temperature. We fit the SED obtained from the 5 Swift/UVOT observations with all six filters used with blackbody curves. We showed the $T_{\rm BB}$ and $L_{\rm BB}$ in Fig. 6. The temperature keeps nearly constant as the luminosity declines and is still $>20,000$ K at a phase of $+220$ day. Such a late-time high temperature is inconsistent with SNe and is consistent with TDEs (e.g., van Velzen et al. 2020; Zabludoff et al. 2021). In addition, as can be seen in Fig. 5, the LCs in the dropping phase ($>$60 days) are roughly in line with $t^{-5/3}$ curves, also consistent with TDEs.

To further examine the similarity between the flare and TDEs and measure the peak luminosity and the energy budget, we simultaneously fit the multi-band LC data with the TDE model of van Velzen et al. (2021). The model assumes blackbody spectra with luminosity and temperature vary with time as:

where $t_{\rm peak}$ and $L_{\rm peak}$ are the peak time and peak bolometric luminosity, $\sigma$ describes the how the pre-peak luminosity rises, and $t_{0}$ and $p$ describes the how the post-peak luminosity declines. And the temperature $T(t)$ is assumed as:

We fit the data using the Markov Chain Monte Carlo (MCMC) code of emcee, and set the priors following van Velzen et al. (2021). The MCMC realizations of bolometric LCs are shown in Fig. 7, and the parameters are listed in Table 2. The peak time obtained using LCs in all bands agrees with the $o$-band peak time. The peak bolometric luminosity is $4.4\pm 0.4\times 10^{44}$ erg s^-1, and the energy released up to $+$700 day reaches $5.2\pm 0.4\times 10^{51}$ erg. We compared the parameters with those of the 17 ZTF TDEs of van Velzen et al. (2021), and found that most of the parameters of the 2021 flare fall within the range of typical values of TDEs, and the peak luminosity and the rise time are close to the maximum value of TDEs.

The UV LCs show a bump in $+220$ to $+260$ days on the basis of a power-law decline (Figure 5 and Figure 7). Such late-time bumps, appearing as re-brightening or flattening on the LCs, are also seen in some TDEs and candidates, such as ASASSN-15lh (Leloudas et al. 2016; Brown et al. 2016) and some recent TDEs (Yao et al. 2023). Possible explanations include a transition from super-Eddington to sub-Eddington fallbacks (Strubbe & Quataert 2009), tidal disruptions of double stars (Mandel & Levin 2015), and delayed formation of accretion disk (Liu et al. 2022). Unfortunately, there are not enough observational data to test these possible explanations.

There are four spectra taken between $+$95 and $+$109 days (Fig. 2). We first investigated the spectrum taken by Magellan on $+$108 day, which has the best data quality and the broadest wavelength coverage. We modelled the spectrum with the sum of a stellar model and a third-order polynomial, the latter representing the flare’s continuum emission, by masking the wavelength ranges affected by narrow emission lines from the host galaxy and telluric absorption bands. We describe the construction of the stellar model and the identification of narrow emission lines in detail in Appendix B. We fit using MCMC code emcee assuming flat priors for all the parameters. Although this model reproduces most parts of the spectrum, it leaves significant residuals in three wavelength ranges around 6550 Å, 4830 Å, and 5890 Å, implying the presence of emission features (Fig. 8(a)). We identified the 6550 Å feature as H$\alpha$ and fitted it by adding a single Gaussian component to the model. It is very broad with a Full Width at Half Maximum (FWHM) of $19200\pm 400$ km s^-1 and is slightly blueshifted with a velocity of $550\pm 150$ km s^-1. We fit the 4830 Å feature by adding another single Gaussian and found that the central wavelength is $4830\pm 8$ Å, the FWHM is $22400\pm 1200$ km s^-1, and the flux is $0.62\pm 0.03$ times that of H$\alpha$. If identified as H$\beta$, the width and shifted velocity differ a lot from those of H$\alpha$, and the ratio of H$\alpha$/H$\beta$ of $1.62\pm 0.08$ is less than the Case B value of 2.70¹²¹²12Here we adopted the value in $T_{e}=20,000$ K and $n_{e}=10^{10}$ cm^-3 (Storey & Hummer 1995).. These are not reasonable; hence, we considered that the 4830 Å feature is a mixing of H$\beta$ and He II $\lambda$4686 emission lines. We deblended them by adding two Gaussians for H$\beta$ and He II each and requiring that the FWHM and shifted velocity of H$\beta$ to be within $\pm$3,000 km s^-1 of those of H$\alpha$, following Charalampopoulos et al. (2022). We also added another single Gaussian for the 5890 Å feature, identified as He I $\lambda$5876. The model fits the data well, as seen in Fig. 8(a), and the parameters are listed in Table C1. The inferred H$\alpha$/H$\beta$ ratio of $3.5^{+2.4}_{-0.6}$ is consistent with the Case B value. The He II FWHM is $25100\pm 3800$ km s^-1 and the He II/H$\alpha$ ratio is $0.35\pm 0.08$.

To explore how much the assumption of the flare’s continuum model would affect the analyses of BELs, we tried two other models, double blackbody and double power-law (they fit much better than either single blackbody or single power-law). Through analyses similar to those described in the previous paragraph, we found that all four lines are significantly detected regardless of which model is used. The best-fitting He II FWHM is $\sim$60,000 km s^-1. We considered this value unreasonable because such broad He II emission lines have never been seen in TDEs or AGNs. Therefore, we set the He II FWHM to be less than 25,000 km s^-1 (the value assuming polynomial continuum) in the MCMC. The resultant BEL parameters are present in Table C1. Combining the results of the three models, the H$\alpha$ FWHM is 16,000–19,000 km s^-1, the He II FWHM is $>21,000$ km s^-1, and He II/H$\alpha$ ratio is 0.3–0.7.

The three spectra observed between $+$95 and $+$109 days by SALT and LJT resemble the Magellan spectrum (Fig. 2), and VBELs can be seen visually. We attempted to analyze them in the same way as we did for the Magellan spectrum. For each of the three spectra, adding VBEL components significantly improved the goodness of fit. This indicates that the VBEL persisted over this period of time. However, the parameters of the VBELs have large uncertainties when different flare’s continuum models are considered. The reason may be the short spectral coverage, especially the absence of the red wing of H$\alpha$. We will not use the results of these three spectra for further discussion.

The spectra observed by Gemini-S and P200 (Fig. 8(b) and  8(c)) on $+$271 and $+$325 days still show BELs though the continuum flux decreases significantly. We analyzed the two spectra similarly, except that we did not add a He I $\lambda$5876 component, which is not visible on the spectrum, and list the BEL parameters in Table C1. H$\alpha$ and He II are still significantly detected, while the H$\beta$ is too weak to be detected. In the Gemini-S spectrum, the width of H$\alpha$ BEL is 18,000–19,000 km s^-1, similar with that in the Magellan spectrum. While in the P200 spectrum, the H$\alpha$ BEL FWHM is 7,000–11,000 km s^-1, much narrower. The inferred He II/H$\alpha$ ratio is in the ranges of 0.3–0.6 and 0.9–2.3 for the two spectra, respectively.

Our analyses show that newly formed very broad emission lines (VBEL) of Balmer and He II $\lambda$4686 appear between $+$95 and $+$325 days after the 2021 flare. The FWHMs of the VBELs exceed 7,000 km s^-1 consistently. The width of H$\alpha$ VBEL can reach 16,000–19,000 km s^-1 at maximum, while that of He II can even reach $\gtrsim$21,000 km s^-1.

The VBELs in F01004-2237 seen 43–273 days after the optical maximum disagree with VBELs seen in SLSNe, which occur only near the optical peak (Könyves-Tóth & Vinkó 2021). Combined with the analysis of late-time photospheric temperature in section 3.2, we ruled out the SLSNe interpretation for the 2021 flare. In addition, the VBELs seen in F01004-2237 can be explained by TDE because such broad lines have been seen in other TDEs (Charalampopoulos et al. 2022).

We had obtained a UV spectrum in the rest-frame wavelength range of 1042–1530 Å of the 2021 flare on July 30, 2022 ($+$294 day) using the HST/STIS spectra with G140L grating taken in 2022 and 2000 (section 2.4). The spectrum shows a broad emission feature around $\sim$1230 Å on the continuum (Fig. 9). It also shows some narrower emission lines, located at 1216, 1241 and 1304 Å, which were identified as Ly$\alpha$, N V $\lambda$1240 and O I $\lambda$1304, respectively.

We fit the spectrum with a model including a blackbody for the continuum, a Gaussian for the broad emission feature, and a Gaussian for each narrower emission line. We masked the wavelength ranges affected by the narrow absorption lines (see details in Appendix D) and the Milky Way’s emission lines. The best-fitting model is shown in Fig. 9. The broad emission feature has a central wavelength of $1231.5\pm 0.8$ Å, a FWHM of $20700\pm 500$ km s^-1, and a luminosity of $4.5\pm 0.1\times 10^{42}$ erg s^-1. If identified as Ly$\alpha$, it is redshifted by $3900\pm 200$ km s^-1, and the FWHM is close to that of H$\alpha$ (18,000–19,000 km s^-1) in the GMOS spectrum with a time interval of 23 days. The difference in the shifted velocity between Ly$\alpha$ and H$\alpha$ (around zero velocity) may be caused by intrinsic variation of the emission line region over the 23-day interval, the difference in conditions of blueshifted and redshifted emitting gas, or mixing of a possible N V $\lambda$1240 VBEL. The flux ratio of Ly$\alpha$/H$\alpha$ is $10.4\pm 0.4$, consistent with the Case B value of 13.0 considering the difference in UV and optical spectra apertures and a possible dust reddening. Our analysis of the UV spectrum again confirms the presence of VBELs and suggests that UV and optical VBELs may have the same origin.

The narrower emission lines of Ly$\alpha$, N V $\lambda$1240 and O I $\lambda$1304 have a FWHM of $1500\pm 80$ km s^-1, and is unshifted relative to the systematic redshift (velocity $20\pm 30$ km s^-1). This emission line system is also seen in C IV $\lambda$1550, He II $\lambda$1640, N III] $\lambda$1750, Si III] $\lambda$1892, C III] $\lambda$1907 and Mg II $\lambda$2800 in the G230L spectrum (Fig. 3). It may be the light echo from the nuclear ISM of the 2010 flare, the 2021 flare, or both of them. A detailed analysis is beyond the scope of this paper and will be described in detail in the following works.

Prior to its 2021 flare, F01004-2237 was observed and detected by VLA at 1.4, 5.5, 9.0, and 14.0 GHz between 2015 Oct and 2016 Aug (Hayashi et al. 2021). As shown in Fig. 12, the radio SED can be fitted with a synthrotron emission model using an MCMC fitting technique (e.g., Goodwin et al. 2022; Zhang et al. 2024). In order to probe whether the radio emission is enhanced following the 2021 flare of F01004-2237, i.e., a nascent jet or outflow was launched, we proposed high frequency (Ku- and K-band) observations with VLA (Program ID 22A-507; PI, Shu) on 2022 July 21. The data were reduced following standard procedures with the CASA package (CASA Team et al. 2022), and the source is clearly detected at both bands. We used the IMFIT task in CASA to fit the radio emission component with a two-dimensional elliptical Gaussian model to determine the position, and the integrated and peak flux density. A comparison of the integrated flux to the peak flux density indicates that the radio emission is compact, and no extended component is detected. By plotting the integrated flux density in Fig. 12, we found no excess in the radio emission in comparison with the flux densities extrapolated from the best-fit synthrotron emission model. Therefore, there is no obvious radio brightening since the 2021 optical flare, at least at the time of our VLA radio observations.

IRAS F01004-2237 was also observed by the Very Long Baseline Array (VLBA) on 2022 May 6 (Hayashi et al. 2024) although the observation was not designed to study this flare. A compact source was detected at 8.4 GHz on a 1 pc scale. Its position is consistent with the galaxy center measured by Gaia with an offset of only 3 mas, and its peak flux is 0.6 mJy beam^-1. The compact source could be related to the long-standing AGN, or it could be related to the optical flare. With the available data, we cannot distinguish between these two possibilities. New high-resolution radio observations in the future may detect changes in the flux or position of this compact source, thereby solving the mystery of its origin.

We checked the level of possible nuclear activity before the flare occurred. This is not easy, as F01004-2237 shows vigorous starburst activity, which could overwhelm possible AGN signals. To make matters worse, there are $\sim 10^{5}$ Wolf-Rayet (WR) stars in the nucleus (Armus et al. 1988; Farrah et al. 2005; Tadhunter et al. 2017), producing some features similar to AGNs.

No broad emission lines (BELs) were detected in any of the optical spectra before the 2010 flare, including CTIO spectrum in 1985 (Armus et al. 1988), KPNO spectrum in 1995 (Veilleux et al. 1999; Lípari et al. 2003), HST/STIS spectrum in 2000 (Farrah et al. 2005), and WHT spectrum in 2005 (Rodríguez Zaurín et al. 2013). There is neither Pa$\alpha$ BEL in the NIR spectrum taken in 1995 (Veilleux et al. 1997). The non-detection of BELs, even in the NIR band, cannot be explained by the obscuration of extremely thick dust. This is because the vicinity of the SMBH must be seen directly as VBELs with FWHMs $>16,000$ km s^-1 were detected after the 2021 flare, with fluxes varying within months, and the ratio of Ly$\alpha$/H$\alpha\sim 10$ indicates little dust reddening. Therefore, the non-detection of BELs is because they are intrinsically weak.

In addition, the X-ray spectra before the 2010 flare can be explained by starburst alone, without needing a power-law AGN component (Teng & Veilleux 2010; Nardini & Risaliti 2011). This cannot be explained by a Compton thick AGN, as proposed by (Nardini & Risaliti 2011), because X-ray emission with a steep spectrum varying in weeks was detected after the 2021 flare, indicating little gas absorption in the line of sight. Thus, we concluded that the nuclear activity must be weak before the 2010 flare. We had measured a 2–10 keV luminosity of $2.7\pm 0.8\times 10^{41}$ erg s^-1 using the 2009 spectrum, consistent with that measured by Nardini & Risaliti (2011). We adopted this luminosity as the upper limit of X-ray luminosity of AGN in 2009 and estimated an upper limit of bolometric luminosity of $\sim 3\times 10^{42}$ erg s^-1 using a correction factor of 12 from Hopkins et al. (2007).

Some literature claimed to have found evidence for luminous AGN (e.g., Veilleux et al. 2009; Yuan et al. 2010). We have re-examined their results and will discuss them in section 6.1.

We discuss the two possible interpretations of the 2021 flare in terms of shapes of LC, UV/optical SED and luminosities, parameters of VBELs, and X-ray to UV/optical luminosity ratios.

TDEs typically have smooth LCs rising quickly and falling slowly, in line with the observation of the 2021 flare, which rose in $\sim$50 days and dropped by 2–3 magnitudes in one year. However, the LCs of AGN flares, taking observations of changing-look AGNs as examples (e.g., Ricci & Trakhtenbrot 2023), showed rich diversity without a uniform pattern. Some AGN flares have LCs resembling those of TDEs, such as AT 2018dyk (Frederick et al. 2019). Thus, the LCs support a TDE interpretation but do not exclude the possibility of AGN flare.

The UV/optical SEDs of the 2021 flare can be well described by blackbody curves with temperatures $\sim$22,000 K, consistent with optical TDEs (e.g., van Velzen et al. 2021; Hammerstein et al. 2023). Meanwhile, the UV-optical SEDs of AGNs are generally close to power-law curves. However, it is challenging to distinguish a blackbody with a temperature of $\sim$22,000 K and a slightly extinct power-law, as they have little difference in observable wavelength ranges. We tried to model the observed SEDs from Swift/UVOT with power-law curves, yielding similarly good fits as blackbody curves with no statistically significant $\chi^{2}$ differences. We did not make a similar distinction using spectral data considering the interference from observed VBELs, possible other VBELs from He and metal elements, and possible Balmer continuum and Fe II emission. If assuming a SED similar to AGNs, we estimated a peak bolometric luminosity of $\sim 9\times 10^{44}$ erg s^-1 using the measured peak $\nu L_{\nu}$ at 4400 Å of $9\times 10^{43}$ erg s^-1 and adopting a bolometric correction of $\sim$10 from Hopkins et al. (2007). We had demonstrated in section 4.1 that the pre-flare nuclear activity, if present, had a bolometric luminosity of $<3\times 10^{42}$ erg s^-1. Then, the variation has an amplitude of more than two orders of magnitude, which is rare among AGN flares. Thus, the UV/optical SED and luminosities support the TDE interpretation but do not exclude AGN flare.

The 2021 flare shows VBELs in H$\alpha$, H$\beta$, He II $\lambda$4686, He I $\lambda$5876 and Ly$\alpha$ between $+$105 and $+$325 days. The FWHM of H$\alpha$ can reach 16,000–19,000 km s^-1, and that of He II can even reach $>$21,000 km s^-1, although the latter relies on the decomposition of He II with H$\beta$. VBELs and strong He II emission distinguish TDEs from AGNs (Zabludoff et al. 2021). As can be seen from Fig. 13, the H$\alpha$ FWHMs and He II/H$\alpha$ ratios in F01004-2237 are consistent with those of TDEs (Charalampopoulos et al. 2022). To understand the properties of BELs in AGNs, we selected 3119 SDSS DR7 quasars with a signal-to-noise ratio of H$\alpha$ BEL $>$10 using measurements by Shen et al. (2011). Among them, 99.87% have H$\alpha$ FWHM $<$16,000 km s^-1. Shen et al. (2011) did not measure the fluxes of He II, nor could we find any other literature that systematically investigated the He II fluxes in AGNs. This is not surprising, as in most quasars and Seyfert 1 galaxies, He II BELs can not be detected because they are weak and overwhelmed in the Fe II bump. Vanden Berk et al. (2001) measured the He II and H$\alpha$ fluxes in the composite spectrum of SDSS quasars. With their measurements, we estimated an average He II/H$\alpha$ ratio in quasars of 0.0045, much smaller than in F01004-2237. We found four SDSS DR7 quasars (0.13%) with H$\alpha$ with FWHM $>$16,000 km s^-1, and show their spectra in Fig. 14. One of them is a well-studied quasar 3C 332. It shows BELs with clear double-peak structures, which are thought to be emission from an accretion disk viewed from an inclined line of sight (Halpern 1990). Of the other three, J1027+6050 and J1605+2309 also show H$\alpha$ with a double-peak structure, and J1334-0138 has an H$\alpha$ emission line profile that deviates from a single Gaussian. None of the four show clear He II BELs, as their weak emission features around 5000 Å can be explained by H$\beta$. All of these are inconsistent with VBELs in F01004-2237. Therefore, the properties of VBELs strongly favour the TDE interpretation other than the AGN flare interpretation.

Most of the time after the 2021 flare, the X-ray emission is weak. During two X-ray flares around $+$280 and $+$350 day, there are significant enhancements in X-ray emission lasting no more than weeks. On one hand, weak and highly variable X-ray emission is consistent with TDEs. Only 4 out of 17 ZTF TDEs in the sample of van Velzen et al. (2021) show clear X-ray detection. Their optical-to-X-ray ratios, expressed as ratios between UV/optical blackbody luminosity and 0.3–10 keV X-ray luminosity, are between 1 and 2000. The same ratio during XF1 in F01004-2237 is 2–6, which is calculated assuming a power-law or a blackbody X-ray spectrum, and is consistent with X-ray bright TDEs. On the other hand, X-ray emission is prevalent in AGNs. The relative strength of X-ray emission in AGNs is described using $\alpha_{\rm OX}$, calculated with the fluxes at 2 keV and 2500 Å (Tananbaum et al. 1979). We computed $\alpha_{\rm OX}$ of the 2021 flare using XMM-Newton data on $+$247 day and Swift/XRT data on $+$281 day, and list the results in Table 4. The former was considered to be a lower limit because the X-ray emission on $+$247 day does not necessarily come from the 2021 flare. The latter was calculated assuming a power-law or a blackbody X-ray spectrum, respectively. Using the correlation between $\alpha_{\rm OX}$ and $L(\nu)_{2500}$ in AGNs (Lusso et al. 2010), the predicted $\alpha_{\rm OX}$ is only $\sim 1.2$. The observed values are larger than the predicted values after considering the scatter of the correlation of 0.15, indicating that the X-ray emission is weaker relative to the UV emission during the 2021 flare than in AGNs. Therefore, the X-ray to UV/optical luminosity ratios favor the TDE interpretations.

In summary, we concluded that the 2021 flare is a TDE.

As mentioned in the introduction section, Tadhunter et al. (2017) considered the 2010 flare as a TDE other than an AGN flare because of the He II BEL with large He II/H$\beta$ ratio in the spectrum $\sim$5 years after the flare occurred. In recent years, a new class of BF AGN flares has been identified, and examples are AT 2017bgt, OGLE17aaj, AT 2021loi, and AT 2019aalc (Trakhtenbrot et al. 2019; Gromadzki et al. 2019; Makrygianni et al. 2023; Milán Veres et al. 2024). These flares are characterized by a strong N III $\lambda$4640 BF emission line with similar intensity to He II $\lambda$4686, both commonly seen in TDEs. Trakhtenbrot et al. (2019) and subsequent works demonstrated that these flares are not TDEs because their emission lines are narrower than typical TDEs and their LCs show double-peak structures, and suggested that they originate from a sudden enhancement of the long-existing accretion flow. They also placed the 2010 flare in F01004-2237 into this class as they claimed similarity between the spectra of F01004-2237 and BF AGN flares.

However, the claimed similarity is questionable because there are two crucial differences that need to be explained¹³¹³13We caution readers that it may be misleading to compare the 2015 spectrum directly to those of BF flares without removing the long-standing WR features and NELs from the host galaxy as some works did.. One is that the BELs in F01004-2237 are broader and show a significant change in width. The FWHMs of BELs are $\sim$6,000 km s^-1 in the 2015 spectrum (Tadhunter et al. 2017), and can reach 16,000–19,000 km s^-1 after the 2021 flare. For comparison, the BEL widths of the BF flares are between 2,000 and 5,000 km s^-1 and do not change dramatically over time. The other is that He II in F01004-2237 is stronger relative to H$\beta$. In the 2015 spectrum, the He II/H$\beta$ ratio is $1.82\pm 0.09$ if emission from all sources are included, and is $5.3\pm 0.9$ if contributions from WR stars and ISM are removed (Tadhunter et al. 2017). For comparison, the BF flares show He II/H$\beta$ ratios of 0.4–0.7. The large and variable width and the large He II/H$\beta$ ratio of BELs are prevalent in TDEs, and BF lines are more common in TDEs than in AGNs (e.g. van Velzen et al. 2020). Thus, we considered that the 2010 flare in F01004-2237 is more similar to TDEs than AGN flares.

We re-examined the nature of the 2010 flare in light of new observations of the 2021 flare. If the 2010 flare were indeed an AGN flare as proposed by Trakhtenbrot et al. (2019), the He II emission lines in the 2015 spectrum would indicate a broad line region with bulk velocities of several $10^{3}$ km s^-1. The broad line region would still exist after the 2021 flare as the dynamic time scale is typically $>10$ yr, and would produce an emission component with a similar width as in 2015. In fact, no such component in He II or other emission lines has been detected in early-time spectra taken after the 2021 flare. The difference in the ionization continuum can hardly explain the no detection because the 5100 Å luminosities at the time when these spectra were taken are $0.7-4\times 10^{43}$ erg s^-1, similar to that in 2015. Thus, the new observations do not support the AGN flare interpretation, and we preferred that a TDE causes the 2010 flare. Note that the questions raised by Trakhtenbrot et al. (2019) about the TDE interpretation have been answered by Tadhunter et al. (2021) and Cannizzaro et al. (2021).

If related, the recurring flares in F01004 may be caused by repeating partial TDEs or double TDEs.

After a partial TDE occurs, the stellar remnant may return, producing the next TDE. The partial TDE of a star intruding in a parabolic orbit could transfer the stellar remnant into an elliptical orbit due to the weak asymmetry of the disrupted material. However, the period would be $>400$ yr (Ryu et al. 2020), which cannot explain the 10.3 yr interval between the two flares. Therefore, repeating TDEs with an interval of 10.3 yr requires that the star was already in an elliptical orbit around the SMBH with a period of $\sim 10$ years before the disruption.

The observed period corresponds to an orbital semi-major axis $a\sim 1400$ AU for an SMBH mass of $\sim 2.5\times 10^{7}$ $M_{\odot}$. For a partial TDE to occur, the pericenter $r_{p}$ must be similar to $r_{t}$, which is $\sim 1.4$ AU for sun-like stars. Then, the orbit must be highly eccentric with $1-e\sim 10^{-3}$. Such an orbit can be explained by the Hills mechanism. If so, the orbital semi-major axis of the captured star is related to the initial orbital parameters of close binary stars (Pfahl 2005) as:

where $a_{b}$ and $m_{b}$ are the binary semi-major axis and total mass, and $m_{\rm ej}$ is the mass of the ejected star. Assuming both the binary stars have a mass of 1 $M_{\odot}$, the observed period corresponds to $a_{b}\approx 8r_{\odot}$. The inferred binary semi-major axis is reasonable as it is far enough to avoid a common envelope, and close enough to avoid disintegration in collisions with fellow stars (Hills 1988).

Simulations show that the flare produced by a partial TDE has a steeper descent in the late-time, obeying $L\propto t^{-9/4}$ rather than $L\propto t^{-5/3}$ as in full TDEs (e.g., Coughlin & Nixon 2019; Miles et al. 2020). We tested this with the observations of the 2021 flare in F01004-2237. We fit the ATLAS/o, ATLAS/c and Swift UVW2 LCs after $+$60 day with t^-9/4 curves. The procedure was similar to the fitting with t^-5/3 curves (Section 3.2), except that we replaced the powerlaw index in the equation (3) with $-9/4$. This slightly improves the fit as the $\chi^{2}$/d.o.f. decreases from 3.57 to 3.31, with d.o.f.$=324$. Therefore, t^-9/4 curves agree slightly better with the data than t^-5/3 curves, which supports the partial TDE model.

There is a rich diversity of possible outcomes for the binary-SMBH encounter, with one possibility being that both stars are tidally disrupted (Mandel & Levin 2015; Mainetti et al. 2016). This phenomenon, known as double TDE, occurs when the pericenter of the binary’s centre of mass is close enough to the black hole. However, the time interval between the two TDEs is generally less than $\sim 10^{2}$ days, resulting in only one flare being observed, which is inconsistent with the situation in F01004-2237.

In the double TDE regime, two possible scenarios could explain the interval of 10.3 yr. One scenario was proposed by Mandel & Levin (2015) to interpret the two flares 20 years apart in IC 3599. A possible outcome of the binary-SMBH encounter is that one star is directly disrupted, causing the first TDE, and then the other star is captured into an elliptical orbit and produces the second TDE when it returns to the pericenter. According to the simulations of Mandel & Levin (2015), the orbital period of the captured star can range from 6 months to several decades, explaining the observations in F01004-2237.

The other scenario involves the encounter between a stellar binary and a mpc-scale SMBHB binary. As shown in the simulations by Wu & Yuan (2018) and Coughlin et al. (2018), there is a significant probability of delayed double TDEs, among which the time interval of two TDEs can be more than ten years. A mpc-scale secondary SMBH may pass through the stellar debris stream of the TDE as it orbits the primary SMBH and cause periodic gaps on the X-ray LC (Liu et al. 2014; Shu et al. 2020). Unfortunately, the X-ray emission from the flare in F01004-2237 is weak and cannot be used to test the presence of mpc-scale secondary SMBH.

We considered the possibility that the two flares are not associated physically. The probability of a TDE (the 2021 flare) occurring in a time interval $\Delta t$ after a specific event (the 2010 flare) is $p=r_{\rm TDE}\times\Delta t$, where $r_{\rm TDE}$ is the incidence rate of TDE per galaxy in a unit of yr^-1. Adopting the rate of optical TDEs of $3.2^{+0.8}_{-0.6}\times 10^{-5}$ yr^-1 (Yao et al. 2023) and $\Delta t=10.3$ yr, we estimated a small probability of $p\sim 3.3\times 10^{-4}$.

For $p$ to increase to $>0.1$, the TDE rate needs to be $>10^{-2}$ yr^-1, $>300$ times higher than in normal galaxies. Thus, assuming the two flares are independent requires an extremely high TDE rate in F01004-2237. On the contrary, assuming that the two flares are repeating partial TDEs or double TDEs does not require such a high TDE rate because they should be treated as “a single event” when calculating the probability.

As there are many theories to explain the high TDE rate in post-starburst galaxies, we examined if they also apply to F01004-2237. First, a starburst can cause an unusually high concentration of stars in the nucleus, enhancing the TDE rate (Stone & Metzger 2016; Stone & van Velzen 2016). In theory, this effect subsides over time (Stone et al. 2018). The starburst in F01004-2237 started only 3–6 Myr ago, as estimated by the age of the WR stars (Tadhunter et al. 2017). Based on the results of Stone et al. (2018), we inferred that the TDE rate in F01004-2237 may be 1–2 orders of magnitude higher than that in post-starburst galaxies with stellar ages of $\sim$100 Myr, and hence 2–3 orders of magnitude higher than in normal galaxies.

Second, a nuclear star cluster (NSC) can enhance the TDE rate by orders of magnitude in galaxies with different masses and various types (Pfister et al. 2020). The HST images show a point-like source in the nucleus of F01004-2237, with SED consistent with a young stellar population with an age of $<10$ Myr and a mass of $10^{8}$ $M_{\odot}$ (Surace et al. 1998). The source is barely resolved, indicating a $<100$ pc size. This supports the existence of an NSC in F01004-2237, with mass and age unusual among NSCs in galaxies. How much such an NSC can raise the TDE rate remains to be answered by theoretical studies.

Third, the presence of a companion SMBH with a mass ratio of 0.01–0.1 can significantly boost the TDE rate (e.g., Ivanov et al. 2005). According to the simulations (e.g., Chen et al. 2009; Li et al. 2017), in a phase lasting for $\sim$Myr, during which the SMBHs are $\sim$pc apart, the TDE rate could be as high as $10^{-2}\sim 1$ yr^-1. Although there is no direct evidence, F01004-2237 may also host a pc-scale SMBH binary, given that starbursts are generally associated with mergers, which can lead to SMBH binaries. There are no obvious merger signals in the HST image (Surace et al. 1998). However, this cannot negate the scenario of SMBH binary, because non-equal mass SMBH binaries may be related to minor mergers, and fall to the pc-scale long after mergers.

Therefore, a high TDE rate of $>10^{-2}$ yr^-1 in F01004-2237 is possible in theory. As a result, we cannot rule out the possibility that the two flares are independent.

Some literature claimed to have found evidence for luminous AGN based on Seyfert-like optical narrow emission line (NELs) ratios (e.g., Allen et al. 1991; Yuan et al. 2010). We noticed that the classification based on narrow emission lines is controversial in the literature, as Veilleux et al. (1999) classified it as star-forming. The optical NELs show a narrow component with FWHM of $\sim$90 km s^-1 and broad components blueshifted by $\sim 1000$ km s^-1, which is considered to be related with outflow (Lípari et al. 2003; Rodríguez Zaurín et al. 2013; Tadhunter et al. 2021). Rodríguez Zaurín et al. (2013) shows that while the narrow component has starforming-like line ratios, the broadest component has Seyfert-like line ratios. This may be why previous literature made contradictory classifications, as pointed out by Tadhunter et al. (2017). The outflow component that shows Seyfert-like line ratios has a radial extent to $\sim$3.5 kpc, as measured by Rodríguez Zaurín et al. (2013) using long-slit spectroscopy. Meanwhile, it is not detected in the STIS 2000 spectrum extracted within an aperture corresponding to 350 pc in physical (Farrah et al. 2005). Thus, its distance to the centre is $>350$ pc and thus traces AGN activity $\sim 10^{3}-10^{4}$ years ago. However, the levels of nuclear activity may vary dramatically on such a time scale, an example being Arp 187 (Ichikawa et al. 2019). We are more concerned about recent activity, for which NELs in the nuclear spectrum are needed for diagnosis. As measured by Farrah et al. (2005), the NEL line ratios in the nuclear region fall around the classification lines between AGN and starforming (Kewley et al. 2001, 2006), and classifications based on [S II]/H$\alpha$ and [O I]/H$\alpha$ are different. Note that the result depends on the choice of classification lines; for example, the line ratios will be the explicit starforming type if the classification lines of Veilleux & Osterbrock (1987) are used. The classification lines of Kewley et al. were calculated assuming the Salpeter’s initial mass function, while the presence of a large population of WR stars in the nuclear region of F01004-2237 may result in line ratios closer to Seyfert galaxies than in normal star-forming galaxies. Thus, spectral analysis in the nuclear region does not show evidence of strong AGN in recent $10^{3}$ years. This is consistent with the absence of BELs and the weak hard X-ray emission.

The MIR spectrum of F01004-2237 shows only week PAH emission features with equivalent widths lower than expected for starburst galaxies (Tran et al. 2001; Imanishi et al. 2007; Veilleux et al. 2009). This was considered to be evidence of AGN, whose intense X-ray and EUV emissions can destroy PAHs (e.g., Voit 1992). However, a weak AGN at present and a strong AGN in the past can also explain the destruction of PAH, and a strong AGN at present is not necessary.

F01004-2237 shows a hotter IR SED than starforming galaxies. The IR SED indexes, such as flux ratios between 30 and 15 $\mu$m, and between MIR to FIR, are all consistent with ULIRGs with clear AGNs (Veilleux et al. 2009). However, the MIR continuum of F01004-2237 shows only a weak 9.7 $\mu$m silicate absorption feature, and does not show strong dust temperature gradients (Imanishi et al. 2007; Veilleux et al. 2009), both of which are expected to be signatures of a buried AGN. This leads us to suspect that the warm dust radiating the MIR continuum may have illumination sources other than the AGN. As previously mentioned, F01004-2237 hosts an NSC with an age of $<10$ Myr, a mass of $10^{8}$ $M_{\odot}$, and a size of $<100$ pc. Such an NSC may create a dense radiation field, producing a higher dust temperature than normal star-forming galaxies. Moreover, considering that two energetic flares occurred in recent 14 years in F01004-2237, other flares possibly occurred in decades prior to the Spitzer observations in 2004, although no observations were available to test them. The 2010 flare was accompanied by a bright dust echo lasting for at least ten years, resulting in increases of monochromatic luminosities at 3.4 and 4.6 $\mu$m of $2-3\times 10^{44}$ erg s^-1 (Dou et al. 2017). If the hypothetic earlier flares also had similar total energies, they might significantly increase the temperature of the nuclear dust, producing an IR SED mimicking those of AGNs. Therefore, it is not reliable to estimate the level of AGN with MIR luminosity.

In this section, we compared the recurring flares in F01004-2237 and other galaxies. We check whether all of them can be interpreted under the framework of repeating partial TDEs, as this is how most of the literature explains them. The information on these flares is listed in Table 6. X-ray QPEs are not listed here because their time scales, luminosities and energy budgets are quite different from flares in F01004-2237.

In the scenario of repeating partial TDEs formed by the Hills mechanism, the outcome of the captured star depends on parameter $\beta\equiv r_{t}/r_{p}$. The star is partially disrupted if $\beta$ is greater than the critical value and is not disrupted otherwise. Depending on $\beta$, repeating partial TDE can be divided into two groups. In one group, $\beta$ is only slightly above the critical value, causing only a tiny fraction of stellar mass to be disrupted, leaving flares with low energies and short durations (e.g., Ryu et al. 2020; Nixon et al. 2021). Since a partial TDE has little effect on $r_{p}$ and the stellar mass decreases only a little, $\beta$ is almost constant and similar partial TDEs can repeat many times. In the other group, $\beta$ is significantly larger than the critical value, causing a large fraction to be disrupted, resulting in energetic flares. Only several flares may be detected because the remaining flares are too weak due to small residual stellar mass.

This scenario can explain the diversity of the observational features of recurring flares. We divided all recurring flares into two groups based on the above analysis. The first group include ASASSN-14ko and X-ray QPEs, where $>9$ flares were detected, and a single flare has a UV energy $<10^{51}$ erg or an X-ray energy $<10^{50}$ erg. The second group include the other recurring flares. We preferred that the difference in $\beta$ is decisive, while the other differences are not essential. The two groups appear to have different periods $P$ ($P<114$ days for the first group and $P>223$ days for the second group). However, this might be a selection effect. Long-period, low-energy repeating partial TDEs might exist but be extremely difficult to identify. Short-period, energetic partial TDEs might also exist but be misclassified as a single TDE because the flares overlap. In addition, most flares in the first group were detected in the X-ray band, while flares in the second group could be detected in both the UV/optical and X-ray bands. The possible reason is that a larger debris mass leads to a higher chance of forming an optically thick envelop through debris collision or outflow, producing strong UV/optical emission (e.g., Dai et al. 2018).

In the scenario of repeating partial TDEs, will the debris mass of the next TDE be greater than, less than, or roughly equal to that of the previous one? It is challenging to give an exact answer because the remnant fraction has a complex relationship with the stellar mass and $\beta$ parameters, as shown in Fig. 4 of Ryu et al. (2020). To complicate matters further, stellar remnants that survive a partial TDE do not have enough time to reach equilibrium before returning (Ryu et al. 2020), and hence, when calculating the remnant fraction of the next TDE, the results of stars of the same mass cannot be simply reproduced. What is clear, however, is that the debris in the next TDE is unlikely to be much more massive than in the previous one. This is because if the stellar core survives a partial TDE as only (part of) the stellar envelope lost, its structure changes little before returning to the pericenter since the time scale of thermal relaxation is much longer than the orbital period (Ryu et al. 2020), and thus it is likely to survive the next partial TDE.

We tested this with observational data. Most recurring flares meet theoretical expectations as the next flare is weaker than or similar to the previous one in total energy and peak luminosity. The only exception is AT 2020vdq, where the second flare is 15 times brighter than the first and has total energy 30 times higher (Somalwar et al. 2023). This is difficult to understand in the scenario of repeating partial TDEs. Meanwhile, in the double TDEs scenario, two flares are produced by two different stars in the binary system, and there is no direct correlation between their luminosities and energies. In addition, only two flares have been detected so far for AT 2020vdq, not contradicting a double TDEs scenario. Thus, AT 2020vdq may be a case of double TDEs, rather than repeating partial TDEs.

If AT 2020vdq is double TDEs, then the other recurring flares appearing only twice may also be double TDEs, because situations where the second flare is stronger than, weaker than or similar to the first one can all occur in the double TDEs scenario. However, it is more likely that ASASSN-14ko and eRASSt J045650.3-203750 are repeating partial TDEs because they show multiple flares. Therefore, we proposed that the recurring flares detected so far could not be explained with only one scenario, and repeating partial TDEs and double TDEs may both exist.

For sources that have flared only twice currently, the best way to distinguish between the repeating partial TDEs and the double TDEs scenarios is to wait and check if there is an expected third flare. We list the timing of these expected flares in the repeating partial TDEs scenario in Table 6. The next flare expected in F01004-2237 will peak between December 2032 and August 2033. However, it may come earlier, considering that the flare interval in eRASSt J045650.3-203750 is gradually shortening (Liu et al. 2024). A flare that arrives at the scheduled time and is weaker than the previous one supports an interpretation of repeating partial TDEs, while the nondetection of the expected flare supports interpretations of double TDEs or independent TDEs. Before that, AT 2018fyk, AT 2020vdq and AT 2022dbl will go to trials in 2025 and 2026.

In section 5.4, we demonstrated that a high TDE rate of $>10^{-2}$ in F01004-2237 is required if the two flares are independent. Theoretically, a galaxy with a young starburst and massive NSC, like F01004-2237, may have a high TDE rate. This can be tested by future monitoring of a large number of WR galaxies (e.g., Brinchmann et al. 2008) with similar properties. With the advent of new generation sky surveys such as the 2.5-meter Wide Field Survey Telescope (Wang et al. 2023) and the Large Synoptic Survey Telescope (Ivezić et al. 2019), TDEs with magnitudes $<23$ can be detected, meaning that the monitoring we mentioned earlier can be performed in $z<1$. Within this large volume, there are at least a few hundred galaxies with properties similar to F01004-2237, taking into account the cosmological evolution of the star formation rate. Monitoring them over several years is sufficient to prove or disprove the possible high TDE rate.