Light-curve Evolution of the Nearest Tidal Disruption Event: A Late-time, Radio-only Flare

We present a number of VLA observations of NGC 4845 (Table 1), collected here for the first time. Some are new data that we have obtained, some are re-imaged from previous data in the VLA archive, and some are reduced images from previous surveys. Together with previously known values, we now have 20 data points from 15 separate dates. This fairly extensive light curve is one of the most extensive and long-lasting of any TDE. The only one with similar numbers of long-term points, extending for several years is that of Sw J1644+57, which has more data points but does not extend for quite as many years (Alexander et al. 2020).

In Table 1, the first two rows correspond to historical observations of which we include here as a check on the pre-TDE flux level. Following these in the table are the first four epochs of data for NGC 4845 corresponding to the actual TDE event (including epoch T1). These observations were taken in 2011 and 2012 as a part of the CHANG-ES project (CHANG-ES, observation ID: 10C-119, PI: J. Irwin; Irwin et al. 2012). These include L-band and C-band (centering at 1.6 and 6.0 GHz, respectively) data in three array configurations, B, C, and D. The data reductions are described in Irwin et al. (2015). In 2015, the source was observed in the L and C bands as part of VLA project 15A–357 (PI: E. Meyer). These data were previously described and published in Perlman et al. (2017).

The remaining entries in the table correspond to additional subsequent observations of NGC 4845, observed in multiple array configurations and bands (L, S, C, and X). The separate observations all share a number of steps in reduction procedure, which we describe here. We have attempted to measure each new data point in a consistent manner; however, an unusual number of coauthors were involved in individually reducing different data sets. Minor differences in, e.g., settings during deconvolution imaging, and the coauthor responsible for the analysis are noted below and are not expected to impact any of the results.

2.1.1. General VLA Data Processing

2.1.2. 2015 High-resolution Observations (ID: 15A-400)

2.1.3. 2016 Swift-concurrent Observations (ID: 16A-420)

2.1.4. 2019 June L-band Observation (ID: 19A-425)

2.1.5. Archival Data from Project AB0879

The X-ray data comes from the Neil Gehrels Swift Observatory (Swift hereafter) and XMM-Newton telescopes with detailed information in Table 2. Some early XMM-Newton and Swift observations were published in Nikolajuk & Walter (2013) but are reanalyzed here. The 2016 March and 2016 May observations were meant to coincide closely with two of our VLA observations. Before calibration, we combined any two Swift/XRT observations that are only minutes or hours apart because sudden gamma-ray bursts can divide a pre-planned Swift observation into two. For XMM, we only keep data from the European Photon Imaging Camera (EPIC) because it is the only instrument on board XMM-Newton covering soft X-rays with 2D imaging. Then, we use the xrtpipeline command for Swift and the emchain/epchain command in SAS (Science Analysis System), respectively, for MOS/PN CCD data of XMM/EPIC to conduct necessary calibrations including calibrating bias, removing bad pixels, and removing flares.

With Swift and XMM-Newton event files that were cleaned of background flares, we choose source and background radii according to the central PSFs of the telescopes, and extracted source and background counts using typical procedures. ¹⁰ We generated Swift response files by selecting the correct RMFs in its calibration database based on their observation times and calculate ARFs by xrtmkarf command. For XMM-Newton, the rmfgen and arfgen commands in SAS are used to produce all response files. We use absorbed power law plus absorbed Astrophysical Plasma Emission Code (APEC) models, according to the expression TBabs_G (APEC + TBabs(power law)), where TBabs_G is the Galactic absorption, ¹¹ for the two bright-phase observations in Table 3, accounting for both the nucleus and the surrounding hot gas. We considered the source to be extincted by the Galactic column of ${N}_{{\rm{H}}}^{\mathrm{Gal}}=1.67\times {10}^{20}\,{\mathrm{cm}}^{-2}$ (Dickey & Lockman 1990; Stark et al. 1992), and absorption from the host galaxy was allowed to vary freely. The default versions of the solar abundances were assumed for all APEC models. All spectra are regrouped to contain at least 30 counts per bin, which should be sufficient for using a simple χ² minimization. Although the C-statistic is statistically more accurate, we have not used it because its use would involve complication and uncertainty in having to model the background spectral contribution.

Table 3. X-Ray Source and Background Details

Observation(s)	PSF	Center	Source	Background	Source	Background	0.5–7.0 keV Flux
			Radius	Radius	Counts	Counts	(erg cm−2 s−1)
00031911001+00031911002					2026	283	2.60 ± 0.06 × 10−11
00031911003		decl. =			9	4	${2.4}_{-0.9}^{+1.3}\times {10}^{-13}$
00031911004+00031911005	18''	+1d34m33s,	30''	30''–60''	6	1	< 3.75 × 10−13
00031911006+00031911007					6	2	< 3.75 × 10−13
00031911008					6	3	< 7.35 × 10−13
00031911004-00031911008		R.A. =			18	6	${1.8}_{-0.5}^{+0.7}\times {10}^{-13}$
		12h58m13			MOS1: 14091	MOS 1: 9439
0658400601	6''		10''	10''–30''	MOS2:14930	MOS2: 9902	${2.87}_{-0.03}^{+0.08}\times {10}^{-11}$
					PN: 28475	PN: 16216

Download table as:  ASCIITypeset image

Consistent with Irwin et al. (2015), we define 2011 December 30 as time T1 (t = 0 on the plot), and plot the other data given here and in Perlman et al. (2017) in Figure 1. As can be seen in Figure 1, the nuclear radio flux declined monotonically from the radio peak until approximately 2015. In 2015, we saw a plateau in the radio flux, with 20% temporal variability in the X band between 2015 July 6 and 2015 July 14 (the apparent C-band variability seen in June–July 2015 is more likely a result of the two slightly different radio frequencies). Then, in 2016 (day 1538), a sharp increase began in the radio flux. By day 1604, the flux in the C, S, and L bands had increased by a factor of 2.5–3. This defines the epoch of our new "flare." Unfortunately there is almost no information on the radio variability of the nuclear source of NGC 4845 previous to the TDE, although we note that the two L-band fluxes in 1995 and 1998 are consistent with one another at the 2σ level.

Zoom In Zoom Out Reset image size

Power-law fits to both XMM-Newton and Swift data (Nikolajuk & Walter 2013) found a steep (Γ ≈ 2.35) X-ray spectrum with a column density N_H ≈ 7 × 10²² cm⁻². Our reanalysis confirms that result, but we find that an APEC model is a better fit to the data, with the main improvement being an improved fitting of the soft excess. Figure 2 shows the XMM-Newton and Swift data as well as the best-fit model. A simultaneous fit of all three 2011 data sets yielded a power law Γ = 2.316 ± 0.042, ${kT}={0.238}_{-0.052}^{+0.149}$ keV for the thermal component, and absorbing column N_H = 8.8 ± 0.2 × 10²¹ cm⁻². The decrease in the ${\chi }_{\nu }^{2}$ was highly significant—from 0.973 to 0.410 for the Swift data and from 1.291 to 0.710 for the XMM-Newton observations. The soft X-ray component is also noted by Nikolajuk & Walter (2013), but those authors suggest it is unrelated to the nuclear emission associated with the TDE based on a fit to a blackbody component with temperature kT = 0.33 ± 0.04 keV that is extincted only by the Milky Way's column, which they claim (their Section 3.2) is consistent with diffuse emission from galaxies (Bogdán & Gilfanov 2011; Jia et al. 2012). They do not provide any statistics on that fit; however, it would be one additional degree of freedom beyond our model, which as noted shows a large improvement over a power law. We suggest that the evidence for a smaller column for the thermal component is not persuasive, and moreover we point out that diffuse emission from a galaxy at only 17 Mpc distance would easily be resolved with XMM-Newton. Moreover, the 0.5–7.0 keV flux from the thermal component, 4.3 × 10⁻¹³ erg cm⁻² s⁻¹, differs by more than 3σ from the value observed in 2016–2017 (Table 3). We therefore assess that the soft component is more likely to be nuclear in origin.

Zoom In Zoom Out Reset image size

Table 3 presents the X-ray source and background region information, as well as counts and fluxes associated with each observation. The fluxes for the 2011 XMM-Newton and Swift observations come from the procedure discussed above. For the 2012–2017 Swift observations, the count rates are far too small to allow for any spectral fitting. We thus use the above best-fit spectral model to convert the count rates to the corresponding energy fluxes for comparison. For the two bright-phase (2011) observations plus the 2012 Swift observation, we were able to extract background information from the individual observations. However, for the 2016–2017 Swift observations, where this was not possible with the individual observations, the background levels were taken from the combined 2016–2017 observations. Figure 3 shows the 0.5–7 keV X-ray fluxes in the 0.5–7 keV band (with no absorption correction, thus allowing us to include the INTEGRAL observation) as a function of time. In addition, we show in Figure 3 the X-ray fluxes from Nikolajuk & Walter (2013), extrapolated down to the 0.5–7 keV band using the joint INTEGRAL + XMM-Newton best-fit model in that paper, namely a hard X-ray spectral index of Γ = 2.22 ± 0.03, slightly different from what we find here.

Zoom In Zoom Out Reset image size

Three things are clear from Figure 3. First, IGRJ12580+0134 represented a much larger departure from the "quiescent" level of X-ray flux than was previously realized. Based on the INTEGRAL light curve of Nikolajuk & Walter (2013), one could only say that the event was an increase in X-ray flux of at most a factor 10. Figure 3 makes it clear that IGRJ12580+0134 represented a much larger increase in NGC 4845's X-ray flux—more than two orders of magnitude. Second, as also found by Nikolajuk & Walter (2013), there is no evidence for spectral curvature between 7 keV and the 17.3–80 keV band of INTEGRAL. Not only are the power-law components of the INTEGRAL spectral fit and our XMM-Newton + Swift fit roughly consistent within the errors (they differ by less than 2σ), but in addition, the 2011 Swift/XRT and XMM-Newton fluxes fit smoothly on the curve described by the INTEGRAL points, without any evidence of error due to extrapolation. Finally, there is no evidence that the 2016 radio flare described in Section 3.1 (Figure 1) was accompanied by an X-ray flare. While the two 2016 data points are not detections (and for them as well as the 2017 Swift observation, we show 2σ upper limits), we see that when the three 2016–2017 observations are added together, a significant detection is seen, and the detected flux is within 1σ of the 2012 data point. If indeed there had been an X-ray component to the 2016 flare, the two 2016 Swift observations would dominate that point, and we see no evidence of that.

Recently, Auchettl et al. (2017) have questioned the classification of IGR J12580+0134 as a TDE, based on the fact that the host galaxy, NGC 4845, is a Sey 2/LINER (Ho et al. 1995). They claimed that IGRJ12580+0134 could be explained by active galactic nucleus (AGN) processes. ¹² The nucleus likely has some X-ray variability, so we definitely agree with the general cautions regarding AGN variability expressed by those authors. But we contend that IGR J12580+0134 is still likely to be a TDE, for several reasons.

First, the large X-ray and radio flux increases observed (now seen to be at least two orders of magnitude in X-rays, much larger than could be justified based on the INTEGRAL data in Nikolajuk & Walter 2013) are well beyond what is usual for an AGN, with less than 1% of AGNs having flares of that magnitude (Auchettl et al. 2017). The flares that tend to be similarly extreme are in blazars (e.g., the 1998 flare of PKS 2005–489; Perlman et al. 1999) and narrow-line Sey 1 galaxies (e.g., 1ES 1927+654; Trakhtenbrot et al. 2019, 2019; Boller et al. 2003), and NGC 4845 displays no properties in common with such sources. Also, the nuclear luminosity in quiescence would be ≤ 6.2 × 10³⁹ erg s⁻¹, as observed in the 2016–2017 Swift observations (of which some could be due to stellar or gas processes), much more in line with a quiescent galaxy nucleus or low-luminosity AGNs (e.g., Ptak et al. 1999).

Second, Nikolajuk & Walter (2013) pointed out that the hard X-ray (17.3–80 keV) flux evolution after the peak was consistent with a t^−5/3 law based upon INTEGRAL data. Arguments consistent with Kepler's second law predict that the mass infall rate should approach t^−5/3 at late times for all types of disrupted stars (Phinney 1989; Lodato et al. 2009, note that the evolution close to the peak can depend on the stellar structure). However, if a thermally emitting disk forms, the evolution can be very different, ∝ t^−5/12 if emission sits in the Rayleigh–Jeans tail of a thermal component. Indeed, at much later times, a t^−5/3 evolution is not fully consistent with the X-ray light curve. (Figure 3). Alternative fits are also possible, including an exponential fit and/or power law at earlier time with an exponential decline at later times. Exponential behavior in the X-ray/optical has been observed in some TDEs, including, e.g., iPTF16fnl (Onori et al. 2019), another faint TDE. Such behavior is not unexpected in TDEs. As Lodato & Rossi (2011) point out, the t^−5/3 behavior is expected to continue for only a few months before steepening to become essentially exponential in the X-ray band, where the X-ray emission at later times would be dominated by an outflow.

Figure 3 makes clear that the 2012 Swift flux fits neatly on that same power law as the early INTEGRAL data. Thus the t^−5/3 behavior continued for nearly 500 days. In addition, the quickness of the decline is not consistent with an extinguishing event, which should take place over much longer timescales ( ∼ 10⁴ yr; Shen 2021).

The X-ray spectrum of IGRJ12580+0134 also has several properties that are more in common with TDEs. Its spectral index, Γ = 2.316 ± 0.042, is rare among AGNs (which are dominated by objects with Γ ≈ 1.8; Auchettl et al. 2018; Ricci et al. 2017) but common among TDEs (Auchettl et al. 2018). While the X-ray power-law emission of IGRJ12580+0134 continues to 80 keV (Nikolajuk & Walter 2013), unlike any other TDE, the power-law component of at least one other jetted TDE, Swift J1644+57, continues up to at least ∼30 keV (Bloom et al. 2011; Kara et al. 2016), as evidenced by its Swift/BAT detection. X-ray spectral indices over two are seen in high-frequency spectral peak BL Lac objects (e.g., Perlman et al. 2005; Paliya et al. 2019), where the X-ray spectrum is dominated by the tail of synchrotron emission, but IGRJ12580+0134 shows no evidence for beaming. A thermal excess is almost universally found in TDEs, but the temperature, ${kT}={0.24}_{-0.05}^{+0.15}$ keV, is somewhat higher than seen in most TDEs (typically ∼50 eV, although comparably high temperatures are seen in ∼10% of cases; Auchettl et al. 2018) and AGNs (typically ∼100 eV but a with a handful over 200 eV; Ricci et al. 2017). As discussed in Section 3.2, this thermal component is very likely to be of nuclear origin. Finally, the absorbing column is reasonable both for AGNs and TDEs (Auchettl et al. 2018) and is similar to that seen in the jetted TDE Swift J1644+57.

This re-assessment confirms that IGR J12580+0134 was a TDE in at most a low-luminosity AGN. Two very recent reviews make the point that TDEs and TDE-like events should be common in AGNs (Gezari 2021; Saxton et al. 2021). Saxton et al. (2021) specifically discusses IGR J12580+0134 and mentions this possibility. As Saxton et al. (2021) note, at least a few other possible TDEs fall into this category. One particularly interesting case could even be a recent event in OJ287 where the X-ray light curve shows a clear, t^−5/3 decay after the event (Huang et al. 2021). A much more doubtful case is Arp 299, where Mattila et al. (2018) suggested the presence of a TDE, attributing the lack of signatures in other bands to a large obscuration. However, that event was more than an order of magnitude smaller in amplitude than IGR J12580 and lacked the characteristic light-curve variability. A recent paper by Zabludoff et al. (2021) discussed a more sophisticated way to discriminate TDEs from interloper events, which will be critical for future surveys such as the Vera G. Rubin Large Synoptic Survey Telescope (LSST).

The clarification of IGR J12580+0134's nature as a TDE also allows us to discuss other aspects of its X-ray spectral behavior. The flow of matter into the near vicinity of an active black hole might make TDEs more common in AGNs than in quiescent galaxies, although this hypothesis needs to be tested by LSST observations. Stolc & Karas (2019) have pointed out that in such systems one might expect an evolving, relativistic Fe Kα line. Our X-ray spectral data do not show such a line, but the signal to noise at 6.4 keV is low.

As discussed in Section 3, the radio flux of IGR12580+0134 has continued to evolve, as shown in Figure 1. Up to and including 2015, VLA data show a descent to a plateau in the radio flux. During that plateau, ∼20% variability was seen in the X band, which is common in radio-loud AGNs and also more quiescent galaxies. For an example of the latter, see the light curve of Sgr A* shown in Subroweit et al. (2017). Very little data exist on the long-term radio light curves of TDE (Alexander et al. 2020), so it is difficult to comment on the frequency of such small-scale radio variability in that population. The plateau in the radio emission seen in 2015 is reminiscent of events seen in the X-ray light curve of ASASSN-14li, 1–2 yr after the TDE (Bright et al. 2018). In that object, Pasham & van Velzen (2018) have detected a time-lag between the soft X-ray and radio emission, and interpreted it as linear disk-jet coupling. Unfortunately, we do not have sufficient radio or X-ray data in IGR J12580+0134 to perform such a test. In addition, we have found a major radio flare in 2016. Our sole radio point after 2016, taken in 2019, shows a large decline in radio flux; however, due to the lack of monitoring between 2016–2019, we cannot speculate on the speed of the decline.

It is interesting to look at the evolution of the radio spectrum during these events, shown in Figure 4. In that plot, two of the spectra are truly simultaneous (March 2016 and May 2016), while the others simply show all data during the indicated periods. In 2015, we saw a C–L band spectral index of α_CL = 0.89, while the X–C band spectral index in July 2015 was α_XC = 1.67 (the slightly different 2015 June and July C–band flux points are clearly a result of observing at two different central frequencies, as Figure 4 makes clear). The spectral index at lower frequencies (L and S bands) was at its flattest during the flare's brightest stage, in May 2016, when a value of α_SL = 0.64 (F_ν ∝ ν^−α ) was seen, as compared to a value of α_SL = 1.51 in 2016 March when the flare was in an earlier stage. However, in 2016 May, the higher-frequency (C to S band) spectral index was significantly steeper (α_CS = 1.45) than it was in 2016 March (α_CS = 1.15). Thus the radio spectrum in 2016 March steepened moderately as one went to lower frequencies, while in 2016 May, it flattened drastically as one went to lower frequencies. The new flare has not been monitored sufficiently to show the decline behavior between 2016–2019. However, in 2019, the radio spectrum between the S and L bands was α_SL = 1.25.

Zoom In Zoom Out Reset image size

IGR J12580+0134 joins ASASSN-15oi (Horesh et al. 2021) as one of only two TDEs with late-time radio flares observed. Comparing the two, ASASSN-15oi during the peak of its late-time radio flare was two orders of magnitude fainter in flux than IGRJ12580+0134, but due to its much greater distance (216 Mpc), its peak luminosity was a factor of a few higher, actually higher than its radio peak during the primary TDE event. The increase in radio luminosity observed in ASASSN-15oi was actually larger than that observed in IGRJ12580+0134—a full order of magnitude, as compared to a factor of three here. The two occur at similarly late times after the TDE (4 yr after the TDE in ASASSN-15oi as compared to 5.5 yr after the TDE in IGRJ1258+0134). However, apparently no radio monitoring of ASASSN-15oi exists after the late-time rise in 2019, at least up until the publication this year of Horesh et al. (2021).

There are, to our knowledge, five other TDEs that show late-time X-ray flares. These are PS10adi (Jiang et al. 2017), OGLE16aaa (Kajava et al. 2020), ASASSN-15oi (at the same time as its radio flare; Gezari et al. 2017; Holoien et al. 2018), ASASSN-18jd (Neustadt et al. 2020), and AT2019azh (Liu et al. 2019). All five have been modeled as possibly due to a second fallback event due to circularization of the orbits of accreting material (van Velzen et al. 2019; Gezari 2021) or a sub-relativistic shockwave (Yalinewich et al. 2019; Alexander et al. 2020). As discussed below, the radio-only flare in IGR J12580+0134 cannot be explained in this context, similar to ASASSN-15oi where a similar conclusion was reached (Horesh et al. 2021).

The 2016 radio flare seems more likely to be the result of the jet intercepting a circum-nuclear cloud and accelerating particles to relativistic energies. Following the modeling of an off-axis jet propagating in the CNM employed to explain the multiwavelength evolution of radio emission from IGR J12580+0134 (Lei et al. 2016; Perlman et al. 2017), we fit the fluxes assuming that a denser gas cloud exists in the path of the propagating jet. The jet dynamics is described by a set of hydrodynamical equations, including the evolution of the radius and Lorentz factor of the jet, and the swept mass (Huang et al. 2000; Liu et al. 2020). The magnetic field and accelerated electron population share parts of the thermal energy of the downstream medium of the shock, with fractions _B and _e , respectively. Electrons were accelerated to form a power-law spectrum with an index of − p. Synchrotron radiation was then calculated given the magnetic field and electron distribution at each time.

Figure 5 shows the light curves from the model calculation, compared with the data. The main model parameters include: the total kinetic energy of the ejecta E_k = 2.5 × 10⁵² erg, the initial ejecta Lorentz factor Γ_i = 8.5, _e = 0.27, _B = 0.09, the jet opening angle θ_j = 54, the viewing angle θ_obs = 36°, the electron spectral index p = 2.80, and the background medium density n = 5.2 cm⁻³. A gas cloud with an average density of n_cl = 10.3 cm⁻³ spanning from 0.95 to 0.97 pc from the central black hole, as indicated by the bottom panel of Figure 5, is found to be able to fit the radio flare. The cloud mass swept by the jet (with opening angle of θ_j ) is estimated to be about 1.8 ×10⁻⁴ M_⊙, and the jet velocity is about 0.17c when it hits the cloud. Given the soft spectrum of electrons and the cooling effect, we estimate that the flux in the X-ray band (at 1 keV) during the radio flare time is about 3.4 × 10⁻¹⁷ erg cm⁻² s⁻¹ keV⁻¹, which is several orders of magnitude lower than the upper limits shown in Figure 3.

Zoom In Zoom Out Reset image size

Note that there is degeneracy among the model parameters, and the above parameters are one example of those with reasonable values. Nevertheless, some of these parameters can be constrained by the data. While the CNM density degenerates with the total kinetic energy as well as the energy partition fractions, the density contrast between the gas cloud and the CNM background can be constrained relatively well (about two, although the exact value is affected somehow by the thickness of the shell). The swept mass, derived by the gas density and cloud thickness, is constrained by the declining profile and post-flare fluxes.

This simplified model fits the observational data, including the radio flare, fairly well. The flattening of the radio spectral index at the brightest phase of the late-time radio flare (Figure 4) is consistent with this model, as such an interaction would cause a shock in the jet. The rapid steepening of the radio spectral index in the later 2016 epoch is consistent with a small size for the cloud. The late-time L-band flux in 2019 is a little bit higher than the model expected flux. The X-band fluxes also show a quick decrease that was not observed in other bands. Further refinement of the model may improve the match between the model and data. If this scenario is correct, the radio monitoring of jetted TDEs could allow us to probe the surrounding medium of the galactic nuclei. The density profile of the background CNM can in principle be probed with the light curves (Berger et al. 2012). It is, however, difficult for IGR J12580+0134 given the sparse data points during and immediately after the flare, since, as can be seen, additional observations would have been helpful.

Horesh et al. (2021) reached a similar conclusion regarding the nature of the late-time radio flare of ASASSN-15oi. They modeled it as due to the collision of the radio jet with a circumstellar shell of material, after rejecting fallback events and shocks within a structured jet due to the steepness of the radio rise. We cannot distinguish that possibility from the model we propose in this paper, and in fact such a shell would be consistent with the possibility we model, although resulting in a considerably greater mass for the disturbance.