The rise and fall of the iron-strong nuclear transient PS16dtm

Context. Thanks to the advent of large-scale optical surveys, a diverse set of ﬂares from the nuclear regions of galaxies has recently been discovered. These include the disruption of stars by supermassive black holes at the centers of galaxies – nuclear transients known as tidal disruption events (TDEs). Active galactic nuclei (AGN) can show extreme changes in the brightness and emission line intensities, often referred to as changing-look AGN (CLAGN). Given the physical and observational similarities, the interpretation and distinction of nuclear transients as CLAGN or TDEs remains di ﬃ cult. One of the obstacles of making progress in the ﬁeld is the lack of well-sampled data of


Introduction
In recent years, the detection of a variety of nuclear transients has become attainable due to the advent of large-scale optical robotic surveys such as the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry 2011), the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014), the Pan-STARRS Survey for Transients (PSST; Chambers et al. 2016b), the Gaia Photometric Science Alerts program (Hodgkin et al. 2021), the intermediate Palomar Transient Factory (iPTF; Kulkarni 2013), and, its successor, the Zwicky Transient Facility (ZTF; Bellm et al. 2019).It has long been known that galaxies with an active galactic nucleus (AGN), where matter is accreted onto the central supermassive black hole (SMBH), show variability in their brightness (e.g., Fitch et al. 1967).Apart from regular, low-level stochastic variability, some AGN occasionally show exceptionally large changes in the luminosity, spectral shape, and/or Xray absorption (e.g., Frederick et al. 2019;MacLeod et al. 2019;Sánchez-Sáez et al. 2018;Zhang 2021;Hon et al. 2022;Green et al. 2022;Ren et al. 2022).The most notable of these changes is when a Seyfert 1 type AGN transforms into a Seyfert 2 AGN and vice versa, the so-called changing-look AGN (CLAGN).The physical explanation of this phenomenon is still debated.It has been argued that the variability arises in the intrinsic changes in the accretion state of the SMBH, which include several scenarios (Noda & Done 2018;Sniegowska et al. 2020;Stern et al. 2018;Mehdipour et al. 2022).A number of works have argued that there might be more than one mechanism that explains CLAGN, given the diversity in the timescale of changes, and their amplitude (e.g., Wang et al. 2012;Campana et al. 2015;Komossa et al. 2017;Trakhtenbrot et al. 2019;Śniegowska et al. 2022).
One important class of optical transients occurring in the nuclear regions of the host galaxies are the flares from tidally disrupted stars by SMBHs (Rees 1988;Phinney 1989).They are particularly interesting since there are several suggestions of using tidal disruption events (TDEs) to study SMBH properties such as its mass (Mockler et al. 2019) and perhaps even more elusive, its spin (Leloudas et al. 2016;Gafton & Rosswog 2019).Optically discovered TDEs differ between each other in the spectroscopic features, shape, and timescale of the light curve, and have been detected in both quiescent (e.g., Gezari 2021) and active galaxies (e.g., Wyrzykowski et al. 2017;Freder-

Summary of previous work on PS16dtm
In this section, we summarize previously published work on PS16dtm. 1 It was discovered by PSST on 12 August 2016 (MJD 57612) (Chambers et al. 2016a), and it was also observed independently by ATLAS and ASSA-SN.The host galaxy of PS16dtm, SDSSJ015804.75-005221.8, is a narrow-line Seyfert 1 (NLSy1) galaxy at z = 0.0804 with a SMBH with ∼ 10 6 M , measured by the velocity dispersion method (Xiao et al. 2011) and from the AGN luminosity and the radius of the broad-line region (BLR) (B17).B17 reported the optical and UV light curves of PS16dtm which brightened approximately two magnitudes above the archival host brightness in ∼ 50 days, with two recognizable peaks distanced ∼ 50 days between them.Furthermore, the light curves showed little color evolution, and stayed approximately at the Eddington luminosity of the SMBH.
B17 considered scenarios of a variable AGN, a Type IIn supernova, and a TDE as possible explanations for PS16dtm.They concluded that the small evolution in brightness and color seen in PS16dtm, is not typical for supernovae where cooling is expected due to expansion and radiative losses, but it is similar to that of TDEs.Another argument in favor of the TDE scenario put forward by B17 is the short timescale of the rising part of the PS16dtm light curve (∼ 50 days) and the amplitude (2 magnitudes rise over the host level).The typical variability amplitudes of AGN are low in large AGN samples and they vary on timescales of years (∼ 0.1 − 0.2 magnitudes per year, see e.g.Sánchez-Sáez et al. 2018).Another interesting aspect of PS16dtm is that after the outburst, X-ray emission from the nucleus decreased by at least an order of magnitude, compared to archival measurements (Pons & Watson 2014).B17 argued that neither the supernova scenario nor AGN variability can explain the X-ray dimming, but that it could be explained by the accreted stellar debris that obscures the X-ray emitting region of the AGN accretion disk.They proposed a simple model which assumes a face-on orientation for the preexisting, X-ray emitting disk, and a nearly edge-on orientation for the disk newly formed in the disruption, which also blocks the X-ray emission.They argued that such a geometrical configuration can also explain the spectral properties of both the debris disk and the host galaxy.
The spectra shown in B17 are dominated by hydrogen Balmer and Fe II emission lines.B17 compared the PS16dtm spectra to those of Type IIn supernovae, and found no reasonable matches.B17 interpreted the spectra as being more similar to those of some NLSy1 galaxies, as optical Fe II (4000-5400 Å) emission are common features in NLSy1 spectra.
PS16dtm also flared in the mid-infrared (MIR) and it was detected by WISE as part of the NEOWISE survey (Jiang et al. 2017, J17;hereafter).J17 reported three NEOWISE detections starting 11 days before the optical ASSA-SN survey and extending to 327 days after.The authors conclude that the MIR flare is consistent with a dust echo, despite the fact that the MIR data shows no delay with respect to the optical.They estimated that the inner radius of the preexisting dust torus increased from ∼ 10 light days to ∼ 70 light days, where the emptied region was replaced with a gas torus.They also argued that the detected Fe II emission lines are produced in the gas coming from the evaporated dust due to the strong radiation field.
The peculiarity of PS16dtm has prompted other authors to examine what powers its luminosity.Frederick et al. (2021), who looked at NLSy1 transients from the literature, put PS16dtm in the class of TDE with strong Balmer and Fe II complexes, despite that it does not satisfy all their defined requirements for TDE classification.Moriya et al. (2017) argued that PS16dtm can be explained by changes related to the AGN, and not as a result of tidal disruption of a star by the SMBH.They explain it as an interaction between accretion disk winds and clouds in the BLRs surrounding them.They also argue that the observed broad (∼ 10, 000 km s −1 ) Mg II absorption in the UV (B17) could be related to the fast SMBH disk wind.Nevertheless, their model predicts that the emission timescale is ∼ 110 days, while PS16dtm has stayed bright above the host-galaxy baseline for much longer, as will be shown in the next sections.

Observations
In this section, we present new observations of PS16dtm out to ∼ 2000 days after discovery.We also present ATLAS and PSST data of PS16dtm for the first time.Notably, ATLAS detected PS16dtm before the PSST and ASSA-SN surveys, and this enabled us to roughly estimate the time of outburst, as it will be shown in Sec.3.2.1.The evolution of the host-subtracted photometry corrected for Milky Way extinction is shown in Figs. 1  and 2. The gaps in the data (from February to August) are when PS16dtm was behind the Sun.We note that at ∼ 2000 days after the discovery, PS16dtm is fading and the emission is approaching preoutburst, host levels (see Hinkle et al. 2021, for measured and synthetically computed PS16dtm host magnitudes).All our photometry will be available through the WISeREP archive2 (Yaron & Gal-Yam 2012).

Swift observations
PS16dtm was monitored by the Neil Gehrels Swift Observatory (Gehrels et al. 2004) with the UV/Optical Telescope (UVOT Roming et al. 2005) in six filters; three optical (V at 5468 Å, B at 4392 Å, U at 3465 Å), and three near-UV (UVW1 at 2600 Å, UV M2 at 2246 Å, and UVW2 at 1928 Å).There are 83 epochs of UVOT/Swift observations spanning from MJD 57632 to MJD 59493.The Swift/UVOT observations were reduced following the standard pipeline from HEAsoft3 .The photometry was extracted using the tool UVOTSOURCE and a source extraction region with a radius of 5".We corrected for Galactic extinction and subtracted the contribution from the host galaxy, for which we used the PS16dtm host galaxy magnitudes from Hinkle et al. (2021).Hinkle et al. (2021) provided corrections to the NUV photometry of TDEs published after 2015 in the literature, which were using Swift/UVOT.Their work was motivated by an update by the Swift team to the UVOT calibration to correct for the loss of sensitivity over time.Hinkle et al. (2021) fitted archival multiwavelength photometry from GALEX, 2MASS, SDSS and WISE of the host galaxy and modeled the spectral energy distribution (SED) of the host galaxy, from which they updated the magnitudes in the NUV UVOT filters.Therefore, the NUV Swift/UVOT magnitudes of the PS16dtm host galaxy are different to those used in B17, leading to an average difference in the host-subtracted photometry of PS16dtm compared to the ones in B17 by 0.58 magnitudes for UVW2 and by 0.02 − 0.07 magnitudes for U, B, V, UVW1 and UV M2.
Simultaneous with the Swift/UVOT observations, PS16dtm was observed with the Swift X-ray Telescope (XRT).After building the XRT light-curve using the online tool4 we conclude that there was no confident detection by Swift/XRT and we can only place upper-limits to the X-ray emission, similar to those presented in B17 (F X < 1×10 −14 erg s −1 cm −2 and L X < 1.7×10 41 erg s −1 ).

ATLAS
ATLAS is a robotic survey that uses at least two identical 50centimeter telescopes, located at Haleakala and Mauna Loa observatories in Hawaii.PS16dtm was observed by ATLAS in two bands, c (cyan) and o (orange), that cover the wavelength range of 4200−6500 and 5600−8200 Å, respectively.The ATLAS image processing is done with a fully automated pipeline that performs flat fielding, astrometric calibration, and photometric calibration.The photometry of PS16dtm is publicly available online, of which we used the reduced photometry from the latest ATLAS data release (Tonry et al. 2018;Smith et al. 2020).Since the reference images contain the flux from PS16dtm, we calculated the mean flux between MJD 57250 and 57500 (prediscovery epochs) and subtract the flux from all the data subsequently.As visible from Fig. 2, the first ATLAS detection of PS16dtm was made at MJD 57591.6, which is ∼ 6 days before the first detection reported by ASSA-SN and ∼ 20 days before Pan-STARRS.We used a second order polynomial to fit the points in the rising part of the light curve and found the intersection with the host level.We found that the outburst likely happened around MJD 57577, with an uncertainty of ∼ 10 days.Throughout the paper we use this date as a reference epoch.We note that this is slightly different than the reference epoch used by B17, which was arbitrarily set to MJD 57600.

Pan-STARRS
Pan-STARRS1 uses a 1.8 m telescope and it is located at Haleakala, Hawaii.After Pan-STARRS1 discovered PS16dtm on 12 August 2016, it also provided photometry from the follow-up observations in two broad filters, the w PS1 and i PS1 bands, which have wavelength ranges 4000-8300 Å and 6800-8300 Å, respectively (Tonry et al. 2012).Images obtained by the Pan-STARRS1 system are processed automatically with the Image Processing Pipeline and transient sources are identified through analysis of difference images, created by subtracting a template from the observed image taken as part of the search for the counterpart (see details in Magnier et al. 2013;Huber et al. 2015).

Liverpool telescope
We collected ten epochs of multiband (griz) imaging of PS16dtm with the IO:O instrument at the robotic 2-m Liverpool telescope at the Roque de los Muchachos Observatory, Spain (Steele et al. 2004).The images were reduced with the IO:O5 pipeline and were subtracted against PSST (Tonry et al. 2012) reference imaging, leaving only the transient light.Then PSF photometry was done on the source if it is detected after the subtraction and the photometry was calibrated relative to PSST photometric standards.
Photometric evolution of PS16dtm.The photometry has been corrected for the Galactic extinction and the host contribution has been subtracted.The reference MJD for the outburst is 57577 obtained from a second order polynomial fit to the ATLAS o photometry in the rising part of the light curve.The photometry has been arbitrarily shifted in the y axis for easier viewing, as indicated in the legend.

NEOWISE MIR photometry
NEOWISE is a project by the Wide-field Infrared Survey Explorer, WISE (Wright et al. 2010), which surveys the sky in 3.4 (W1) and 4.6 µm (W2) (Mainzer et al. 2014).We retrieved the photometry from the IRSA public data archive 6 .We first computed the variance-weighted average for 1-day bins and filter for all epochs after PS16dtm was detected.To remove the host contribution from the transient light curve, we computed the variance-weighted average of all preoutburst data and subtracted the host contribution from the transient light curve.In Fig. 2 we show the host-subtracted MIR light curve together with the optical data from ATLAS and PSST.ner with the aid of the PESSTO pipeline (Smartt et al. 2015), to apply bias and flat-field corrections, determine a wavelength solution, and calibrate the relative flux with a standard star observed in the same setup.Following common practice, the intraday spectra of grism 11 and grism 16 were combined, since they are characterized by almost the same spectral resolution.The spectra were corrected for the Galactic extinction and the absolute flux calibration was improved by scaling the spectra with the aid of photometry.A spectroscopic log can be found in Table 1 and the spectral series is shown in Fig. 3.All reduced spectra will be available through the WISeREP archive7 (Yaron & Gal-Yam 2012).

VLT/X-shooter spectrum
We also observed PS16dtm with X-shooter (Vernet et al. 2011) on the Very Large Telescope (VLT), on MJD 59595, that is at +1868 rest-frame days.X-shooter is a medium resolution spectrograph covering the wavelength range from 3000 to 24800 Å in three spectroscopic arms.We used slit widths of 1.0", 0.9" and 0.9" for the UVB, VIS and NIR arms respectively, resulting in nominal resolutions of R=5400, 8900 and 5600.The data were reduced by employing the X-shooter pipeline in the EsoReflex GUI environment (Freudling et al. 2013), as implemented in Selsing et al. (2019).

X-ray observations
The Chandra X-ray Observatory (CXO) observed PS16dtm on three epochs on 2019 January 10 (MJD 58493), 2019 November 11 (MJD 58798), and 2020 September 29 (MJD 59121) 8 for a net exposure of 10, 10, and 20 ks, respectively.We reduced these IDs with the CIAO v4.14 software (Fruscione et al. 2006) in presence of the latest calibration files using the chandra_repro pipeline.Standard filtering and analysis methods are further applied to Advanced CCD Imaging Spectrometer mode data in our study.We found weak X-ray emission in the 0.5-8 keV band at the location of PS16dtm.The source is moderately present in a 1 arcsec region at a significance level of 2.6, 1.6, and 2.5σ level in the first, second, and third Chandra epochs with a count rate of about 4×10 −4 , 2.3×10 −4 , and 2.7×10 −4 counts s −1 , respectively.No significant emission, however, is detected in the 0.5-1 keV range.
We subsequently performed the spectral analysis to examine the long-term X-ray flux evolution of the host galaxy and its connection to the transient.As the number of detected source photons is limited to <5 in each observation, all three Chandra spectra were stacked together to carry out the analysis using the Cash statistics in XSPEC (Arnaud 1996).Using a redshifted power-law model with a column density fixed at 2.5×10 20 cm −2 (Pons & Watson 2014, B17), the photon index is found to be 0.8 +8.5  −0.1 along with a fit null hypothesis probability of 0.9 at four degrees of freedom.The errors are calculated for a 90% confidence interval using the Markov chain Monte Carlo simulation method in XSPEC.The 0.3-10 keV unabsorbed source flux is (1.1±0.8)×10−14 erg s −1 cm −2 that is obtained by the cflux convolution model.Similarly, in place of the power-law component, a red-shifted blackbody model at the above-fixed column density is also tested on the Chandra data.We found a blackbody temperature of 1.5 +2.0 −0.5 keV along with a null hypothesis probability of 0.8 at four degrees of freedom.The 0.3-10 keV unabsorbed model flux is estimated to be (8±3)×10 −15 erg s −1 cm −2 in the latter case.Based on our spectroscopy, we conclude that the source flux remains close to the last measurement (upper limit of 1×10 −14 erg s −1 cm −2 ) provided by XRT in 2016 and 2017 (B17).In 2019 and 2020 the system was still present at an order of magnitude lower than the 2005 XMM-Newton detection (Pons & Watson 2014) before the transient.

Light curve and SED analysis
The data presented here, make PS16dtm one of the few nuclear outbursts that have been observed with a photometric and spectroscopic campaign that spans over six years of the transient activity.After the first peak, at which PS16dtm reached the absolute magnitude of M V = −22.0mag in the UVOT V band, the luminosity dropped for ∼ 50 days, after which it rose again to a second peak and stayed at almost plateau level for ∼ 100 days.Then the flux decreased steadily with time at all wavelengths, as seen in Fig. 1.However, the decline is most pronounced up to +270 days in the NUV/UVOT filters, afterward the decline is slower.Interestingly, there have been few other TDEs which have shown double peaks in the light curve, such as ASSASN-15lh (Leloudas et al. 2016) and AT 2018fyk (Wevers et al. 2019), but the physical process behind this feature is still not clear.
In order to explore the physical parameters of the transient, we attempted to fit a blackbody to SEDs constructed from Milky Way extinction-corrected and host-subtracted photometry.We used the Markov chain Monte Carlo (MCMC) method with a publicly available code9 that uses the Python package emcee (v 3.0.2) (Foreman-Mackey et al. 2013).It was already pointed out in B17 that a simple blackbody is not a good fit to the SED, which they constructed with the first 200 days of PS16dtm data.Therefore, B17 fitted a blackbody by excluding the photometry from the three NUV/UVOT filters by arguing that there must be a significant obscuring material in the line of sight.This implies that true luminosity could be higher, but they argued that it should not be more than a factor 2 (that is not higher than 5 × 10 44 erg s −1 ).We also note the updated Swift/UVOT calibration correction from Hinkle et al. (2021) which only affected the NUV filters, making them fainter compared to B17 measurements.This means that this problem persists and it is even more pronounced than from what B17 estimated.
We repeated the same exercise of fitting a blackbody to our photometry without the NUV/UVOT bands, requiring at least four photometric measurements in each epoch.The results are shown in Fig. 4. Keeping in mind this problem with likely underestimated temperature due to the omission of NUV/UVOT bands, the blackbody fits would indicate an initial rise in the temperature to ∼ 15000 K contemporaneous with the rise in the luminosity.After that, the temperature drops to ∼ 10000 K, and is held relatively constant.After ∼ 500 rest-frame days, the blackbody fit to the SED appears increasingly erratic, as at later epochs when the host-subtracted photometry of PS16dtm becomes fainter, thus noisier.
We explored the possibility that the inability to fit a single blackbody to the NUV/optical data is caused by the extinction in the line of sight by the circumnuclear dust near the SMBH, despite that there is no clear consensus regarding the wavelength dependence of the extinction in the line of sight to AGN.However, some studies based on individual reddened AGN have suggested they have an extinction similar to that measured in the Small Magellanic Cloud (SMC).10SMC extinction strongly affects the UV part of the SED (Willott 2005) and compared to the )  Milky Way or the one used for starburst galaxies (Calzetti 2001), SMC average extinction curve has essentially no 2175 Å bump and it has a strong far-UV rise (see e.g., Salim & Narayanan 2020, for a recent review).We attempted to fit blackbody SEDs reddened with the SMC-like extinction curve (Gordon et al. 2003) where we add E(B − V) as a free parameter.It was pos-sible to fit the SED to a blackbody with T ∼ 2 − 4 • 10 5 K and a bolometric luminosity of L ∼ 10 48 − 10 49 erg s −1 , at a corresponding blackbody radius of R ∼ 1 • 10 15 cm, and SMC-like extinction curve with E(B − V) ∼ 0.5 mag.This best fit model is clearly nonphysical, with such high temperatures, which would also imply an extremely bright X-ray source.In fact, this solution for sufficiently large (that is above 20000 K) temperatures corresponds to a degeneracy between temperature and radius, where all the measurement points lie on the high-wavelength slope of the distribution, and for any arbitrarily large temperature there exists a radius that produces the observed light curve.In conclusion, a single blackbody does not provide a satisfactory fit to the SED, even when we included extinction.However, the fact that the SED cannot be described with a single blackbody does not exclude the interpretation of PS16dtm as a TDE.According to the radiative transfer calculations by Roth et al. (2016) and Thomsen et al. (2022), the optical/UV continuum of TDEs is not necessarily described by a single blackbody.In addition, it is possible that the SEDs of TDEs in AGN are different compared to those in inactive galaxies due to the preexisting AGN disk, as suggested by the simulations in Chan et al. (2021).
We also fitted the fading light curve with a power-law profile, where the power-law index is fit freely, Our best fit parameter, α ∼ 5/7 for the NUV/UVOT light curves for the initial (158 < t < 270) days, while ∼ 270 days after the outburst, it follows a α ∼ 1/6 decay.For the bolometric luminosity, the best-fit power-law has the index α ∼ 0.98 (see Fig. 4).Interestingly, a similar decline trend is also seen in the bolometric light curve of the TDE AT 2017gge (L ∝ t −1 , Onori et al. 2022), which also shows a MIR echo and coronal emission lines in the late-time spectra, although more numerous and more intense than what it has been detected in PS16dtm, as shown in the next sections.Nevertheless, the luminosity evolution of PS16dtm cannot be simply explained by the t −5/3 decline as expected from simulations of the fallback rate for the stellar debris onto the SMBH (Rees 1988;Gezari et al. 2009).At the moment, the origin of the optical emission from TDEs remains a puzzle, but it might be powered by an inner accretion disk or by shocks of the intersecting debris streams (see van Velzen et al. 2020, for a recent review).In the first months, optical and UV light curves of TDEs in quiescent galaxies often decline as t −5/3 , but after a few months, models that include reprocessing of the disk emission by outer debris predict weaker temperature evolution, are expected to show weaker temperature evolution and to follow a flatter decline as t −5/12 (Lodato & Rossi 2011).Furthermore, if the disruption of the star is only partial, the expected fallback rate could also be shallower (Guillochon & Ramirez-Ruiz 2013).Interestingly, van Velzen et al. (2019), by observing TDEs at late times in the far-UV (FUV), found that for the light curves of TDEs from low-mass SMBHs (M BH < 10 6.5 M ), the early-time decay follow a t −5/3 power-law decline, but the later-time evolution is much shallower.They conclude that this could be the sign of different disk emission mechanisms operating at early and late times.They argue that unobscured accretion disk models, rather than reprocessing and circularization paradigms, can explain the late-time FUV emission.

Blackbody fits with the NEOWISE MIR photometry.
As mentioned in the previous sections, PS16dtm exhibits MIR emission almost simultaneously with the optical one.In Fig. 2 we show that the first NEOWISE detection is ∼ 4 days before the first ATLAS detection.In the first 500 days, PS16dtm shows a steep rise in the MIR light curve.After that, it keeps rising, albeit at a much slower rate.Another noticeable trait from Fig. 2 is that not only the rise in brightness, but also the change in color that happens in the first epochs, that is at +9 days W1 − W2 has a negative value of −0.37, then at +165 days it is −0.07, and in the epochs afterward to the last epoch, from +349 to +1700 days, W1 − W2 remains steady with a positive value of 0.24 − 0.35 mag.
First, we attempted to fit a single blackbody to the data by including the optical and the MIR photometry at quasi simultaneous epochs.As visible from Fig. 2, the first MIR point has a very different color than the rest of the MIR light curve, so we examined the possibility that it is compatible with the Rayleigh-Jeans tail of a single blackbody.As shown in the example in Fig. 5, the fit would indicate a temperature of 8900 K at the earliest epochs after the outburst, however the single blackbody is a poor fit to the data.Second, we attempted to fit a double-blackbody model to the optical and the MIR NEOWISE photometry, as it was done in J17, for the first three NEOWISE epochs.The first NEOWISE epoch at +9 days, yields ∼ 2300 K for the dust blackbody temperature, while the second epoch at +165 days, yields ∼ 1300 K, and the epochs afterward (from +349 to +1700 days) settle at ∼ 900 − 1000 K. J17 interpreted the MIR emission as light coming from the dust thermal emission heated by the TDE; in the first epoch, the temperature is above the sublimation temperature; after the second epoch, the dust temperature dropped below the sublimation temperature, the MIR emission is from larger distances from the SMBH, so the dust in the inner region must be optically thin.Before we proceed with our spectroscopic analysis, we summarize here the findings presented in B17 in which they used ten spectra that spanned from +54 to +198 days: -The main spectral features of PS16dtm are the multicomponent hydrogen Balmer lines and Fe II lines, visible at all epochs.-The Hα exhibited a complex asymmetric profile consisting of a broad component and a superimposed narrower component that has a slight shoulder on the blue side of the peak.-Several broad features appear near the [O II] λλ7320, 7330, [Ca II] λλ7291, 7324, which they interpret as a blend of several Fe II lines.The Fe II lines became stronger with time.-In the NIR part of the spectra, they noticed several broad emission features with asymmetric profiles, mainly Fe II lines, out to 1.2 µm, beyond which the spectrum is relatively featureless, -Paschen γ and other lines of the Paschen series were detected.Paschen α and β were not detected because they fall in the NIR telluric bands.-The Ca II triplet shows an unusual shape where the central line in the triplet is much stronger than the others.-Blueward of ∼ 4500 Å the spectrum shows a very complex combination of narrow features superimposed on broad features.They identify [O II] λλ3727, Balmer lines, and additional Fe II lines.-In their UV spectrum obtained with HST, they found evidence for absorption, in particular, broad Mg II λλ2800 absorption lines were detected, with ∼ 10, 000 km s −1 , that can perhaps indicate the presence of an outflow.
B17 also studied the temporal evolution of the hydrogen Balmer lines Hα and Hβ.They fitted the spectral lines with the sum of three Gaussians.They reported that the widths of the intermediate component of Hα is around 750 km s −1 that narrows as a function of time, going from 900 km s −1 to 600 km s −1 .For the broad component, they found 3500 km s −1 which increases to 4000 km s −1 in the decline phase of the light curve.They found that the narrow component is 100 km s −1 from their spectrum with the highest resolution.For their low-resolution spectra the narrow component is unresolved, so they fix it at the instrumental resolution.Furthermore, they found that the width in the earliest epoch of their intermediate component was similar to the width of the preexisting broad component in the archival host spectra.For this reason, B17 argues that the intermediate component that they measure in PS16dtm might be associated with the BLR of the NLSy1 host galaxy.The flux of the narrow component of Hα has remained unchanged relative to the host spectrum.

Spectral evolution to 1868 rest-frame days after the outburst
Now we turn to the analysis of our spectra, 16 of which were taken with the low-resolution NTT/EFOSC spectrograph at +155 to +1743 rest-frame days, and one medium-resolution spectrum with the VLT/X-shooter at +1868 days.The spectra are shown in Fig. 3. Similar to B17 (see their Fig.7), we also found that the main spectral characteristics are the hydrogen Balmer lines and the Fe II complexes.The most noticeable variability in our spectra is that the blue continuum becomes weaker and almost disappears in the final X-shooter spectrum at +1868 restframe days after the outburst.Perhaps even more striking are the Fe II optical lines, in the blue (around Hβ) and red (around Hα) line, that become weaker as time progresses (see the zoomedin view around the Hβ line in Fig. 6).This is a first time such strong Fe II emission is seen in a TDE, even when compared to other nuclear transients for which Fe II emission has been claimed, such as J123359.12+084211.We show the spectroscopic comparison of the iron-rich TDE candidates in Fig. 7. Some TDE and flaring AGN spectra show O III and N III lines which have been explained as due to the mechanism of Bowen fluorescence (Blagorodnova et al. 2019;Leloudas et al. 2019;Trakhtenbrot et al. 2019;Onori et al. 2019).
We also searched for these lines, but in PS16dtm these are difficult to resolve due to strong blending with Fe II emission.An extremely weak feature is resolved on the place of a N III λ4640 line in the X-shooter spectrum (see later discussion).We also detected strong features in the red part of the spectrum, in the range 7000-8000 Å, which is fading with time, and completely disappearing in the X-shooter spectrum.These are most likely blended Fe II emission, as already noted by B17, with some probable contribution from the other nearby lines, such as O I λ8446, is clearly identified in the +1514 spectrum (see Fig. 3).

Iron emission-line model for the spectral fitting
Next, we proceed to identify the emission lines and understand their temporal evolution, despite that, given the richness of features and their strong variability, their behavior is difficult to follow and disentangle since all emission lines are strongly affected by blending (see Fig. 6).The complex Fe II ion can produce thousands of line transitions, thus these lines are typically blended and difficult to identify in the spectra.With this caveat in mind, we fit all spectra using a python-based tool called ) Fully Automated pythoN tool for AGN Spectra analYsis (FAN-TASY) 11 , that was initially developed for fitting AGN optical spectra (Ilić et al. 2020;Rakić 2022, Ilić et al. 2022, in prep.).
The code fits the underlying broken power-law continuum and sets of emission lines, with the Fe II model based on the atomic parameters of Fe II.In contrast to widely used fully empirical iron templates (for recent discussion on different Fe II templates, see Park et al. 2022 and references therein), such as that by Boroson & Green (1992), our approach is slightly different and it was first developed by Kovačević et al. (2010).The semi-empirical model of Kovačević et al. (2010) relies not only on the observed properties of AGN spectra, but also on atomic properties of the transitions, so the Fe II line sets are grouped according to the same lower energy level in the transition with line ratios con-  (Ilić et al., in prep).This approach can be useful to understand newly discovered transients, which show strong and complex Fe II emission.We assume the following constraints for the emission line ratios, widths and shifts: (i) the broad Balmer emission lines (Hα, Hβ, Hγ)  In Fig. 8 we show an example of the results of the multicomponent spectral fitting in the 4100 − 7000 Å rest-frame wavelength range.Notably, our model was able to reconstruct the observed spectra, but only when assuming a broad component of Fe II multiplets.The presence of the broad Fe II component has been previously detected in AGN (e.g., Dong et al. 2011;Park et al. 2022); however, since this component is typically very weak, it is often merged with the continuum emission, especially in poor spectral resolution data.We also detected narrow features on top of the broad Fe II components, which our code usually identified as narrow Fe II emission, despite the limited instrumental resolution.We were also able to identify the broad features blueward from Hα region with the Fe II model (see Fig. 8).We performed a series of tests with other emission-line models (such as no broad Fe II, no forbidden Fe II lines, etc.), but none were able to produce reasonable fit based on the fitting residuals.The goodness of the fit was evaluated using the χ 2 parameter.In Fig. 9 we show the time evolution of the FWHM of the most prominent emission lines, broad Hα and Hβ, broad Fe II, and all narrow lines.It is evident that the strength of the lines slowly decreases with time.The width of the narrow lines remains the same, but this is constrained by the instrumental resolution.
In Fig. 10 we show the time evolution of the integrated emission-line features, that is, for the broad Hα and Hβ line, total broad Fe II, total narrow and forbidden Fe II lines.The strength of broad line-blends (Hα, Hβ, total Fe II broad) are showing similar trends, slight increase followed by long decrease with time.A similar evolution is seen in the narrow Fe II blend.The exact time of peak is difficult to identify due the fluctuations in line fluxes.Furthermore, the extraction of single lines Hα and Hβ from the low-resolution spectra is difficult; it depends on the possibility for the code to identify the narrow [O III] λ5007 line.Nevertheless, the general similar trend is evident in all light curves, especially of line-blends.
In Fig. 10 we also show the ratio of the bolometric luminosity and the Eddington luminosity, L bol /L Edd .The Eddington luminosity was estimated as L Edd = 1.26 × 10 38 (M BH /M ) erg s −1 with the SMBH mass of M BH ∼ 10 6 M , which we determine in Sect.4.3.8.For the bolometric luminosity L bol we used a different approach than the one used in Sect.4.1.Here, to obtain the bolometric luminosity we used a standard procedure for quasars, where L bol = k bol λL λ by applying the mean quasar bolometric correction k bol ≈ 10 (e.g., Richards et al. 2006;Runnoe et al. 2012) to the continuum luminosity at 5100 Å extracted from the multicomponent fitting.This assumes that the underlying continuum originates from the accretion disk of an AGN, and it remains to be investigated if this is a valid approximation for an AGN hosting a TDE.We subtracted the host-galaxy contribution to the L 5100 luminosity before calculating the Edding-ton ratios (see Fig. 11).The host-galaxy spectrum was extracted from the SDSS spectrum, using the principal component analysis (Vanden Berk et al. 2006), as explained in Ilić et al. ( 2020) (see Sect. 5.2 for further discussion).A word of caution is given for the flux at 5100 Å measured from the fits, since it is sensitive to the spectral multicomponent fitting, but even more to the absolute spectral calibration.Here all spectra are calibrated using photometry, however, for more accurate absolute calibration in the AGN monitoring, one usually applies more precise intercalibration based on the constant narrow-line flux, such as [O III] or [S II] lines (van Groningen & Wanders 1992;Fausnaugh 2017).
In order to assess the source of ionization of the broad lines and blends, in Fig. 12 we plot the correlation between the continuum flux at 5800 Å (F cnt ), and the flux of Hα (left), Hβ (middle), and total Fe II blend (right).F cnt is measured from the observed spectra as the median of 5790-5810 Å range, since this part is identified as free of emission lines.Apart from the first three epochs (+155, +179, +196) the line fluxes correlate with the continuum flux, supporting that photoionization by the central continuum source is the main heating source of the line emitting gas (see e.g., Netzer 2006).The outliers are seen in each plot, indicating that in the first period the emission may be induced by other mechanisms, such as shock ionization.This is also supported by the higher gas velocities measured from the line widths.It is important to note that for consistency, all epochs were uniformly fitted with the same initial constraints, listed above.

Spectral features in the 7000 − 9000 Å range
In the spectral region redward of the Hα line, there are broad features visible in the spectra shown in Fig. 3.These are most notably He I lines, namely He I 6680 (which is attached to the Hα red wing), He I λ7065, He I λ7281.There are also Ca II triplet (λ8498, λ8542, λ8662), well known oxygen lines O I λ7774 and O I λ8446, as well as the Mg II λ7892 line present.All of these broad features are of similar width as the broad Hα component, being at the maximum at a similar time, and decreasing in intensity and width across time in the similar fashion.A broad feature around 7300 Å was identified as [O II] λ7320, λ7330, [Ca II] λλ7291, 7324 by B17, but the authors suggested that these fea- ture could be associated to Fe II emission lines.However, we believe that this a mixture of Fe II blend and He I λ7281, as well as the feature centered around 8200 Å.As it will be shown in the next section, all broad features disappeared in the last VLT/Xshooter spectrum, with possibly only weak Mg II and O I being present.
4.3.6.Analysis of the medium-resolution VLT/X-shooter spectrum Given that the resolution of the VLT/X-shooter spectrum is much higher compared to those taken with NTT/EFOSC2, we were able to clearly resolve the broad component in the hydrogen Balmer Hα and Hβ lines (Fig. 13), whereas the broad component of the Balmer Hγ line is too weak to be detected).In the X-shooter spectrum, taken at +1868 days, it is noticeable that the broad Fe II emission has fully faded, and that the galaxy is most likely returning to its preoutburst state.Nevertheless, we were able to identify the strongest narrow Fe II lines coming from the a 6 S, a 4 G, and b 4 F multiplets, as well as forbidden [Fe II] lines, displayed in the zoom-in plots in Fig. 13.The narrow emission lines can be now unambiguously identified and used for the diagnostics of physical conditions in the ionized gas.We have extracted the following narrow lines: Hα, Hβ We fit the Hα and Hβ line region with the multicomponent model in the same manner as described in Sect.4.3.3 to measure the narrow emission line fluxes, as well as to extract the broad component of the Balmer lines.The strongest narrow emission lines, two [O III] lines and Hα are fitted with two Gaussians, one for the core and the other for the line wings (Fig. 14).These are all constrained to have the same width and shift, and two a 6 S Fe II lines ratios are fixed according to our model (see Sect. 4.3.3).
The width of Fe II lines is ∼220 km s −1 , larger than the narrow lines, which have widths that are closer to ∼100 km s −1 .This is evident for all iron lines, including the coronal [Fe VII] λ6087 line.The extracted broad Hα and Hβ components (obtained after subtracting narrow lines) have similar profiles (Fig. 14, bottom panel), with Hα being slightly broader with the width of ∼1900 km s −1 .Both broad lines show slightly asymmetric profiles, with possibility that the line peak is a slightly shifted redward.

The presence of coronal lines in the X-shooter spectrum
The spectral lines coming from forbidden transitions of highly ionized ions (ionization potential larger than 54.4 eV) are typically referred to as coronal lines, and are known to exist in the optical spectra of Seyfert galaxies (see for example Korista & Ferland 1989;Gelbord et al. 2009).Komossa et al. (2008) and Wang et al. (2012) discovered several coronal line emitters which can be interpreted as light echos from past TDEs.The high ionization potential of coronal lines indicates that soft Xrays are required for their production.These lines are typically of somewhat larger width than the narrow lines, with higher critical densities (∼ 10 7.5 cm − 3 Osterbrock & Ferland 2006).It is assumed that they come from either the inner narrow-line region or from the inner edge of the torus (e.g., Rose et al. 2015).Their strengthening has been associated with some sort of transition event, such as the awaking of a dormant AGN (Ilić et al. 2020) or a TDE in a gas rich environment (Wang et al. 2012).Recently, Onori et al. (2022) reported the detection of coronal lines in the TDE AT 2017gge 1700 days after the initial outburst.
In the X-shooter spectrum, we identified a weak [Fe VII] λ6087 line (see bottom mid-panel in Fig. 13).Moreover, two other [Fe VII] lines at 5159 Å and 5276 Å could be also misidentified as Fe II lines.Another well-known coronal line [Fe X] 6374 seemed to be present as well, but some telluric disturbance makes its detection quite hard.We examined all previous spectra and the [Fe VII] λ6087 is either not present or too weak to be resolved in the low-resolution spectra.We also re- trieved the spectrum with higher resolution from B17 at +192 taken with the MagE spectrograph to check if the aforementioned line was present, but it is also either too weak to be seen or not present.We note that there could be many other coronal lines (for a review of coronal lines, see e.g., Rose et al. 2015) hidden in the strong iron blends throughout the evolution of transient events such as PS16dtm.Future monitoring campaigns of nearby transient events with high resolution instruments and high S/N spectra could shed light on the variability of these lines and on their connection with the occurrence of such transient events.It remains unclear what could be the origin of coronal lines in PS16dtm, since the X-ray emission remains weak in this object (see Sect. 3.6) contrary to the case for AT 2017gge and AT 2019qiz which experienced considerable late-time X-ray brightening.

SMBH mass and AGN bolometric luminosity estimation
To be able to measure the SMBH mass from the X-shooter spectrum, we first needed to check if the transition to the preoutburst state, has occurred.Therefore, we compared the X-shooter with the archival, preoutburst SDSS spectrum.We attempted to scale the X-shooter spectrum to have the same spectral resolution as the SDSS spectrum, for which we used the [O I] 6300 Å that is clearly seen in both spectra and is not contaminated with satellite lines.However, since the SDSS and X-shooter spectra are taken with different apertures (2 and 0.9 − 1 , respectively), the absolute re-scaling of the spectra using [O III] or [S II] emission line fluxes was not possible.This was due to the fact that the narrow emission lines are predominantly originating from the host galaxy (see Fig. 11), and the host-galaxy contribution was significantly different for the two used apertures.Therefore, we compared the X-shooter spectrum smoothed to the SDSS spectral resolution, to the pure AGN component extracted from the SDSS spectrum, which was scaled-up so that the Hα or Hβ have the same broad line intensity, as shown in Fig. 15.We then concluded that the broad line profiles of Hα or Hβ in the archival SDSS spectrum and in the X-shooter spectrum at +1868 days, are similar.Using the approach from Sect.4.3.4where the bolometric luminosity is expected to scale with the luminosity at 5100 Å , we estimated the AGN bolometric luminosity from the X-shooter spectrum to be L bol = 3.7 × 10 43 erg s −1 , which is almost reaching the preoutburst value, which we extracted from the SDSS spectrum to be L bol = 1.6 × 10 43 erg s −1 .Both values are host-corrected using the same constant host contribution estimated from the SDSS spectrum.
Given that the broad Hβ line can be clearly extracted from the X-shooter spectrum, we calculated the mass of the SMBH using the single-epoch method (see e.g., Gaskell 2009, and references therein) and the FWHM of the Hβ line of 1160 ± 190 km s −1 .The radius of the BLR is estimated using the radius-luminosity relation from Bentz et al. (2013) applied to the AGN continuum luminosity at 5100 Å.The M BH is then calculated using the virial theorem assumption, with the virial factor of f = 0.75 (the same as in Xiao et al. 2011, who estimated the mass for this object).The obtained mass is M BH = 10 6.07±0.18M , similar to the previous estimates by B17 and Xiao et al. (2011).We note that the obtained value is dependent on the assumption for the virial factor f (see discussion in e.g., Dalla Bontà et al. 2020, and references therein).

Spectral classification of the host galaxy
Previous works concluded that the host of PS16dtm is a NLSy1 galaxy (see Sect. 2 in B17), based on the detection of a weak broad-line component in the Hα line (Greene & Ho 2007;Xiao et al. 2011) and high X-ray luminosity of L 2−10keV ∼ 10 42 erg s −1 (Pons & Watson 2014).The later could not be easily accounted for by star formation only, as it would require star formation rates of at least 200 M yr −1 (Ranalli et al. 2003).Some of the most common features of NLSy1 galaxies are strong and rich Fe II multiplets, however we note that in the quiet phase, the archival optical SDSS spectrum of the host shows no (or very weak) Fe II emission, which classifies this object as a rare type of NLSy1 (Zhou et al. 2006).It is interesting to note that the spectrum of the host galaxy of the nuclear transient CSS100217 (SDSS J102912.58+404219.7)shows more typical NLSy1 features, such as stronger broad hydrogen Balmer lines, He I, and very prominent Fe II emission around Hβ (Drake et al. 2011).However, both host galaxies are not typical NLSy1 objects.Their location on the Baldwin-Phillips-Terlevich (BPT) diagrams (Baldwin et al. 1981;Kewley et al. 2006) is either left to or at the border of the AGN-starburst composite objects location, suggesting the presence of significant star formation.
B17 placed the host AGN more on the border of the starburst-AGN division based on the line-ratios measured from the archival SDSS spectrum.Given that the X-shooter spectrum is showing that is returning to its preoutburst state, we can make use of its higher resolution to perform the diagnostics using the emission line ratios of [O III] λ5007/Hβ, [N II]    (Baldwin et al. 1981;Kewley et al. 2006).In Fig. 16 we plot the diagnostics diagrams using the narrow emission lines measured from the X-shooter spectrum.
The location on the BPT diagrams, place the host galaxy in the area of AGN-starburst composite objects.This means that the narrow-line emission has significant contribution to the ionization from stellar sources, much more than previously shown in B17.

PS16dtm in terms of normal AGN variability
B17 discarded the possibility that PS16dtm can be explained in terms of normal AGN variability.After six years of monitoring this nuclear transient, we can re-open the same question.The variability of an AGN can be assessed through the mean fractional variability parameter F var (Peterson 2001), that can be approximated with the ratio of the variance and mean value.NLSy1 are known to be low-level variability AGN (Ai et al. 2010) also on longer time-scales (Shapovalova et al. 2012), except for the radio-loud NLSy1 which exhibit blazar like variability (Berton et al. 2015).For PS16dtm, the F var is 0.43, 0.26, 0.58 for Hβ, Hα and total broad Fe II, respectively.Again, we can confirm, with data that spans over six years that such high variability up to 50% indicates that the outburst is not due to regular AGN activity.
In Fig. 17, we investigated the location of the PS16dtm outburst on the AGN main-sequence plane, defined by two measured quantities: the FWHM of Hβ and R FeII (Sulentic et al.During 1600 days of observations, R FeII and FWHM(Hβ) are slowly decreasing toward the preoutburst state, however the location on the main-sequence remains within the extreme accreators, far from the majority of AGN (indicated by SDSS quasars from Shen et al. (2011) catalog, gray dots) and extreme end of NLSy1 (data from Rakshit et al. (2017) catalog, blue dots).Such strong R FeII indicates high-accretion rates, most likely super-Eddington accretion (Marziani et al. 2018;Panda et al. 2018).This is further supported with our estimated L bol /L Edd ratio, shown in the bottom panel of Fig. 10.The strength of Fe II emission in PS16dtm is remarkably high, even when we consider that there are other TDE candidates with Fe II emission lines such as AT 2018fyk (Wevers et al. 2019) and PS1-10adi (Kankare et al. 2017;He et al. 2021), with their spectroscopic comparison is shown in Fig. 7. What is striking is that the R FeII is ∼5 times larger than most AGN from the literature.We note that fitted spectra are photometrically calibrated, thus, the absolute fluxes should be taken with caution, however, the R FeII is representative of the line ratio, therefore it should be more reliable.

Dust echo from MIR emission and the time lags with respect to the optical peak
There are several optically discovered TDE candidates with MIR flares (see e.g., Jiang et al. 2021).They are detected with a delay with respect to the optical ones, so they are interpreted as dust echoes of the TDE optical emission (van Velzen et al. 2016).
The "reverberation time lag," due to light-travel times, between the optical/UV/X-ray variations and the IR response is a commonly used technique in AGN science.Such a time delay can be directly inferred by the light curves and it is presumed to provide an estimate of the distance to the outer radius of the dust torus.However, many parameters influence the IR reverberation response, including the convolution of the UV/optical light curve with a transfer function that contains information about the geometry and structure of the torus (see Almeyda et al. 2017, for a detailed discussion).
For galaxies with AGN, it is possible to estimate the inner radius of the dusty torus from the sublimation radius as (see e. where L 45 is the bolometric luminosity of the AGN in units of 10 45 erg s −1 , a 0.1 is the grain size in units of 0.1 µm, and T 1800 is the temperature of the dust, normalized to the expected sublimation temperature of 1800 K.For PS16dtm, it is possible to make this calculation using both the pre-and post-outburst bolometric luminosity, which will yield different values, corresponding to PS16dtm clearing out a larger radius by sublimating the dust closer to the nucleus (J17).Using our estimates of the bolometric luminosity (1.6 × 10 43 erg s −1 and ∼ 10 45 erg s −1 pre-and postoutburst), we calculated the inner torus radius to change from ∼18 and ∼144 light days, respectively.It is worth noting that several works have found that the innermost torus radii based on dust reverberation were systematically smaller than the theoretical prediction of Equation 1 by a factor of few (see the recent review van Velzen et al. 2021).
The MIR light curve was steeply rising already at +9 days with respect to the estimated optical outburst (Fig. 2), although the W1 − W2 color at this early phase is bluer and the MIR could be consistent with the Rayleigh-Jeans tail of the UV/optical blackbody (Fig. 5).The MIR light curve is still rising at the end of our monitoring campaign, albeit slowly (see Fig. 2) suggesting that the time delay with respect to the optical must be at least 1500 rest-frame days, implying for the outer radius of the dust medium to be R 4 × 10 18 cm, or 1.3 pc.Some authors have used the duration of the MIR echo as a discriminator between the CLAGN and TDE scenario (see e.g., Kool et al. 2020), suggesting that CLAGN MIR echos evolve on much longer timescales of several years (see the work by Sheng et al. 2017Sheng et al. , 2020)), compared to TDE echos that evolve on timescales of months.In more recent works that gathered even more TDE and CLAGN candidates, this distinction is not clear anymore.This is because there are CLAGN candidates that have shorter dust echos and vice versa, there are TDE candidates with longer MIR echos.Jiang et al. (2021) looked at 23 TDE candidates from the literature and found 11 with NEOWISE emission; the time of the MIR peak ranges from 0 to 800 days after the optical peak, and their duration are from tens to more than 1000 days.In another TDE in a luminous infrared galaxy, a MIR echo lasting for at least 13 years and still going today, is observed (Mattila et al. 2018).Lyu et al. (2022) looked at the WISE data of 13 CLAGN and found the time lags of the variation in the mid-IR behind that in the optical band for 13 CLAGN with strong mid-IR variability, from tens to hundreds of days.We have also inspected the NEOWISE light curve of AT 2018dyk, which was classified as CLAGN by Frederick et al. (2019), despite that initially was suspected as TDE candidate (Arcavi et al. 2018).The AT 2018dyk MIR light curve has rosen more that a magnitude above the host level and it is still declining at more than 4 years after the peak.The nuclear transient Gaia16aax which could be explained by both the TDE scenario or some variation in the acretion rate of the active SMBH, showed 140 days of time lag of the NEOWISE light curve compared to the optical (Cannizzaro et al. 2020).In conclusion, at the moment the duration of the MIR echo cannot be used as a discriminator between the physical cause of the sudden optical flare.
Another proposed discriminator between TDEs and CLAGN is the covering factor of the dust, which can be determined from the ratio of the total energy radiated by the MIR by the energy absorbed by the dust (Gezari 2021).This has been justified by the work of Jiang et al. (2021), who found values on the order of 1% for TDEs in quiescent galaxies13 , while the characteristic values of CLAGN covering factors are two orders of magnitude higher.The dust covering factor of PS1-10adi is ∼ 40% (Jiang et al. 2019), while for PS16dtm this is estimated to be > 10% (J17, Jiang et al. 2021) as the MIR peak has not been reached yet.We note that, since it is possible that TDEs in AGN have similar covering factors with CLAGN, the dust covering factor is not a discriminator for TDE vs. CLAGN scenario, but perhaps can indicate the presence of preexisting dusty torus in the cases where the AGN has not been already known.
In the case of PS16dtm we measure a time lag between the peak of the broad Fe II component and the continuum light curves (Fig. 10).The nominal value of this lag is about ∼300 days using the nominal peak of the Fe II (∼400 days) and the continuum light curves (70 days for the first peak), although this is uncertain and could could be smaller as the timescales between 200 -400 days were not well probed by our spectroscopy.This could be consistent with the inner radius of the dust torus after PS16dtm sublimated the dust closest to the black hole, and the idea proposed by He et al. (2021) for PS1-10adi that Fe II emission is related to the torus inner radius and to Fe that was initially locked in dust, which was evaporated by the TDE.Given the large uncertainties in the time scales, it is not possible to directly confirm their hypothesis.We therefore investigated among TDEs with reported Fe II emission (see Fig. 7) and we found that nearly all of them have strong associated MIR flares (Fig. 18), with the exception of AT 2018fyk, which is also the one with the weakest Fe II.This is seems to strengthen the potential relation between Fe II and dust, which needs to be investigated further in the future.We note, however, that the Fe II and the MIR emissions in PS16dtm evolve on very different time scales (dropping after 400 days and increasing steadily for 1500 days), so it is not obvious how to establish a straightforward relation between them.
We finally note that a time lag is also seen between the continuum and Hα (see Fig. 10), similar to what has been observed for other TDEs, but a factor of 10 longer (Charalampopoulos et al. 2022).This time lag is more evident than for Fe II and might or might not be of similar value, with both light curves peaking ∼ 400 days after the outburst.

AGN extreme variability and the TDE scenario
The nuclear outburst PS16dtm may shed light on an important question on what triggers or maintains the AGN, especially in the cases where extreme variability has been observed.In fact, TDEs is one of the suggested processes that can fuel gas to the SMBH (Hills 1975;Frank & Rees 1976), and it could work for the low-luminosity end of AGNs where the average accretion rate from tidal disruptions is enough to account for the luminosity (see also Combes 2001, for a review on this subject).After the publication of B17, ∼ 200 CLAGN have been identified, with the detection of the event after it has already occurred, so detailed studies have not been possible (Yang et al. 2018;Frederick et al. 2019;MacLeod et al. 2019;Graham et al. 2020;Sánchez-Sáez et al. 2021;López-Navas et al. 2022;Green et al. 2022;Hon et al. 2022).The classification into a CLAGN is based on the extreme changes in magnitude or emission line intensities, with a disappearance or appearance of the broad component in emission lines (see e.g., Lyutyj et al. 1984;Kollatschny & Fricke 1985;Oknyansky et al. 2019;Ilić et al. 2020).For example, MacLeod et al. ( 2019) searched for |∆g| > 1 mag, |∆r| > 0.5 mag variability over any of the available time baselines probed by the SDSS and Pan-STARRS1 surveys.The timescales of the variability of these CLAGN cannot be pinpointed precisely since they do not have good coverage, but typical timescales ranges from a roughly a year to 13 years in the rest-frame.
There is still a vivid debate about what could cause the extreme variability in AGN, with many authors leaning toward the change in the accretion rate as the most plausible scenario (Noda & Done 2018;Sniegowska et al. 2020), despite that it is difficult to explain the cases with rapid (less than a year) rise in the light curve, as the time scales that govern the dynamics of an accretion disk are much longer in the classical theoretical framework of accretion disk physics (Cannizzaro et al. 2020).Other possible explanations are the obscuration of the broad-line region (e.g., Elitzur 2012), though questioned by many works (e.g., Stern et al. 2018;Mehdipour et al. 2022), close binary supermassive black holes (Wang & Bon 2020), and supernova explosions (Campana et al. 2015).Growing number of studies are suggesting that some CLAGN are due to change in the accretion state in the SMBH triggered by a TDE (Komossa et al. 2017;Trakhtenbrot et al. 2019;Zhang 2021).In fact, the observational signatures of CLAGN can be similar to TDEs that happen in galaxies with AGN, in sense that there is a transition phase that can be accompanied by a drastic change in the AGN continuum flux, so the optical and MIR become brighter when the CLAGN turn on.
PS16dtm could have been classified as an extreme CLAGN with the change of brightness of several magnitudes (see Figs. 1  and 2) according to the criteria in MacLeod et al. ( 2019), but it would have been clearly distinguished from an ordinary AGN variability (e.g., Vanden Berk et al. 2004).Given that PS16dtm was discovered with an all-sky survey which is biased toward the live selection of the most variable object, it is natural that PS16dtm may not be representative of a large population.However, the rise of PS16dtm luminosity on such short time-scales (< 50 days), makes it difficult to explain in terms of the intrinsic change to the AGN (e.g., the changes in the accretion disk structure, see e.g., Stern et al. 2018;Cannizzaro et al. 2020).From the spectroscopic point of view, the boost in the intensities of broad emission lines, as well as the width of the broad lines, may have also classified as CLAGN.The most striking change is the boost of Fe II emission, especially the appearance of the broad Fe II emission, that completely disappears in time (or gets so weak that it blends with the underlying continuum emission).Interestingly, the nuclear transient SDSS J123359.12+084211.5 in NLSy1 reported by MacLeod et al. (2019) and classified as a CLAGN, also shows a remarkable change in the Fe II emission.However, that event is a noticeable outlier in the CLAGN sample, especially based on the large Eddington ratio compared to the other CLAGN which are typically seen in lower Eddington ratio AGN (MacLeod et al. 2019).Despite that the data of SDSS J123359.12+084211.5 are very limited, it is quite possible that PS16dtm and SDSS J123359.12+084211.5 are powered by the same physical mechanism.We conclude that PS16dtm, without the extensive follow-up campaign shown here, might easily have been easily classified as a CLAGN, especially if the rising part of the light curve was missed.This implies that some CLAGN, such as SDSS J123359.12+084211.5, may actually be powered by a TDE.

Conclusions
Here, we studied PS16dtm, which is a nuclear outburst in a NLSy1 galaxy, with photometric and spectroscopic data that spans up to six years after its discovery.These are the main conclusions from these observations: -The NUV/optical light curve (see Fig. 1), after an initial rise of ∼ 50 days, displays another peak after ∼ 50 days.After a plateau phase, it started to decay slowly.In the first ∼ 270 days, the NUV light curve scales with a t −5/7 decay, then afterward with a shallower t −1/6 decline.The NUV/optical photometry shows very little color evolution.A single blackbody model is not a good fit to the NUV/optical data.We also attempted to correct for extinction with a SMC-like extinction law, for which some authors have suggested that can be applied to the line of sight to AGN, but this did not yield satisfactory results.-The MIR light curve which precedes the first ATLAS optical detection by ∼ 4 days, has steeply risen in the first ∼ 500 days and it is still slowly rising ∼ 1700 days after the first outburst, but the maximum value has not been reached already.The MIR color in the first epochs compared to the later epochs are quite different.A blackbody fit to the MIR data indicates that the temperature during the rising part of the light curve was ∼ 2300 K, it then dropped to ∼ 1200 K, and then settled to ∼ 900 K. -Our spectra taken between +155 and +1868 days past the outburst, show the following main properties.The blue continuum has dropped over time, the spectral evolution over the years proceeded with a very slow rate and almost returned to preoutburst level in our last spectrum.The most striking spectroscopic features, except the broad Balmer lines, are the Fe II emission lines.We highlight that this broad Fe II emission is transient, in the sense that it was not present in the archival preoutburst spectrum and almost completely disappeared in our last spectrum at +1868 days after the outburst.
We made a spectroscopic comparison with other TDEs with reported broad Fe II emission in the optical, and conclude that PS16dtm is the strongest iron emitter so far.We have also investigated the NEOWISE light curves of these TDE candidates and found that three out of four of them have also shown exceptional MIR flaring.-Thanks to our semi-empirical model which we developed here, we were able to identify and study the evolution of the emission lines.The flux of the broad Hα, Hβ, broad Fe II and narrow Fe II lines, after a rise in the first ∼ 200 days after the outburst, their strength dropped over time after ∼ 400 days after the outburst.We also found that the majority of the flux of the broad Balmer and Fe II lines is consistent with being produced by photoionization at least starting 200 days past the outburst.-We have also studied the MIR light curve and the light curve of the broad Fe II emission lines in order to explore the surrounding region.Given that the quasi-simultaneous rise of the MIR, optical and Fe II emission, it is tempting to assume that the regions from which they originate are nearby and that even the Fe II emitting region might be in the inner radius of the AGN torus.On one hand, there are the broad Hα, Fe II emission line that reach their peaks at ∼ 350 days after the optical data, even though this is uncertain given the gaps in the data.On the other hand, we estimated that the sublimation radius of the host AGN is ∼ 18 light days.In addition, after the rise, the MIR, optical and Fe II light curve evolve on different timescales.Therefore, it was not possible to establish a direct link between the dust and the iron emitting regions.-We found a weak coronal [Fe VII] 6087 line in our X-shooter spectrum at +1868 days.There are possibly other coronal lines, but their detection is not certain.-B17 predicted that the X-ray emission, which dimmed after the PS16dtm outburst, will reappear when the TDE accretion rate declines, that should be approximately a decade after the outburst.We found only weak X-ray emission in the 0.5-8 keV band at the location of PS16dtm, at +848, +1130 and +1429 days after the outburst.This is consistent with the levels observed by B17, so this reappearance has still not occurred yet, despite that the optical photometry and spectroscopy indicate that the emission is returning to the preoutburst level.-We have reexamined the host galaxy properties and found that the host galaxy, beyond the AGN emission, has also an important contribution from the ongoing star formation.
We have further discussed whether it is possible to use the duration of the dust echo or the dust covering factor as a indicator of the CLAGN or TDE scenario, and concluded that at this point this is not feasible given the overlap of the values of the two.We have reopened the question of whether PS16dtm can fit within the category of what is considered a "normal" AGN variability, and found that the answer is no.However, the extensive PS16dtm data shown here, suggests that perhaps some AGN extreme variability can be explained as a change in the accretion state in the SMBH triggered by a TDE.
The Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), will provide a sizeable sample of nuclear outbursts, such as TDEs, supernovae, and variable AGN, during its 10 years of surveying the sky (see e.g., Bricman & Gomboc 2020;Roth et al. 2021).Therefore, a large-scale systematic study in real time from the newly discovered nuclear candidates with the Vera Rubin Observatory will be needed in order to make progress in understanding CLAGN and TDEs.

Fig. 2 .
Fig. 2. Optical (ATLAS and PSST) and mid-infrared NEOWISE light curves of PS16dtm.The photometry has been extinction corrected for Milky Way dust and host subtracted.The limits from ASSASN in V band are also shown.The WISE error-bars are smaller than the symbols.The inset shows a zoom-in at the epochs around first detection.

Fig. 3 .
Fig.3.Temporal evolution of the spectra taken with NTT/EFOSC2 (and the final X-shooter spectrum at +1868 days).The spectra are dominated by the Balmer series and a plethora of Fe II lines.The spectra are colored as indicated in the legend and the phases refer to rest-frame days after MJD 57577.Vertical solid lines indicate the position of hydrogen Balmer lines whereas the shaded pink areas mark the location of the strongest Fe II templates.

Fig. 4 .
Fig. 4. Evolution of the bolometric luminosity, temperature and radius for PS16dtm, by fitting a blackbody to the Galactic extinction and host-corrected photometry, excluding the photometry in the three NUV UVOT filters and requiring at least 4 photometric measurements in each epoch.The purple line in the upper insert represents the best fit of the fading light curve with a power-law profile, where the power-law index is fit freely, L = L o (t − t o ) −α .

Fig. 5 .
Fig.5.Examples of three epochs of blackbody fits to the photometry that includes the NEOWISE data.In the left panel, the best fit by using one blackbody is shown.For the first epoch of the MIR detection, together with the quasi-simultaneous ATLAS photometry, the best fit temperature is ∼ 8900 K.A better fit is obtained by using a double blackbody model, and the resulting blackbody temperatures are shown.In the right panel, the fits of the two-blackbody model are shown to the photometry of two epochs, +165 and +835 days where the NUV bands are plotted, but have not been included in the fitting procedure.

Fig. 6 .Fig. 7 .
Fig. 6.Zoomed-in view of the spectra containing the H β and the Fe II multiplets.The phase indicated in the legend refer to the rest-frame days after MJD 57577.Vertical solid lines indicate the position of hydrogen Balmer lines whereas the shaded pink areas mark the location of the strongest Fe II features.

Fig. 8 .Fig. 9 .
Fig.8.Examples of multicomponent spectral fitting in the 4100 − 7000 Å rest-frame range for the +179 days spectrum (solid gray line).The underlying continuum emission and all emission-line components are plotted with different colors, as indicated in the legend.For details on the assumed emission-line components see text.The final model (solid red line) which is the sum of all components is also shown, whereas the residual spectrum, is plotted below.

Fig. 10 .Fig. 11 .
Fig. 10.Temporal evolution of the extracted emission lines and lineblends during the campaign.The bottom panel shows the evolution of L bol /L Edd ratio, in which the bolometric luminosity is estimated from the spectral continuum at 5100 Å.The vertical dashed line indicates +400 days after the estimated outburst time to guide the eye.

Fig. 12 .
Fig. 12. Correlations between the continuum flux at 5800 Å, and the flux of Hα (left), Hβ (middle), and total broad Fe II blend (right) in units of 10 −15 erg s −1 cm −2 .The color-bar indicates the day after MJD 57577.The clear correlation of the continuum and line fluxes (apart from the first three epochs) supports the photoionization as the main heating source (see text for details).

Fig. 13 .
Fig. 13.The upper panels show the Balmer Hβ, Hα, and Hγ lines from the medium-spectral resolution VLT/X-shooter spectrum.The lower panels zoom in at the strongest identified Fe II and [Fe II] lines, as well as He II, [O I] and [S II].The weak coronal line [Fe VII] λ6087 is also detected.

Fig. 14 .
Fig.14.Fits of the H β (top) and H α (middle) line regions in the Xshooter spectrum.The bottom panel compares the broad components (obtained after subtracting the narrow lines) of the H β and H α in the velocity space.

Fig. 15 .
Fig.15.X-shooter spectrum, scaled to the SDSS spectrum of poorer spectral resolution, compared to the SDSS AGN component spectra (from Fig.11), which was scaled-up so that the broad line wings fit to the broad line of Hα (left) and Hβ (right).

Table 1 .
Spectroscopic observations of PS16dtm Notes.* In the rest frame with respect to the estimated time of outburst MJD 57577.
Kovačević et al. (2010) transition oscillatory strengths.Kovačević et al. (2010)produced a multicomponent template covering 4000-5500 Å. Building up on that work, we extended the Fe II model to include the wavelength range up to 7000 Å to cover the area around Hα line where Fe II emission is also very strong, using the atomic data from Kurucz database12 11https://fantasy-agn.readthedocs.io/en/latest/ are fixed to have the same width and shift; from other strong broad lines in this range, we included He I λ5876, He II λ4686, the Na I doublet λλ5890, 5896, and O I λ6046; (ii) the model of broad Fe II, with all lines having the same widths and shifts, and line ratios connected as described above.The broad Fe II lines are assumed to have the same profiles as broad Balmer lines (see e.g., Dong et al. 2011, and references therein), so this component is set to have the same width and shift as other broad lines; (iii) the AGN narrow emission lines, such as Hα, Hβ, Hδ, [O III], [N II], [S II], Ti II, Cr II (see Véron-Cetty et al. 2006; Park et al. 2022, for the list of significant narrow lines) are constrained to have the same width and shift.In addition to fixing the ratio of [O III] and [N II] doublets to the theoretical values of ≈ 3 (Dimitrijević et al. 2007); (iv) the model of narrow Fe II, plus the forbidden Fe II lines, have the same width as narrow lines (Véron-Cetty et al. 2006; Park et al. 2022).
4.3.4.Results of the modeling of the spectra in the 4100 − 7000 Å range