Live to Die Another Day: The Rebrightening of AT 2018fyk as a Repeating Partial Tidal Disruption Event

The classic prediction for the mass fallback rate generated by a star being tidally disrupted by a supermassive black hole (SMBH) is an asymptotic, t^−5/3 decay (Rees 1988; Phinney 1989). While some of the tidal disruption events (TDEs) identified so far have displayed such long-term behavior, a significant fraction show different light curve evolution, which in some cases is completely decoupled from the mass fallback rate (e.g., Gezari et al. 2017; Kajava et al. 2020) and may be expected (e.g., Guillochon & Ramirez-Ruiz 2013; Hayasaki & Jonker 2021). Auchettl et al. (2017) found that the X-ray light curves can be well described by power-law indices ranging from –0.5 to –2. Hammerstein et al. (2022) defined three types of behavior for the UV/optical light curves, labeling them power-law decay (with indices ranging from –1 to –3 and a sizable fraction that decay consistent with a t^−5/3 law), plateau, and structured light curves. Deviations from the late-time t^−5/3 decay have also been suggested theoretically; Hayasaki et al. (2013) and Cufari et al. (2022a) found that stars on eccentric orbits can lead to a prompt shutoff in the light curve, while Guillochon & Ramirez-Ruiz (2013) found that partial TDEs—in which the dense stellar core survives the tidal encounter with the SMBH—can lead to significant deviations from t^−5/3. Coughlin & Nixon (2019) predicted that partial TDEs should generically exhibit a t^−9/4 decay. More dramatic deviations, including truncation, order-of-magnitude dips, and reflaring, can be induced by TDEs in SMBH binaries (e.g., Liu et al. 2009; Ricarte et al. 2016; Coughlin et al. 2017). Recently, the source ASASSN-14ko was interpreted to be a repeating partial TDE, such that the star is on a bound orbit about the SMBH and partially stripped of its mass—thus feeding a new accretion flare—each pericenter passage (Payne et al. 2021).

In this work, we report on the renewed X-ray and UV activity of the transient AT 2018fyk, a proposed TDE originally described in Wevers et al. (2019), ∼1200 days after discovery. We compare the observational properties to those of the previously observed accretion flow properties in Section 2, after which we explore a repeating partial TDE scenario in Section 3 to explain the long-term properties. We present the implications and predictions of our model in Section 4 before summarizing and concluding in Section 5.

2.1. A Brief History and Basic Properties of AT 2018fyk

2.2. Rebrightening at Very Late Times

Following the dramatic dimming after ∼600 days seen at X-ray (by a factor of >6000) and UV (by a factor of ∼15) wavelengths, SRG/eROSITA (Predehl et al. 2021) scanned the position of AT 2018fyk four times at phases of 611, 796, 978, and 1163 days after discovery. It was not detected in these epochs (see Figure 1), providing 3σ upper limits (in the 0.3–2 keV band) of L_X ≈ 1–3 × 10⁴² erg s⁻¹.

Zoom In Zoom Out Reset image size

Fifty-three days after the last eROSITA nondetection, the source was detected again by Neil Gehrels Swift (hereafter Swift) monitoring observations obtained 1216 days after discovery with a luminosity of 8 × 10⁴² erg s⁻¹. This implies a relatively quick reappearance of the X-ray emission. The X-ray brightness has increased by a factor of at least 100 compared to the deepest upper limit 700 days before and a factor of 2–3 compared to the last eROSITA upper limit; the UV emission (0.03–3 μm) has also brightened by a factor of ≈10 to L_UV = 7 × 10⁴² erg s⁻¹. Such behavior is both unprecedented and unexpected in the classical scenario of a star being fully disrupted by the SMBH. The data reduction for all new observations used in this work is described in the Supplementary Materials.

2.2.1. SED and X-Ray Spectrum

By modeling the available (host galaxy–subtracted) Swift UV data with a blackbody function, we find a blackbody temperature of ∼25,000–35,000 K, similar to the temperature at early times. We use this temperature to convert the UVW1 luminosity into the 0.03–3 μm emission ¹⁰ (representing the total UV/optical emission, L_UV).

We use XSPEC (Arnaud 1996) to model a new XMM-Newton (EPIC/PN) X-ray spectrum (obtained through director's discretionary time) with a phenomenological model (TBAbs×(diskbb + powerlaw)) consisting of a thermal component and a power law absorbed by a Galactic column of n_H = 1.15 × 10²⁰ cm⁻² (i.e., the same model used in Wevers et al. 2021). We find a temperature of 113 ± 31 eV for the thermal component, a power-law index of Γ = 2.2 ± 0.2, and a power-law fraction of emission (defined as the ratio of the power-law flux to the total X-ray flux in the 0.3–10 keV band) of 80% ± 10% (full details are provided in Table A1). These values are all consistent with the previous hard state properties. Given these parameters, we calculate the conversion factor to first translate the count rate to 0.3–10 keV luminosity and then convert this luminosity into the 0.01–10 keV luminosity (L_X). We then calculate the Eddington fraction of emission as ${{\rm{f}}}_{\mathrm{Edd}}=\tfrac{{{L}}_{\mathrm{UV}}+{{L}}_{X}}{{{L}}_{\mathrm{Edd}}}$; that is, we assume that L_bol = L_UV + L_X.

Table A1. Best-fit Parameters Obtained from X-Ray Spectral Modeling of the (Stacked) Swift and XMM Data

Spectrum	Count Rate	State	${t}_{\exp }$	kT	Norm(kT)	Γ	log10(norm (Γ))	log10(flux)	PL Frac.	χ2 (dof)
XRT	0.011	F	45,650	175 ± 60	${8}_{-6}^{+42}$	2.15 ± 0.4	–4.2 ± 0.2	–12.35 ± 0.04	79 ± 10	23 (24)
PN (0601)	0.25	F	9200	113 ± 31	${102}_{-71}^{+458}$	2.16 ± 0.2	–4.10 ± 0.08	–12.37 ± 0.03	80 ± 10	114 (113)
PN (1401)	0.30	F	3500	152 ± 80	${24}_{-20}^{+900}$	2.4 ± 1.4	–4.07 ± 0.5	–12.24 ± 0.2	80 ± 10	203 (170)

Note. The mean count rate for each spectrum is given in the second column. The effective exposure time ${t}_{\exp }$ is given in kiloseconds. The spectral model used is TBabs*zashift*(diskbb + powerlaw); kT is the temperature of the thermal component, while Γ denotes the power-law spectral index. Normalizations for the thermal and power-law components are listed in the norm(kT) and norm(Γ) columns. The flux is integrated from 0.3–10 keV. "PL Frac." denotes the fractional contribution of the power-law component to the total X-ray flux. The final column lists the reduced χ² and degrees of freedom (dof).

Download table as:  ASCIITypeset image

Combining this power-law fraction and index (consistent with the values obtained from Swift/XRT data) with the host-subtracted and extinction-corrected Swift/UVW1 fluxes, we calculate the UV–to–X-ray slope α_OX of the late-time emission. The full light curve and α_OX as a function of bolometric Eddington ratio f_Edd are compared to the earlier evolution in Figure 1 (top left and top right panels, respectively), while the late-time SED (including the XMM4 and Swift/UVOT data) is shown in Figure 2.

Zoom In Zoom Out Reset image size

We find that the spectral properties of AT 2018fyk are very similar to those observed in the previous hard state observations (states C and D), just before the source became faint around day 500. The EPIC/PN light curve (Figure 1, bottom right panel) does not show statistically significant variability on timescales of 100–1000 s, although the uncertainties are large due to the relatively low count rate. ¹¹ When NICER restarted monitoring observations with a roughly twice daily cadence, several X-ray flaring episodes were observed (Figure 1, bottom left panel), which is not evident from the (lower-cadence) Swift/XRT light curve. Significant variability on timescales of 6–12 hr is present throughout the NICER observations.

2.3. Previous Models for the X-Ray and UV Dimming

2.4. The Host Galaxy Is Not an AGN

In order to investigate the presence of an AGN, we inspect publicly available Multi Unit Spectroscopic Explorer (MUSE) and X-shooter data (see Supplementary Materials) for emission lines. After modeling and subtracting the stellar continuum (see Figure A2), both the MUSE and X-shooter data show very weak (EW ∼ 0.5–1 Å) emission lines. We measured the line strengths and ratios to produce the BPT (Baldwin et al. 1981) diagnostic diagram. The host galaxy is located in the low-ionization nuclear emission-line region (LINER; Figure 3). First thought to be exclusively produced by weak AGNs (e.g., Heckman 1980; Kewley et al. 2006), other ionization mechanisms (unrelated to accretion) can also produce consistent line ratios (Stasińska et al. 2008). The WHAN diagram (Cid Fernandes et al. 2011) can differentiate a weak AGN from a retired galaxy ¹² by substituting the [O iii]/Hβ ratio for the equivalent width (EW) of Hα. Figure 3 shows the classification using AT 2018fyk's host as a retired galaxy by the WHAN diagram. The lack of an increase in the Hα EW toward the nucleus (see Figure A1) also supports the non-AGN scenario.

Zoom In Zoom Out Reset image size

To further investigate this interpretation, we also look at the infrared data: (i) prior to 2018, the WISE (Wright et al. 2010) IR W1–W2 color of the host galaxy is ≈0.05, inconsistent with IR AGN selection criteria (e.g., Stern et al. 2012); and (ii) Jiang et al. (2021) have analyzed the NEOWISE IR light curve of AT 2018fyk and measured a covering factor (f_c ) ¹³ ∼ 0.01, indicating a dust-/gas-poor circumnuclear environment unlike those found in AGNs (f_c > 0.3; Roseboom et al. 2013). These results, consistent between data sets and wavelengths, provide the most robust evidence to date for the absence of an AGN in AT 2018fyk. The implication is that AGN variability is strongly disfavored to explain the dramatic UV and X-ray variability seen in AT 2018fyk.

For a 10^7.7 M_⊙ SMBH, at most (for maximal spin), ≲5% of stars with masses and radii comparable to those of the Sun (or smaller) will enter within the tidal radius, be destroyed completely, and not swallowed whole (Kesden 2012; Ryu et al. 2020; Coughlin & Nixon 2022). The tidal radius is also highly relativistic, suggesting that—even for partial TDEs—disk formation will be prompt, which is consistent with the observed properties of AT 2018fyk (e.g., the presence of low-ionization Fe ii lines in the optical spectrum, the persistent X-ray brightness at UV/optical peak, the thermal X-ray spectrum at early times, and its short-timescale variability in the X-ray; Wevers et al. 2019, 2021). These arguments suggest that the star that initially fueled the outburst from AT 2018fyk by virtue of producing an observable flare was partially disrupted (most TDEs will result in unobservable direct captures for the high black hole mass; see also Coughlin & Nixon 2022). Typically, tidally disrupted stars are on approximately parabolic orbits (e.g., Merritt 2013), which begs the question of how a partial TDE could yield a rebrightening because, as noted by Cufari et al. (2022a), tidal dissipation within the partially disrupted star yields a minimum orbital period of a few × 10³ yr for a 10^7.7 M_⊙ SMBH (see their Equation (1)). One can bind the partially disrupted star more tightly if the star was initially part of a binary system that was destroyed through Hills capture (Hills 1988). In this case, the orbital period one would expect for the captured star is (Cufari et al. 2022b)

$$\begin{eqnarray}&&{T}_{\mathrm{orb}}\simeq \displaystyle \frac{2\pi }{\sqrt{{{GM}}_{\star }}}{\left(\displaystyle \frac{{a}_{\star }}{2}\right)}^{3/2}{\left(\displaystyle \frac{M}{{M}_{\star }}\right)}^{1/2},\end{eqnarray} \tag{ 1 }$$

where a_⋆ is the binary semimajor axis, and M_⋆ is the mass of the primary. A schematic of the different phases of the repeated partial disruption scenario and the timescales involved is shown in Figure 4.

Zoom In Zoom Out Reset image size

With a host galaxy velocity dispersion of σ = 158 km s⁻¹, the maximum separation that a binary can have and still survive in the galactic nucleus is a_⋆ ≲ GM_⋆/(4σ²) ≃ 0.01 au (e.g., Hills 1975; Gould 1991; Quinlan 1996; Yu 2002). With a_⋆ = 0.01 au, M_⋆ = 1 M_⊙, and M = 10^7.7 M_⊙, Equation (1) gives T_orb ≃ 2.5 yr. A dynamical exchange can therefore produce a star on an orbit about the SMBH with a period as short as ∼a few years. For separations ≲0.01 au, the tidal disruption radius of the binary is comparable to the tidal disruption radius of the star (increased by stellar rotation and relativistic effects; Gafton et al. 2015; Golightly et al. 2019; Gafton & Rosswog 2019), and a partial TDE will occur (Figures 4(a) and (b)). The tight required separation of the initial binary provides constraints on the maximum size of the stars, in this case, ≲2 R_⊙. Such systems would require either two low-mass stars or a main sequence–compact object binary; the latter (with the main-sequence star captured) is favored in order to reproduce the overall energetics and timescales of the TDE, as we now discuss (see Section 4 for additional motivation for this type of binary).

Upon being partially disrupted, the material returns to the SMBH on a timescale that is approximately (Lacy et al. 1982; Rees 1988)

$$\begin{eqnarray}&&{T}_{\mathrm{acc}}\simeq \displaystyle \frac{2\pi }{\sqrt{{{GM}}_{\star }}}{\left(\displaystyle \frac{{R}_{\star }}{2}\right)}^{3/2}{\left(\displaystyle \frac{M}{{M}_{\star }}\right)}^{1/2}\end{eqnarray} \tag{ 2 }$$

and has a peak magnitude

$$\begin{eqnarray}&&\begin{array}{l}{{L}}_{\mathrm{acc}}\simeq 1.5\times {10}^{45}\left(\displaystyle \frac{\eta }{0.1}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-3/2}\\ \times {\left(\displaystyle \frac{M}{{10}^{7}{M}_{\odot }}\right)}^{-1/2}{{\rm{}}\,{\rm{erg}}\,{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 3 }$$

though these are generally longer and lower, respectively, for partial disruptions (e.g., Guillochon & Ramirez-Ruiz 2013; Miles et al. 2020; Nixon et al. 2021; see also Section 4 below). The proportionality coefficient in Equation (3) matches simulations that yield a peak accretion rate equal to Eddington for M = 10⁷ M_⊙, M_⋆ = 1 M_⊙, and a radiative efficiency η = 0.1 (Wu et al. 2018). Setting M = 10^7.7 M_⊙ and taking solar-like values gives T_acc ≃ 0.8 yr and L_acc ≃ 6.7 × 10⁴⁴ erg s⁻¹ (Figure 4(b)). Partial TDEs typically rise, peak, and decay as ∝ t^−9/4 (Coughlin & Nixon 2019; Miles et al. 2020; Nixon et al. 2021), but for a star on a bound orbit, the fallback rate plummets as the star returns to pericenter (Liu et al. 2022). The reason for this sharp decline in the fallback rate is that the stellar core has a Hill sphere—an approximately spherical region within which the star's gravitational field dominates over that of the SMBH—near which the stream density is much smaller than that of the bulk of the stream (Figure 4(e)). This feature of the fallback rate can physically explain the rapid shutoff displayed in Figure 1 at ∼600 days.

While it likely does not inhibit the formation of a disk, nodal precession—assuming the SMBH has a modest spin—is probably important for its subsequent evolution; over many orbits of the material in the innermost regions of the disk, nodal and apsidal precession, coupled to the (likely) large misalignment angle between the spin axis of the SMBH and the angular momentum of the gas, will cause fluid annuli to precess independently instead of conforming to a smooth, warped disk (Nixon et al. 2012; Liska et al. 2021). The orbit of the returning star also precesses and leads to a time-dependent feeding angle of the flow; thus, the gas is likely morphologically complex and, we suggest, closer to spherically symmetric than in the form of a traditional disk (see also Patra et al. 2022). If we assume that the returning debris stream is cylindrical with a cross-sectional radius of ∼R_⊙ and length $a\simeq \tfrac{{a}_{\star }}{2}{(M/{M}_{\odot })}^{2/3}$ (Cufari et al. 2022b), then taking a_⋆ = 0.01 au, M = 10^7.7 M_⊙, and with 0.05 M_⋆ contained in the stream (see Section 4),

$$\begin{eqnarray}&&\rho \simeq \displaystyle \frac{0.05{M}_{\odot }}{\pi {R}_{\odot }^{2}{a}_{\star }{\left(M/{M}_{\odot }\right)}^{2/3}}\simeq 3.2\times {10}^{-7}{{\rm{}}\,{\rm{g}}\,{\rm{cm}}}^{-3}.\end{eqnarray} \tag{ 4 }$$

Taking v = 0.1c as the speed of the material as it shocks (recall that the pericenter is highly relativistic), the shocked-gas pressure is

$$\begin{eqnarray}&&p\simeq \rho {v}^{2}\simeq 2.9\times {10}^{12}{{\rm{}}\,{\rm{erg}}\,{\rm{cm}}}^{-3}.\end{eqnarray} \tag{ 5 }$$

The fluid is radiation pressure–dominated with a temperature

$$\begin{eqnarray}&&T\simeq {\left(\displaystyle \frac{3p}{a}\right)}^{1/4}\simeq 5.8\times {10}^{6}{\rm{}}\,{\rm{K}}\simeq 0.5{\rm{}}\,{\rm{keV}}.\end{eqnarray} \tag{ 6 }$$

Equation (6) represents the self-intersection temperature near the horizon. The gas expands roughly adiabatically from the self-intersection point (e.g., Jiang et al. 2016), which reduces the temperature and density. At early times, the gas will be optically thick, the photosphere at large radii, and the peak emission at temperatures below Equation (6). However, as time advances, the fallback rate declines, the density drops due to the continued expansion of the gas, and the flow becomes more optically thin to reveal the hot, inner regions, thus providing a plausible interpretation of the late-time dominance of the X-ray emission.

From Wevers et al. (2019) and the additional data obtained since then, the total amount of energy radiated is equivalent to ≈9 × 10⁵¹ erg. This energy could be up to a factor of ∼5 smaller ¹⁴ if the majority of the energy is not radiated at UV/optical wavelengths, as we have assumed in calculating the bolometric luminosity (and thus the total radiated energy; see Figure 2). If we adopt a radiative efficiency of η = 0.1, the radiated energy amounts to ≈0.05 M_⊙ of accreted mass. In normal TDEs (i.e., where the center of mass is on a parabolic orbit), approximately half of the stellar mass is accreted, and this implies that the star lost at most ∼0.1 M_⊙ during the tidal encounter. We note that for typical binaries, the ratio of the binding energy of the binary to that of its stellar constituents is very small (on the order of the ratio of the stellar radius to the binary separation); hence, the approximation that only half of the material is accreted is usually warranted. Here, however, the binary must be very tight to reproduce the observed timescales, meaning that the binding energy of the binary is not substantially smaller than that of the star itself, and the "unbound debris" featured in panel (b) of Figure 4 may actually remain bound to the SMBH. If this is the case, we expect the luminosity of the "lesser bound" tail to be significantly lower than that of the more tightly bound tail owing to the longer return time. This material may be of sufficiently low density that it is substantially affected/destroyed by interactions with circumnuclear gas (Bonnerot et al. 2016) and the surviving core as it passes through pericenter a second time. Additional and more detailed investigations are required to constrain the energetics of the unbound/less bound tail.

Because the surviving core is spun up to near its breakup velocity, the tidal radius moves out (Golightly et al. 2019), and it is possible that the star was completely destroyed on its second pericenter passage (Figure 4(f)). If the mass lost from the star is closer to the upper limit of ∼0.1 M_⊙ that is inferred from the bolometric luminosity, then it could be that the star was completely destroyed on the second passage, and we would expect the accretion rate to monotonically decline with time. On the other hand, if the bolometric inference significantly overestimates the energy radiated and the mass accreted is closer to ∼0.01 M_⊙, it is likely that the star survived and will return to cause another dimming and future flare. Future observations will show if the star survived the second encounter to generate a third flare.

From the eROSITA nondetections between days ∼600 and ∼1200 (although note that the last eROSITA upper limit is only a factor of ∼2–3 below the observed flux level), the fallback time of the material tidally stripped from the star during its second pericenter passage is, from Figure 1, ∼600 days (note that ordinarily, this timescale is virtually impossible to constrain from observations of single TDEs, but it was possible because of fortuitous eROSITA data points). As noted above, the canonical timescale for a TDE between a solar-like star and a 10^7.7 M_⊙ SMBH is T_acc ≃ 0.8 yr ≃ 300 days (see Equation (2)), which is a factor of ∼2 shorter than the observed fallback time.

However, the return time of the most bound debris from a partial TDE can be significantly longer than the canonical value because of the gravitational influence of the surviving core, which is obvious from the fact that the fallback time becomes infinitely long in the limit that no mass is lost. From Figure 4 of Nixon et al. (2021), the return time of the most bound debris from a 1 M_⊙ zero-age main-sequence star increases by a factor of ∼ 2–3 in going from β ≃ 2 (where the disruption is full) to β ≃ 1, where the star loses ∼10% of its mass (see Figure 4 of Guillochon & Ramirez-Ruiz 2013 and Figure 5 of Mainetti et al. 2017). Nixon et al. (2021) showed (top panel of their Figure 2) that if the peak in the return time is extended to ∼0.2 yr, implying a return time of ∼0.2 yr × 10^1.7/2 ∼ 520 days for a 10^7.7 M_⊙ SMBH (comparable to what is observed for AT 2018fyk), we would need β ≲ 0.7 if the star is somewhat evolved and conceivably smaller if the star is near zero-age (Figure 1 of the same paper).

For β ≃ 0.7, Figure 4 of Nixon et al. (2021) predicts a peak luminosity of ∼4 × 10⁴³ erg s⁻¹ (adopting a radiative efficiency η = 0.1) for a 10^7.7 M_⊙ SMBH, which is slightly less than but still in rough agreement with the X-ray luminosity in the top left panel of Figure 1. At this value of β, the amount of mass lost from the star is also predicted from Newtonian simulations to be ≲0.01 M_⊙ (e.g., Guillochon & Ramirez-Ruiz 2013; Law-Smith et al. 2020), which is in tension with the estimates from the bolometric luminosity that give a value that is closer to ∼0.1 M_⊙. Nonetheless, as we noted, the bolometric luminosity (and thus the total energy radiated) is uncertain for this system, as it is based on a classic accretion disk model where the bulk of the energy is emitted at wavelengths for which we have no data, and the picture outlined here and shown in Figure 4 is clearly quite distinct from a standard disk; the total energy radiated could thus be a factor of ∼5 smaller than the value used to infer the estimate of 0.05 M_⊙ accreted (see the discussion at the beginning of this section). Furthermore, general relativistic simulations indicate that more mass is lost from the star for the same β as compared to Newtonian estimates; for example, Figure 3 of Gafton et al. (2015) shows that a β = 0.7 encounter between a solar-like 5/3 polytrope and a 4 × 10⁷ M_⊙ SMBH—for which the pericenter distance is ∼5.8 GM/c²—strips ∼70% of the stellar mass, while a 10⁶ M_⊙ SMBH removes only ∼25% for the same β. Figures 8–12 of Gafton & Rosswog (2019) show, nonetheless, that the timescales of the TDE remain similar. For our case, in which a Sun-like star is disrupted by a 10^7.7 M_⊙ SMBH, the tidal radius is r_t ≃ 3.5GM/c² and thus highly relativistic. Hence, even for β ≲ 0.7, we would expect a larger fraction of the mass to be lost than would be predicted in the Newtonian limit. Thus, while more detailed modeling is required to more accurately constrain the properties of, e.g., the disrupted star, we find that the overall duration and energetics of the flare are consistent with the partial disruption of a near-solar star. On the other hand, increasing the mass and size of the star would increase the timescale, luminosity, and accreted mass and thus reduce these tensions, but the small separation of the binary restricts the size of the star to ≲1–2 R_⊙ to avoid a common envelope phase (see also the last paragraph of this section).

Assuming that the first detection was approximately coincident with the time of the initial outburst, which is consistent with the lack of optical variability (e.g., from the pre-peak ASAS-SN light curve), we infer that the orbital period of the star is ∼1200 days, or ∼3.3 yr (i.e., the star's first pericenter passage was at day ∼−600 relative to discovery). We therefore predict that—if the star was not destroyed on its second pericenter passage—the source will abruptly decline in luminosity again around day 1800 (2023 August) before flaring for a third time (presuming the star is not destroyed on its third pericenter passage) around ¹⁵ day ∼2400 (2025 March).

Finally, if the orbital period of the captured star is ≈1200 days, then Equation (1) with M_⋆ = 1 M_⊙ and M/M_⋆ = 10^7.7 suggests that the separation of the initial binary—which was ripped apart to yield the captured star—had a separation of ∼0.012 au. As noted above, the M–σ relationship with a black hole mass of 10^7.7 M_⊙ implies that binaries must have a separation of less than ∼0.01 au to survive; hence, this binary separation is consistent with the high velocity dispersion in the nucleus of the galaxy. The distributions of observed binaries that are near solar are roughly uniform in semimajor axis or, for higher-mass stars, in ${Log}(a)$ (Opik's law) and thus peaked toward small separations (Offner et al. 2022). Since the hardening rate is roughly constant once the binary has reached a hardened separation (Quinlan 1996), from a probabilistic standpoint, we would also expect those with the widest (but hardened) initial separations to survive long enough to be fed into the galactic nucleus and tidally destroyed.

From the timescales and energetics arguments above (see the discussion around Equations (2) and (3)), the captured star that is repeatedly partially disrupted likely must be near solar in terms of its mass and size. With a separation a ≲ 0.01 au ∼ 2 R_⊙, the companion object—which was ejected during the separation of the binary (see Figures 4(a) and (b))—is therefore likely required to be a compact object to avoid being in a common envelope phase (as also argued in Cufari et al. 2022b in the context of the event ASASSN-14ko). If the companion was a white dwarf, which is most likely from a statistical standpoint, then the small binary separation appears consistent with the substantial population of detached white dwarf–main-sequence binaries with semimajor axes ∼1 R_⊙ (likely as a consequence of a previous common envelope phase; e.g., Willems & Kolb 2004; Parsons et al. 2015; Mu et al. 2021; Hernandez et al. 2021; Zheng et al. 2022). Thus, in addition to being required from a survivability standpoint and to reproduce the orbital period of the captured star, the small separation of the binary is consistent if the companion is a white dwarf.

After ∼600 days of quiescence, the TDE AT 2018fyk showed an anomalous rebrightening in both the UV and X-ray bands to luminosities to within a factor 10 of their peak values, a behavior that is unprecedented in observations of TDEs. The model we propose to explain this behavior is that the initial flare was caused by the partial disruption of a star that was part of a binary system. The partially disrupted star was captured onto a relatively tight orbit through the destruction of the binary (i.e., Hills capture), thus generating a repeating, partial TDE (as well as a high-velocity star flung out from the system) and the late-time flare. This model is not only consistent with the observations but also predicts that (1) the fallback time of the tidally stripped debris is ∼600 days (a timescale that is, we note, ordinarily very hard to constrain from observations of full TDEs), (2) the orbital time of the captured star is ∼1200 days, and (3) the source should once again dim at day ∼1800 (when the core is expected to return again) and brighten a third time at day ∼2400 if the star was not completely destroyed on its second pericenter passage; on the other hand, if it was completely destroyed, we would expect—as it is then an ordinary TDE—a roughly power-law decay in its luminosity (although, if the star is on a bound orbit, it may exhibit a double-peaked light curve, depending on the eccentricity; Cufari et al. 2022a).

We briefly remark that qualitatively similar behavior, including a late-time rebrightening in the background X-ray emission to ∼60% of its peak magnitude around day ∼3600 postdiscovery, has recently been observed in a source exhibiting quasiperiodic X-ray eruptions (QPEs; Miniutti et al. 2022). To explain their properties, QPEs have also been hypothesized to be the result of repeated tidal stripping, particularly of white dwarfs by low-mass SMBHs (e.g., Arcodia et al. 2021; King 2020; Miniutti et al. 2022). Their host galaxies share several peculiar properties with those of TDEs, including low-mass black holes and a preference for poststarburst galaxies (Wevers et al. 2022).

We finish by highlighting the importance of X-ray and UV monitoring observations of TDEs at late times. Almost all TDEs identified so far lack long-term (yearslong) follow-up. This leaves significant uncertainty as to whether similar behavior has occurred in other TDEs. For example, van Velzen et al. (2019) reported a deep UV upper limit for the source SDSS-TDE1, but no other meaningful constraints exist in the ∼6 yr prior to that observation. Similarly, the majority of TDEs have either no or extremely sparse UV and X-ray constraints at late times. One exception to this is the recently reported observations of AT 2021ehb, a TDE that similarly shows accretion state transitions at late times (Yao et al. 2022), although a partial TDE scenario is not necessary to explain that behavior. Long-term monitoring observations of TDEs—particularly for those with high-mass SMBHs, where partial TDEs are very likely—may provide more evidence for partial TDEs in the future. Indeed, highly periodic flaring may be among the most unambiguous signatures of a partial TDE in general.

We thank the anonymous referee for thoughtful comments and suggestions that helped to improve the manuscript. We thank the XMM, Swift, and NICER PIs (Norbert Schartel, Bradley Cenko, and Keith Gendreau) and their operations teams for approving and promptly scheduling the requested observations. We warmly thank M.I. Saladino for help in creating Figure 4. E.R.C. thanks Chris Nixon for useful discussions and acknowledges support from the National Science Foundation through grant AST-2006684 and the Oakridge Associated Universities through a Ralph E. Powe Junior Faculty Enhancement Award. Raw optical/UV/X-ray observations are available in the NASA/Swift archive (http://heasarc.nasa.gov/docs/swift/archive; target names: AT 2018fyk, ASASSN-18UL) and XMM-Newton Science Archive (http://nxsa.esac.esa.int; ObsID: 0911790601, 0911791601); NICER data are publicly available through the HEASARC: https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl. Based on observations collected at the European Southern Observatory under ESO programs 0103.D-0440(B) and 0106.21SS.

This work is based on data from eROSITA, the soft X-ray instrument on board SRG, a joint Russian–German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE).

The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA.

The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium.

A.1. Observations and Data Reduction

A.2. Swift XRT and UVOT

A.3. XMM-Newton

A.4. NICER/XTI

NICER is a nonimaging detector with 52 coaligned concentrators that focus X-rays onto silicon drift detectors at their respective foci. It has a field of view of 31 in radius and a nominal bandpass of 0.2–12 keV. But, depending on the source brightness and background, the usable bandpass can vary. NICER's large effective area of >1700 cm² at 1 keV enabled by its 52 focal plane modules (FPMs) and ability to steer rapidly to any part of the sky and monitor sources for extended periods of months and years makes it an excellent telescope for tracking long-term transients like TDEs.

Following the Swift/XRT detection of AT 2018fyk, NICER started a high-cadence monitoring program as part of an approved guest observer program (ID: 5070; PI: Pasham). NICER data are organized in the form of ObsIDs that represent a collection of short exposures or GTIs varying from 100 s to up to 2000 s over the time span of a day. While NICER monitoring of AT 2018fyk continues at the time of writing of this paper, we include all data taken prior to 2022 August 22.

We started our NICER data analysis by downloading the raw, unfiltered data from the HEASARC public archive. These were reduced using the standard reduction procedures of running nicerl2 followed by nimaketime. All of the filter parameters except for overonly_range, underonly_range, and overonly_expr were set to the default values as recommended by the data analysis guide: https://heasarc.gsfc.nasa.gov/lheasoft/ftools/headas/nimaketime.html. The reason for not screening on undershoots and overshoots is to ensure we are not throwing away good data in the name of strict default screening values. Instead, we screen each GTI based on the net, i.e., background-subtracted, 0–0.2, 13–15, and 4–12 keV count rates as recommended by Remillard et al. (2022). After a background spectrum is estimated using the 3c50 model, if the absolute value of the net count rate in 0–0.2 keV is more than 2 cps, the absolute value of the net rate in 13–15 keV is more than 0.05 cps, or the absolute value of the net 4–12 keV is more than 0.5 cps, we mark that GTI as bad and omit it from further analysis (see Pasham et al. 2021 for more details).

To improve statistics, we also extracted 18 time-resolved spectra by combining multiple GTIs. Spectra were binned with the optimal binning scheme of Kaastra & Bleeker (2016). To do this, we used the ftool ftgrouppha with an additional requirement to have a minimum of 20 counts per spectral bin. Object AT 2018fyk was above the background in the 0.3–0.7 keV bandpass. Because of this limited bandpass, we fit each spectrum with a simple power law plus a Gaussian model (tbabs*zashift(pow) + Gaussian in XSPEC) and inferred the best-fit power-law index and absorption-corrected 0.3–10 keV luminosities. The Gaussian component was used to model out the variable-strength background oxygen line at 0.54 keV from the Earth's atmosphere. A summary of the spectral modeling is shown in Table A2.

Table A2. Summary of Time-resolved X-Ray Energy Spectral Modeling of AT 2018fyk

Best-fit Parameters from Fitting Time-resolved 0.3–0.7 keV NICER X-Ray Spectra
Start	End	Exposure	FPMs	Phase	Γ	log(Integ. Lum.)	log(Obs. Lum.)	Count Rate	Gaussian
(MJD)	(MJD)	(ks)				(0.3–10 keV)	(0.3–10.0 keV)	(0.3–10.0 keV)	Norm.
59,682.51	59,689.0	0.87	51	L1	${2.64}_{-0.53}^{+0.6}$	${42.8}_{-0.16}^{+0.21}$	${42.7}_{-0.19}^{+0.14}$	0.0055 ± 0.0024	${1.0}_{-1.0}^{+2.3}$
59,693.7	59,698.27	2.6	43	L2	${4.49}_{-0.23}^{+0.24}$	${42.99}_{-0.02}^{+0.02}$	${42.9}_{-0.01}^{+0.01}$	0.0166 ± 0.0012	${23.1}_{-2.4}^{+1.7}$
59,698.27	59,703.0	5.45	47	L3	${3.47}_{-0.21}^{+0.22}$	${42.84}_{-0.03}^{+0.03}$	${42.74}_{-0.02}^{+0.03}$	0.0093 ± 0.0006	${7.7}_{-1.2}^{+1.1}$
59,703.0	59,708.0	8.68	46	L4	${3.26}_{-0.18}^{+0.18}$	${42.88}_{-0.03}^{+0.03}$	${42.77}_{-0.03}^{+0.03}$	0.0095 ± 0.0003	${5.3}_{-0.9}^{+1.0}$
59,708.0	59,718.0	9.1	49	L5	${2.77}_{-0.17}^{+0.17}$	${42.82}_{-0.04}^{+0.05}$	${42.75}_{-0.04}^{+0.04}$	0.0073 ± 0.0003	${5.5}_{-0.8}^{+0.8}$
59,718.0	59,723.0	3.46	51	L6	${2.53}_{-0.35}^{+0.36}$	${42.79}_{-0.11}^{+0.14}$	${42.72}_{-0.12}^{+0.1}$	0.0057 ± 0.0007	${3.9}_{-1.2}^{+1.2}$
59,723.0	59,728.0	5.12	50	L7	${2.39}_{-0.22}^{+0.22}$	${42.84}_{-0.08}^{+0.09}$	${42.77}_{-0.06}^{+0.06}$	0.0055 ± 0.0005	${2.6}_{-0.8}^{+0.8}$
59,728.0	59,733.0	3.69	50	L8	${1.81}_{-0.33}^{+0.34}$	${43.12}_{-0.19}^{+0.23}$	${43.06}_{-0.24}^{+0.27}$	0.0052 ± 0.0007	${0.7}_{-0.7}^{+0.7}$
59,733.0	59,738.0	3.63	52	L9	${2.48}_{-0.33}^{+0.36}$	${42.82}_{-0.11}^{+0.14}$	${42.73}_{-0.11}^{+0.15}$	0.0056 ± 0.0007	${1.9}_{-1.1}^{+1.1}$
59,738.0	59,743.0	3.23	52	L10	${2.33}_{-0.29}^{+0.32}$	${42.91}_{-0.11}^{+0.13}$	${42.83}_{-0.16}^{+0.14}$	0.0055 ± 0.0008	${0.0}_{-0.0}^{+0.3}$
59,743.0	59,748.0	4.76	52	L11	${2.5}_{-0.36}^{+0.38}$	${42.7}_{-0.11}^{+0.15}$	${42.61}_{-0.12}^{+0.16}$	0.0041 ± 0.0005	${0.9}_{-0.9}^{+0.9}$
59,748.0	59,758.0	10.3	52	L12	${2.84}_{-0.22}^{+0.21}$	${42.55}_{-0.04}^{+0.06}$	${42.46}_{-0.05}^{+0.06}$	0.0038 ± 0.0002	${1.9}_{-0.5}^{+0.5}$
59,758.0	59,768.0	2.53	51	L13	${3.44}_{-0.6}^{+0.69}$	${42.59}_{-0.06}^{+0.1}$	${42.51}_{-0.06}^{+0.07}$	0.0056 ± 0.0009	${5.8}_{-1.5}^{+1.5}$
59,768.0	59,778.0	4.14	51	L14	${2.68}_{-0.34}^{+0.35}$	${42.62}_{-0.09}^{+0.11}$	${42.57}_{-0.07}^{+0.09}$	0.0048 ± 0.0006	${5.8}_{-1.0}^{+0.9}$
59,778.0	59,783.0	6.03	52	L15	${2.36}_{-0.21}^{+0.21}$	${42.85}_{-0.07}^{+0.09}$	${42.77}_{-0.07}^{+0.11}$	0.0055 ± 0.0004	${2.1}_{-0.4}^{+0.6}$
59,783.0	59,788.0	5.06	52	L16	${2.38}_{-0.21}^{+0.22}$	${42.92}_{-0.08}^{+0.09}$	${42.83}_{-0.07}^{+0.07}$	0.006 ± 0.0005	${0.6}_{-0.6}^{+0.4}$
59,788.0	59,798.0	7.05	52	L17	${2.37}_{-0.16}^{+0.16}$	${42.96}_{-0.05}^{+0.06}$	${42.87}_{-0.05}^{+0.04}$	0.0066 ± 0.0004	${0.0}_{-0.0}^{+0.0}$
59,798.0	59,820.0	8.22	52	L18	${2.82}_{-0.2}^{+0.2}$	${42.74}_{-0.04}^{+0.05}$	${42.63}_{-0.06}^{+0.06}$	0.0052 ± 0.0003	${0.7}_{-0.5}^{+0.7}$

Note. Here 0.3–0.7 keV NICER spectra are fit with the tbabs*zashift(clumin*pow) + Gaussian model using XSPEC (Arnaud 1996). Start and End represent the start and end times (in units of MJD) of the interval used to extract a combined NICER spectrum. Exposure is the accumulated exposure time during this time interval. The FPMs are the total number of active detectors minus the "hot" detectors. Phase is the name used to identify the epoch, Γ is the photon index of the power-law component, log(Integ. Lum.) is the logarithm of the integrated absorption-corrected power-law luminosity in 0.3–10 keV in units of erg s⁻¹, and log(Obs. Lum.) is the logarithm of the observed, extrapolated 0.3–10.0 keV luminosity in units of erg s⁻¹. Count rate is the background-subtracted NICER count rate in 0.3–0.7 keV in units of counts per second per 50 FPMs. All error bars represent 1σ uncertainties. The total best-fit χ²/degrees of freedom over all spectra is 76.5/65.

Download table as:  ASCIITypeset image

A.5. SRG/eROSITA

A.6. MUSE

On 2019 June 10 (MJD 58,644), AT 2018fyk was observed by MUSE (Bacon et al. 2010) as part of the All-weather MUse Supernova Integral field Nearby Galaxies (AMUSING) survey, ESO ID: 0103.D-0440(B). At this epoch, the transient was already quiescent at optical wavelengths, and the host galaxy emission completely dominated the data. The data cube was analyzed as part of the AMUSING++ Nearby Galaxy Compilation (López-Cobá et al. 2020). However, the authors did not include it in the final sample of the paper due to the lack of strong emission lines, which were the main subject of their study. Nevertheless, we obtained the final products of their stellar population and emission line fitting analyses (López-Cobá, private communication).

A detailed description is presented in López-Cobá et al. (2020). In summary, the following procedure was adopted. First, the raw data cubes were reduced with REFLEX (Freudling et al. 2013) using version 0.18.5 of the MUSE pipeline. Next, the emission lines and stellar population content were analyzed using the PIPE3D pipeline (Sánchez et al. 2016a), a fitting routine adapted to analyze IFS data using the package FIT3D (Sánchez et al. 2016b). The procedure starts by performing a spatial binning on the continuum (V band) to increase the signal-to-noise ratio in each spectrum of the data cube. The stellar population model was derived by performing stellar population synthesis; the PIPE3D implementation adopts the GSD156 stellar library, which comprises 39 ages and four metallicities, extensively described in Cid Fernandes et al. (2013). Then, a model of the stellar continuum in each spaxel was recovered by rescaling the model within each spatial bin to the continuum flux intensity in the corresponding spaxel. The best model for the continuum was then subtracted to create a pure gas data cube. A set of 30 emission lines within the MUSE wavelength range were fit spaxel by spaxel for the pure gas cube by performing a nonparametric method based on a moment analysis. The data products of the pipeline are a set of bidimensional maps of the considered parameters with their corresponding errors.

In Figure A1, we show the sample of these maps with the main parameters of interest for this study. The galaxy shows a centrally concentrated structure, like most TDE hosts (Hammerstein et al. 2022); a very old stellar population (mean age ≥10^9.7 yr); a lack of dust (A_V* ≤ 0.05 in all spaxels); and very faint emission lines (mean EW Hα < 1 Å), without any apparent increase toward the central spaxels.

Zoom In Zoom Out Reset image size

A.7. X-shooter

The host galaxy was observed in long-slit mode with the X-shooter instrument on the Very Large Telescope Unit Telescope 3 on 2020 October 16 (MJD 59,138.08). Slit widths of 10, 09, and 09 were deployed for the UVB, VIS, and NIR arms, respectively, for a total exposure time of 1300 s. The average seeing of 07 during the observations results in a seeing-limited spectral resolution of R = 7700 (UVB), 12,700 (VIS), and 8000 (NIR), equivalent to an FWHM spectral resolution of 40 (at 4000 Å) and 25 (at Hα) km s⁻¹. The data were taken in on-slit nodding mode, but to increase the signal-to-noise ratio of the UVB and VIS arms, we reduce these data using the X-shooter pipeline with recipes designed for stare mode observations.

We modeled the stellar continuum of the X-shooter spectrum with a wavelength range of 4000–7000 Å in the rest frame using the penalized pixel fitting (Cappellari 2017) routine. We masked some emission and absorption lines that are usually significant in galaxy spectra, since they may affect the best fits of stellar continuum models, e.g., Hδ, Hγ, Hβ, Hα, N ii 4640, He ii 4686, [O iii] 4959, 5007, He i 5875, [O i] 6300, [N ii] 6548, 6584, and [S ii] 6717, 6731. We used MILES single stellar population (SSP) models (Vazdekis et al. 2010) as the stellar templates and adopted the SSP model spectra. Given that the initial resolution of the X-shooter spectrum is R ∼ 10,000, much higher than that of the MILES spectra (R ∼ 2000), we convolved the X-shooter spectrum to reduce its resolution to R ∼ 2000. Except for the stellar template, a polynomial with degree = 4 was added to avoid mismatches between galaxy spectra and stellar templates. The residuals were obtained after subtracting the best-fit stellar continuum model. Residual flux errors are the same as the original flux errors. Line ratios and EWs were measured on the residual spectra, and uncertainties were determined by taking into account the flux uncertainties. The resampled galaxy spectrum overlaid with the fit and residuals after template subtraction are shown in Figure A2.

Zoom In Zoom Out Reset image size