11institutetext: Center for Astrophysics and Cosmology, University of Nova
Gorica, Vipavska 11c, 5270 Ajdovščina, Slovenia.
11email<EMAIL_ADDRESS>22institutetext: DTU Space, National Space
Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby,
Denmark.
22email<EMAIL_ADDRESS>33institutetext: Department of Astronomy,
University of Belgrade - Faculty of Mathematics, Studentski trg 16, 11000
Belgrade, Serbia.
33email<EMAIL_ADDRESS>44institutetext: Hamburger Sternwarte,
Universitat Hamburg, Gojenbergsweg 112, 21029 Hamburg, Germany.
55institutetext: Tuorla Observatory, Department of Physics and Astronomy,
University of Turku, FI-20014 Turku, Finland 66institutetext: IAASARS,
National Observatory of Athens, 15236 Penteli, Greece 77institutetext:
Department of Astrophysics, Astronomy & Mechanics, Faculty of Physics,
National and Kapodistrian University of Athens, 15784 Athens, Greece.
88institutetext: Faculty of Natural Sciences and Mathematics, University of
Banjaluka, Mladena Stojanovića 2, 78000 Banjaluka, Bosnia and Herzegovina.
99institutetext: Oskar Klein Centre, Department of Physics, Stockholm
University, SE-10691 Stockholm, Sweden. 1010institutetext: Department of
Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064,
USA. 1111institutetext: European Southern Observatory, Alonso de Córdova 3107,
Casilla 19, Santiago, Chile. 1212institutetext: Institute for Astronomy,
University of Hawaii, 2680 Woodlawn Drive, Honolulu HI 96822, USA.
1313institutetext: Oskar Klein Centre, Department of Astronomy, Stockholm
University, AlbaNova, SE-10691 Stockholm, Sweden. 1414institutetext:
Astronomical Observatory, Volgina 7, 11060 Belgrade, Serbia.
1515institutetext: School of Physics, Trinity College Dublin, The University
of Dublin, Dublin 2, Ireland. 1616institutetext: Astronomical Observatory,
University of Warsaw, Al. Ujazdowskie 4, 00-478 Warszawa, Poland.
1717institutetext: Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer
de Can Magrans, s/n, E-08193 Barcelona, Spain. 1818institutetext: Institut
d’Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain.
1919institutetext: Birmingham Institute for Gravitational Wave Astronomy and
School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT,
UK. 2020institutetext: INAF-Osservatorio Astronomico d’Abruzzo, via M. Maggini
snc, I-64100 Teramo, Italy. 2121institutetext: Cosmic DAWN centre, Niels Bohr
Institute, University of Copenhagen, Rådmandsgade 62-64, 2200, Copenhagen,
Denmark. 2222institutetext: Astrophysics Research Centre, School of
Mathematics and Physics, Queen’s University Belfast, Belfast, BT7 1NN,
Northern Ireland, UK. 2323institutetext: School of Physics, O’Brien Centre for
Science North, University College Dublin, Belfield, Dublin 4, Ireland

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

T. Petrushevska The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm G.
Leloudas The rise and fall of the iron-strong nuclear transient PS16dtmThe
rise and fall of the iron-strong nuclear transient PS16dtm D. Ilić The rise
and fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtm M. Bronikowski The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm P. Charalampopoulos The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtmThe rise
and fall of the iron-strong nuclear transient PS16dtm G. K. Jaisawal The rise
and fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtm E. Paraskeva The rise and fall of the
iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm M. Pursiainen The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm N. Rakić
The rise and fall of the iron-strong nuclear transient PS16dtmThe rise and
fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtm S. Schulze The rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtm K. Taggart The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm C. K. Wedderkopp The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm J. P. Anderson The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm T. de
Boer The rise and fall of the iron-strong nuclear transient PS16dtmThe rise
and fall of the iron-strong nuclear transient PS16dtm K. Chambers The rise and
fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtm T. W. Chen The rise and fall of the
iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtm G. Damljanović The rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtm M. Fraser The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm H. Gao The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm A.
Gomboc The rise and fall of the iron-strong nuclear transient PS16dtmThe rise
and fall of the iron-strong nuclear transient PS16dtm M. Gromadzki The rise
and fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtm N. Ihanec The rise and fall of the iron-
strong nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm K.
Maguire The rise and fall of the iron-strong nuclear transient PS16dtmThe rise
and fall of the iron-strong nuclear transient PS16dtm B. Marčun The rise and
fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtm T. E. Müller-Bravo The rise and fall of
the iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm M. Nicholl The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm F. Onori
The rise and fall of the iron-strong nuclear transient PS16dtmThe rise and
fall of the iron-strong nuclear transient PS16dtm T. M. Reynolds The rise and
fall of the iron-strong nuclear transient PS16dtmThe rise and fall of the
iron-strong nuclear transient PS16dtm S. J. Smartt The rise and fall of the
iron-strong nuclear transient PS16dtmThe rise and fall of the iron-strong
nuclear transient PS16dtm J. Sollerman The rise and fall of the iron-strong
nuclear transient PS16dtmThe rise and fall of the iron-strong nuclear
transient PS16dtm K. W. Smith The rise and fall of the iron-strong nuclear
transient PS16dtmThe rise and fall of the iron-strong nuclear transient
PS16dtm T. Wevers The rise and fall of the iron-strong nuclear transient
PS16dtmThe rise and fall of the iron-strong nuclear transient PS16dtm Ł.
Wyrzykowski The rise and fall of the iron-strong nuclear transient PS16dtmThe
rise and fall of the iron-strong nuclear transient PS16dtm

(Received XX; accepted XX)

###### Abstract

Context. Thanks to the advent of large-scale optical surveys, a diverse set of
flares 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
difficult. One of the obstacles of making progress in the field is the lack of
well-sampled data of long-lived nuclear outbursts in AGN.

Aims. Here, we study PS16dtm, a nuclear transient in a Narrow Line Seyfert 1
(NLSy1) galaxy, which has been proposed to be a TDE candidate. Our aim is to
study the spectroscopic and photometric properties of PS16dtm, in order to
better understand the outbursts originating in NLSy1 galaxies.

Methods. Our extensive multiwavelength follow-up that spans around $2000$ days
includes photometry and spectroscopy in the UV/optical, as well as mid-
infrared (MIR) and X-ray observations. Furthermore, we improved an existing
semiempirical model in order to reproduce the spectra and study the evolution
of the spectral lines.

Results. The UV/optical light curve shows a double peak at $\sim 50$ and $\sim
100$ days after the first detection, and it declines and flattens afterward,
reaching preoutburst levels after 2000 days of monitoring. The MIR light curve
rises almost simultaneously with the optical, but unlike the UV/optical which
is approaching the preoutburst levels in the last epochs of our observations,
the MIR emission is still rising at the time of writing. The optical spectra
show broad Balmer features and the strongest broad Fe II emission ever
detected in a nuclear transient. This broad Fe II emission was not present in
the archival preoutburst spectrum and almost completely disappeared +1868 days
after the outburst. We found that the majority of the flux of the broad Balmer
and Fe II lines is produced by photoionization. We detect only weak X-ray
emission in the 0.5-8 keV band at the location of PS16dtm, at +848, +1130, and
+1429 days past the outburst. This means that the X-ray emission continues to
be lower by at least an order of magnitude, compared to archival, preoutburst
measurements.

Conclusions. We confirm that the observed properties of PS16dtm are difficult
to reconcile with normal AGN variability. The TDE scenario continues to be a
plausible explanation for the observed properties, even though PS16dtm shows
differences compared to TDE in quiescent galaxies. We suggest that this event
is part of a growing sample of TDEs that show broad Balmer line profiles and
Fe II complexes. We argue that the extreme variability seen in the AGN host
due to PS16dtm may have easily been misclassified as a CLAGN, especially if
the rising part of the light curve had been missed. This implies that some
changing look episodes in AGN may be triggered by TDEs. Imaging and
spectroscopic data of AGN with good sampling are needed to enable testing of
possible physical mechanisms behind the extreme variability in AGN.

##  1 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 X-ray 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; Frederick et al., 2021). The discoveries of TDEs in AGN were initiated
later compared to TDEs in quiescent galaxies, partly because, initially,
transients in galaxies with known AGN have been excluded in optical surveys
from follow-up investigation, so as to avoid being overwhelmed by spurious
candidates. Furthermore, TDEs arising in a galaxy with an existing AGN are
more difficult to identify simply because the nucleus itself can be a variable
source.

In recent years, transient phenomena in galaxies with AGN have been observed,
with extraordinary changes in photometric and spectroscopic properties on
short timescales (e.g., Frederick et al., 2021). However, the physical
mechanism to explain these phenomena is unclear. One example is the transient
AT 2018dyk, initially classified as a TDE (Arcavi et al., 2018); however, a
subsequent study argued that it is a CLAGN (Frederick et al., 2019). Another
example is the transient CSS100217 which Drake et al. (2011) argued was a
supernova, while Frederick et al. (2021) favored the AGN scenario.
Furthermore, Zhang et al. (2022) and Cannizzaro et al. (2022) suggested that
it could be explained as a TDE, after gathering data for more than 10 years.
The reason for this ambiguity is that, at present, there is no single
observational signature that disentangles, unambiguously, the physical
mechanism that causes the luminous phenomena (see Zabludoff et al., 2021, for
a recent review). Since a TDE happening in AGN can also cause intrinsic change
in the accretion disk, it can also provide one of the channels to explain
CLAGN (see e.g., Kool et al., 2020; Cannizzaro et al., 2020; Zhang, 2021;
Hinkle et al., 2022). An ongoing issue in understanding long-lived nuclear
outbursts in AGN is, on the one hand, the low numbers of such events, and on
the other hand, often the data are sparsely sampled from the long-lived
nuclear outbursts, so it does not allow for one to study them in detail.

Here, we focus on the nuclear transient, PS16dtm, for which we gathered an
extensive dataset over a time scale of almost six years. PS16dtm was first
classified as supernova Type IIn (Terreran et al., 2016; Dong et al., 2016),
but further in-depth analysis based on follow-up observations concluded that
it is more consistent with a TDE interpretation (Blanchard et al., 2017,
hereafter B17). B17 presented observations up to 200 days after it was
detected. In this paper, we present photometric and spectroscopic data for a
period of $\sim 2000$ days after the first detection of PS16dtm. As we
subsequently show in this paper, PS16dtm is perhaps the nuclear transient with
the most dramatic increase of Fe II emission after the outburst, although
there are other examples, albeit not with such strong Fe II emission, such as
J123359.12+084211.5 (MacLeod et al., 2019), AT 2019dsg (Cannizzaro et al.,
2021), AT 2018fyk (Wevers et al., 2019), CSS100217 (Drake et al., 2011), and
PS1-10adi (Kankare et al., 2017). We also show that PS16dtm exhibits strong
MIR emission similar to other nuclear flares in hosts with AGN. Intriguingly,
the X-ray emission continues to be dimmed compared the archival preoutburst
state, despite that UV/optical photometry and spectroscopy are approaching the
preoutburst levels.

The paper is organized as follows. In Sect. 2, we summarize the findings on
PS16dtm that have been previously published. In Sect. 3 we present our
observations, while in Sect. 4 we perform the analysis. The discussion is
found in Sect. 5 and finally, in Sect. 6 we present our conclusions. We use
the AB magnitude system throughout this work. Furthermore, we assume a
cosmology with $H_{0}=67$ km s-1 Mpc-1, $\Omega_{m}=0.32$ and
$\Omega_{\Lambda}=0.68$ as in B17, which implies a luminosity distance of 381
Mpc to PS16dtm. The Galactic extinction along the line of sight to PS16dtm is
$\rm A_{V}=0.069\pm 0.001$ mag (Schlafly & Finkbeiner, 2011).

##  2 Summary of previous work on PS16dtm

In this section, we summarize previously published work on PS16dtm.111PS16dtm
is registered as SN 2016ezh on the Transient Name Serverhttps://www.wis-
tns.org/object/2016ezh. We use the discovery name throughout this paper as the
supernova interpretation has been disfavored. 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 $\sim 10^{6}M_{\odot}$, 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 $\sim 50$ days, with two recognizable peaks
distanced $\sim 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
($\sim 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 ($\sim 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 $\sim 10$ light days to $\sim 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 ($\sim 10,000$ $\rm 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 $\sim 110$ days, while PS16dtm has
stayed bright above the host-galaxy baseline for much longer, as will be shown
in the next sections.

##  3 Observations

In this section, we present new observations of PS16dtm out to $\sim 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 $\sim 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
archive222https://www.wiserep.org (Yaron & Gal-Yam, 2012).

###  3.1 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 $\AA$, $B$ at 4392 $\AA$, $U$ at 3465 $\AA$), and
three near-UV ($UVW1$ at 2600 $\AA$, $UVM2$ at 2246 $\AA$, and $UVW2$ at 1928
$\AA$). 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 HEAsoft333https://heasarc.gsfc.nasa.gov/docs/software/lheasoft/.
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 $UVM2$.

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 tool 444https://www.swift.ac.uk/user_objects/ 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\times$10-14
erg s-1 cm-2 and $L_{X}<1.7\times$1041 erg s-1).

Figure 1: 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.

###  3.2 Optical ground-based photometry

####  3.2.1 ATLAS

ATLAS is a robotic survey that uses at least two identical 50-centimeter
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$ $\AA$, 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 $\sim
6$ days before the first detection reported by ASSA-SN and $\sim 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 $\sim 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.

####  3.2.2 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_{\rm
PS1}$ and $i_{\rm PS1}$ bands, which have wavelength ranges 4000-8300 $\AA$
and 6800-8300 $\AA$, 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).

####  3.2.3 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:O555https://telescope.livjm.ac.uk/TelInst/Inst/IOO/ 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.

###  3.3 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 $\mu m$ ($W2$)
(Mainzer et al., 2014). We retrieved the photometry from the IRSA public data
archive666https://irsa.ipac.caltech.edu/cgi-bin/Gator/nph-
scan?mission=irsa&submit=Select&projshort=WISE. 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. The first detection by NEOWISE is on MJD 57587 (2016 July 18),
it is $\sim 4$ days before the first ATLAS detection of PS16dtm which was made
at MJD 57591.6 and it is the earliest optical detection, to our knowledge.
With the ATLAS data reported here, the gap between the first PS16dtm optical
and MIR outburst are closer that previously shown in J17, which presented the
ASSA-SN data only.

Figure 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.

###  3.4 NTT/EFOSC2 spectroscopy

We carried out observations of PS16dtm within the Public ESO Spectroscopic
Survey for Transient Objects (PESSTO; Smartt et al., 2015), and its
continuations, ePESSTO and ePESSTO+. PESSTO uses the low-resolution ESO Faint
Object Spectrograph and Camera v.2 (EFOSC2) on the New Technology Telescope
(NTT) in La Silla Observatory, Chile. We have collected NTT/EFOSC2 spectra of
PS16dtm at 16 different epochs, using various settings including the standard
grims 11, 13 and 16 (Smartt et al., 2015). The spectra were reduced in a
standard manner 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
archive777https://www.wiserep.org (Yaron & Gal-Yam, 2012).

Figure 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. Table 1: Spectroscopic observations of PS16dtm

Date | MJD | Phase∗ | Wavelength range | Exposure time | Grism
---|---|---|---|---|---
2016-12-22 | 57744 | +155 | 3380–10320 | 1500 | gr11 & gr16
2017-01-17 | 57770 | +179 | 3380–10320 | 900 | g11 & gr16
2017-02-05 | 57789 | +196 | 3380–10320 | 900 | g11 & gr16
2017-08-20 | 57985 | +378 | 3380–10320 | 1500 | g11 & gr16
2017-09-11 | 58007 | +398 | 3380–10320 | 1500 | g11 & gr16
2017-10-19 | 58045 | +433 | 3380–10320 | 1500 | g11 & gr16
2017-11-09 | 58065 | +452 | 3380–10320 | 1500 | g11 & gr16
2017-12-13 | 58099 | +483 | 3380–10320 | 1500 | g11 & gr16
2018-08-17 | 58347 | +713 | 3380–10320 | 1500 | g11 & gr16
2018-09-16 | 58377 | +740 | 3685–9315 | 1500 | g13
2019-01-25 | 58508 | +862 | 3380–10320 | 1500 | g11 & gr16
2019-07-29 | 58693 | +1033 | 3380–7520 | 1800 | g11
2020-10-23 | 59145 | +1451 | 3380–10320 | 1800 | g11 & gr16
2020-12-30 | 59213 | +1514 | 3380–10320 | 1800 | g11 & gr16
2021-09-04 | 59460 | +1743 | 3380-10320 | 2700/1800 | g11 & 2xgr16
2022-01-16 | 59595 | +1868 | 3000-24800 |  | UVB, VIS, NIR
888∗ In the rest frame with respect to the estimated time of outburst MJD
57577.

###  3.5 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 $\AA$ 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).

###  3.6 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)999P.I. Blanchard P., ID 21460, 21461 and 22618 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$\sigma$ level in the first, second, and third Chandra
epochs with a count rate of about 4$\times$10-4, 2.3$\times$10-4, and
2.7$\times$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$\times$1020 cm-2 (Pons & Watson,
2014, B17), the photon index is found to be 0.8${}_{-0.1}^{+8.5}$ 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$\pm$0.8)$\times$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${}_{-0.5}^{+2.0}$
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$\pm$3)$\times$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$\times$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.

##  4 Analysis

###  4.1 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.0$ mag in the UVOT $V$
band, the luminosity dropped for $\sim 50$ days, after which it rose again to
a second peak and stayed at almost plateau level for $\sim 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 code101010https://github.com/nblago/utils 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\times 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 $\sim
15000$ K contemporaneous with the rise in the luminosity. After that, the
temperature drops to $\sim 10000$ K, and is held relatively constant. After
$\sim 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).111111We note that this implies the presence of small dust grains, in
contrast to what other authors have argued, that the largest grains are those
able to survive in the vicinity of the SMBH. SMC 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 $\AA$ 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 possible to fit the SED to a
blackbody with $T\sim 2-4\cdot 10^{5}$ K and a bolometric luminosity of $L\sim
10^{48}-10^{49}$ erg $s^{-1}$, at a corresponding blackbody radius of $R\sim
1\cdot 10^{15}$ cm, and SMC-like extinction curve with $E(B-V)$ $\sim 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, $L=L_{o}(t-t_{o})^{-\alpha}$. Our best fit
parameter, $\alpha\sim{5/7}$ for the NUV/UVOT light curves for the initial
($158<t<270$) days, while $\sim 270$ days after the outburst, it follows a
$\alpha\sim{1/6}$ decay. For the bolometric luminosity, the best-fit power-law
has the index $\alpha\sim{0.98}$ (see Fig. 4). Interestingly, a similar
decline trend is also seen in the bolometric light curve of the TDE AT 2017gge
($L\propto 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
($\rm M_{BH}<10^{6.5}M_{\odot}$), 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.

Figure 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})^{-\alpha}$.

###  4.2 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 $\sim 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 $\sim 2300$ K for the dust
blackbody temperature, while the second epoch at +165 days, yields $\sim 1300$
K, and the epochs afterward (from +349 to +1700 days) settle at $\sim
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.

Figure 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 $\sim
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.

###  4.3 Spectroscopic analysis

####  4.3.1 Summary of the spectral features presented in B17

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 $\rm H\alpha$ 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] $\lambda\lambda 7320$, 7330, [Ca
II] $\lambda\lambda 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 $\mu m$, beyond which
the spectrum is relatively featureless,

  * •

Paschen $\gamma$ and other lines of the Paschen series were detected. Paschen
$\alpha$ and $\beta$ 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 $\sim 4500$ $\AA$ the spectrum shows a very complex combination of
narrow features superimposed on broad features. They identify [O II]
$\lambda\lambda 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 $\lambda\lambda 2800$ absorption lines were detected,
with $\sim 10,000$ $\rm km\,s^{-1}$, that can perhaps indicate the presence of
an outflow.

B17 also studied the temporal evolution of the hydrogen Balmer lines $\rm
H\alpha$ and $\rm H\beta$. They fitted the spectral lines with the sum of
three Gaussians. They reported that the widths of the intermediate component
of $\rm H\alpha$ is around 750 $\rm km\,s^{-1}$ that narrows as a function of
time, going from 900 $\rm km\,s^{-1}$ to 600 $\rm km\,s^{-1}$. For the broad
component, they found 3500 $\rm km\,s^{-1}$ which increases to 4000 $\rm
km\,s^{-1}$ in the decline phase of the light curve. They found that the
narrow component is 100 $\rm 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 $\rm H\alpha$ has remained unchanged
relative to the host spectrum.

####  4.3.2 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 rest-frame days after the outburst. Perhaps even more
striking are the Fe II optical lines, in the blue (around H$\beta$) and red
(around H$\alpha$) line, that become weaker as time progresses (see the
zoomed-in view around the H$\beta$ 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.5 in MacLeod et al. (2019) and PS1-10adi (Jiang et al.,
2017; He et al., 2021), AT 2018fyk (Wevers et al., 2019) and AT 2019dsg
(Cannizzaro et al., 2021). 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
$\lambda$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 $\lambda$8446, is clearly identified in the +1514 spectrum (see
Fig. 3).

Figure 6: Zoomed-in view of the spectra containing the $\rm H_{\beta}$ 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. Figure 7: Spectroscopic comparison of PS16dtm with TDEs from
the literature where Fe II emission has been identified: PS1-10adi (Kankare et
al., 2017), AT 2018fyk (Wevers et al., 2019), and AT 2019dsg (Cannizzaro et
al., 2021). Vertical solid lines indicate the position of hydrogen Balmer
lines, the dashed line denotes He II $\lambda$4686, whereas the shaded pink
areas mark the location of the strongest Fe II complexes. The spectra of the
three comparison TDEs have been scaled up and offset for display.

####  4.3.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 (FANTASY)121212https://fantasy-
agn.readthedocs.io/en/latest/, 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 connected through the line 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 $\AA$ to cover the area around H$\alpha$ line
where Fe II emission is also very strong, using the atomic data from Kurucz
database131313https://lweb.cfa.harvard.edu/amp/ampdata/kurucz23/sekur.html
(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:

  1. (i)

the broad Balmer emission lines (H$\alpha$, H$\beta$, H$\gamma$) are fixed to
have the same width and shift; from other strong broad lines in this range, we
included He I $\lambda$5876, He II $\lambda$4686, the Na I doublet
$\lambda\lambda$5890, 5896, and O I $\lambda$6046;

  2. (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;

  3. (iii)

the AGN narrow emission lines, such as H$\alpha$, H$\beta$, H$\delta$, [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 $\approx 3$ (Dimitrijević et al., 2007);

  4. (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$ $\AA$
range

In Fig. 8 we show an example of the results of the multicomponent spectral
fitting in the $4100-7000$ $\AA$ 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$\alpha$ 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 $\chi^{2}$ parameter.

Figure 8: Examples of multicomponent spectral fitting in the $4100-7000$ $\AA$
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.
Figure 9: Evolution of the FWHM of different line components (denoted in the
upper right corner) during the spectroscopic campaign. Figure 10: Temporal
evolution of the extracted emission lines and line-blends during the campaign.
The bottom panel shows the evolution of $L_{\rm bol}/L_{\rm 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.

In Fig. 9 we show the time evolution of the FWHM of the most prominent
emission lines, broad H$\alpha$ and H$\beta$, 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$\alpha$ and H$\beta$ line, total broad Fe
II, total narrow and forbidden Fe II lines. The strength of broad line-blends
(H$\alpha$, H$\beta$, 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$\alpha$ and H$\beta$ from the low-resolution spectra is difficult; it
depends on the possibility for the code to identify the narrow [O III]
$\lambda$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_{\rm bol}/L_{\rm Edd}$. The Eddington luminosity was
estimated as $L_{\rm Edd}=1.26\times 10^{38}(M_{\rm BH}/M_{\odot})$
$\rm{erg\,s}^{-1}$ with the SMBH mass of $M_{\rm BH}\sim 10^{6}M_{\odot}$,
which we determine in Sect. 4.3.8. For the bolometric luminosity $L_{\rm 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_{\rm bol}=k_{\rm bol}\lambda L_{\lambda}$ by applying the mean quasar
bolometric correction $k_{\rm bol}\approx 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
Eddington 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 inter-calibration
based on the constant narrow-line flux, such as [O III] or [S II] lines (van
Groningen & Wanders, 1992; Fausnaugh, 2017).

Figure 11: The pure AGN components for H$\alpha$ (left) and H$\beta$ (right)
region are shown, obtained when the host-galaxy contribution is subtracted
from the observed total spectrum in the archival, preoutburst, SDSS spectrum.

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_{\rm
cnt}$), and the flux of H$\alpha$ (left), H$\beta$ (middle), and total Fe II
blend (right). $F_{\rm 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.

Figure 12: Correlations between the continuum flux at 5800 $\AA$, and the flux
of H$\alpha$ (left), H$\beta$ (middle), and total broad Fe II blend (right) in
units of $10^{-15}$ $\rm{erg\,s}^{-1}\rm{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).

####  4.3.5 Spectral features in the $7000-9000$ Å range

In the spectral region redward of the H$\alpha$ 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$\alpha$ red wing), He I
$\lambda$7065, He I $\lambda$7281\. There are also Ca II triplet
($\lambda$8498, $\lambda$8542, $\lambda$8662), well known oxygen lines O I
$\lambda$7774 and O I $\lambda$8446, as well as the Mg II $\lambda$7892 line
present. All of these broad features are of similar width as the broad
H$\alpha$ 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] $\lambda$7320, $\lambda$7330, [Ca II]
$\lambda\lambda$7291, 7324 by B17, but the authors suggested that these
feature could be associated to Fe II emission lines. However, we believe that
this a mixture of Fe II blend and He I $\lambda$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/X-shooter spectrum, with possibly only
weak Mg II and O I being present.

####  4.3.6 Analysis of the medium-resolution VLT/X-shooter spectrum

Figure 13: The upper panels show the Balmer H$\beta$, H$\alpha$, and $H\gamma$
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] $\lambda$6087 is also
detected.

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$\alpha$ and H$\beta$ lines (Fig. 13),
whereas the broad component of the Balmer H$\gamma$ 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 a6S, a4G, and b4F
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$\alpha$, H$\beta$, [O I] 6300, [S II]
6716, 6731, and [O II] 3727, 3729.

We fit the H$\alpha$ and H$\beta$ 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$\alpha$ 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
a6S Fe II lines ratios are fixed according to our model (see Sect. 4.3.3).

Figure 14: Fits of the $\rm H_{\beta}$ (top) and $\rm H_{\alpha}$ (middle)
line regions in the X-shooter spectrum. The bottom panel compares the broad
components (obtained after subtracting the narrow lines) of the $\rm
H_{\beta}$ and $\rm H_{\alpha}$ in the velocity space.

The width of Fe II lines is $\sim$220 $\rm km\,s^{-1}$, larger than the narrow
lines, which have widths that are closer to $\sim$100 $\rm km\,s^{-1}$. This
is evident for all iron lines, including the coronal [Fe VII] $\lambda$6087
line. The extracted broad H$\alpha$ and H$\beta$ components (obtained after
subtracting narrow lines) have similar profiles (Fig. 14, bottom panel), with
H$\alpha$ being slightly broader with the width of $\sim$1900 $\rm
km\,s^{-1}$. Both broad lines show slightly asymmetric profiles, with
possibility that the line peak is a slightly shifted redward.

####  4.3.7 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 X-rays are required
for their production. These lines are typically of somewhat larger width than
the narrow lines, with higher critical densities ($\sim 10^{7.5}{\rm 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] $\lambda$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] $\lambda$6087 is either not present or too weak to be
resolved in the low-resolution spectra. We also retrieved 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 near-by 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.

####  4.3.8 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$\alpha$ or H$\beta$ have the same broad line
intensity, as shown in Fig. 15. We then concluded that the broad line profiles
of H$\alpha$ or H$\beta$ in the archival SDSS spectrum and in the X-shooter
spectrum at +1868 days, are similar.

Using the approach from Sect. 4.3.4 where the bolometric luminosity is
expected to scale with the luminosity at 5100$\,\AA\,$, we estimated the AGN
bolometric luminosity from the X-shooter spectrum to be $L_{\rm bol}=3.7\times
10^{43}$ $\rm{erg\,s}^{-1}$, which is almost reaching the preoutburst value,
which we extracted from the SDSS spectrum to be $L_{\rm bol}=1.6\times
10^{43}$ $\rm{erg\,s}^{-1}$. Both values are host-corrected using the same
constant host contribution estimated from the SDSS spectrum.

Given that the broad H$\beta$ 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$\beta$
line of $1160\pm 190$ $\rm 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\,\AA$. The MBH 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_{\rm BH}=10^{6.07\pm 0.18}M_{\odot}$, 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).

###  4.4 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$\alpha$ line (Greene & Ho, 2007; Xiao et al., 2011) and high X-ray
luminosity of $L_{\rm 2-10keV}\sim 10^{42}$ $\rm{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 $\rm M_{\odot}$ 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$\beta$ (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] $\lambda$5007/H$\beta$, [N II]
$\lambda$6583/H$\alpha$, [O I] $\lambda$6300/H$\alpha$, [S II]
$\lambda\lambda$6717,31/H$\alpha$, and [O III] $\lambda$5007/[O II]
$\lambda\lambda$3726,3729 (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.

Figure 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$\alpha$
(left) and H$\beta$ (right). Figure 16: Diagnostic diagrams based on the line
ratios of [O III] $\lambda$5007/H$\beta$, [N II] $\lambda$6583/H$\alpha$, [O
I] $\lambda$6300/H$\alpha$, [S II] $\lambda\lambda$6717,31/H$\alpha$, and [O
III] $\lambda$5007/[O II] $\lambda\lambda$3726,3729 (Baldwin et al., 1981;
Kewley et al., 2006). The blue filled circle denotes position of the host
galaxy of PS16dtm measured from the X-shooter spectrum.

###  4.5 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_{\rm 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_{\rm var}$ is 0.43, 0.26, 0.58 for H$\beta$, H$\alpha$ 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$\beta$
and $R_{\rm FeII}$ (Sulentic et al., 2009; Marziani & Sulentic, 2014; Shen &
Ho, 2014), where $R_{\rm FeII}$ is the ratio of the equivalent width of the Fe
II (we used here only the broad component since it is dominant, measured in
the 4434–4684 Å wavelength range) to H$\beta$. Typical AGN have $R_{\rm
FeII}\sim 0.4$ with $\sim 90\%$ of objects in the range 0.1 to 1. The PS16dtm
outburst is located in the extreme right of the main sequence, in the location
of highly accreting population A objects with $R_{\rm FeII}>1$ (see e.g.,
Marziani & Sulentic, 2014).

During 1600 days of observations, $R_{\rm FeII}$ and FWHM(H$\beta$) 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_{\rm 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_{\rm bol}/L_{\rm 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_{\rm FeII}$ is $\sim$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_{\rm
FeII}$ is representative of the line ratio, therefore it should be more
reliable.

Figure 17: Location of PS16dtm (full circles) on the AGN main-sequence plane,
that is the FWHM H$\beta$ line vs. $R_{\rm FeII}$, where $R_{\rm FeII}$ is the
ratio of the equivalent width of the FeII, measured in the 4434–4684 Å
wavelength range to H$\beta$. The SDSS quasars from Shen et al. (2011) catalog
(gray dots) as well as the NLSy1 from Rakshit et al. (2017) catalog (blue
dots) are also shown. Solid lines at 4000 $\rm km\,s^{-1}$ and $R_{\rm
FeII}=1$ indicate separations between different populations of AGN (Marziani
et al., 2018). The color-bar indicates the day after MJD 57577.

##  5 Discussion

###  5.1 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.g., Namekata & Umemura,
2016; Jiang et al., 2019):

$R_{sub}=0.121\;L_{45}^{0.5}\;T_{1800}^{-2.804}\;a_{0.1}^{-0.51}\;(pc)$ (1)

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 $\mu$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\times 10^{43}$ $\rm{erg\,s}^{-1}$ and $\sim
10^{45}$ $\rm{erg\,s}^{-1}$ pre- and post-outburst), we calculated the inner
torus radius to change from $\sim$18 and $\sim$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 $\gtrsim 1500$ rest-
frame days, implying for the outer radius of the dust medium to be $R\gtrsim
4\times 10^{18}$ cm, or $\gtrsim 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., 2017, 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 galaxies141414with the exception of the
superluminous TDE ASSASN-15lh., while the characteristic values of CLAGN
covering factors are two orders of magnitude higher. The dust covering factor
of PS1-10adi is $\sim 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 $\sim$300 days using the nominal peak of the Fe II
($\sim$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$\alpha$ (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 $\sim 400$ days after the outburst.

Figure 18: Host-subtracted MIR light curves of TDEs from the literature where
the Fe II have been identified, PS1-10adi (Kankare et al., 2017; Jiang et al.,
2019; He et al., 2021), AT 2018fyk Wevers et al. (2019), and AT 2019dsg
(Cannizzaro et al., 2021).

###  5.2 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, $\sim 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 $|\Delta g|>1$ mag, $|\Delta
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 super-massive 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.

##  6 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 $\sim$ 50
days, displays another peak after $\sim 50$ days. After a plateau phase, it
started to decay slowly. In the first $\sim 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
$\leavevmode\nobreak\ \sim 4$ days, has steeply risen in the first $\sim 500$
days and it is still slowly rising $\sim 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 $\sim 2300$ K, it then dropped to $\sim 1200$ K, and then
settled to $\sim 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$\alpha$, H$\beta$, broad Fe II and narrow Fe II lines, after a rise in the
first $\sim 200$ days after the outburst, their strength dropped over time
after $\sim 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$\alpha$, Fe II emission line that reach their
peaks at $\sim 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 $\sim 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.

###### Acknowledgements.

TP acknowledges the financial support from the Slovenian Research Agency
(grants I0-0033, P1-0031, J1-8136 and Z1-1853). This work was supported with
travel grants by the Royal Swedish Academy of Sciences and the COST Action
CA16104 GWverse. GL, PC and MP were supported by a research grant (19054) from
VILLUM FONDEN. DI acknowledges funding from the grant 451-03-68/2022-14/200104
of the Ministry of Education, Science, and Technological Development of the
Republic of Serbia, and the support of the Alexander von Humboldt Foundation.
FO acknowledges support from MIUR, PRIN 2017 (grant 20179ZF5KS) ”The new
frontier of the Multi-Messenger Astrophysics: follow-up of electromagnetic
transient counterparts of gravitational wave sources” and the support of
AHEAD2020 grant agreement n.871158. TEMB acknowledges financial support from
the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de
Investigación (AEI) 10.13039/501100011033 under the PID2020-115253GA-I00
HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC)
under the PIE project 20215AT016 and the I-LINK 2021 LINKA20409, and the
program Unidad de Excelencia María de Maeztu CEX2020-001058-M. MF is supported
by a Royal Society - Science Foundation Ireland University Research
Fellowship. MG is supported by the EU Horizon 2020 research and innovation
programme under grant agreement No 101004719. TMR acknowledges the financial
support of the Vilho, Yrjö and Kalle Väisälä Foundation of the Finnish academy
of Science and Letters. KM is funded by the EU H2020 ERC grant no. 758638. AG
acknowledges the financial support from the Slovenian Research Agency
(research core funding P1-0031, infrastructure program I0-0033, project grants
J1-8136, J1-2460). NI was partially supported by Polish NCN DAINA grant No.
2017/27/L/ST9/03221. GD was supported by the Ministry of Education, Science
and Technological Development of the Republic of Serbia (contract No
451-03-68/2022-14/200002) and by the observing grant support from the
Institute of Astronomy and Rozhen NAO BAS through the bilateral joint research
project ”Gaia Celestial Reference Frame (CRF) and fast variable astronomical
objects.” ŁW was partially supported from the Polish NCN grants: Harmonia No.
2018/30/M/ST9/00311, Daina No. 2017/27/L/ST9/03221, MNiSW grant DIR/WK/2018/12
and European Commission’s H2020 OPTICON RadioNet Pilot grant No. 101004719.
The Liverpool Telescope is operated on the island of La Palma by Liverpool
John Moores University in the Spanish Observatorio del Roque de los Muchachos
of the Instituto de Astrofisica de Canarias with financial support from the UK
Science and Technology Facilities Council. This research made use of
Astropy151515http://www.astropy.org, a community-developed core Python package
for Astronomy (Astropy Collaboration et al. 2013, 2018).

## References

  * Ai et al. (2010) Ai, Y. L., Yuan, W., Zhou, H. Y., et al. 2010, ApJ, 716, L31
  * Almeyda et al. (2017) Almeyda, T., Robinson, A., Richmond, M., Vazquez, B., & Nikutta, R. 2017, ApJ, 843, 3
  * Arcavi et al. (2018) Arcavi, I., Burke, J., French, K. D., et al. 2018, The Astronomer’s Telegram, 11953, 1
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Barlow, T., et al. 2019, PASP, 131, 068003
  * Bentz et al. (2013) Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ, 767, 149
  * Berton et al. (2015) Berton, M., Foschini, L., Ciroi, S., et al. 2015, A&A, 578, A28
  * Blagorodnova et al. (2019) Blagorodnova, N., Cenko, S. B., Kulkarni, S. R., et al. 2019, ApJ, 873, 92
  * Blanchard et al. (2017) Blanchard, P. K., Nicholl, M., Berger, E., et al. 2017, ApJ, 843, 106
  * Boroson & Green (1992) Boroson, T. A. & Green, R. F. 1992, ApJS, 80, 109
  * Bricman & Gomboc (2020) Bricman, K. & Gomboc, A. 2020, ApJ, 890, 73
  * Calzetti (2001) Calzetti, D. 2001, PASP, 113, 1449
  * Campana et al. (2015) Campana, S., Mainetti, D., Colpi, M., et al. 2015, A&A, 581, A17
  * Cannizzaro et al. (2020) Cannizzaro, G., Fraser, M., Jonker, P. G., et al. 2020, MNRAS, 493, 477
  * Cannizzaro et al. (2022) Cannizzaro, G., Levan, A. J., van Velzen, S., & Brown, G. 2022, MNRAS, 516, 529
  * Cannizzaro et al. (2021) Cannizzaro, G., Wevers, T., Jonker, P. G., et al. 2021, MNRAS, 504, 792
  * Chambers et al. (2016a) Chambers, K. C., Huber, M. E., Flewelling, H., et al. 2016a, Transient Name Server Discovery Report, 2016-562, 1
  * Chambers et al. (2016b) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016b, arXiv e-prints, arXiv:1612.05560
  * Chan et al. (2021) Chan, C.-H., Piran, T., & Krolik, J. H. 2021, ApJ, 914, 107
  * Charalampopoulos et al. (2022) Charalampopoulos, P., Leloudas, G., Malesani, D. B., et al. 2022, A&A, 659, A34
  * Combes (2001) Combes, F. 2001, in Advanced Lectures on the Starburst-AGN, ed. I. Aretxaga, D. Kunth, & R. Mújica, 223
  * Dalla Bontà et al. (2020) Dalla Bontà, E., Peterson, B. M., Bentz, M. C., et al. 2020, ApJ, 903, 112
  * Dimitrijević et al. (2007) Dimitrijević, M. S., Popović, L. Č., Kovačević, J., Dačić, M., & Ilić, D. 2007, MNRAS, 374, 1181
  * Dong et al. (2016) Dong, S., Chen, P., Bose, S., et al. 2016, The Astronomer’s Telegram, 9843, 1
  * Dong et al. (2011) Dong, X.-B., Wang, J.-G., Ho, L. C., et al. 2011, ApJ, 736, 86
  * Drake et al. (2011) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2011, ApJ, 735, 106
  * Elitzur (2012) Elitzur, M. 2012, ApJ, 747, L33
  * Fausnaugh (2017) Fausnaugh, M. M. 2017, PASP, 129, 024007
  * Fitch et al. (1967) Fitch, W. S., Pacholczyk, A. G., & Weymann, R. J. 1967, ApJ, 150, L67
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Frank & Rees (1976) Frank, J. & Rees, M. J. 1976, MNRAS, 176, 633
  * Frederick et al. (2019) Frederick, S., Gezari, S., Graham, M. J., et al. 2019, The Astrophysical Journal, 883, 31
  * Frederick et al. (2021) Frederick, S., Gezari, S., Graham, M. J., et al. 2021, ApJ, 920, 56
  * Freudling et al. (2013) Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6270, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. D. R. Silva & R. E. Doxsey, 62701V
  * Gafton & Rosswog (2019) Gafton, E. & Rosswog, S. 2019, MNRAS, 487, 4790
  * Gaskell (2009) Gaskell, C. M. 2009, New A Rev., 53, 140
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005
  * Gelbord et al. (2009) Gelbord, J. M., Mullaney, J. R., & Ward, M. J. 2009, Monthly Notices of the Royal Astronomical Society, 397, 172
  * Gezari (2021) Gezari, S. 2021, ARA&A, 59 [arXiv:2104.14580]
  * Gezari et al. (2009) Gezari, S., Heckman, T., Cenko, S. B., et al. 2009, ApJ, 698, 1367
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279
  * Graham et al. (2020) Graham, M. J., Ross, N. P., Stern, D., et al. 2020, MNRAS, 491, 4925
  * Green et al. (2022) Green, P. J., Pulgarin-Duque, L., Anderson, S. F., et al. 2022, ApJ, 933, 180
  * Greene & Ho (2007) Greene, J. E. & Ho, L. C. 2007, ApJ, 670, 92
  * Guillochon & Ramirez-Ruiz (2013) Guillochon, J. & Ramirez-Ruiz, E. 2013, ApJ, 767, 25
  * He et al. (2021) He, Z., Jiang, N., Wang, T., et al. 2021, ApJ, 907, L29
  * Hills (1975) Hills, J. G. 1975, Nature, 254, 295
  * Hinkle et al. (2021) Hinkle, J. T., Holoien, T. W. S., Shappee, B. J., & Auchettl, K. 2021, ApJ, 910, 83
  * Hinkle et al. (2022) Hinkle, J. T., Holoien, T. W. S., Shappee, B. J., et al. 2022, ApJ, 930, 12
  * Hodgkin et al. (2021) Hodgkin, S. T., Harrison, D. L., Breedt, E., et al. 2021, A&A, 652, A76
  * Hon et al. (2022) Hon, W. J., Wolf, C., Onken, C. A., Webster, R., & Auchettl, K. 2022, MNRAS, 511, 54
  * Huber et al. (2015) Huber, M., Chambers, K. C., Flewelling, H., et al. 2015, The Astronomer’s Telegram, 7153, 1
  * Ilić et al. (2020) Ilić, D., Oknyansky, V., Popović, L. Č., et al. 2020, A&A, 638, A13
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Jiang et al. (2021) Jiang, N., Wang, T., Hu, X., et al. 2021, ApJ, 911, 31
  * Jiang et al. (2019) Jiang, N., Wang, T., Mou, G., et al. 2019, ApJ, 871, 15
  * Jiang et al. (2017) Jiang, N., Wang, T., Yan, L., et al. 2017, ApJ, 850, 63
  * Kankare et al. (2017) Kankare, E., Kotak, R., Mattila, S., et al. 2017, Nature Astronomy, 1, 865
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961
  * Kollatschny & Fricke (1985) Kollatschny, W. & Fricke, K. J. 1985, A&A, 146, L11
  * Komossa et al. (2017) Komossa, S., Grupe, D., Schartel, N., et al. 2017, in New Frontiers in Black Hole Astrophysics, ed. A. Gomboc, Vol. 324, 168–171
  * Komossa et al. (2008) Komossa, S., Zhou, H., Wang, T., et al. 2008, ApJ, 678, L13
  * Kool et al. (2020) Kool, E. C., Reynolds, T. M., Mattila, S., et al. 2020, MNRAS, 498, 2167
  * Korista & Ferland (1989) Korista, K. T. & Ferland, G. J. 1989, ApJ, 343, 678
  * Kovačević et al. (2010) Kovačević, J., Popović, L. Č., & Dimitrijević, M. S. 2010, ApJS, 189, 15
  * Kulkarni (2013) Kulkarni, S. R. 2013, The Astronomer’s Telegram, 4807, 1
  * Leloudas et al. (2019) Leloudas, G., Dai, L., Arcavi, I., et al. 2019, ApJ, 887, 218
  * Leloudas et al. (2016) Leloudas, G., Fraser, M., Stone, N. C., et al. 2016, Nature Astronomy, 1, 0002
  * Lodato & Rossi (2011) Lodato, G. & Rossi, E. M. 2011, MNRAS, 410, 359
  * López-Navas et al. (2022) López-Navas, E., Martínez-Aldama, M. L., Bernal, S., et al. 2022, MNRAS, 513, L57
  * Lyu et al. (2022) Lyu, B., Wu, Q., Yan, Z., Yu, W., & Liu, H. 2022, ApJ, 927, 227
  * Lyutyj et al. (1984) Lyutyj, V. M., Oknyanskij, V. L., & Chuvaev, K. K. 1984, Soviet Astronomy Letters, 10, 335
  * MacLeod et al. (2019) MacLeod, C. L., Green, P. J., Anderson, S. F., et al. 2019, ApJ, 874, 8
  * Magnier et al. (2013) Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30
  * Marziani et al. (2018) Marziani, P., Dultzin, D., Sulentic, J. W., et al. 2018, Frontiers in Astronomy and Space Sciences, 5, 6
  * Marziani & Sulentic (2014) Marziani, P. & Sulentic, J. W. 2014, MNRAS, 442, 1211
  * Mattila et al. (2018) Mattila, S., Pérez-Torres, M., Efstathiou, A., et al. 2018, Science, 361, 482
  * Mehdipour et al. (2022) Mehdipour, M., Kriss, G. A., Brenneman, L. W., et al. 2022, ApJ, 925, 84
  * Mockler et al. (2019) Mockler, B., Guillochon, J., & Ramirez-Ruiz, E. 2019, ApJ, 872, 151
  * Moriya et al. (2017) Moriya, T. J., Tanaka, M., Morokuma, T., & Ohsuga, K. 2017, ApJ, 843, L19
  * Namekata & Umemura (2016) Namekata, D. & Umemura, M. 2016, MNRAS, 460, 980
  * Netzer (2006) Netzer, H. 2006, in Physics of Active Galactic Nuclei at all Scales, ed. D. Alloin, Vol. 693, 1
  * Noda & Done (2018) Noda, H. & Done, C. 2018, MNRAS, 480, 3898
  * Oknyansky et al. (2019) Oknyansky, V. L., Winkler, H., Tsygankov, S. S., et al. 2019, MNRAS, 483, 558
  * Onori et al. (2019) Onori, F., Cannizzaro, G., Jonker, P. G., et al. 2019, MNRAS, 489, 1463
  * Onori et al. (2022) Onori, F., Cannizzaro, G., Jonker, P. G., et al. 2022, MNRAS, 517, 76
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Panda et al. (2018) Panda, S., Czerny, B., Adhikari, T. P., et al. 2018, ApJ, 866, 115
  * Park et al. (2022) Park, D., Barth, A. J., Ho, L. C., & Laor, A. 2022, ApJS, 258, 38
  * Peterson (2001) Peterson, B. M. 2001, in Advanced Lectures on the Starburst-AGN, ed. I. Aretxaga, D. Kunth, & R. Mújica, 3
  * Phinney (1989) Phinney, E. S. 1989, in The Center of the Galaxy, ed. M. Morris, Vol. 136, 543
  * Pons & Watson (2014) Pons, E. & Watson, M. G. 2014, A&A, 568, A108
  * Rakić (2022) Rakić, N. 2022, MNRAS, 516, 1624
  * Rakshit et al. (2017) Rakshit, S., Stalin, C. S., Chand, H., & Zhang, X.-G. 2017, ApJS, 229, 39
  * Ranalli et al. (2003) Ranalli, P., Comastri, A., & Setti, G. 2003, A&A, 399, 39
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523
  * Ren et al. (2022) Ren, W., Wang, J., Cai, Z., & Guo, H. 2022, ApJ, 925, 50
  * Richards et al. (2006) Richards, G. T., Lacy, M., Storrie-Lombardi, L. J., et al. 2006, ApJS, 166, 470
  * Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95
  * Rose et al. (2015) Rose, M., Elvis, M., & Tadhunter, C. N. 2015, MNRAS, 448, 2900
  * Roth et al. (2016) Roth, N., Kasen, D., Guillochon, J., & Ramirez-Ruiz, E. 2016, ApJ, 827, 3
  * Roth et al. (2021) Roth, N., van Velzen, S., Cenko, S. B., & Mushotzky, R. F. 2021, ApJ, 910, 93
  * Runnoe et al. (2012) Runnoe, J. C., Brotherton, M. S., & Shang, Z. 2012, MNRAS, 422, 478
  * Salim & Narayanan (2020) Salim, S. & Narayanan, D. 2020, ARA&A, 58, 529
  * Sánchez-Sáez et al. (2021) Sánchez-Sáez, P., Lira, H., Martí, L., et al. 2021, AJ, 162, 206
  * Sánchez-Sáez et al. (2018) Sánchez-Sáez, P., Lira, P., Mejía-Restrepo, J., et al. 2018, ApJ, 864, 87
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Selsing et al. (2019) Selsing, J., Malesani, D., Goldoni, P., et al. 2019, A&A, 623, A92
  * Shapovalova et al. (2012) Shapovalova, A. I., Popović, L. Č., Burenkov, A. N., et al. 2012, ApJS, 202, 10
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48
  * Shen & Ho (2014) Shen, Y. & Ho, L. C. 2014, Nature, 513, 210
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45
  * Sheng et al. (2020) Sheng, Z., Wang, T., Jiang, N., et al. 2020, ApJ, 889, 46
  * Sheng et al. (2017) Sheng, Z., Wang, T., Jiang, N., et al. 2017, ApJ, 846, L7
  * Smartt et al. (2015) Smartt, S. J., Valenti, S., Fraser, M., et al. 2015, A&A, 579, A40
  * Smith et al. (2020) Smith, K. W., Smartt, S. J., Young, D. R., et al. 2020, PASP, 132, 085002
  * Sniegowska et al. (2020) Sniegowska, M., Czerny, B., Bon, E., & Bon, N. 2020, A&A, 641, A167
  * Śniegowska et al. (2022) Śniegowska, M., Grzedzielski, M., Czerny, B., & Janiuk, A. 2022, Astronomische Nachrichten, 343, e210065
  * Steele et al. (2004) Steele, I. A., Smith, R. J., Rees, P. C., et al. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5489, Proc. SPIE, ed. J. Oschmann, Jacobus M., 679–692
  * Stern et al. (2018) Stern, D., McKernan, B., Graham, M. J., et al. 2018, ApJ, 864, 27
  * Sulentic et al. (2009) Sulentic, J. W., Marziani, P., & Zamfir, S. 2009, New A Rev., 53, 198
  * Terreran et al. (2016) Terreran, G., Berton, M., Benetti, S., et al. 2016, The Astronomer’s Telegram, 9417, 1
  * Thomsen et al. (2022) Thomsen, L. L., Kwan, T. M., Dai, L., et al. 2022, ApJ, 937, L28
  * Tonry (2011) Tonry, J. L. 2011, PASP, 123, 58
  * Tonry et al. (2018) Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, PASP, 130, 064505
  * Tonry et al. (2012) Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99
  * Trakhtenbrot et al. (2019) Trakhtenbrot, B., Arcavi, I., MacLeod, C. L., et al. 2019, ApJ, 883, 94
  * Trakhtenbrot et al. (2019) Trakhtenbrot, B., Arcavi, I., Ricci, C., et al. 2019, Nature Astronomy, 3, 242–250
  * van Groningen & Wanders (1992) van Groningen, E. & Wanders, I. 1992, PASP, 104, 700
  * van Velzen et al. (2020) van Velzen, S., Holoien, T. W.-S., Onori, F., Hung, T., & Arcavi, I. 2020, Space Science Reviews, 216
  * van Velzen et al. (2016) van Velzen, S., Mendez, A. J., Krolik, J. H., & Gorjian, V. 2016, The Astrophysical Journal, 829, 19
  * van Velzen et al. (2021) van Velzen, S., Pasham, D. R., Komossa, S., Yan, L., & Kara, E. A. 2021, Space Sci. Rev., 217, 63
  * van Velzen et al. (2019) van Velzen, S., Stone, N. C., Metzger, B. D., et al. 2019, ApJ, 878, 82
  * Vanden Berk et al. (2006) Vanden Berk, D. E., Shen, J., Yip, C.-W., et al. 2006, AJ, 131, 84
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105
  * Véron-Cetty et al. (2006) Véron-Cetty, M. P., Joly, M., Véron, P., et al. 2006, A&A, 451, 851
  * Wang & Bon (2020) Wang, J.-M. & Bon, E. 2020, A&A, 643, L9
  * Wang et al. (2012) Wang, T.-G., Zhou, H.-Y., Komossa, S., et al. 2012, ApJ, 749, 115
  * Wevers et al. (2019) Wevers, T., Pasham, D. R., van Velzen, S., et al. 2019, MNRAS, 488, 4816
  * Willott (2005) Willott, C. J. 2005, ApJ, 627, L101
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Wyrzykowski et al. (2017) Wyrzykowski, Ł., Zieliński, M., Kostrzewa-Rutkowska, Z., et al. 2017, MNRAS, 465, L114
  * Xiao et al. (2011) Xiao, T., Barth, A. J., Greene, J. E., et al. 2011, ApJ, 739, 28
  * Yang et al. (2018) Yang, Q., Wu, X.-B., Fan, X., et al. 2018, ApJ, 862, 109
  * Yaron & Gal-Yam (2012) Yaron, O. & Gal-Yam, A. 2012, PASP, 124, 668
  * Zabludoff et al. (2021) Zabludoff, A., Arcavi, I., La Massa, S., et al. 2021, Space Sci. Rev., 217, 54
  * Zhang et al. (2022) Zhang, W. J., Shu, X. W., Sheng, Z. F., et al. 2022, A&A, 660, A119
  * Zhang (2021) Zhang, X. 2021, ApJ, 919, 13
  * Zhou et al. (2006) Zhou, H., Wang, T., Yuan, W., et al. 2006, ApJS, 166, 128