Extreme variability in an active galactic nucleus: Gaia16aax
G. Cannizzaro, M. Fraser, P. G. Jonker, J. E. Pringle, S. Mattila, P. C. Hewett, T. Wevers, E. Kankare, Z. Kostrzewa-Rutkowska, Ł. Wyrzykowski, F. Onori, J. Harmanen, K. E. S. Ford, B. McKernan, C. J. Nixon
MMNRAS , 1–21 (2019) Preprint 22 January 2020 Compiled using MNRAS L A TEX style file v3.0
Extreme variability in an active galactic nucleus:Gaia16aax
G. Cannizzaro , (cid:63) , M. Fraser , P. G. Jonker , , J. E. Pringle , S. Mattila ,P. C. Hewett , T. Wevers , E. Kankare , Z. Kostrzewa-Rutkowska , , , (cid:32)L. Wyrzykowski ,F. Onori , J. Harmanen , K.E.S. Ford , , , B. McKernan , , , C. J. Nixon SRON, Netherlands Institute for Space Research, Sorbonnelaan, 2, NL-3584CA Utrecht, the Netherlands Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, the Netherlands School of Physics, O’Brien Centre for Science North, University College Dublin, Belfield, Dublin 4, Ireland Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK Tuorla Observatory, Department of Physics and Astronomy, FI-20014 University of Turku, Finland Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands Warsaw University Astronomical Observatory, Al. Ujazdowskie 4, PL-00-478 Warszawa, Poland Istituto di Astrofisica e Planetologia Spaziali (INAF), Via Fosso del Cavaliere 100, Roma, I-00133, Italy Department of Astrophysics, American Museum of Natural History, Central Park West at 79th Street, New York, NY 10024 Department of Science, Borough of Manhattan Community College, City University of New York, New York, NY 10007 Physics Program, The Graduate Center, City University of New York, New York, NY 10016 Department of Physics and Astronomy, University of Leicester, Leicester LE1 7RH, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present the results of a multi-wavelength follow up campaign for the luminousnuclear transient Gaia16aax, which was first identified in January 2016. The transientis spatially consistent with the nucleus of an active galaxy at z=0.25, hosting a blackhole of mass ∼ × M (cid:12) . The nucleus brightened by more than 1 magnitude in theGaia G-band over a timescale of less than one year, before fading back to its pre-outburst state over the following three years. The optical spectra of the source showbroad Balmer lines similar to the ones present in a pre-outburst spectrum. During theoutburst, the H α and H β emission lines develop a secondary peak. We also report on thediscovery of two transients with similar light curve evolution and spectra: Gaia16akaand Gaia16ajq. We consider possible scenarios to explain the observed outbursts. Weexclude that the transient event could be caused by a microlensing event, variabledust absorption or a tidal encounter between a neutron star and a stellar mass blackhole in the accretion disk. We consider variability in the accretion flow in the innerpart of the disk, or a tidal disruption event of a star ≥ M (cid:12) by a rapidly spinningsupermassive black hole as the most plausible scenarios. We note that the similaritybetween the light curves of the three Gaia transients may be a function of the Gaiaalerts selection criteria. Key words: galaxies: active – quasars: supermassive black holes – galaxies: nuclei– accretion, accretion disks – transients: tidal disruption events – quasars: individual:Gaia16aax
Lying at the center of galaxies, supermassive black holes(SMBHs) have an enormous impact on the region surround-ing them. Once circumnuclear matter is close enough to theSMBH to fall into the gravitational potential, the black hole (cid:63)
E-mail: [email protected] will start accreting this matter through the formation ofan accretion disk if the material has angular momentum.These Active Galactic Nuclei (AGNs) typically emit overthe whole electromagnetic spectrum. Broad (1 000–20 000km s − ) and/or narrow (300–1 000 km s − ) highly ionisedemission lines are typically present in the optical spectra ofAGNs. These lines come from clouds of matter in differentregions and at different distances from the central black hole: © a r X i v : . [ a s t r o - ph . H E ] J a n G. Cannizzaro et al. the Broad Line Region (BLR) and the more distant Nar-row Line Region (NLR). The presence or absence of broadand/or narrow lines lead to the classification of AGNs in“types”: Type 1 when both kinds of lines are present, Type2 when only narrow lines are present. Intermediate types,depending on the presence of single Balmer lines, exist, e.g.type 1.9 characterised by the presence of a broad H α emis-sion line but with an absent broad H β (Osterbrock 1981).The unified model explains this dichotomy in terms of ourviewing angle to the AGN, where a parsec-scale dusty torusmay partially or totally obscure our view of the central en-gine and the BLR (Antonucci 1993; Urry & Padovani 1995).AGN are intrinsically variable objects, with stochasticvariations in brightness of order of 20% over time-scale ofmonths to years and larger variations happening on longertimescales (MacLeod et al. 2010, 2012). In recent years, moreextreme examples of Quasar variability have been discov-ered, e.g. Graham et al. (2017) and Rumbaugh et al. (2018)found large samples of objects that showed variability ofmore than 1 magnitude over a timescale of several yearsin an archival search in the Catalina Real-Time TransientSurvey, Sloan Digital Sky Survey and Dark Energy Survey.Lawrence et al. (2016) investigated the nature of a sam-ple large amplitude nuclear transients discovered during thePan-STARSS π survey. The majority of these objects areclassified as hypervariable AGNs. In some cases, AGNs willshow slow and significant optical variability of more than 1magnitude over several years, accompanied by spectral vari-ability (e.g Matt et al. 2003; LaMassa et al. 2015; MacLeodet al. 2016). These so-called ”changing look AGNs” havebeen observed transitioning between different AGN types,with the (dis)appearance of broad emission lines (typicallythe Balmer series). While various physical mechanisms havebeen proposed to explain this extreme variability, such asvariable obscuration across the line of sight, microlensing oraccretion disk instabilities, no single explanation has beenfound. Merloni et al. (2015) proposes as an explanation theoccurrence of a Tidal Disruption Event (TDE) in the nucleusof the galaxy. A TDE happens when a star passes so close toa SMBH that the tidal forces of the SMBH are stronger thanthe star’s self-gravity. This leads to the (partial) disruptionof the star (Hills 1975; Rees 1988; Evans & Kochanek 1989).The disruption of the star and the subsequent accretion ofroughly half of the stellar material (the other half will beexpelled on unbound orbits) gives rise to a short and lu-minous flare that typically peaks in the UV or soft X-rays.The bolometric light curve decay of the transient event isexpected to follow the fallback rate of the material with apowerlaw decline ∝ t − / on a timescale of some months upto a year (Evans & Kochanek 1989; Cannizzo et al. 1990;Lodato et al. 2009). TDEs have been invoked as a possibleexplanation also for numerous objects in the large samplesdescribed in Graham et al. (2017); Lawrence et al. (2016);Rumbaugh et al. (2018).Here, we report the follow up observations of Gaia16aax,a transient event happening at the center of a galaxy host-ing a known QSO (SDSS J143418.47+491236.5, Abolfathi We use the term Quasar or QSO to mean the most luminousfraction of the AGN population, with bolometric luminosity L (cid:38) erg s − et al. 2018) at z=0.25, discovered by the Gaia Science Alertsproject (GSA, Hodgkin et al. 2013). The transient was giventhe designation AT2016dbt in the Transient Name Server(TNS). The center of the galaxy brightened by ∼ α and H β ), already present in thepre-outburst spectrum, during the outburst show a signif-icantly different morphology, with two peaks of differentintensity and separated by ∼
100 ˚A. The initial classifica-tion spectrum showed a strong blue continuum with broad,double-peaked Balmer lines. After the initial classification asan AGN outburst, we started a multiwavelength follow-upcampaign with optical and NIR photometry, optical spec-troscopy, and X-ray observations.In Section 2 we describe both the available, pre-outburstdata and our own follow-up data, in Section 3 we describeour analysis of the photometric data, the Spectral EnergyDistribution (SED) and the fit to the emission lines, finallyin Section 4 we discuss our results in the broader frameworkof extreme AGN variability and we discuss some possiblescenarios to explain the outburst.Throughout the paper we assume a cosmology with H =67.7 km s − Mpc − , Ω M =0.309, Ω Λ =0.691, consistentwith Planck Collaboration et al. (2016). Gaia16aax was first alerted on by the Gaia Photometric Sci-ence Alerts program on 2016 January 26. The transient is ∼ . (cid:48)(cid:48) from the position of the host galaxy in the GaiaData Release 2 (DR2, Gaia Collaboration et al. 2018) and ∼ . (cid:48)(cid:48) from the position of the host galaxy in the Sloan Dig-ital Sky Survey (SDSS, Abolfathi et al. 2018). As shown inKostrzewa-Rutkowska et al. (2018), where nuclear transientsdetected by Gaia were matched astrometrically with SDSSsources, the astrometry provided by the Gaia PhotometricScience Alerts system is accurate to ∼ . (cid:48)(cid:48) . We hence regardGaia16aax as consistent with having no offset from the hostnucleus. The Gaia lightcurve is shown in Fig.1: it shows asteady, smooth rise up to its maximum in less then a yearand a slower decay on a timescale of more than two years toits pre-outburst, quiescent level. The apparent magnitude atpeak in the Gaia broad G filter (which corresponds to whitelight, see Jordi et al. 2010) is G = . , which corresponds toan absolute magnitude M G = − . , given a source redshiftz (cid:39) D L = . Gpc, assumingno extinction. The date of the alert issue corresponds withthe date at which, according to the Gaia lightcurve, the flarereached its peak emission. https://wis-tns.weizmann.ac.il/object/2016dbt http://gsaweb.ast.cam.ac.uk/alerts/alert/Gaia16aaxMNRAS , 1–21 (2019) aia16aax Figure 1.
Light curve from the
Gaia satellite for Gaia16aax. Thearrow indicates the date the transient was alerted on. Magnitudesare in the
Gaia G broad band filter (Jordi et al. 2010). The blackcircles at the top indicate the epochs at which an optical spec-trum was taken. The black triangles and squares at the bottomshow the epochs at which optical and NIR images were taken,respectively. The two y-axes are apparent (left axis) and absolute(right axis) magnitudes.
Figure 2.
Archival spectrum from the SDSS legacy survey. Theobject shows broad Balmer emission lines and narrow forbiddenoxygen and nitrogen lines.
The source is present in SDSS Data Release 14 as SDSSJ143418.47+491236.5. It was imaged on 2002 May 8 andmorphologically classified as a galaxy with r=17.76 ± ± g,r,i,z ) are similar to, though not formally con-sistent with the ones reported in SDSS. The different valuesare likely due to the differences in the PSF modeling usedin the two surveys. from here onwards all times are in UTC. Figure 3.
Historical lightcurve from the Catalina Sky Survey. Inthe period covered by the survey the source has not exhibited anoutburst similar to Gaia16aax.
The source has been observed by the Catalina Sky Survey(CSS) (Drake et al. 2009) over a period of eight years (fromJune 2005 to June 2013, source ID CSS J143418.6+491236).As shown in Fig. 3, the object has not shown an enhancedemission state like Gaia16aax during the eight years of CSSmonitoring. The source is present in the Two Micron All-Sky Survey(2MASS) as 2MASS J14341848+4912366. It has been ob-served once on 2000 May 07 in the three Near Infra-Redbands J , H and K s , with magnitude values J = . ± . , H = . ± . and K s = . ± . . The source has been observed in the infrared band by boththe WISE and the NEO-WISE (Near-Earth Objects Wide-field Infrared Survey Explorer, Mainzer et al. 2011) phasesof the
WISE satellite between January 2014 and June 2017.This period has overlap with both the Gaia data (see Sec.2.1)and our own follow-up data. The WISE satellite provides ob-servations in four filters W1, W2, W3 and W4 (3.4, 4.6, 11.6and 22.1 µ m respectively), while the reactivation missionNEO-WISE only provides values for the W1 and W2 filters,in which it takes 10-20 single exposures for each epoch ofobservation. We selected only data points with photometricflags A or B (for which the minimum flux signal-to-noiseratio is above 10 and 3, respectively). The NEO-WISE mag-nitudes in the W1 and W2 bands in the quiescent phase(before MJD 57200) are compatible with the ones measuredduring the original WISE mission (Wright et al. 2010): W1= 13.783 ± ± The CSS survey is funded by the National Aeronautics andSpace Administration under Grant No. NNG05GF22G issuedthrough the Science Mission Directorate Near-Earth Objects Ob-servations Program.MNRAS000
WISE satellite between January 2014 and June 2017.This period has overlap with both the Gaia data (see Sec.2.1)and our own follow-up data. The WISE satellite provides ob-servations in four filters W1, W2, W3 and W4 (3.4, 4.6, 11.6and 22.1 µ m respectively), while the reactivation missionNEO-WISE only provides values for the W1 and W2 filters,in which it takes 10-20 single exposures for each epoch ofobservation. We selected only data points with photometricflags A or B (for which the minimum flux signal-to-noiseratio is above 10 and 3, respectively). The NEO-WISE mag-nitudes in the W1 and W2 bands in the quiescent phase(before MJD 57200) are compatible with the ones measuredduring the original WISE mission (Wright et al. 2010): W1= 13.783 ± ± The CSS survey is funded by the National Aeronautics andSpace Administration under Grant No. NNG05GF22G issuedthrough the Science Mission Directorate Near-Earth Objects Ob-servations Program.MNRAS000 , 1–21 (2019)
G. Cannizzaro et al. available during the NEO-WISE phase of the mission) areW3 = 10.93 ± ± The source was observed with all the instruments on boardthe Neils Gehrels Swift Observatory (Gehrels et al. 2004,
Swift from here on) on 9 and 14 September 2006. UVdata were analysed with the
Swift task
UVOTSOURCE , us-ing a (cid:48)(cid:48) aperture to estimate the source brightness anda (cid:48)(cid:48) aperture, centered on a nearby empty region ofsky, to estimate the background levels. Combining the twoepochs of observations we derived that the source hasUV magnitudes UVW1=21.5 ± ± ± . ± . × − erg cm − s − , obtained modeling the datawith a power-law with photon index α ∼ . . This, at theluminosity distance reported previously, corresponds to a lu-minosity L X = . ± . × erg s − .We observed the source again with Swift on 18 July2019. The UV and X-ray data was reduced and analysed inthe same way as the archival data. We derived UV mag-nitudes UVW1 = . ± . , UVM2= . ± . andUVW2= . ± . . The 0.3 – 10.0 keV X-ray luminosity is L X = ( . ± . ) × erg s − . The object is therefore stillbright in the UV and X-rays, with respect to the archivalvalues.We observed Gaia16aax also with XMM-Newton, forwhich we were awarded a Director’s Discretionary Time ob-servation. The observation started on June 30, 2016. Dur-ing the observations, the pn camera was operated in FullWindow mode, providing a time resolution of 73.4 ms. Boththe MOS detectors were operated in Partial Window mode,where the Small Window option was used, implying thatonly the central 100 ×
100 pixels of 1.1 (cid:48)(cid:48) each were read out.This yields a time resolution of 0.3 seconds. The source fluxis too low to yield meaningful Reflection Grating Spectrom-eter data.We run the default
SAS v17 (20180620) tool xmmex-tractor under the HEASOFT ftools software version 6.24to extract source light curves and spectra. After filtering forperiods of high background, we are left with a total expo-sure for the pn and MOS1 and MOS2 of 46.1, 59.7 and 59.4ksec, respectively. We use the
SAS command epatplot toassess if photon pile–up is important during our observation.We conclude pile–up is not important, as expected given theoverall number of photons detected and given that the lightcurve does not show large flares.The pn source spectrum was extracted from a circu-lar region with radius of 30 (cid:48)(cid:48) centered on the known op-tical source position. Background photons were extractedfrom a source-free circle with a radius of 59.59 (cid:48)(cid:48) centered onRight Ascension 218.527 and Declination 49.2065. The MOSsource spectra were extracted from circular regions with ra-dius of 30 (cid:48)(cid:48) centered on the known optical source position.The MOS background spectra are obtained from an annu-lus centered on the source position, with inner and outer radius of 330 and 660 (cid:48)(cid:48) , respectively. The extracted spectraare rebinned to yield a minimum of 30 counts per bin.We fitted the spectra with a powerlaw modified byGalactic foreground extinction. The reduced χ of the fit is1.24 for 250 degrees of freedom. The normalisation for the pnand MOS power laws are within a few percent of each other,and we tie the other parameters (N H and the powerlaw in-dex) to require them to be the same. The best–fitting powerlaw index is 1.85 ± × cm − . The 0.3–10 keV absorbed and unabsorbed fluxis ( . ± . ) × − , ( . ± . ) × − erg cm − s − , respec-tively, corresponding to a luminosity ( . ± . )× erg s − .Simultaneously to the X-ray observation, we obtainedphotometry in the UV filters UVM2 and
UVW1 and opticalfilter Johnson U with the Optical Monitor (OM) on boardof XMM-Newton . The resulting magnitudes are
UVM2 = . ± . , UVW1 = . ± . and U = . ± . , inthe AB magnitude system. We started a follow-up campaign of Gaia16aax using the Al-hambra Faint Object Spectrograph (ALFOSC) using grism ∼ (cid:48)(cid:48) ) at the NordicOptical Telescope (NOT), in the framework of the NOT Un-biased Transient Survey . For all observations the slit wasplaced at parallactic angle. We obtained 7 epochs of opticalspectroscopy between 2016 February 09 and 2017 July 01(see Table 1). All observations were reduced using foscgui which is a pipeline based on standard iraf data reductionsprocedures: flat field and bias correction, cosmic-ray clean-ing, wavelength and flux calibration with arc lamps andstandard stars and telluric line correction (Tody 1986). Thepipeline also performs second order contamination removalfor the spectra following Stanishev (2007). The sequence ofALFOSC spectra is shown in Fig.4.On 2018 July 08 we took one last spectrum to verify thatthe object had returned to quiescence as suggested by itsGaia lightcurve (see Fig.1). Indeed, this spectrum appears tobe similar to the pre-outburst SDSS spectrum and dividingthe two spectra results in an almost featureless spectrum butwith a slight slope. The slope and the presence of residualH α and [O iii ] lines can both be introduced due to the useof different spectrographs by SDSS (optical fibers) and NOT(long-slit): an optical fiber will collect more light from thewhole galaxy, while the slit is centered on the nuclear region,hence the SDSS spectrum will have an excess of light comingfrom the circumnuclear region (see Fig. 5). Besides our spectroscopic follow-up, we started a series ofphotometric observations to monitor the decay of the tran-sient to complement the lightcurve from Gaia. These pho-tometric observations were also performed in the framework http://csp2.lco.cl/not/ Foscgui is a graphic user interface aimed at extracting SN spec-troscopy and photometry obtained with FOSC-like instruments.It was developed by E. Cappellaro. A package description can befound at http://sngroup.oapd.inaf.it/foscgui.htmlMNRAS , 1–21 (2019) aia16aax Figure 4.
Eight optical spectra taken with ALFOSC at the NOT, shifted to the rest-frame of the host galaxy. For each spectrum weannotate to the right the number of days passed since the publishing of the Gaia Science Alert (2016 January 26, see Sec.2.1). Note thelarge gap in time between the first and second spectroscopic observation. The dotted lines indicate the wavelengths of the main emissionlines, the dashed lines indicate the secondary components (red wings) of the H α and H β lines. The grey bands represent wavelengthsaffected by the telluric absorption, one of which affects the H α line region. Figure 5.
Comparison between the SDSS pre-outburst spectrum(in red) and the last spectrum of our follow-up (July 2018, inblue). The division of the two spectra results in a nearly fea-tureless spectrum, shown in the bottom panel. In the residualspectrum there is still a trace of the H α and [O iii ] emission lines(marked in the bottom panel). Table 1.
ALFOSC/NOT spectroscopic observationsMJD ( ) UTC Date exposure time slit airmass[d] [s] [ (cid:48)(cid:48) ]57 428.27 2016 Feb 10 2400 1.0 1.0757 571.06 2016 Jul 02 1800 1.0 1.5657 600.97 2016 Jul 31 2700 1.0 1.4957 732.25 2016 Dec 10 1800 1.0 1.7657 780.26 2017 Jan 27 2400 1.0 1.1157 845.15 2017 Apr 02 2700 1.0 1.0757 935.94 2017 Jul 01 2000 1.3 1.1058 307.01 2018 Jul 08 2700 1.0 1.32
Note. (1) Modified Julian Day of observations. of the NOT Unbiased Transient Survey using NOTCam(Nordic Optical Telescope near-infrared Camera and spec-trograph) at the NOT for the NIR bands J, H and K s . Inaddition we employed the Optical Wide Field camera (IO:O)at the Liverpool Telescope (LT) using the Sloan ugriz andBessel B and V filters. After the source went back to its pre- MNRAS000
Note. (1) Modified Julian Day of observations. of the NOT Unbiased Transient Survey using NOTCam(Nordic Optical Telescope near-infrared Camera and spec-trograph) at the NOT for the NIR bands J, H and K s . Inaddition we employed the Optical Wide Field camera (IO:O)at the Liverpool Telescope (LT) using the Sloan ugriz andBessel B and V filters. After the source went back to its pre- MNRAS000 , 1–21 (2019)
G. Cannizzaro et al. outburst emission level, we obtained a final set of observa-tions in all necessary filters to subtract the host-galaxy/pre-outburst contribution to the flux detected in outburst. NOT-CAM images were reduced using a modified version of theexternal NOTCAM iraf package (version 2.5) and theirzero points calibrated using 2MASS (Two Micron All-SkySurvey, Skrutskie et al. 2006) catalogue. Optical images werereduced through the Liverpool Telescope pipelines and theirzero points were calculated using stars in SDSS. For the cal-culation of the zero points in Bessel B and V, we applied thefilter transformation equations found in Jester et al. (2005).For each epoch of observation we performed differentialphotometry to measure the magnitude of the transient. Forthis we subtracted a template image, taken when the objectwas back at its pre-outburst emission level, from our scien-tific images, using HOTPANTS (Becker 2015). This programuses the algorithm from Alard & Lupton (1998) for the sub-traction of images taken under different seeing conditions:it applies a spatially varying kernel to one of the two im-ages to match the PSF of the other one and then subtractsthe template image to obtain a final image in which onlythe transient is visible, while all constant luminosity sourceswill have been subtracted.On these subtracted images we performed aperture pho-tometry using the iraf task APPHOT with variable aperturesdepending on the seeing conditions. Uncertainties on themagnitudes are obtained by adding in quadrature the photo-metric error from the aperture photometry and the standarddeviation due to the scatter in the sources used to calculatethe zero points. We did not apply band-pass corrections asour uncertainties are dominated by systematic errors due tothe image-subtraction process.
In the first months of 2016, other two objects with verysimilar characteristics to Gaia16aax were discovered bythe Gaia Photometric Science Alerts pipeline: Gaia16ajq (AT2016dvz on the TNS) and Gaia16aka (AT2016dwk onthe TNS).Gaia16ajq was discovered on 2016 March 31, as an increasein brightness by ∼ ∼ ∼ ∼ https://github.com/acbecker/hotpants https://gsaweb.ast.cam.ac.uk/alerts/alert/Gaia16ajq/ https://gsaweb.ast.cam.ac.uk/alerts/alert/Gaia16aka/ Figure 6.
Light curves of the three similar sources found by theGaia Alerts system: Gaia16aax (black squares), Gaia16ajq (bluediamonds) and Gaia16aka (red circles). For the sake of compari-son, all three lightcurves have been shifted vertically and horizon-tally. The 0 on the x-axis is the date at which each object reachesits peak emission.
Figure 7.
Spectra of the three similar Gaia objects: Gaia16aka(at the top, in red), Gaia16ajq (in the middle, in magenta) andGaia16aax (at the bottom, in black). With ⊕ we indicate telluricabsorption bands. decay part of the lightcurve, with Gaia16ajq showing twovery prominent ones at ∼
200 and ∼
400 days after peak.The three objects show similarities also in their spectra.For Gaia16ajq and Gaia16aka we only have a classificationspectrum, taken with the NOT a few days after discovery.Both spectra and the first spectrum of Gaia16aax, for com-parison, are shown in Fig. 7. All three objects show a bluebump - with Gaia16aka showing the most prominent one- broad Balmer lines and narrow forbidden lines typical ofAGNs. In all three spectra the continuum between H γ andH β is high, this could be due to the presence of unresolvedlines due to the Bowen fluorescence mechanism or to a forestof unresolved Fe lines.While the H α and H β lines in Gaia16aax showed avery distinct double-peaked shape during the outburst, bothGaia16ajq and Gaia16aka do not show the same morphologyin their spectra. MNRAS , 1–21 (2019) aia16aax Figure 8.
Lightcurve from the NEO-WISE mission, comparedwith the Gaia lightcurve. Plotted are the magnitudes in the twofilters W1 (blue) and W2 (green), shifted along the y-axis to allowto compare with the Gaia data. Solid blue squares (W1) and greencircles (W2) are the mean values of the exposures taken for eachepoch. The NEO-WISE lightcurve reaches its peak with a delayof ∼
140 days with respect to the Gaia lightcurve.
The Gaia lightcurve shown in Fig. 1 shows a decay from thepeak of the emission to the pre-outburst level of emissionover more than 2 years. The lightcurve shows two bumpsat around 57700 and 58000 MJD (around 300 and 600 daysfrom peak, respectively). The decay rate is well fit by an ex-ponential decay ∼ t − . ± . (with χ ν = . ). As shown inFigure 8, where we compare the NEO-WISE lightcurve withthe Gaia one, Gaia16aax reaches the peak of its MIR emis-sion with a delay of ∼
140 days with respect to the opticaldata. We assumed that the epoch at which the NEO-WISEreaches its peak to be the one corresponding to the max-imum observed value (MJD ∼ ± days. Assuming that the MIR emission comes from a dustytorus surrounding the Broad Line Region (BLR), this delaywould imply a distance of this torus from the central engineof . ± . pc.The lightcurve of NIR and optical bands is shownin Figure 9, and the apparent magnitudes after the host-subtraction are listed in Table A1. The lightcurve shows asmooth decay, consistent with the behaviour observed in theGaia lightcurve. There is a bump in the H band around57800 MJD, this is around 100 days after the bump in theGaia lightcurve at 57700 MJD, a value that is consistentwith the delay present in the NEO-WISE lightcurve. Thecolor evolution is shown in Fig. 10: all colors have been fit-ted both with a constant value and with a slope. The reduced χ in all cases is well below 1, meaning that the uncertain-ties on the color are too large to determine an evolutionarytrend of the colors. We also performed a KS test to assessthe goodness of fit and found that the colors B-V and r - i are better fit with a slope than with a constant. Therefore whilethe g - r and u - g remain constant over the period of obser-vation, the color B-V decreases from − . ± . to − . ± . and r - i increases from . ± . to . ± . .The host-subtracted and extinction corrected magni-tudes were used to model the Spectral Energy Distribution(SED) and calculate the bolometric luminosity of the objectusing a python code adapted from the program superbol (Nicholl 2018). For all filters the extinction correction wasdone using the values in Schlafly & Finkbeiner (2011) whoassume a reddening law with R v =3.1. Our NIR and opticaldata are not coincident in time, therefore we focused on theepochs in which we had optical data, adding the NIR pointsto better constrain the fit where possible, extrapolating orinterpolating the NIR data to match the epochs of opticalimaging. However, since our first NIR observations (MJD57 509.04) is almost coincident with our fourth epoch of op-tical data (MJD 57 505.03), we decided to add NIR dataonly from MJD 57 491.01 (the third epoch of optical imag-ing) onward, to avoid large extrapolation of the NIR data.The total bolometric luminosity is then calculated by con-necting the points in the optical wavelenghts with straightlines and integrating the area under the resulting curve, plusa black-body extrapolation in the UV and NIR regions. Forevery epoch we applied a K-correction to the data and thenfit the fluxes derived from each band with two black-bodies.An example of the fit is shown in Figure 11The absence of UV data-points hinders our ability toconstrain the peak of the black-body emission. Using theUV data from our single XMM -Newton observations we cancheck our choice of using two black-bodies instead of onefor the optical and NIR observations closest in time to the
XMM-Newton observation: subtracting the pre-outburst UVflux measured by
Swift from the values measured from the
XMM -Newton observation for the UVM2 and UVW1 fil-ters (central wavelength 2310 ˚A and 2910 ˚A, respectively),we obtain a luminosity density of ∼ × erg s − ˚A − and ∼ × erg s − ˚A − , respectively. These values are higherthan the peak of our fit, suggesting that by using two black-bodies we are not overestimating the emission from the ob-ject.Modeling the SED with two black-body curves does notyield a good fit (reduced χ ∼ − for all epochs). On the onehand, our NIR and optical magnitudes are not coincident intime, thus introducing some systematic uncertainty in ourSED. On the other hand, there is probably a non-thermalcomponent that we are not considering in our fit. This is alsohinted at by the high values of the UV and X-ray luminosi-ties. The absence of UV data, though, hinders our ability toconstrain the parameters of a power-law that would prob-ably describe the non-thermal component. On top of this,we see emission from more than one component: the opticaland the (delayed) infrared emission (see Fig.8). All in all, weassume that our fit using two black-bodies is a reasonablefirst order description of SED and allows for an estimate ofthe bolometric luminosity.The bolometric luminosity resulting from the black-body fits is plotted in Fig. 12. The bolometric luminosityof the object decreases by a factor of 2 over ∼
300 daysand then remains constant within uncertainties at around ∼ × erg s − . The temperature and radius for both black-bodies remain constant during all the period of observation. MNRAS000
300 daysand then remains constant within uncertainties at around ∼ × erg s − . The temperature and radius for both black-bodies remain constant during all the period of observation. MNRAS000 , 1–21 (2019)
G. Cannizzaro et al.
The peak of the two black-bodies curves is at ∼ ∼ µ m . The black-body temperature of the second black-body component (associated with the IR emission) is high,above 2000K, indicating the presence of dust possibly abovethe evaporation temperature.We calculated the bolometric correction to the abso-lute magnitude from Gaia with the values of the bolometricluminosity obtained from the photometric measurements.We then calculated the average of the correction over theepochs of observations and, under the assumption that thiscorrection remains constant, used it to estimate the bolo-metric luminosity from the Gaia absolute magnitudes overthe duration of the whole outburst. With this method wewere able to calculate the estimated total energy irradiated: E tot = ( . ± . ) × erg . The uncertainty on this measurehas been calculated by propagating the uncertainty on thebolometric luminosity from our photometric measurements.As the Gaia magnitudes are given without uncertainty, theerror on the total energy radiated during the outburst de-rives only from the error on the bolometric luminosity wecalculated.The bolometric luminosity evolution as a function oftime calculated from the Gaia magnitudes mimics the shapeof the Gaia lightcurve by design, therefore the decay trend ofthe bolometric luminosity is also well fit by an exponentialdecay ∝ t − . ± . .We use the luminosity in the X-ray, calculated fromthe XMM observation, in the band 2–10 keV, L −
10 keV = ( . ± . ) × erg / s , to have another estimate of the bolo-metric luminosity. To do this, we calculate the values of thebolometric correction ( k bol ) using the methods describedin Marconi et al. (2004) and Netzer (2019). We obtain k bol (cid:39) , using the equation of Marconi et al. (2004) and k bol (cid:39) . from Netzer (2019). The resulting luminosityvalues are L = . × erg s − and L = . × erg s − , re-spectively. The uncertainties in these values are dominatedby the uncertainties in the evaluation of the bolometric cor-rection, but they are not easy to calculate. The range inthe extreme values of k bol can be as large as an order ofmagnitude (Netzer 2019) and depends strongly on the con-straints used to calculate k bol . We note that the value of thebolometric luminosity estimated from the X-ray luminosityis higher than the one calculated from the SED fit. This isin agreement with the statement that our SED fit probablyunderestimates the bolometric luminosity, given the absenceof UV and X-ray data points. If we, instead, compare the L − keV with the bolometric luminosity obtained from theSED fit, we would get a value of k bol ∼ . . This value ismuch lower than what reported in the literature (Brightmanet al. 2017; Netzer 2019; Marconi et al. 2004). We thereforeconsider that the value of the total energy radiated over theoutburst, obtained from the SED fit, is underestimated. We used the first spectrum (MJD 57 428.27) to classify thetransient event. This was done by cross-correlating our clas-sification spectrum with a library of spectra using the
SNID code (Blondin & Tonry 2007), resulting in a good matchwith an AGN spectrum (Mattila et al. 2016). The redshiftcalculated from the transient spectrum is consistent with theone reported in SDSS.
Figure 9.
Light curve of the NOT and LT observations in theoptical and NIR bands, compared with the Gaia lightcurve. Thevertical dashed line at MJD 57 413 corresponds to the date thealert was published while the vertical dashed line at MJD 57 555and the gray area indicate the time frame in which the NEO-WISE lightcurve reaches its peak. The magnitudes plotted arethe ones resulting from the host-subtraction process.
Figure 10.
Color evolution of r - i (black diamonds), g - r (magentacircles), B-V (green squares) and u - g (red triangles). For all thecolors both a fit with a constant line (black dotted lines) and a fitwith a slope (blue dashed lines). The colors are calculated fromthe host-subtracted magnitudes. In Fig.4 we plot all the spectra of our follow-up cam-paign. Looking at the follow-up spectra by eye, the first op-tical spectrum shows a blue component, while during thedecay of the outburst the spectra become redder. BroadBalmer lines (H α , H β , H γ ) are present in all spectra, as wellas narrow [O iii ] and [N ii ] lines. These lines are present inthe SDSS quiescent spectrum, but both their intensity andstructure appear to be changed in the outburst spectra. Thisis likely due to the difference of spectral and spatial resolu-tions of the spectrographs employed by SDSS and NOT. Theregion of the H β broad emission line is complex with vari-ous broad components and narrow [O iii ] ( λλ α line is contaminated by the telluricabsorption around 8250 ˚A; this hindered our ability to fit indetail the properties of the emission lines in this wavelengthregion. At the first epoch there is a bump in the contin-uum between H γ and H β . This is possibly due to a complex MNRAS , 1–21 (2019) aia16aax Figure 11.
Fit of two black-body curves to the luminosity den-sity. The black squares are the values of the luminosity for eachfilter (optical - MJD 57505 and NIR - MJD 57509), the black cir-cles are the two values of the monochromatic luminosity for thetwo UV filters UVM2 and UVW1 (MJD 57569, not considered forthe fit), the solid magenta line represent the fit using two blackbody components, these two components are represented by theblue and red dotted lines. On the x-axis there is the rest-framewavelength. of unresolved lines such as that caused by the Bowen fluo-rescence mechanism (Bowen 1934, 1935) blended with He ii at 4686 ˚A. Metal lines originating from the Bowen fluores-cence have been recently reported to be present in TDEs(Leloudas et al. 2019; Blagorodnova et al. 2019; Onori et al.2019). In the subsequent spectra, this blend is not presentanymore.The absence of other strong lines typically associ-ated with the Bowen fluorescence mechanism, such as N iii λλ iii ] λ pPXF (Cappellari & Emsellem 2004; Cappel-lari 2017). The method approximates the galaxy spectrumby convolving a series of N template spectra T ( x ) with an ini-tial guess of f ( v ) , the Line Of Sight Velocity Dispersion func-tion (LOSVD) to the observed spectrum. The galaxy modelis obtained through this parametrisation (in pixel space x ): G mod ( x ) = (cid:205) Nn = w n {[ T n ( x ) ∗ f n ( cx )] (cid:205) Kk = a k P k ( x )} + (cid:205) Ll = b P l ( x ) + (cid:205) Jj = c j S j ( x ) , (1)where the w n are the spectral weights, the P k and P l aremultiplicative and additive orthogonal polynomials and S j are the sky spectra. The polynomials and sky spectra areoptional components of the parametrisation. The LOSVD Figure 12.
Results from the black-body fits to the flux densityat each epoch. In the top panel we plot the bolometric luminosity,in the second and third panel we plot the temperature for the twoblack-bodies and in the fourth and last panel the radius for thetwo components (blue squares for the first black-body and redcircles for the second one, in all relevant panels. The first black-body is the one that peaks at shorter wavelengths, the secondone is the one that peaks at larger wavelengths, see Fig. 11).The large uncertainties on the black-body parameters in the firsttwo epochs (especially for the radius of the second black-bodycomponent) are due to the absence of NIR data. f ( cx ) = f ( v ) is parametrised by a series of Gauss-Hermitepolynomials as: f ( v ) = σ √ π exp (cid:32) ( v − V ) σ (cid:33) (cid:20) + (cid:205) Mm = h m H m (cid:18) v − V σ (cid:19)(cid:21) , (2)where V is the mean velocity along the line of sight, σ is thevelocity dispersion, H m are the Hermite polynomials and h m their coefficients. The best fitting template is then foundby χ minimisation. To compute our synthetic galactic spec-trum, we chose the The Indo-US Library of Coud´e Feed stel-lar spectra (Valdes et al. 2004), a library of 1273 stellar spec-tra covering a broad range of parameters (effective tempera-ture, metallicity, surface gravity) with a wavelength range of3460 – 9464 ˚A and a spectral resolution of 1.35 ˚A (FWHM), σ ∼
30 km s − , R ∼ pPXF method convolvingthe library to our last spectrum, taken after the object wentback to its pre-outburst state. We excluded from the fit theregions in which the AGN lines are present, the regions oftelluric absorption and the edges of the spectrum. The tem- MNRAS000
30 km s − , R ∼ pPXF method convolvingthe library to our last spectrum, taken after the object wentback to its pre-outburst state. We excluded from the fit theregions in which the AGN lines are present, the regions oftelluric absorption and the edges of the spectrum. The tem- MNRAS000 , 1–21 (2019) G. Cannizzaro et al.
Figure 13.
Fit of the stellar libraries to our last spectrum (2018July 08, 894 days after peak) using the pPXF method. In black theobserved spectrum, in red, the best fit of the templates from thestellar library. The grey bands represent the areas excluded fromthe fit. In green and blue the subtracted spectrum. plates fit is shown in Fig.13. We then subtract the galaxyspectrum obtained with this method from the spectra atall other epoch. The analysis presented from here onward isperformed on the galaxy-subtracted spectra.To analyse the complex emission line structure we fitthe spectra with a combination of Gaussian function us-ing a python code employing the lmfit package (Newvilleet al. 2014). The results of the line fits are listed in TableA2. During the outburst, except for the last two epochs,the H β emission line complex can be well described by twoGaussian components: component A and a red wing, com-ponent B (H β A and H β B from here on), see Fig. A1. H β A is blue-shifted with respect to the rest-frame wavelength ofH β (the shift is between ∼ ∼
20 ˚A, corresponding toa velocity between ∼
300 and ∼ − ) while the redwing is red-shifted with respect to the restframe wavelengthby a factor that varies between ∼
50 and ∼
80 ˚A (correspond-ing to a velocity between ∼ ∼ − ). In theH β region there are also two [O iii ] emission lines at theirrestframe wavelength (4959 ˚A and 5007 ˚A). The separationbetween the [O iii ] lines and their FWHM have been keptfixed during the fit.The presence of the telluric absorption that falls on topof the H α region at the redshift of the object made a preciseanalysis of the H α region difficult. Overall, the H α emis-sion line complex displays a morphology similar to the oneof H β (see Fig. A2): two components of which one (com-ponent A, H α A from here on) is blue-shifted with respectto the rest-frame wavelength (shift between ∼ ∼
15 ˚A,corresponding to a velocity between ∼
300 and ∼
700 km s − )and a red wing (component B, H α B from here on) whichis shifted by around 100 ˚A (corresponding to a velocity of ∼ − ).In addition to these two components, we fit also the two[N ii ] lines (6550 and 6585 ˚A) and [S ii ] (6718 and 6733 ˚A, https://lmfit.github.io/lmfit-py/ unresolved in our spectra). To reduce the numbers of freeparameters in the fit of the H α region, we constrained theFWHM of the narrow lines to be the same as the value of[O iii ] in the fit to the H β region at the same epoch, underthe assumption that the lines coming from the Narrow LineRegion do not change on the timescale of the outburst we areanalysing in this manuscript and the assumption that linesof different metals have the same FWHM. The separationbetween the two [N ii ] lines has also been kept fixed duringthe fitting procedure. On the spectrum of 2016 December09 (MJD 57 732) a cosmic ray hit the detector exactly onthe H β A component, therefore we did not display the resultsfrom the line fitting at this epoch for both H β components.During the time covered by our observations, the line param-eters show variations above the statistical noise, without aclear evolution with time.At the first two epochs (see Figures A2a and A2b) thereis a narrow component (FWHM ∼ km s − ) on top of thered wing of H α , centered at ∼ i i ii α line in all epochs areshown. It is worth noting that the H α line retains the doublepeaked nature until our last spectrum, taken after the Gaialightcurve reached the pre-outburst level of emission, whilethe H β is well described by only one component in the lasttwo epochs. This could be due to the lower signal to noiseratio of these spectra and to the lower contrast between theH β and the continuum, with respect to the H α line.In Figure 14 we show the position of the source ina Baldwin, Phillips & Terlevich (BPT) diagram (Baldwinet al. 1981). For the calculation of the ratios of the fluxes weintroduced a narrow H α and H β component. To do this, wefirst fitted this additional component in the first spectrum(the one with the highest signal to noise ratio), constrain-ing its FWHM to be the same of the other narrow lines atthe same epoch. In the subsequent spectra, we kept the pa-rameters of the narrow H α and H β component fixed, withinuncertainties, under the assumption that the narrow lines donot change over the timescale of our follow-up. We measuredthe ratios for the NOT spectra taken during the outburst de-cay and our last spectrum taken when the source was backin its pre-outburst state. At all epochs in which the narrowH α and H β components could be constrained, the posi-tion in the BPT diagram is consistent with the AGN area ofthe diagram, meaning that the NRL is dominated by AGNionisation. We use the M − σ ∗ relation (Ferrarese & Ford 2005) to ob-tain an estimate of the BH mass of Gaia16aax, using thewidth of the Ca H+K lines to estimate the velocity disper-sion of the galaxy. For this calculation, we use the secondspectrum of our follow-up, as it is the one with the high-est signal to noise ratio (SNR), after the first classificationspectrum, but less contaminated by the outburst light. Af-ter correcting for the instrumental broadening (16.2 ˚A for a1.0 (cid:48)(cid:48) slit), we calculate σ = ± km s − and a BH mass of MNRAS , 1–21 (2019) aia16aax Figure 14.
Plot of the position of the source on a Baldwin,Phillips & Terlevich (BPT) diagram (Baldwin et al. 1981) us-ing the narrow emission lines detected in the host-spectrum ofGaia16aax, during the outburst (blue squares) and after it wentback to the pre-outburst state (magenta triangle). In all cases theposition is consistent with the AGN region of the diagram. The lines that separate the different activity regions come fromKewley et al. (2001) (dotted line) and Kauffmann et al. (2003)(solid line). M bh = ( . ± . )× M (cid:12) . The uncertainty on this measure islarge, considering the modest SNR of the spectrum and theresolution of the instrument. The scatter in the M − σ ∗ rela-tion, which is 0.34 dex (Ferrarese & Ford 2005), contributesto the uncertainty of this mass estimate. The resulting Ed-dington luminosity is L Edd = ( . ± . ) × er g / s .We are able to use the same method also for Gaia16ajq,where the Ca H+K lines are visible in the outburst spec-trum. We obtain a BH mass of M bh ∼ . × M (cid:12) . In thecase of Gaia16aka, the Ca H+K lines are not visible, there-fore we use two Single Epoch (SE) scaling relations betweenthe continuum luminosity, the width of the H α (Greene et al.2010) or H β (Vestergaard & Peterson 2006) line and the BHmass. To obtain the continuum luminosity of the AGN andthe width of the line, we first subtract the galactic compo-nent from the spectrum of Gaia16aka, using the pPPXF proce-dure described before. We then fit the subtracted spectrumwith multiple Gaussian components. From the SE relationdescribed in Greene et al. (2010), using the width of theH α line, we obtain a BH mass of M bh ∼ . × M (cid:12) , whilefrom the relation between the BH mass and the width ofthe H β line we obtain M bh ∼ . × M (cid:12) . It is important tonotice that these single epoch relations are intrinsically un-certain and may only yield an order of magnitude estimate(Vestergaard & Peterson 2006). Before discussing in detail the possible interpretations ofGaia16aax, it is worth summarising the main properties ofthe event: the transient is coincident with the nucleus of agalaxy that hosts a radio-quiet QSO ( ∼ . (cid:48)(cid:48) separation),with an inferred M bh of ∼ × M (cid:12) . The source bright-ened in the optical by 1 magnitude over a timescale of ∼ M G =-22.17 andstarted a smooth decay with t − . ± . , going back to itspre-outburst level of emission over 2 years. Flaring activityof this magnitude was detected both in the MIR (with a timelag of ∼
140 days) and in the X-rays. The first spectra showeda strong blue continuum, while in the subsequent epochs theobject became redder. In the first spectrum a blend of linesbetween H γ and H β is present, this could be due to a blend ofHe ii and Bowen fluorescence emission lines. During the out-burst, the Balmer lines underwent a dramatic change in theirmorphology, showing a clearly separated double-peaked pro-file. The narrow lines did not change over the period of ourobservations. In the first epochs of observations the spectrashowed a strong blue continuum that disappeared with time.The red wings of both the H α and H β lines are offset fromtheir rest-frame wavelength by several thousand km s − andthis shift does not vary significantly in time. The width ofthe lines and their equivalent widths do not show significantevolution with time. The SED cannot be satisfactorily fitusing a single black body, and a reasonable fit, although notperfect, is obtained using two black body components. The enhanced emission and change of appearance of theemission lines could have been caused by variability in theaccretion disk. While typical AGN variability is of tenths ofmagnitudes over timescales of some years (van Velzen et al.2011), more extreme cases of variability have been found. InRumbaugh et al. (2018) a big sample ( ∼ > Edd < ≤ days. It is of interest to compare the timescales found inGaia16aax, (i.e. a rise timescale of hundreds of days and adecay timescale of years) with the characteristic timescalesassociated with different aspects of accretion disks in AGN.Lawrence (2018) has recently refocused attention ona long-standing problem (see also, Antonucci 2018) that
MNRAS , 1–21 (2019) G. Cannizzaro et al. quasar variability is not compatible with standard accre-tion disk theory. The problems are succinctly summarised byDexter & Begelman (2019). In agreement with Lawrence’sconclusion that ”the optical output we see comes entirelyfrom reprocessing of a central source”, Kynoch et al. (2019)argue that the most favoured explanation of the observationsis that the variability at all wavelengths can be accountedfor by an intrinsic change in the luminosity of the centralobject (and the brightening observed in Gaia16aax in UVand X-rays seems to confirm this) – presumably the centralregions of the black hole accretion disk. Assuming that thisis so, then it appears that the most severe problem is thetimescale of the variability (see also Stern et al. 2018).Such large amplitude variations in luminosity must becaused by variations in the accretion rate in the inner ac-cretion disk (e.g. within R < − R g , where R g = GM / c is the gravitational radius). At such radii for the more lu-minous quasars, standard accretion disk theory (Shakura &Sunyaev 1973) finds that the disks are radiation pressuredominated. In that case the standard viscous timescale forradial inflow t ν ≈ R / ν , where ν is the effective kinematicviscosity, can be written (e.g. Pringle 1981) t ν ≈ ( Ω α ) − ( H / R ) , (3)where Ω is the local disk angular velocity, and H / R thedisk aspect ratio, or, equivalently, in this case (cf. Dexter& Begelman 2019) t ν ≈ (cid:16) α . (cid:17) − (cid:18) κκ T (cid:19) − (cid:18) M M (cid:12) (cid:19) (cid:18) (cid:219) m . (cid:19) − (cid:18) R R g (cid:19) / d . (4)Here α is the usual viscosity parameter, and we haveadopted the observed value of α ≈ . for fully ionised disks(Martin et al. 2019; King et al. 2007), κ is the opacity, κ T the electron scattering opacity, M the black hole mass, and (cid:219) m = . (cid:219) Mc / L Edd is the dimensionless accretion rate, where (cid:219) m = gives a luminosity at the Eddington limit, L Edd for anassumed radiative efficiency of (cid:15) = . .Thus we see that while the inflow timescale can be asshort as months at very small radii, the timescale dependsstrongly on radius. Thus, for example, in order to accountfor changing look behaviour in which to whole of the innerdisk regions (say, out to R ≈ − R g ) need to be removedand/or added would require, according to Equation (4) atimescale of some 30 – 370 years. In addition to the timescale problem with standard disks,there is a more fundamental problem, which is that the rea-son for such disks to be variable at all has not been iden-tified. Most of the ideas on what might cause fluctuationsis the disk mass flow centre on local variability within thedisk. Kelly et al. (2009) and Dexter & Agol (2011) pro-pose ad hoc thermal fluctuations within the disk (see also,Nowak & Wagoner 1995; Ruan et al. 2014). Ingram & Done(2011) propose ad hoc local fluctuations in the accretion rate, (cid:219) M . More physically, fluctuations in local magnetic processeshave been proposed by Poutanen & Fabian (1999) (flares inthe corona), Hawley & Krolik (2001), (instabilities in theplunging region of the inner disk) Hogg & Reynolds (2016),(fluctuations caused by local dynamo processes) and Riols et al. (2016) (cyclic dynamo and wind activity). As pointedout by King et al. (2004), all these ideas have the samefundamental drawbacks: the timescales for the fluctuationsare too short (typically a few local dynamical timescales),in that they are much less than the local inflow timescaleand so do not propagate far radially. Thus, being essentiallylocalised, they all produce low amplitude fluctuations. The problem of producing large amplitude fluctuations inthe luminosity of the disk requires a large amplitude fluctu-ation at small radii (where most of the accretion energy isreleased) but on a timescale much longer than the dynami-cal timescales at those radii (typically hours to days). Thisproblem has been addressed by King et al. (2004), thoughin a different context. King et al. drew on an earlier idea byLyubarskii (1997) who pointed out that if some mechanismcould be found for varying the accretion rate at large radius, on the timescale for inflow at that large radius , then the re-sulting accretion rate fluctuation would be able to propagateto small radii, and so produce a large amplitude luminosityfluctuation on that timescale. This is the underlying basisof the ad hoc fluctuation analysis of Ingram & Done (2011,2012).King et al. (2004) proposed a physical model which wasable to generate accretion rate fluctuations in the disk ontimescales much longer than the local dynamical timescale.The basis of the proposal is that from time to time at leastsome of the angular momentum of disk material is removedby a magnetic wind (cf. Blandford & Payne 1982) and notjust by disk viscosity. In order for such a wind to be effectiveat a given radius, it is necessary that there is a strong enoughpoloidal field component threading the disk at that radius.Lubow et al. (1994a) demonstrated that such a field cannotbe transported inwards by the accretion flow itself, because,for an effective magnetic Prandtl number of unity (likely forMagneto Hydro-Dynamic turbulence), the field can moveradially through the disk on a timescale t B ≈ ( H / R ) t ν , thatis with a speed v B ∼ ( R / H ) v R , where v R ≈ ν / R is the usualviscous flow speed. Thus it is more likely that such a field, ifit exists, must be generated locally by local MHD processesinvolving an inverse cascade in order to produce a field withspatial scale ∼ R necessary for a wind, from the dynamodisk scale ∼ H . That such a process can occur in disks hasbeen suggested by Tout & Pringle (1996) and Uzdensky &Goodman (2008).Lubow et al. (1994b) pointed out that such a locally-driven wind can in principle produce an inflow velocity thatis faster than the usual viscous flow speed v R . They pro-posed that if the effect of the wind became locally strongenough that the wind-driven inflow speed exceeded v B , thenthe poloidal field responsible for the wind could be draggedinwards. In that case, they suggested, there is the possibil-ity for a wind-driven avalanche which could sweep the innerdisk regions inwards. Cao & Spruit (2002) and Campbell(2009) concur with this conclusion. An example of how suchan avalanche might operate is to be found in the numericalsimulations by Lovelace et al. (1994).From the point of view of timescales, this mechanism (ifand when it occurs) has the advantage that it occurs on atimescale shorter than the usual viscous timescale (Equation MNRAS , 1–21 (2019) aia16aax
4) by a factor of H / R . that is, on a timescale t B ≈ (cid:16) α . (cid:17) − (cid:18) κκ T (cid:19) − (cid:18) M M (cid:12) (cid:19) (cid:18) (cid:219) m . (cid:19) − (cid:18) R R g (cid:19) / d . (5)At a radius of R ≈ R g this corresponds to a time-scale ofabout a year [335 days].What remain, of course, very uncertain are the time-scales on which such avalanche episodes might recur, andhow these might depend on specific disk properties. Kinget al. (2004) proposed a specific (speculative and, of course, ad hoc ) model, which has some physical basis. They proposethe basic idea that each disk annulus of width ∼ H acts asan independent producer of vertical field B z which variesstochastically on the local dynamo timescale ∼ − Ω − .They suggest, further, that occasionally enough ( ∼ R / H )neighbouring annuli have B z aligned to produce a poloidalfield with length-scale ∼ R , and that when that happensangular momentum can be lost to a locally produced mag-netically driven wind. They find that most of the time suchprocesses give rise to frequent small amplitude luminosityfluctuations, but they do give one example of one such large-scale fluctuation which began at a few hundred R g , whichswept rapidly inwards, reducing the overall disk surface den-sity by almost an order of magnitude. To comprehensively survey all other plausible scenarios, wealso consider if the observed outburst could be explained bya supernova explosion in the nuclear region of the galaxy. Asupernova explosion in the proximity of the nucleus (or su-perimposed spatially on our line of sight) would not explainthe shape change of the Balmer lines. On top of this, the ab-solute magnitude of the peak of the outburst ( M G (cid:39) − ) isuncomfortably high to be caused by a SN and SNe are rarelyseen emitting at X-ray wavelengths (Dwarkadas & Gruszko2012) and would not explain the observed enhanced X-rayluminosity. Finally, even with efficient conversion of kineticenergy to radiation, the total radiated energy of Gaia16aaxstrains most reasonable supernova scenarios. Nonetheless,It is worth noting that Kankare et al. (2017) discuss a SNorigin for the nuclear transient PS1-10adi, that happened inthe nucleus of a galaxy hosting an AGN and radiated a totalenergy of ∼ . × erg. Microlensing events have been invoked to explain highlyvariable AGNs: in Graham et al. (2017), 9 of the 51 ob-jects analysed are well described by a single-lens model andLawrence et al. (2016) have proposed that microlensing pro-vides a good description for many of the objects in theirsample.A microlensing event could in principle explain the in-tensity of the outburst, but it would not be possible to ex-plain the change in shape of the Hydrogen lines and theirtime evolution via this phenomenon alone. On top of this,if a microlensing event would indeed be the cause of the en-hanced emission, we would expect the lightcurve to be sym-metric, while in our case the decay is much shallower thanthe rise to peak. The presence of a secondary bump at late times in the Gaia lightcurve also disfavours a microlensingevent.
For completeness, we also investigate whether the changeof appearance in Gaia16aax could be due to variable ab-sorption in our line of sight. In the case of Gaia16aax it iseasy to see why this is not a viable explanation with: us-ing the method described in MacLeod et al. (2016) from themonochromatic luminosity of the QSO at 5100 ˚A, we can es-timate the size of the BLR to be R
BLR ∼
10 light days (with λ L λ ( ) (cid:39) × erg s − , from the SDSS spectrum) or ∼ . × cm ∼ ∼
100 days (the rise timeof the outburst), that means that the cloud must travel at(10/100)c or v (cid:39) × km/s. To obtain the required red-dening of A v (cid:39) N H (cid:39) × cm − ,assuming a spherical cloud this would lead to a volume den-sity n (cid:39) N H / R BLR (cid:39) . × cm − . Such an object wouldcorrespond to a dense core within a molecular cloud (Blitz& Williams 1999). Assuming that the object is in virialequilibrium, it must have an internal dispersion velocity of σ (cid:39) ( GM cloud R − ) / (cid:39) ( Gm p N H R BLR ) / (cid:39)
30 km / s . Thedifficulty of this scenario would be to have this cloud moveat the required speed of ∼ f D L , with 0 < f < D L theluminosity distance (1.26 Gpc), then the required size of thecloud is R= f R BLR and the required velocity is v = f × km s − . The column density N H must be kept fixed becauseof the required reddening value, therefore the volume den-sity (assuming a spherical cloud) must increase accordingly: n (cid:39) . f − × cm − . To maintain the cloud internal virialsupport we will have σ ∝ f / .If we consider a cloud in the outskirts of an interveninggalaxy with a velocity of v ∼ km s − , we would need f = × − and a cloud of size R = . × cm (cid:39) anda density of n (cid:39) × cm − . Considering instead a cloudin the local group we would need f = − to have d = kpc. Then we would have v (cid:39) km s − and n (cid:39) . × cm − .If the cloud is instead inside the Milky Way, we would con-sider f = − , d = . kpc and n (cid:39) . × cm − . Inall these cases the requirements on the velocity and/or thevolume density of the intervening clouds make this scenarioimprobable. The fact that the outburst was observed withsimilar amplitude in the NIR and X-rays also plays againstan absorption scenario: we would expect some form of repro-cessing of the radiation at different wavelengths. In additionto this, the source is undergoing an outburst rather than adimming episode and it is the first time it has been observedin such a state. This would mean that we would need a cloudwith the described properties constantly obscuring the cen-tral engine, until a ’hole’ with the right density gradientwould expose the source to us. On top of this, it is unclearhow a variable absorption scenario would explain the ap-pearance of the double peaks in the H α and H β emissionlines. The delay observed between the peak of the Gaia andWISE lightcurves also disfavours this variable absorptionscenario. In fact, if we suppose that the enhanced emission MNRAS000
30 km / s . Thedifficulty of this scenario would be to have this cloud moveat the required speed of ∼ f D L , with 0 < f < D L theluminosity distance (1.26 Gpc), then the required size of thecloud is R= f R BLR and the required velocity is v = f × km s − . The column density N H must be kept fixed becauseof the required reddening value, therefore the volume den-sity (assuming a spherical cloud) must increase accordingly: n (cid:39) . f − × cm − . To maintain the cloud internal virialsupport we will have σ ∝ f / .If we consider a cloud in the outskirts of an interveninggalaxy with a velocity of v ∼ km s − , we would need f = × − and a cloud of size R = . × cm (cid:39) anda density of n (cid:39) × cm − . Considering instead a cloudin the local group we would need f = − to have d = kpc. Then we would have v (cid:39) km s − and n (cid:39) . × cm − .If the cloud is instead inside the Milky Way, we would con-sider f = − , d = . kpc and n (cid:39) . × cm − . Inall these cases the requirements on the velocity and/or thevolume density of the intervening clouds make this scenarioimprobable. The fact that the outburst was observed withsimilar amplitude in the NIR and X-rays also plays againstan absorption scenario: we would expect some form of repro-cessing of the radiation at different wavelengths. In additionto this, the source is undergoing an outburst rather than adimming episode and it is the first time it has been observedin such a state. This would mean that we would need a cloudwith the described properties constantly obscuring the cen-tral engine, until a ’hole’ with the right density gradientwould expose the source to us. On top of this, it is unclearhow a variable absorption scenario would explain the ap-pearance of the double peaks in the H α and H β emissionlines. The delay observed between the peak of the Gaia andWISE lightcurves also disfavours this variable absorptionscenario. In fact, if we suppose that the enhanced emission MNRAS000 , 1–21 (2019) G. Cannizzaro et al. comes from de-obscuration rather than an accretion-relatedevent, we would not expect the IR emission to change simi-larly to the optical emission. This is because in this scenario,the IR emission would emerge, in an isotropic fashion, fromthe dust heated by the central engine. If the central emissionis not varying, even if the obscuration in our line of sight ischanging, we would expect the IR emission from the dust toremain unchanged.
Tidal disruption events typically show a fast rise to a peakluminosity around erg s − , followed by a decay thatfollows t − / . The decay time of our outburst is too longto fall into the canonical picture of TDEs, but there havebeen examples of long TDEs, the most extreme cases being adecade-long TDE candidate, see Lin et al. (2017) and a longlived TDE in a merging galaxy pair Arp 299 (Mattila et al.2018). TDEs are known to show broad (up to 10 000 km s − )Hydrogen and/or Helium lines. The absence of color andblack body temperature evolution of Gaia16aax, albeit onthe limited time span of our photometric and SED coverage,is compatible with a TDE scenario (Arcavi et al. 2014).The high mass inferred for the SMBH, though, com-plicates the TDE scenario. In fact, a sun-like star will beswallowed whole by the SMBH if its mass is above M (cid:12) .To explain Gaia16aax as a TDE, the disrupted star mustbe more massive and/or the SMBH must be rapidly spin-ning (Hills 1978, see also the TDE candidate ASASSN-15lh,Leloudas et al. 2016). For a maximally spinning BH, thelimit on the BH mass for a TDE to take place can increaseby an order of magnitude (Kesden 2012). The high spin, ontop of the high BH mass, introduces new arguments in favourof a TDE interpretation. Depending on the BH spin and ifthe to-be-destroyed star’s orbit is prograde or retrograde, theradiation efficiency can vary by a factor of 10 (Bardeen et al.1972), therefore the presence of a spinning BH could justifythe high energy output of the event. On top of this, sim-ulations performed by Guillochon & Ramirez-Ruiz (2015)found that for massive BHs ( (cid:38) M (cid:12) ) the strong generalrelativistic effects produce a rapid circularisation of the de-bris giving rise to a prompt flare, explaining the fast rise topeak of Gaia16aax.If the BH is maximally spinning, the inner disk temper-ature will be higher than in the case of a non rotating BH.We tried to constrain the upper temperature of the hotterof the two black bodies that describe the SED by fittingthe combined UV and optical magnitudes (up to ∼ ∼ K) black body,as the emitting region associated with such a black bodyhas a size of ∼ r g where r g is the gravitational radius.Therefore we cannot constrain the spin of the SMBH.TDEs are usually found and studied in inactive galaxies,although one of the canonical TDEs ASASSN 14li occurredin a low-luminosity AGN (Prieto et al. 2016). The theoreticalpredictions and observational properties of TDEs in galaxiesthat harbor an AGN are less well-constrained. In Chan et al. (2019) it is shown that a stream of debris from a disruptedstar will impact the accretion disk and drain the accretiondisk from the point of impact inwards on timescales shorterthan the inflow time from a Shakura & Sunyaev (1973) disk,with timescales that depend on the orbital parameters ofthe disrupted star. The disk will then replenish itself onthe (longer) inflow timescale. These two different timescalescould explain the difference between the observed rise anddecay times. This scenario could also help to explain theshape of the H α and H β emission lines and the secondarypeaks we observe in the light curve decay. Initially, the de-bris stream coming from the disrupted star slams into theaccretion disk, starting the drainage of the disk material, re-sulting in the enhanced luminosity. Due to the impact, someof the material will splash out at an angle. If we considerthis outflowing material as part of the source of the H α andH β emission lines (the other one being the BLR), the incli-nation of the outflow angle with respect to our line of sightwould explain the double peaked shape of the lines. If theinitial debris stream is dense enough, the stream will piercethrough the accretion disk, and slam into the disk again ata subsequent passage. This multiple encounters between thestream and the disk would help explain the presence of thetwo bumps in the light curve decay at ∼
300 and ∼
600 daysfrom peak.It is worth noting that in the case of a TDE inter-acting with an AGN disk, the interaction between the dis-rupted material and the accretion disk - and thus the ob-served properties - will highly depend on the angle of inci-dence. One could argue that the discovery of Gaia16ajq andGaia16aka could play against a TDE interpretation, giventhe almost identical light curves and spectra of the threeobjects. Nonetheless, there are differences between the threeobjects that could reflect a diversity of encounters betweenthe tidal debris and the accretion disk. Only Gaia16aaxshows two components in the H α and H β emission lines,which in our picture are associated with splashing materialafter the encounter between the debris stream and the ac-cretion disk. Moreover, the light curve of Gaia16aka decayssmoothly, without secondary peaks, while Gaia16ajq showsa very strong bump ∼ days after peak. As stated before,the bumps in the light curve are associated with eventualmultiple interactions between the debris stream and the ac-cretion disk.If the three events are due to TDEs around spinningSMBH, Fig. 6 requires the characteristic timescales involved(e.g. mass return time and cooling time) as well as the totalenergy output to be very similar in each case. The SMBHmasses in these sources vary by a factor of ∼ , implyingcooling and mass return times that vary by a factor ∼ − for similar stellar mass and pericenter (see Chan et al. (2019)for the relevant equations). For the lightcurves to lie on topof each other as in Fig. 6, we would need fine tuning be-tween the stellar mass and radius, the SMBH mass and thepericenter distance of the encounter. It is also possible thatthe Gaia selection criteria mean that we are biased towardsflares of this form. Implying that there could be a large num-ber of flares due to TDEs in AGN that are not being de-tected efficiently by the Gaia Alerts system (e.g. Kostrzewa-Rutkowska et al. 2018).We also note that Chan et al. (2019) suggest thatparabolic TDEs should cause a dip in the AGN X-ray flux MNRAS , 1–21 (2019) aia16aax as the corona is obscured by the tidal tail. In Gaia16aax, theX-ray emission increases by almost one order of magnitudewith respect to the pre-outburst value. If the X-ray flux ob-served is associated with the TDE, it implies the TDE mayhave occurred close to the plane of the AGN disk. The factthat the object is still bright in the UV and X-rays yearsafter the peak is also in line with what seen in other TDEs(e.g. Brown et al. (2017); Jonker et al. (2019)).It is interesting to consider also the case of the TDE Arp299-B-AT1, where Mattila et al. (2018) observed a high totalradiated energy above . × erg. The TDE happened ina nucleus hosting a known AGN and the transient showeda very significant, slowly evolving IR emission coming fromthe dust surrounding the AGN. Another possible scenario that could potentially explain thesimilarity between the flares observed in the three Gaia ob-jects is that of a neutron star (NS) tidally disrupted by astellar mass BH in the AGN accretion disk. The mass of theBH required for the tidal disruption depends on the NS equa-tion of state and, to large extent, on the BH spin, see Lat-timer & Schramm 1976; Shibata & Taniguchi 2011. Whilethe electromagnetic signal coming from the NS-BH mergerwould be too weak to be visible against the high luminosityof the AGN, part of the NS material would be outflowing atrelativistic speeds and hit the Hydrogen rich accretion disk,creating the observed flare. In this scenario, the flare is notdirectly related to the SMBH - as in the case of an accretionevent - therefore it can more naturally explain the similarbolometric luminosity of the three events. A tidal encounterbetween a NS and a stellar mass BH in the proximity ofan AGN accretion disk could be relatively common. Con-sidering a distribution of NS and BH around the SMBH, afraction of these objects will either have orbits on the diskplane, or will have orbits that intersect the disk and thatwill evolve to the orbital plane during the disk lifetime. Dif-ferent objects within the disk will suffer different dynamicalfriction and will encounter each other at low relative veloc-ities, favouring the creation of binaries (Stone et al. 2017;McKernan et al. 2012).In Fig. 15 we show a simple calculation of the energyreleased in such a scenario. On the Y-axis we have the massin the outflow (a fraction of the total NS mass) and on theX-axis the velocity of the outflow. The solid lines representthe kinetic energy ( / M v ) of the outflowing material (with M its mass and v its velocity).In this scenario, to reach an energy release of the or-der of erg, we would need a fraction of the NS ma-terial, around 10%, to be traveling at ∼ M (cid:12) and . c are the maximum values forthe ejecta mass and velocity, respectively.Gaia16aax released a total amount of energy of ∼ × erg, therefore the energy provided by the dynamical ejectais not enough to explain the emitted energy. Barbieri et al.2019 suggests that a jet could also carry O( erg ) and inDeaton et al. (2013), a short-lived neutron-rich accretion Figure 15.
Energy ( / Mv ) as curves for erg (blue solidcurve) and erg (red dashed curve) plotted in mass of theejecta and its velocity. As shown in Barbieri et al. (2019), wedon’t expect more than . M (cid:12) to be available for the NS ejecta,this limit is represented by the gray shaded area in the plot. disk around a rapidly spinning stellar mass black hole canhave thermal energy also of O( erg ). Thus, we regard ∼ × erg as a reasonable upper limit for the energy outputfrom a NS TDE, requiring a jet, fast and relatively massiveand high-velocity ejecta and a rapidly spinning BH.As said, the appeal of this picture lies in that it natu-rally explains the similar bolometric luminosity of the threeevents, something which is more difficult to do if it is relatedto accretion onto the central BH. The energy released dur-ing the Gaia16aax outburst (and, by extension, of the otherGaia objects), though, stretches this scenario, since its en-ergy release is consistent with the upper limit calculated forthis scenario. However, the total energy output of Gaia16aaxcould be underestimated, as our calculation of the bolomet-ric luminosity was done mainly using the optical photometryfrom Gaia. If this is the case, the scenario discussed in thissection must be ruled out. We present multiwavelength follow up of Gaia16aax, a nu-clear transient discovered by the Gaia science alerts pro-gram. The transient is coincident with a quasar at z = . (SDSS J143418.47+491236). The target brightened by morethan 1 magnitude in the optical over less than a year andwent back to its pre-outburst level over more than two years.Variability of similar amplitude has been detected also inthe NIR (with a time delay of ∼
140 days), UV and X-rays.The most striking property of this object is the shape of theH α and H β broad emission lines. These lines were presentin the SDSS pre-outburst spectrum, but during the outburstthey show two distinct peaks, at different shifts with respectto the rest-frame wavelength. Gaia16aax is part of a groupof three nuclear transients discovered by the Gaia sciencealerts program that showed nearly identical light curves.The other two are Gaia16aka and Gaia16ajq, for which adetailed follow-up is not available. The three objects showsimilarities also in their spectra. MNRAS000
140 days), UV and X-rays.The most striking property of this object is the shape of theH α and H β broad emission lines. These lines were presentin the SDSS pre-outburst spectrum, but during the outburstthey show two distinct peaks, at different shifts with respectto the rest-frame wavelength. Gaia16aax is part of a groupof three nuclear transients discovered by the Gaia sciencealerts program that showed nearly identical light curves.The other two are Gaia16aka and Gaia16ajq, for which adetailed follow-up is not available. The three objects showsimilarities also in their spectra. MNRAS000 , 1–21 (2019) G. Cannizzaro et al.
We discuss various scenarios present in the literatureto explain large amplitude flares in AGNs to try to explainthe observed properties of the outburst of Gaia16aax. Theshort timescale of the rise to peak and the total energyoutput are the most challenging properties to explain. Ofthese scenarios, some are easily ruled out: a microlensingevent, a SN explosion, and variable absorption in the lineof sight cannot explain the observed properties. Other phe-nomena are more promising, although none can explain allthe observed properties in a straightforward manner. Theoutburst of Gaia16aax can be explained by a change in theaccretion flow onto the central BH. The low Eddington ra-tio of Gaia16aax is a characteristic shared by many otherquasars that showed variability of similar amplitude. But,in the theoretical framework of accretion disk physics, it isdifficult to explain the rapid rise shown by Gaia16aax, asthe time scales that govern the dynamics of an accretiondisk are much longer than those at play in Gaia16aax. Wereview some proposed mechanisms in the literature for vari-ability in the inner part of the accretion disk, finding that,with the aid of some magnetic wind-driven loss of angularmomentum, a high amplitude variability on timescales of ∼ ≥ M (cid:12) . The presence of star formation may allow forthis. Certainly for the stars of around a solar mass the BHmust be rapidly spinning. This would help explain both thehigh energy release of the event and the short timescale.The interaction between the debris stream and the accre-tion disk could help explain the shape of the light curve andthe shape of the emission lines. The encounter between thedebris stream and the accretion disk should give rise to dif-ferent observable properties in the three Gaia objects, giventhe different SMBH masses and the dependence of such en-counters on multiple parameters (e.g. density of the stream,incidence angle). The similar timescale and peak brightnessof the three transients question the validity of this scenario.We cannot exclude that the detection of the three similarlyshaped flares is a product of the Gaia selection criteria andthat we are missing more TDEs in AGNs. It is important tonote that transient candidates detected by Gaia are vettedby eye before being published (Kostrzewa-Rutkowska et al.2018) and it is therefore very difficult to gauge the selectionfunction of Gaia.We also explore the possibility that Gaia16aax is dueto a NS – BH merger happening in the AGN disk. Thisscenario would help explain the three different transientsindependently of the central SMBH properties, but the hightotal energy output of Gaia16aax is on the high side of whatis conceivably produced in this picture. ACKNOWLEDGEMENTS
GC, PGJ and ZKR acknowledge support from European Re-search Council Consolidator Grant 647208. MF is supportedby a Royal Society - Science Foundation Ireland Univer-sity Research Fellowship. FO acknowledges the support ofthe H2020 European Hemera program, grant agreement No 730970. JH acknowledges financial support from the FinnishCultural Foundation and the Vilho, Yrj¨o and Kalle V¨ais¨al¨aFoundation (of the Finnish Academy of Science and Letters).BM & KESF are supported by NSF grant 1831412. CJN issupported by the Science and Technology Facilities Coun-cil (STFC) (grant number ST/M005917/1). We thank thereferee for their useful comments and suggestions. We ac-knowledge ESA Gaia, DPAC and the Photometric ScienceAlerts Team (http://gsaweb.ast.cam.ac.uk/alerts). NUTS issupported in part by IDA (The Instrument Centre for Dan-ish Astronomy). The data presented here were obtained [inpart] with ALFOSC, which is provided by the Instituto deAstrofisica de Andalucia (IAA) under a joint agreement withthe University of Copenhagen and NOTSA. Based [in part]on observations made with the Nordic Optical Telescope, op-erated by the Nordic Optical Telescope Scientific Associationat the Observatorio del Roque de los Muchachos, La Palma,Spain, of the Instituto de Astrofisica de Canarias. The Liver-pool Telescope is operated on the island of La Palma by Liv-erpool John Moores University in the Spanish Observatoriodel Roque de los Muchachos of the Instituto de Astrofisicade Canarias with financial support from the UK Science andTechnology Facilities Council. The CSS survey is funded bythe National Aeronautics and Space Administration underGrant No. NNG05GF22G issued through the Science Mis-sion Directorate Near-Earth Objects Observations Program.The CRTS survey is supported by the U.S. National ScienceFoundation under grants AST-0909182 and AST-1313422.Based on observations obtained with XMM-Newton, an ESAscience mission with instruments and contributions directlyfunded by ESA Member States and NASA.
REFERENCES
Abolfathi B., et al., 2018, ApJS, 235, 42Alard C., Lupton R. H., 1998, ApJ, 503, 325Antonucci R., 1993, ARA&A, 31, 473Antonucci R., 2018, Nature Astronomy, 2, 504Arcavi I., et al., 2014, ApJ, 793, 38Baldwin J. A., Phillips M. M., Terlevich R., 1981, PASP, 93, 5Barbieri C., Salafia O. S., Perego A., Colpi M., Ghirlanda G.,2019, arXiv e-prints, p. arXiv:1908.08822Bardeen J. M., Press W. H., Teukolsky S. A., 1972, ApJ, 178, 347Becker A., 2015, HOTPANTS: High Order Transform of PSFANd Template Subtraction, Astrophysics Source Code Li-brary (ascl:1504.004)Blagorodnova N., et al., 2019, ApJ, 873, 92Blandford R. D., Payne D. G., 1982, MNRAS, 199, 883Blitz L., Williams J. P., 1999, in Lada C. J., Kylafis N. D., eds,NATO Advanced Science Institutes (ASI) Series C Vol. 540,NATO Advanced Science Institutes (ASI) Series C. p. 3Blondin S., Tonry J. L., 2007, ApJ, 666, 1024Bowen I. S., 1934, PASP, 46, 146Bowen I. S., 1935, ApJ, 81, 1Brightman M., et al., 2017, ApJ, 844, 10Brown J. S., Holoien T. W. S., Auchettl K., Stanek K. Z.,Kochanek C. S., Shappee B. J., Prieto J. L., Grupe D., 2017,MNRAS, 466, 4904Campbell C. G., 2009, MNRAS, 392, 271Cannizzo J. K., Lee H. M., Goodman J., 1990, ApJ, 351, 38Cao X., Spruit H. C., 2002, A&A, 385, 289Cappellari M., 2017, MNRAS, 466, 798Cappellari M., Emsellem E., 2004, PASP, 116, 138MNRAS , 1–21 (2019) aia16aax Chan C.-H., Piran T., Krolik J. H., Saban D., 2019, arXiv e-prints,Deaton M. B., et al., 2013, ApJ, 776, 47Dexter J., Agol E., 2011, ApJ, 727, L24Dexter J., Begelman M. C., 2019, MNRAS, 483, L17Drake A. J., et al., 2009, ApJ, 696, 870Dwarkadas V. V., Gruszko J., 2012, MNRAS, 419, 1515Evans C. R., Kochanek C. S., 1989, ApJ, 346, L13Evans P. A., et al., 2009, MNRAS, 397, 1177Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523Gaia Collaboration et al., 2018, A&A, 616, A1Gehrels N., et al., 2004, ApJ, 611, 1005Graham M. J., Djorgovski S. G., Drake A. J., Stern D., MahabalA. A., Glikman E., Larson S., Christensen E., 2017, MNRAS,470, 4112Greene J. E., Peng C. Y., Ludwig R. R., 2010, ApJ, 709, 937Guillochon J., Ramirez-Ruiz E., 2015, The Astrophysical Journal,809, 166Hawley J. F., Krolik J. H., 2001, ApJ, 548, 348Hills J. G., 1975, Nature, 254, 295Hills J. G., 1978, MNRAS, 182, 517Hodgkin S. T., Wyrzykowski L., Blagorodnova N., Koposov S.,2013, Philosophical Transactions of the Royal Society of Lon-don Series A, 371, 20120239Hogg J. D., Reynolds C. S., 2016, ApJ, 826, 40Ingram A., Done C., 2011, MNRAS, 415, 2323Ingram A., Done C., 2012, MNRAS, 419, 2369Jester S., et al., 2005, AJ, 130, 873Jonker P. G., Stone N. C., Generozov A., van Velzen S., MetzgerB., 2019, arXiv e-prints, p. arXiv:1906.12236Jordi C., et al., 2010, A&A, 523, A48Kankare E., et al., 2017, Nature Astronomy, 1, 865Kauffmann G., et al., 2003, MNRAS, 346, 1055Kelly B. C., Bechtold J., Siemiginowska A., 2009, ApJ, 698, 895Kesden M., 2012, Phys. Rev. D, 85, 024037Kewley L. J., Heisler C. A., Dopita M. A., Lumsden S., 2001,ApJS, 132, 37King A. R., Pringle J. E., West R. G., Livio M., 2004, MNRAS,348, 111King A. R., Pringle J. E., Livio M., 2007, MNRAS, 376, 1740Kostrzewa-Rutkowska Z., et al., 2018, MNRAS, 481, 307Kynoch D., Ward M. J., Lawrence A., Bruce A. G., Landt H.,MacLeod C. L., 2019, MNRAS, 485, 2573LaMassa S. M., et al., 2015, ApJ, 800, 144Lattimer J. M., Schramm D. N., 1976, ApJ, 210, 549Lawrence A., 2018, Nature Astronomy, 2, 102Lawrence A., et al., 2016, MNRAS, 463, 296Leloudas G., et al., 2016, Nature Astronomy, 1, 0002Leloudas G., et al., 2019, ApJ, 887, 218Lin D., et al., 2017, Nature Astronomy, 1, 0033Lodato G., King A. R., Pringle J. E., 2009, MNRAS, 392, 332Lovelace R. V. E., Romanova M. M., Newman W. I., 1994, ApJ,437, 136Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994a, MNRAS,267, 235Lubow S. H., Papaloizou J. C. B., Pringle J. E., 1994b, MNRAS,268, 1010Lyubarskii Y. E., 1997, MNRAS, 292, 679MacLeod C. L., et al., 2010, ApJ, 721, 1014MacLeod C. L., et al., 2012, ApJ, 753, 106MacLeod C. L., et al., 2016, MNRAS, 457, 389Mainzer A., et al., 2011, ApJ, 731, 53Marconi A., Risaliti G., Gilli R., Hunt L. K., Maiolino R., SalvatiM., 2004, MNRAS, 351, 169Martin R. G., Nixon C. J., Pringle J. E., Livio M., 2019, New As-tron., 70, 7Matt G., Guainazzi M., Maiolino R., 2003, MNRAS, 342, 422Mattila S., Harmanen J., Reynolds T., Fraser M., Hodgkin S., Blagorodnova N., Wyrzykowski L., Kankare E., 2016, TheAstronomer’s Telegram, 8669Mattila S., et al., 2018, Science, 361, 482McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MN-RAS, 425, 460Merloni A., et al., 2015, MNRAS, 452, 69Netzer H., 2019, MNRAS, 488, 5185Newville M., Stensitzki T., Allen D. B., Ingargiola A., 2014,LMFIT: Non-Linear Least-Square Minimization and Curve-Fitting for Python, doi:10.5281/zenodo.11813, https://doi.org/10.5281/zenodo.11813
Nicholl M., 2018, Research Notes of the AAS, 2, 230Nowak M. A., Wagoner R. V., 1995, MNRAS, 274, 37Onori F., et al., 2019, MNRAS, 489, 1463Osterbrock D. E., 1981, ApJ, 249, 462Planck Collaboration et al., 2016, A&A, 594, A13Poutanen J., Fabian A. C., 1999, MNRAS, 306, L31Prieto J. L., et al., 2016, ApJ, 830, L32Pringle J. E., 1981, ARA&A, 19, 137Rees M. J., 1988, Nature, 333, 523Riols A., Ogilvie G. I., Latter H., Ross J. P., 2016, MNRAS, 463,3096Ruan J. J., Anderson S. F., Plotkin R. M., Brand t W. N., BurnettT. H., Myers A. D., Schneider D. P., 2014, ApJ, 797, 19Rumbaugh N., et al., 2018, ApJ, 854, 160Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Shibata M., Taniguchi K., 2011, Living Reviews in Relativity, 14,6Skrutskie M. F., et al., 2006, AJ, 131, 1163Stanishev V., 2007, Astronomische Nachrichten, 328, 948Stern D., et al., 2018, ApJ, 864, 27Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946Tody D., 1986, in Crawford D. L., ed., Proc. SPIEVol. 627, Instru-mentation in astronomy VI. p. 733, doi:10.1117/12.968154Tout C. A., Pringle J. E., 1996, MNRAS, 281, 219Urry C. M., Padovani P., 1995, PASP, 107, 803Uzdensky D. A., Goodman J., 2008, ApJ, 682, 608Valdes F., Gupta R., Rose J. A., Singh H. P., Bell D. J., 2004,ApJS, 152, 251Vestergaard M., Peterson B. M., 2006, ApJ, 641, 689Wright E. L., et al., 2010, AJ, 140, 1868van Velzen S., et al., 2011, ApJ, 741, 73
APPENDIX A: ADDITIONAL MATERIAL
MNRAS000
MNRAS000 , 1–21 (2019) G. Cannizzaro et al. (a) 20160209 (b) 20160701 (c) 20160731(d) 20161209 (e) 20170401 (f) 20170701(g) 20180719 (h) SDSS - 20020705
Figure A1.
Evolution in time of the H β line region with the fit to the various components with multiple Gaussian curves. In all thepanels the blue dashed lines are the two H β components, the magenta dashed line is the narrow H β component and the red dashed linesare the two [O iii ] emission lines. The spectra of 2017 January 26 is not shown since a cosmic ray hit the CCD on top of the H β line,making the results of the fit unreliable even after several attempts at correcting the data. Panel g is the H β line region at an epoch whenthe Gaia lightcurve has reached the pre-outburst level. The double-peaked shape of H β shown in outburst is not clearly present anymoreat the last two epochs (panels f and g). The last panel is the SDSS spectrum, obtained well before the outburst start.MNRAS , 1–21 (2019) aia16aax (a) 20160209 (b) 20160701 (c) 20160731(d) 20161209 (e) 20170126 (f) 20170401(g) 20170701 (h) 20180719 (i) SDSS - 20020705 Figure A2.
Evolution in time of the H α line region with the fit to the various components with multiple Gaussian curves. In all panelsthe red dashed lines are the two [N ii ] lines, the blue dashed lines are the two H α components, the magenta dashed line is the H α narrowcomponent, the brown dashed line is the He i line (only present in the first two panels, a and b) and the orange dashed line is the blendof the two [S ii ] lines. Panel h is the spectrum taken after the object returned to its pre-outburst state, according to the Gaia lightcurve.MNRAS000
Evolution in time of the H α line region with the fit to the various components with multiple Gaussian curves. In all panelsthe red dashed lines are the two [N ii ] lines, the blue dashed lines are the two H α components, the magenta dashed line is the H α narrowcomponent, the brown dashed line is the He i line (only present in the first two panels, a and b) and the orange dashed line is the blendof the two [S ii ] lines. Panel h is the spectrum taken after the object returned to its pre-outburst state, according to the Gaia lightcurve.MNRAS000 , 1–21 (2019) G . C a nn i zz a r o e t a l . Table A1.
Magnitude values in the Optical and NIR.MJD u B V g r i z y
J H K s [ d ] (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12)51 671.30 (13) · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± ± · · · · · · ± ± ± ± · · · · · · · · · · · ·
56 938.05 (15) · · · · · · · · · ± ± ± ± ± · · · · · · · · ·
57 462.24 18.94 ± ± ± ± ± ± · · · · · · · · · · · · · · ·
57 476.24 19.12 ± ± ± ± ± ± ± · · · · · · · · · · · ·
57 491.01 19.01 ± ± ± ± ± ± ± · · · · · · · · · · · ·
57 505.03 19.21 ± ± ± ± ± ± ± · · · · · · · · · · · ·
57 509.04 · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± ± ± ± ± ± ± ± · · · · · · · · · · · ·
57 588.95 19.32 ± ± ± ± ± ± ± · · · · · · · · · · · ·
57 611.95 · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± · · · · · · · · · · · · · · · · · · · · · · · · · · · ± ± · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± · · · · · · · · · · · · · · · · · · · · · · · · ± ± ± Notes: (1) MJD date of observations; (2), (5), (6), (7), (8), (9) apparent magnitudes with 1- σ uncertainties in the optical bands u , g , r , i , z and y , in the AB system; (3) and (4) apparentmagnitudes with 1- σ uncertainties in the optical bands B and V, in the Vega system; (10), (11) and (12) extinction-corrected apparent magnitudes with 1- σ uncertainties in the NearInfra-Red filters J, H and K s , in the Vega system. The first three lines are the archival magnitudes from 2MASS (13), SDSS (14) and PanSTARRS (15), all the other values are theresult of the image subtraction procedure. The symbol · · · indicates an epoch in which an observation for the specific filter was not available. Table A2.
Results from the line fitting procedure on H α and H β .MJD ( ) H β A FWHM ( ) H β B FWHM ( ) H β A flux ( ) H β B flux ( ) H β A shift ( ) H β B shift ( ) H α A FWHM ( ) H α B FWHM ( ) H α A flux ( ) H α B flux ( ) H α A shift ( ) H α B shift ( ) [km s − ] [km s − ] [erg cm − s − ] [erg cm − s − ] [km s − ] [km s − ] [km s − ] [km s − ] [erg cm − s − ] [erg cm − s − ] [km s − ] [km s − ]57428.27 3134 ±
206 4775 ±
466 197 ±
19 210 ±
22 -833 ±
104 3593 ±
201 4168 ±
152 4464 ±
541 850 ±
51 318 ±
46 -520 ±
114 4464 ± ±
238 3724 ±
495 316 ±
20 96 ±
15 -394 ±
105 4492 ±
292 3778 ±
120 3856 ±
391 1180 ±
70 414 ±
47 -560 ±
88 4186 ± ±
482 5914 ± ±
34 107 ±
35 -1064 ±
302 2905 ± ±
178 3482 ±
493 526 ±
42 168 ±
26 -667 ±
117 4322 ± ±
445 6192 ±
884 113 ±
15 93 ±
25 -795 ±
133 3514 ±
636 4069 ±
212 6400 ±
704 468 ±
47 338 ±
40 -537 ±
122 5265 ± ( ) · · · · · · · · · · · · · · · · · · ±
141 5294 ±
161 637 ±
34 342 ±
11 -671 ±
110 4766 ± ±
281 5450 ± ±
13 77 ±
18 -655 ±
133 5403 ±
374 4050 ±
114 5301 ±
370 686 ±
51 383 ±
31 -936 ±
97 4371 ± ( ) ± · · · ± · · · ± · · · ±
267 6234 ±
620 297 ±
32 295 ±
32 -1130 ±
114 4955 ± Notes: (2) Full Width Half Maximum; (3) Equivalent Width; (4) shift of the central wavelength with respect to the laboratory wavelength of H α and H β . A negative value correspondsto a blueshift, a positive value to a redshift. (5) In the spectrum taken on 2016 December 10 a cosmic ray hit the CCD at the location of the H β A component, therefore the results ofthe fit on the H β line is not reliable for this epoch. (6) At this epoch, the H β B component is not present anymore. The last observation (MJD 58 307.01) was not analysed as the objectwas back to its ”quiescent”, pre-outburst state at this epoch. M N R A S , ( ) aia16aax This paper has been typeset from a TEX/L A TEX file prepared bythe author.MNRAS000