There are numerous physical processes that lead to variability in the nuclei of galaxies. Every massive galaxy likely houses a supermassive black hole (SMBH; Kormendy & Richstone 1995 , Kormendy & Ho 2013 , and the past several decades have been spent unraveling their accretion and variability processes (for reviews, see, e.g., Ulrich et al. 1997 , Ho 2008 , Heckman & Best 2014 , Yuan & Narayan 2014 , Padovani et al. 2017 , Hickox & Alexander 2018 , Komossa 2018 , Blandford et al. 2019 . Without the ability to spatially resolve the immediate vicinity of the SMBH, other methods must be used to probe accretion physics.
Variable accretion in active galactic nuclei (AGN) is the primary driver of nuclear variability. Most quasars appear to vary in brightness stochastically with statistical properties that can be modeled relatively well by a Damped Random Walk (DRW; e.g., Kelly et al. 2008 , Koz lowski et al. 2010 , MacLeod et al. 2010 , Zu et al. 2013 . Signatures that deviate from DRW behavior, namely periodic or semi-periodic features, have been suggested as possible indicators for a binary system at the galaxy's core (e.g., Komossa 2006) .
This has led to searches for periodic signals in AGN light curves to identify SMBH binaries. For example, Graham et al. (2015a) used the Catalina Realtime Transient Survey (CRTS) to search for sub-parsec SMBH binaries. They reported 111 candidates that showed evidence of periodicity associated with a Keplerian orbit. Another study by Charisi et al. (2016) used the Palomar Transient Factory to identify 33 candidates with evidence for periodic variability. Liu et al. (2015; searched for periodicity in Pan-STARRS1's Medium Deep Survey, ultimately finding one candidate, PSO J185. Other candidates with quasiperiodic/periodic variability include NGC 4151 with an estimated period of ∼16 years (e.g., Oknyanskij 1978 , Pacholczyk et al. 1983 , Guo et al. 2006 , Oknyanskij & Lyuty 2007 , Bon et al. 2012 , and PG 1302-102 with an estimated period of 1,884 days (Graham et al. 2015b) . Simulations have shown that periodic variability is ex- * NASA Fellow † NSF Graduate Research Fellow ‡ Carnegie Fellow § DOE CSGF Fellow pected in the light curves of SMBH binaries at subparsec separations due to a variety of processes, including modulated mass accretion onto the binary (e.g., D'Orazio et al. 2013 , Gold et al. 2014 , Farris et al. 2014 or relativistic Doppler boosting of the minidisks formed as a result of the binary interaction (D'Orazio et al. 2015 ). Aside from low-level variability, AGN can also show outbursts or flares in which the brightness of the AGN increases dramatically for a short period of time before returning to a level of relative quiescence. The best example of quasi-periodic optical flares is the 12-year outburst cycle of OJ 287. These were first reported by Sillanpaa et al. (1988) , who suggested that outbursts are due to perturbations of the primary black hole's accretion disk during pericenter passages of the secondary black hole on a 12 year orbital cycle. Relativistic effects like precession probably alter the orbital geometry so that the events are not strictly periodic (see, Valtonen et al. 2006 , Laine et al. 2020 . The most recent flare detected from OJ 287 brightened in the X-ray, UV, and optical and occurred between April-June 2020, which was consistent with the predictions of the binary black hole model (Komossa et al. 2020 ). Another candidate is IC 3599, which has competing theories for the cause of its X-ray/optical flares. Grupe et al. (2015) proposed accretion disk instabilites as the cause of the recurring ∼9.5 year flares, while Campana et al. (2015) argued for partial tidal disruption events.
Tidal disruption events (TDEs) also lead to variable SMBH activity. In a TDE, a star is torn apart when it passes within the tidal radius of the SMBH at the center of its host galaxy, as the tidal shear forces overwhelm the self-gravity of the star. The disrupted stellar material subsequently produces a luminous transient flare of electromagnetic radiation that we observe as a TDE. Total disruptions from parabolic orbits (Hills 1975; Rees 1988; Evans & Kochanek 1989; Phinney 1989) lead to half of the disrupted star's mass being ejected, while the other half is gravitationally bound and asymptotically returns to pericenter at a fallback rate proportional to t −5/3 .
In the case of a partial TDE, the star survives the encounter with the SMBH and only a fraction of the stellar material is tidally stripped, leaving the stellar core intact. found that fallback rate for partial TDEs at late times becomes steeper than t −5/3 because there is less debris with orbital bind-ing energies close to zero if the stellar core survives the encounter. Coughlin & Nixon (2019) found a fallback rate proportional to t −9/4 , which is effectively independent of the mass of the core that survives the passage close to the black hole. This fallback rate is supported by the hydrodynamical simulations of Miles et al. (2020) . Hydrodynamical simulations also indicate that partial disruptions can repeat, causing episodic mass transfer from the star to the SMBH at every pericenter passage, resulting in a series of low-level flares that repeat on the orbital timescale (MacLeod et al. 2013) . Partial disruptions are most easily achieved for giant stars (e.g., MacLeod et al. 2013 which might also be created as stellar merger products (e.g., Antonini et al. 2011 , MacLeod et al. 2012 .
Most theoretical predictions for TDEs predate any observations of the phenomenon. The first observational claims of TDEs were soft X-ray outbursts from otherwise quiescent galaxies (e.g., Bade et al. 1996 , Komossa & Greiner 1999 , Grupe et al. 1999 , Greiner et al. 2000 , Gezari et al. 2003 , Komossa 2015 . Since then, TDE flares have been discovered at a range of wavelengths, including the hard X-ray (e.g., Bloom et al. 2011; Burrows et al. 2011; Cenko et al. 2012a; Pasham et al. 2015) , soft X-ray (e.g., Komossa & Bade 1999; Donley et al. 2002; Maksym et al. 2010; Saxton et al. 2012) , ultraviolet (e.g., Stern et al. 2004 , Gezari et al. 2006 , Gezari et al. 2008 , Gezari et al. 2009 ), and optical (e.g., van Velzen et al. 2011 , Gezari et al. 2012 , Cenko et al. 2012b , Arcavi et al. 2014 , Holoien et al. 2014a , Vinkó et al. 2015 , Holoien et al. 2016a , Holoien et al. 2016b , Holoien et al. 2018 , Holoien et al. 2019a , Holoien et al. 2019b , Hinkle et al. 2020b , van Velzen et al. 2020 . Due to the intrinsic multi-wavelength properties of both TDEs and AGN, one problem is to identify characteristics that clearly distinguish between the two objects, notably using X-rays . This is becoming more important with the discoveries of more ambiguous transients such as ASASSN-18el (Trakhtenbrot et al. 2019 , Ricci et al. 2020 ) and ASASSN-18jd (Neustadt et al. 2020) .
Here we report the discovery and long-term observation of a series of periodic outbursts from ASASSN-14ko, which is associated with the AGN ESO 253−G003 (z=0.042489, Aguero et al. 1996) . ESO 253−G003 was spectroscopically classified as a Type 2 Seyfert by Lauberts (1982) . In Section 2 we discuss the discovery of ASASSN-14ko and the photometric and spectroscopic data used in this analysis. In Section 3 we discuss the host properties, and in Section 4 we discuss the light curve and the periodic nature of the outbursts. The spectroscopic results are presented in Section 5, and we discuss several theoretical interpretations of these periodic flares in Section 6. For a flat Ω m = 0.3 universe, the luminosity distance is ≈ 188 Mpc and the projected scale is ≈ 0.85 kpc/arcsec. The Galactic extinction is A V = 0.118 mag (Schlafly & Finkbeiner 2011) .
On 2014-11-14.28 UT, the All-Sky Automated Survey for Supernovae (ASAS-SN, Shappee et al. 2014 ) triggered on a nuclear transient associated with ESO 253−G003 at V ∼ 17.0 mag and reported it as ASASSN-14ko (Holoien et al. 2014b (Holoien et al. , 2017 . We will refer to this object as ASASSN-14ko throughout this paper. As reported in Holoien et al. (2014b) , a follow-up spectrum on 2014-11-16 using the Boller and Chivens (B&C, Osip et al. 2004 ) spectrograph on the du Pont 2.5m at Las Campanas Observatory revealed a strong blue continuum and the emission lines were consistent with a Type 2 Seyfert. At the time, the event was considered to be a Type IIn supernova with a blue continuum projected very close to the nucleus of a Type 2 Seyfert, but strong AGN activity was not ruled out as a possibility (Holoien et al. 2014b) .
UVOT and XRT observations by the Neil Gehrels Swift Observatory (Swift hereafter, Gehrels et al. 2004) were taken on UT 2014 -11-16, 2014 -11-19, 2014 -11-21, 2014 -11-23, and 2014 . These observations showed that the central region of the galaxy had significantly brightened in the UV but were consistent with archival magnitudes in the optical. These Swift data also revealed X-ray emission, with fluxes of (2.85±0.8)×10 −13 ergs cm −2 s −1 and (3.1±0.7)×10 −13 ergs cm −2 s −1 on UT 2014-11-16 and UT 2014-11-19, respectively . The X-ray spectrum was consistent with a highly absorbed AGN with a column density of ∼10 23 cm −2 and a luminosity of L X ∼ 3 ×10 42 ergs s −1 (Holoien et al. 2014b) .
As part of ongoing work to examine the long-term behavior of AGN observed by ASAS-SN, a full light curve of ESO 253−G003 was extracted in February 2020. Visual examination of the light curve revealed sixteen flares spaced out roughly evenly over six years, as shown in Figure 1 . The seventeenth outburst in Figure 1 was then predicted and observed. The original ASASSN-14ko trigger corresponds to the second outburst in the series. This initiated the further analysis and photometric and spectroscopic follow-up of ASASSN-14ko which we report here. All photometric data used in this analysis are presented in Table 1 ATLAS o-band data are shown in orange and ATLAS c-band data are shown in cyan. Swift epochs are denoted by dark blue marks, and spectral epochs are shown by orange marks at the bottom of each panel. The fifth-order polynomial fits for each outburst are shown in red along with shaded red regions corresponding to estimates for each peak time and its uncertainty. The magenta vertical lines show the predicted peaks for the model with a period derivative described in Section 4.1. (LCOGT, Brown et al. 2013 ) at five sites around the globe. Each telescope consists of four 14-cm aperture Nikon telephoto lenses with 8. 0 pixels and a 4.5 • × 4.5 • field of view. ASAS-SN's primary objective is to discover supernovae by surveying the entire visible sky every night. The ASAS-SN data shown in Figure 1 includes both V -band and g-band observations. In 2018, the first two ASAS-SN mounts transitioned from Vband to g-band to match the three ASAS-SN telescopes deployed in 2017-2018. The data were reduced using a fully-automated pipeline based on the ISIS image subtraction package (Alard & Lupton 1998; Alard 2000) . Each photometric epoch (usually) combines three dithered 90-second image exposures with a 4.47×4.47 degree field-of-view that is subtracted from a reference image. We then used the IRAF package apphot to perform aperture photometry with a 2-pixel, or approximately 16. 0, radius aperture on each subtracted image, generating a differential light curve. The photometry was calibrated using the AAVSO Photometric All-Sky Survey (Henden et al. 2015) . All low-quality ASAS-SN images of ESO 253−G003 were inspected by-eye, and images with clouds or other systematic problems were removed.
Following the original discovery, we requested Swift UltraViolet/Optical Telescope (UVOT, Roming et al. 2005 ) ToO observations (ToO ID: 33529). Then, after we discovered its periodic nature, we again requested Swift data (ToO IDs: 13836, 13979, 14005) to monitor ASASSN-14ko during quiescence and then during the outburst predicted for UT 2020-05-18.5 (see below). Data were obtained in six filters (Poole et al. 2008) : V (5468Å), B (4392Å), U (3465Å), U V W 1 (2600 A), U V M 2 (2246Å), and U V W 2 (1928Å). We used the HEAsoft (HEASARC 2014) software task uvotsource to extract the source counts using a 16. 0 radius aperture and used a sky region of ∼ 40. 0 radius to estimate and subtract the sky background. This aperture size was chosen to match the ASAS-SN photometry. All fluxes were aperture corrected and converted into magnitudes and fluxes using the most recent UVOT calibration (Poole et al. 2008 , Breeveld et al. 2010 . The UVOT transient, host, and calibration star magnitudes were corrected for Galactic extinction. In order to measure only the transient flux in each epoch, the quiescent host fluxes were subtracted in the same aperture. The host magnitudes are given in Table 2 . We converted the Swift UVOT B and V magnitudes to Johnson B and V magnitudes using the standard color corrections 1 .
There are three additional observations of this galaxy in the Swift data archive under the identification SWIFT J0525.3 − 4600 (Obs ID: 37354, PI: Markwardt). These observations were obtained to investigate the Swift Burst Alert Telescope (BAT) source reported by Baumgartner et al. (2013) in their 70 month catalog and identified with the nearby blazar P KS0524 − 460. However, Lansbury et al. (2017) concluded that ESO 253−G003 was in fact the BAT source rather than PKS 0524 − 460. The three Swift epochs are consistent with the quiescent magnitudes we measure from the later data and there were too few observations to usefully add to the constraints on the times of the outbursts.
The ATLAS survey (Tonry et al. 2018) consists of two 0.5m f/2 Wright Schmidt telescopes on Haleakalā and at the Mauna Loa Observatory. Designed primarily for detecting hazardous asteroids, the telescopes obtain four 30-second exposures of 200-250 fields per night. This corresponds to roughly a quarter of the visible sky. AT-LAS uses two broad-band filters, the 'cyan' (c) filter covering 420 -650 nm and the 'orange' (o) filter covering 560 -820 nm (Tonry et al. 2018) .
The ATLAS pipeline performs flat-field corrections for each image as well as astrometric and photometric calibrations. Reference images of the host galaxy were created by stacking multiple images taken under ideal conditions and this reference was then subtracted from each ATLAS epoch to isolate the flux from the transient. We performed forced photometry on the subtracted ATLAS Table 2 . Host magnitudes of ESO 253−G003 measured using a 16. 0 aperture. The U V W 2-, U V M 2-, U V W 1-and U -band magnitudes were determined from Swift data taken in April 2020 during pre-outburst quiescence. Johnson-Cousins BV R and SDSS gri magnitudes were determined from LCOGT, Swift B and V , and amateur astronomer data also taken in April 2020. The magnitudes were combined by averaging the data weighted by the inverse squares of the uncertainties. All magnitudes are in the AB system.
Magnitude images of ASASSN-14ko as described in Tonry et al. (2018) . ATLAS images taken on the same night were stacked. The resulting ATLAS o-and c-band photometry and 3-sigma limits are included also in Figure 1 .
ASASSN-14ko was observed by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2014 ) during Sectors 4-6, which occurred between 2018-10-18 and 2019-01-07. Similar to the process applied to the ASAS-SN data, we used the ISIS package (Alard & Lupton 1998; Alard 2000) to perform image subtraction on the 30-minute cadence TESS full frame images (FFIs) to obtain high fidelity light curves of this galaxy. This process is fully described in Vallely et al. (2019) .
We construct independent reference images for each sector as opposed to utilizing a single reference image over all sectors to avoid introducing problems created by the field rotations between sectors. The reference images were built using the first 100 good-quality FFIs for each sector. FFIs were considered poor-quality if the sky background levels or PSF widths were above average for the sector. FFIs were also excluded from our analysis if they had data quality flags, or were taken when the spacecraft's pointing was compromised due to instrument anomalies, or when scattered light affected the images.
The measured fluxes were converted into TESS -band magnitudes using an instrumental zero point of 20.44 electrons per second from the TESS Instrument Handbook (Vanderspek et al. 2018) . TESS observes in a single broad-band filter, spanning roughly 6000-10000Å with an effective wavelength of ∼7500Å, and TESS magnitudes are calibrated to the Vega system (Sullivan et al. 2015) . The TESS light curve is also shown in Figure 1 .
We obtained photometric observations from Las Cumbres Observatory Global Telescope (LCOGT, Brown et al. 2013) . The B-, V -, g -, r -, and i -band observations were taken using the 1-meter telescope at Siding Spring Observatory in New South Wales, Australia. The LCOGT photometric observations began on 2020-04-13 in quiescence and continued through the midpoint of the May 2020 outburst when telescope horizon observing limits prevented further observations. Aperture magnitudes were obtained using a 16. 0 radius aperture using the IRAF apphot package using an annulus to estimate and subtract background counts. We used stars with APASS DR 10 magnitudes to calibrate the data. Similar to the process for the UVOT observations, the aperture magnitudes were corrected for Galactic extinction. The host galaxy flux was measured using the April 2020 quiescent data and subtracted to isolate the flux from the transient.
Amateur astronomers at four different observatories observed ASASSN-14ko starting shortly after it was discovered to be periodically flaring. Data were taken at Moondyne Observatory, east of Perth, Australia, using a 0.4-m telescope with AOX adaptive optics between April 28 and June 8, 2020 and on a daily cadence between May 10 and May 28, 2020. B C -, V -, G S -, R S -, and I C -band images were obtained with guided 120 second exposures. The data were reduced and calibrated with standard procedures, and then stacked with 3 and 5 image sets aligned using background stars. Data were also collected using a 41-cm telescope at Savannah Skies Observatory from Queensland, Australia. The bias and dark subtracted data were taken in the B-, V -, R C -, and I C -bands using 180 second exposures. ASASSN-14ko was observed from Bronberg Observatory in South Africa using 14-and 12-inch telescopes in the R band with 15 second exposures. The images were calibrated and then sets of eight images were stacked per night. Finally, ASASSN-14ko was observed from Dogsheaven Observatory in Brazil using the B, V , R C , and I C filters. The images were observed with 120 second exposures using a 14-inch telescope and calibrated.
All these images were then astrometrically calibrated and aperture magnitudes were measured using the IRAF apphot package and a 16. 0 aperture following the same procedures described for the Swift UVOT and LCOGT data.
The first available observation of ESO 253−G003 was obtained by Kewley et al. (2001 Kewley et al. ( ) on 1996 . At that time, ESO 253−G003 was classified as a Seyfert 2 galaxy. When ASASSN-14ko was first discovered, we took a follow-up spectrum using the B&C spectrograph on the du Pont 2.5-m at Las Campanas Observatory on 2014-11-16. Other spectra were taken by PESSTO (Smartt 2015) as part of transient follow-up using the ESO-NTT/EFOSC2-NTT on 2014-11-25, 2014-12-12, 2014-12-28, 2015-01-26, and 2015-01-27 . These spectra are available at WISEREP (Yaron & Gal-Yam 2012) and in the ESO archive 2 . Both spectra taken in November 2014 showed noticeably broadened Balmer emissionlines compared to the archival spectrum from 1996.
We also obtained seven spectra with the LCOGT FLOYDS spectrograph (Sand 2014 ) at the robotic 2-m Faulkes Telescope South located at Siding Spring Observatory (Brown et al. 2013) . These observations were taken on 2020-04-12, 2020-04-15, 2020-04-24, 2020-04-25, 2020-04-27, 2020-05-16, and 2020-05-17 in order to observe any changes in spectral features prior to and during the most recent outburst. All spectra were reduced following standard reduction procedures using IRAF. All observations were taken with an exposure time of 600 seconds and span 4,300 to 10,000Å.
We used the analysis tool mapspec 3 (MCMC Algorithm for Parameters of Spectra; Fausnaugh et al. 2016) to calibrate the longslit spectra onto the same absolute flux scale using the [OIII]λ 5007 flux to scale the spectra. We assume that the [OIII] narrow-line flux is constant because it originates in an extended region too large to vary on these short time scales. mapspec uses MCMC methods to adjust the flux, wavelength shift, and resolution of each individual spectrum to match that of the reference spectrum. The reference spectrum was defined by an average of the spectra sample.
In addition to the Swift UVOT observations, we also obtained simultaneous Swift X-Ray Telescope (XRT, Burrows et al. 2005) photon-counting observations of ASASSN-14ko. All observations were reprocessed from level one XRT data using the Swift xrtpipeline version 0.13.2, producing cleaned event files and exposure maps. Standard filter and screening criteria 4 were used, as well as the most up-to date calibration files.
To extract background-subtracted count rates, we used a source region with a radius of 50 centered on the optical position of ASASSN-14ko. To define the background, we used a 150. 0 radius source free region centered at (α,δ)=(05 h 25 m 18.08 s , −46 • 00 21.0 ). All count rates are aperture corrected.
To improve the signal to noise, we merged the most recent Swift XRT observations (ObsIDs: 14005, 13979, 13836) using the HEASOFT tool xselect. From this merged observation, we extracted spectra using the task xrtproducts version 0.4.2 and the same extraction regions. Ancillary response files were obtained by merging the individual exposure maps using XIMAGE version 4.5.1 and the task xrtmkarf. We used the readymade response matrix files that are available with the Swift calibration files. This merged Swift spectrum was grouped to have a minimum of 10 counts per energy bin using the FTOOLS command grppha.
On 2015 August 19, ESO 253−G003 was observed using the MOS and PN detectors onboard XMM-Newton (ObsID: 0762920501, PI: Koss) as part of a program to study heavily obscured AGN. Both of the detectors were operated in full-frame mode using a thin filter. We reduced the data using XMM-Newton science system version 15.0.02 and the most up to date calibration files. Periods of high background/proton flares that could affect the quality of the data were identified by generating a count rate histogram of events between 10 and 12 keV.
The observations were only marginally affected by background flares, leading to effective exposures of 25.6 ks and 24.0 ks for the MOS and PN detectors, respectively. For our analysis, we used standard event screening and flags for both the MOS and PN detectors 5 . All files were corrected for vignetting using EVIGWEIGHT. Spectra were extracted from both detectors using the SAS task EVSELECT and the cleaned event files. We used the same source region used to analyze the Swift observa-tions, and the spectra were grouped using grppha to have a minimum of 20 counts per energy bin.
ESO 253−G003 was also observed using the Nuclear Spectroscopic Telescope Array (NuSTAR) on 2015 August 21 (ObsID: 60101014002, PI: Koss) as part of the same program to observe heavily obscured AGN. We reduced the data using the NuSTAR Data Analysis Software (NuSTARDAS) Version 1.8.0 and NuS-TAR CALDB Version 20170817. We performed the standard pipeline data processing with nupipeline, with the saamode=STRICT to identify the South Atlantic Anomaly passages. Using the nuproducts FTOOL, we extracted source spectra from a 100 -radius region and produced ancillary response files and redistribution matrix files for both the A and B modules. We extracted background spectra from annular regions centered on the source, and we followed the procedure outlined by Wik et al. (2014) to estimate the backgrounds and subtract them from the source spectra using the nuskybgd routines 6 .
To analyze the merged Swift, XMM-Newton and NuS-TAR spectra we used the X-ray spectral fitting package (XSPEC) version 12.10.1f and χ 2 statistics. Finally, to further constrain the X-ray emission, we also analysed the available XMM-Newton slew observations of the region. XMM-Newton slew observations take advantage of the fast read out of the PN detector and its ability to observe the sky without reduction of image quality as XMM-Newton maneuvers between pointed observations. Slew observation can detect X-ray emission down to a 0.2-10.0 keV flux limit of ∼ 10 −12 ergs cm −2 s −1 (Saxton et al. 2008) . We found nine slew observations overlapping the source and analyzed them using the SAS tool eslewchain 7 . To extract count rates, we use the same source and background regions used for the pointed XMM-Newton observation. Due to the low exposure times of each observation no spectra could be extracted.
We fit the quiescent host photometry given in Table  2 using AGNfitter (Calistro Rivera et al. 2016) . In addition to our photometry from Section 2, we also include the J, H, and K s fluxes from the 2MASS All-Sky Point Source Catalog (Skrutskie et al. 2006) , the W 1-W 4 fluxes from the AllWISE Source Catalog (Wright et al. 2010) , and the 12, 25, 60, and 100µm fluxes from the IRAS Faint Source Catalog (Moshir & et al. 1990 ).
AGNfitter uses MCMC methods to model the combined contribution of the AGN accretion disk, dusty torus, stellar population, and cold dust. Based on our fit, ESO 253−G003 has a stellar mass of M * = (3.68 +0.06 −0.05 ) × 10 10 M , a stellar population age of 0.59 ± 0.10 Gyr, and a star formation rate of SFR = 24.4 ± 0.3 M yr −1 .
We used several approaches to estimate the SMBH mass. First, we used the 2014-11-25 PESSTO spectrum to estimate a broad-line region radius (R BLR ) using the scaling relation of Kaspi et al. (2005) . We then fit the broad component of Hβ with a Gaussian to estimate the characteristic velocity of the BLR. Assuming a Keplerian orbit and using a geometry factor of 4.3 (Bentz & Katz 2015) , this yields a SMBH mass log 10 (M BH ) = 8.18 +0.05 −0.07 M . Second, we applied the scaling relation between bulge near-infrared (NIR) luminosity and SMBH mass of Marconi & Hunt (2003) . We individually fit each outburst with a fifth-order polynomial to determine the timings of the peaks as shown in red in Figure 1 . We then measured the errors on the peak times by bootstrap resampling the light curves. These errors are shown by the shaded red regions in Figure 1 . The times and fluxes for each peak in the V -and g-band ASAS-SN light curve are given in Table  3 . Figure 1 visually demonstrates that the outburst peaks recur at consistent intervals. Initially, we analyzed the light curves using the box least squares (BLS) periodogram method (Kovács et al. 2002) to obtain a period of 111.82 ± 1.9 days from the first sixteen observed peaks. We used the Half-Width-at-Half-Maximum (HWHM) as an estimate for the period uncertainty (Mighell & Plavchan 2013) . This period esti- mate was used to predict and plan for the 2020-05-18.5 outburst.
Using this mean period leads to significant residuals between the observed and calculated peak times, as shown in the Observed -Calculated (O-C) diagram of Figure 2 . Such a periodic trend is indicative of a period derivative, so for the present analysis we fit the time of the transient peak arrivals as
where t 0 is a reference time, P 0 is the "mean" period,Ṗ is the period derivative, and n is the peak number starting from the first peak set as n=0. We use MCMC estimates of the 1-σ parameter uncertainties. When fitting without aṖ initially, we found the best-fit parameter P 0 = 111.20 ± 0.06 days but the fit returns a reduced χ 2 of 9.75 for 14 degrees of freedom. The fit noticeably improved when including aṖ with t 0 = 3.4 +2.2 −2.3 days, P 0 = 114.6 ± 0.3 days,Ṗ = −0.0021 ± 0.0002 and a reduced χ 2 of 2.07. This fit is shown in Figure 2 as the solid line. If we expand the time of peak uncertainties in quadrature by 0.63 days to have a χ 2 per degree of freedom of unity, we find best-fit parameters t 0 = 5.9 ± 2.7 days, P 0 = 114.2 ± 0.4 days, andṖ = −0.0017 ± 0.0003. This fit is shown by Figure 2 correspond to 2020 September 7.4 ± 1.1 and 2020 December 26.5 ± 1.4. The V -and g-band ASAS-SN light curves for ASASSN-14ko stacked using the phasing of this model are shown in Figure 3 along with the TESS light curve. Including theṖ significantly reduces the scatter between the stacked light curves. We also binned these stacked light curves. As apparent in the right panel of Figure  3 , the peak times are closely aligned across these three filters and the light curve morphologies are similar between peaks and these filters. The outbursts are characterized by a fast rise to peak followed by a shallower decline.
The Catalina Real-Time Transient Survey (CRTS; Drake et al. 2009 ) data contain observations spanning nine years prior to the start of the ASAS-SN V -band light curve. However, the CRTS data quality made it difficult to identify prior outburst peaks relative to the quiescent baseline. The CRTS photometry includes flux from the entire host galaxy since it uses Source Extractor (Bertin & Arnouts 1996) rather than image subtrac-tion. Since the host galaxy is bright and spatially large, it is difficult to recognize prior outbursts. We also found that data from the All-Sky Automated Survey (ASAS; Pojmanski 1997) were not useful for identifying earlier outbursts and we could find no other earlier data.
The TESS light curve of ASASSN-14ko's November 2018 outburst provides a unique view of the outburst morphology. Because TESS observed ASASSN-14ko over three consecutive sectors, the light curve captures the pre-outburst quiescence, the full rise to peak, and the full decline back to quiescence. We first characterized the early-time rise of ASASSN-14ko with a power law model of f = z when t < t 1 , and
(2)
consisting of residual background z, the time of rise t 1 , a flux scale h, and the power law index α. We use Figure 4 . For this fit we inflated the error bars in quadrature by 0.0086 ergs s −1 cm −2 in order to make the reduced χ 2 of the fit unity for 160 degrees of freedom. The high photometric precision of TESS shows that the early-time rise was smooth, and that the time to peak in the TESS filter is 5.60 ± 0.05 days from t 1 when measuring the peak from the data binned at 8-hour intervals. Three TDEs have estimates of α, all of which are steeper. The first TDE detected by TESS was ASASSN-19bt, which had a power law index of α = 2.10 ± 0.12 (Holoien et al. 2019a ). ASASSN-19dj had α = 1.9 ± 0.4 using the ASAS-SN light curve (Hinkle et al. 2020b ). Values of α ∼ 2 are similar to the fireball model used for the early-time evolution of SNe (Riess et al. 1999 , Nugent et al. 2011 . The third example, AT2019qiz, had a still steeper power-law index rise of α = 2.8 ± 0.3 (Nicholl et al. 2020) .
The photometric precision of the TESS light curve also reveals that the decline was also remarkably smooth, as shown in Figure 5 . Since the rate of the decline in partial TDEs is predicted to be steeper than the canonical t −5/3 model, we fit the TESS light curve starting from five days after peak as
with t 0 being the time of disruption, constrained to be before the start of the rise, t 1 , determined above. We again inflated the error bars in quadrature by 0.01 ergs s −1 cm −2 in order to make the reduced χ 2 unity for 2611 degrees of freedom. The best-fit power-law has the parameters t 0 = 2458429.70 ± 0.45 JD, which corresponds to the start of the rise, z = 0.01 ± 0.001 ergs s −1 cm −2 , h = −5.60 ± 0.64 ergs s −1 cm −2 , and α = −1.24 ± 0.03. As seen in Figure 5 , this smooth power law describes the decline well. The exponent, far from being steeper than the canonical α = −5/3 = −1.66, is actually shallower.
We also fit the decline as an exponential decay of the form f = ae −(t−t peak )/τ + c
which returns the best-fit parameters a = 0.24 ± 0.001 ergs s −1 cm −2 , τ = 12.22 ± 0.11 days, and c = 0.04 ± 0.001 days. We set t peak to the peak of the TESS light curve since its value is degenerate with the other parameters. This model has a reduced χ 2 value of 0.77 for the same uncertainties. The exponential decline is, therefore, a better representation of the decline than the power-law decay model. Since starting the fit five days after peak is somewhat arbitrary, Table 4 gives the results for fits starting 5, 10, 15, and 20 days after peak all for the same error model. The fit parameters are relatively stable and the exponential model is always the better fit.
The May 2020 outburst was the first to be predicted in advance, and it occurred as expected. We initially predicted the outburst to peak on MJD 58990.2 ± 1.9 based on our preliminary period estimate using the periodogram. The ASAS-SN g-band light curve actually peaked on MJD 58987.5 +2.5 −0.8 , as shown in the last panel of Figure 1 , consistent with the prediction.
We requested Swift observations to monitor the outburst along with the LCOGT and amateur ground-based data. The non-host-subtracted light curves along with the X-ray hardness ratios are shown in Figure 6 Figure 5 . The declining phase models of the TESS light curve starting 5 days past peak. The best-fit power-law decline is shown in red and the best-fit exponential decline is shown in blue. The TESS data binned in 8 hour intervals are shown by turquoise squares. The bottom panel shows the residuals of the fits color-coded by model. Magenta and cyan squares show the residuals of the power-law decline and exponential decline compared to the binned TESS data, respectively. Table 4 . Best-fit parameters for the power-law and exponential decline models of the TESS light curve starting the fits for different numbers of days after peak. z for the power-law decline was 0.01 ± 0.001 ergs s −1 cm −2 for all start times post-peak. curve due to weather closures. The combined data set was processed as uniformly as possible, but the scatter prevents unambiguous interpretations of the multi-band light curve. The most interesting feature of Figure 6 is that we clearly see a wavelength dependence to the flux peak. While the rise was fully observed for the B-band and longer wavelengths, the U V and U -band data clearly peaked at still earlier times and the optical bands lagged the U V -bands by several days. This is also apparent in the blackbody fits to the host-subtracted Swift data shown in Figure 7 and Table 5 . The peak blackbody luminosity and temperature occurred approximately four days prior to the ASAS-SN g-band peak. However, the blackbody radius remained roughly consistent over time.
The Swift XRT X-ray light curve did not follow the same trend. At peak optical flux, the X-ray flux had dropped by a factor of ∼4 and the spectrum became harder. The hardening of the emission could mean that the effective radius for the X-ray emission is shrinking, that a harder power-law flux component is increasing, or there is additional obscuration that is attenuating the soft X-ray photons. However, during the outburst peak, the hardness ratio decreased, making it difficult to reconcile the X-ray evolution with the UV/optical.
The UV/optical and X-ray SEDs at two epochs during the May 2020 outburst are shown in Figure 8 . The two epochs correspond to four days before and seven days after the ASAS-SN g-band peak, which approximately corresponds to the U V peak and post-maximum U V quiescence. The peak luminosity is roughly 1% of the Eddington luminosity derived in Section 3.
We first compared each of ASASSN-14ko's outbursts to the peak luminosity-decline rate relation for previously-studied TDEs (Hinkle et al. 2020a) . This relation describes a correlation between the peak luminosity and its decline luminosity over 40 days. First, we bolometrically corrected the ASAS-SN V -and g-band, and TESS photometry using a temperature of 28,800 K, which is the median temperature of the most recent outburst based on the Swift data. Then, following the procedure of Hinkle et al. (2020a) , we calculated a peak luminosity and the decline rate over 40 days for each outburst. In Figure 9 , we compare the peak luminosity and decline rate of these outbursts to known TDEs. Even though the power-law slope of the decline is shallower than t −5/3 (Section 4.2), the actual decline rate is steeper than all TDEs of similar luminosity.
We directly compare the TESS light curves of ASASSN-14ko and ASASSN-19bt (Holoien et al. 2019a) in Figure 10 . ASASSN-14ko clearly evolves much more rapidly than ASASSN-19bt, with both a more rapid rise and a more rapid decline. However, the overall morphologies of the declining light curves are very similar if we compress the time relative to peak of ASASSN-19bt by 30% in order to align the light curve declines as shown in the right panel of Figure 10 .
We used optical spectra observed at different points in time and shown in Figure 11 to first classify ASASSN-14ko using standard BPT diagnostics (Baldwin et al. 1981; Veilleux & Osterbrock 1987; Kauffmann et al. 2003; Kewley et al. 2001; Kewley et al. 2006 ). The 1996 spectrum and a weighted average of the 2014-2015 spectra both have line ratios consistent with an AGN. The measured ratios are given in Table 6 .
Next we examined the evolution of the Hβ and Hα emission-line profiles. The spectra from 2014-2015 were taken after ASASSN-14ko's optical peak, as indicated by orange tick marks in Figure 1 . There is a noticeable change in the the Balmer line strengths and profiles during outburst. Since the line changes associated with the transient are at the redshift of ESO 253−G003, we can be certain that these transients are not Galactic in origin. Directly comparing these spectra should be done with caution because the data were taken with different instruments and long-slit setups. However, it is apparent that the emission-line profile shapes clearly changed during the 2014 November outburst.
We obtained five spectra with the LCOGT FLOYDS spectrograph during quiescence in April 2020 and two spectra during the optical rise of the May 2020 outburst. We compared them to determine any change in the emission-line profiles as shown in Figure 12 . Due to observatory airmass constraints and ASASSN-14ko's low position near the horizon from Siding Springs Observatory at that time, observations were not possible past 2020-05-17. These two observations occur one and two days before the optical peak, respectively, which coincides with the start of the decline in the UV, as shown by the vertical dashed lines in Figure 6 and the orange tick marks in Figure 1 . A noticeable feature of the spectra during outburst is a blue wing around Hβ near 4830Å that is not present in the spectra taken during quiescence. This feature is similar to the broadened wings present in the 2014-11-16 and 2014-25 spectra which were taken during an outburst optical decline. This gives further evidence that the Balmer lines change during the outburst. . Non-host subtracted photometry of the May 2020 outburst. ASASSN-14ko g-band photometry is represented as squares, Swift data as circles, LCOGT data as pentagons, and data from amateur observatories as diamonds. B-band, g-band, and V -band data taken from different telescopes on the same day were averaged. Swift B-and V -band data were converted to Johnson B and V magnitudes before being converted to flux to enable direct comparison with the ground-based data. Error bars are plotted but are frequently smaller than the size of the points. The red shaded region denotes the time of the ASAS-SN g-band peak on JD 2458987.5 +2.4 −0.8 and the magenta vertical solid line shows the predicted peak for our model with a period derivative, as was also shown in Figure 1 Table 5 . The time is relative to the g-band peak on MJD 58987.5. The luminosity and temperature peak occurred several days prior to the g-band peak. Table 5 . Blackbody luminosity, radius, and temperature during the May 2020 outburst derived from the hostsubtracted Swift data. The similarities between the spectroscopic evolution of the two outbursts separated by nearly six years indicates that the spectroscopic change during the flare may be a consistent component of these events. The photometric and spectroscopic evolution appear closely tied, TDEs ASAS-SN V ASAS-SN g TESS Figure 9 . The peak luminosity-decline rate for ASASSN-14ko's outbursts compared to twenty-one previously studied TDEs. Individual ASAS-SN V (green circles), ASAS-SN g (teal circles), and TESS (red circle) epochs are shown along with the TDEs analyzed by Hinkle et al. (2020a) represented as gray squares. The black solid line is the best-fit line for the TDEs and the dashed black lines are the allowed range of uncertainty from the best fit line. (Baldwin et al. 1981 , Veilleux & Osterbrock 1987 , Kewley et al. 2001 , Kauffmann et al. 2003 , Kewley et al. 2006 . and higher temporal resolution of upcoming outbursts may reveal further connections between the rise and decline of the photometric light curves with morphological changes in the Balmer lines.
We compared the two X-ray spectral epochs to characterize the X-ray emission evolution and are shown in Figure 13 . We modeled the 2015 XMM-Newton+NuSTAR spectra as a combination of a soft black body plus a power law. The black body model had a temperature of 0.15 ± 0.01 keV and the power law model had a Γ = 0.87 ± 0.04. There is a strong 6.4 keV Fe line present that we include in the model as a Gaussian. The Swift XRT spectrum was extracted from the merged data taken from ∼30 days prior to ∼9 days after the May 2020 flare, which corresponds to times during both quiescence and outburst. This spectrum was best fit by a power law with Γ = 1.09 +0.20 −0.22 and the fit was not improved by adding an additional black body with a best-fit temperature of 0.13 ± 0.03 keV. The Swift spectrum had combined all the May 2020 observations because the signal-to-noise ratio of the spectrum would otherwise have been too low to assess its characteristics. With only ∼650 counts, poor statistics likely prevent any detection of the Fe line. The fit parameters are summarized in Table 7 .
The XMM-Newton spectrum was taken during a period when the optical flare was quiescent, however the flux derived from the XMM-Newton spectrum is similar to what we observed with Swift. The similarities between the Swift spectrum and the XMM-Newton spectrum may indicate that even during times around the UV peak, the average properties of the X-ray emission during May 2020 may not have changed significantly from the properties of the deep XMM-Newton observation. Future observations with better signal-to-noise will be essential to disentangle potential X-ray emission evolution between quiescence and outburst. There are nine epochs of XMM-Newton slew observations of this galaxy observed between 2007-2019. However, given the sparse sampling and low signal to noise, there is no obvious pattern after phase folding these data.
Here we discuss several scenarios for ASASSN-14ko's periodic outbursts and some of their problems. We consider the possibilities that ASASSN-14ko is a sub-parsec SMBH binary system, a SMBH and perturbing massive star binary system, and, finally, a repeating partial TDE. We do not consider stellar origins for the outbursts. As noted earlier, a Galactic source is ruled out because the emission line changes associated with the outbursts occur at the redshift of the host. At the redshift of the host, the peak luminosities are similar Figure 12 . Spectra of ESO 253−G003 observed with the FLOYDS spectrograph at LCOGT during quiescence in April 2020 and during the May 2020 outburst. The first five spectra were obtained during quiescence and the last two spectra were obtained during the optical rise of May 2020 outburst. The spectra have been scaled using mapspec to put all of them onto a common flux scale. to those of luminous supernovae and so are too high for non-explosive stellar transients. The many quasiperiodic repetitions of similar luminosity rule out explosive possibilities. This appears to leave only phenomena associated with a central SMBH as possible explanations.
As discussed in the introduction, apparently periodic emission from AGN is frequently interpreted as evidence for an SMBH binary. However, there are three main inconsistencies with this scenario for ASASSN-14ko: (1) the period derivative; (2) the short life time of such a binary and (3) the lack of velocity shifts. We will scale the results for black hole masses of M BHa = 5 × 10 7M BHa M , M BHb = 5 × 10 7M BHb M and a period of P = 114P days.
Due to the emission of gravitational waves, the binary will have a period derivative oḟ
where f (e) = (1 + 73e 2 /24 + 37e 4 /96)(1 − e 2 ) −7/2 is the dependence on the orbital ellipticity e, and a time to merger of
where g(e) = (1 − e 2 ) 7/2 (Peters 1964 ). The first problem is that the observed period derivative is over an order of magnitude larger than the scaling in Eqn. 6 for black holes in a circular orbit. This can be solved only by substantially increasing the mass of at least one of the black holes, or by making the orbit significantly elliptical. Both of these solutions to the mismatch in period derivatives will exacerbate the next two problems by reducing the binary lifetime and increasing the binary velocities. If the solution is to change the masses, the lines must be formed around the less massive SMBH, which also exacerbates the velocity problem.
The second problem is that finding a binary so close to merging is very improbable. The first way to phrase this is through what it requires for the properties of SMBH binaries in other galaxies. The number density of L * galaxies is roughly n 0.01h 3 Mpc −3 (e.g., Kochanek et al. 2001) , so there are roughly N g 3 × 10 5 such galaxies inside the distance to ESO 253−G003. To have one system merging in the next t 0 = 1000 years implies that all of these galaxies must contain binaries that will merge in the next N g t 0 300 million years. This in turn implies that they must all have orbital periods shorter than ∼ N 3/8 g P 0 35 years, which would seem to make the problem of finding binary SMBH systems rather trivial rather than being as difficult as it appears to be in practice. The second way to phrase the problem is that if there is one SMBH binary with a period < 1/3 year in this volume, there are 3 8/3 = 19 with periods < 1 year, 350 with periods < 3 years, and 8700 with periods < 10 years.
The third problem for an SMBH binary is the high velocity scale, with for circular orbits. While our phase sampling is poor, we arguably can set a limit on any Hβ line shifts of < 50Å, which corresponds to a velocity shift limit of < 3000 km/s. If our estimate ofṖ is correct, this problem cannot be solved by reducing the masses or invoking a large mass ratio with the emission lines being formed by material associated with the more massive SMBH. Similarly, raising the ellipticity to solve theṖ problem makes this problem worse because of the higher velocities at pericenter compared to a circular orbit of the same period. The velocity problem can only be solved by making the system nearly face on or by relying on the poor spectral sampling to hide the velocity shifts.
Rather than a binary system with two SMBHs, the periodic outbursts could be driven by a star orbiting a single SMBH. Because stars are far more common than SMBHs and the gravitational wave merger lifetimes are now far longer, the probability argument against an SMBH binary is removed. The long gravitational wave merger time does mean that the observedṖ must have a different origin such as viscous interactions between the disk and the star, although estimating this effect is non-trivial.
We now assume a star of mass M * =M * M and radius R * =R * R orbiting a black hole of mass M BH = 5 × 10 7M BH M . There are three length scales of immediate interest, and we show their relative values and their dependence on SMBH mass in Figure 14 . The Schwarzschild radius of the black hole is
the tidal (Roche) radius to disrupt the star is
and the orbital semi-major axis is a = 2.5 × 10 15P 2/3M 1/3 BH cm.
Finally, for a pericentric radius R p , the ellipticity of the orbit is Figure 14 . The Schwarzschild radius (Rs, solid), semimajor axis (a, dotted) and several tidal radii (RT , dashed) as a function of SMBH mass MBH . The tidal radii are for a 0.3M main-sequence star (RT (0.3M )), the Sun (RT (1.0M )), and a star with the mass of the Sun and a radius of 10R (RT (M , 10R ). On these logarithmic scales, the effects of spin on the BH horizon and the tidal limits are modest.
For some star/SMBH scenarios, the true period would be twice the observed period, which would increase the semimajor axis by 0.2 dex. From Figure 14 , we see that there is no problem having a main sequence star orbiting an SMBH in this mass range without serious tidal effects provided the orbital ellipticity is moderate. Particularly over many orbits, there would be some relativistic effects. For example, the orbital precession per orbit iṡ
for a Schwarzschild black hole. This means that we would expect systematic changes with time independent of the mechanism driving the flares. Figure 15 . The stacked ASAS-SN g-band light curves separated by even (blue) and odd (red) outbursts as a function of phase. The light curves are offset for clarity and the data of each phase are given a different color shade. The binned even light curve, shown by star symbols, is superimposed over the odd light curves, and vice versa, with the binned odd light curve shown by diamond symbols.
If the star is not being tidally disrupted, then the flares must be driven by periodic disturbances of the accretion flows as the star passes through the accretion disk. However, a star simply embedded in the disk would represent a continuous perturbation that is unlikely to drive periodic flares. Flaring would seem to require a stellar orbit at a significant inclination angle relative to the disk. The star would then make two passages through the disk per orbit, so the orbital period would likely be twice the flare period.
Each orbit would produce a pair of flares with the spacing dependant on the orbital eccentricity and the argument of periapsis (ω) relative to the accretion disk. The separation between pairs will increase as the orbit becomes more eccentric. While eccentric orbits with ω near 0 or π can have equally spaced encounters, the star will encounter the disk at different radii/temperatures. Given that the flare spacing, profiles, and amplitudes are all essentially constant between cycles, a perturbing star would seem be required to be on an inclined but nearly circular orbit. Even then, the similarity of the flares seems odd because one would expect encounters with the star moving away from the observer and into the disk to differ from the reverse, except for nearly edge on viewing angles. In Figure 15 we see that the even and odd flare profiles are very similar in amplitude, shape, and duration which is difficult to reconcile with this model.
Finally, while an inclined, circular orbit of a star orbiting a SMBH might be able to perturb the accretion disk with the right frequency, there is no obvious time scale in disks to then make the flares so short in duration.
A third possibility to explain ASASSN-14ko's periodic outbursts is as a repeating TDE that is partially disrupted after each passage close to the central SMBH. While it requires fine tuning to have a main sequence star pass close to its tidal limit but remain outside the SMBH horizon, we can see from Figure 14 that it is relatively easy for an evolved star on an elliptical orbit to do so. As discussed in the introduction, giants are also the most likely candidates for partial disruptions that could power periodic flares. MacLeod et al. (2012) and find that the star will begin to lose mass once β = R T /R p < 0.5.
The peak luminosities of the flares, L p 5 × 10 44 erg/s correspond to peak accretion rates ofṀ p 0.1 −1 0.1 M /year where = 0.1 0.1 is the accretion efficiency. If the peaks last ∼ 10 days, the accreted mass of ∆M 0.003 −1 0.1 M is certainly low enough to allow repeated outbursts on this scale. Note, however, that these accretion rates are significantly higher than envisioned by MacLeod et al. (2013) and the time scales are much shorter. Ryu et al. (2020) found that the change in the orbital specific energy of the star in a partial disruption is f ∼ 10 −3 of the specific energy scale GM BH R * /R 2 t of the stripped debris. They find both positive and negative energy changes, so there is no prediction of the sign of the changes. The period derivative measured in Section 4.1 implies a change in the orbital specific energy of GM BHṖ /3a. This means that we should expect a period derivative of
For f ∼ 10 −3 , this implies |Ṗ | ∼ 10 −3 with relatively little sensitivity to the exact values of the parameters and remarkably close to the measured period derivative. The agreement is perhaps more remarkable because the Ryu et al. (2020) simulations were for a single pericentric passage of a main sequence star on an initially parabolic orbit, rather different from the orbit required here. However, the orbit of the puffy stellar merger remnant in Figure 13 (top) of Antonini et al. (2011) has a semi-major axis shrinking as ∆a/a ∼ 10 −3 per orbit, which is the same order of magnitude. The example in the lower panel of this figure shows very little orbital evolution but also shows very little ongoing mass loss. Note that the orbital changes essentially occur with the pericenter fixed because the tidal interactions are only important at pericenter. Because the structure of the star must be changing with each pericentric passage due to the mass loss, torques and heating,Ṗ presumably cannot be constant on longer time scales. The partial TDE hypothesis also seems better able to explain the similarity of the flares since each pericentric passage is almost identical in geometry to the previous and the required mass loss rates appear to be modest. However, they cannot truly be identical since the orbital geometry must slowly change due to precession (relativistic and tidal) and the mass lost over tens of encounters ceases to be modest. Overall, the repeating, partial TDE interpretation seems most consistent with the available observations. 7. SUMMARY Although ASASSN-14ko was first thought to be a supernova, the subsequent six years of ASAS-SN V -and gband data show that the flares occur at regular intervals. The 17 flares observed to date are well modeled using a period of P 0 = 114.2 ± 0.4 days and period derivative oḟ P = −0.0017 ± 0.0003. Adopting this model, the next two flares will peak in the optical on UT 2020-09-7.4±1.1 and UT 2020-12-26.5 ± 1.4.
In addition to ASAS-SN, ATLAS, and Swift multiwavelength photometric data, TESS observed ASASSN-14ko during its November 2018 outburst. The TESS light curve has a decline rate that is dissimilar to previously studied TDEs and a rapid rise to peak occurring over 5.60 ± 0.05 days. The TESS data also show that the rise and decline were smooth and lack shorttimescale variability. The individual outbursts are morphologically very similar over the six year baseline of observations. While a host of problems interfered with studying the May 2020 outburst well, there was clear evidence that the outburst peaks a few days earlier in the UV than in the optical. Spectra taken during and prior to the May 2020 outburst revealed morphological changes around Hβ during the flare which was similar to what occurred during the 2014 outburst. This suggests that morphological changes in the emission lines are consistently associated with the optical outbursts over time. We examined several possible scenarios to explain the cause of this AGN's unusual behavior, including the presence of a SMBH binary, a SMBH plus a perturbing massive star, or a repeating partial TDE. Between these scenarios, we favor a repeating partial TDE. We believe that any stellar transients or explosions whether Galactic or in the host are ruled out. The most important next step is to time and study the flares more closely across the electromagnetic spectrum. The relatively short period and system brightness make this relatively easy. ASASSN-14ko will be observed by TESS again in Sectors 31-33 during the predicted December 2020 outburst. These observations will give further constraints on the nature of these outbursts, and presents a unique opportunity to do a detailed reverberation mapping analysis of the system.
Software: ftools (Blackburn 1995) , NuSTAR Data Analysis Software (v1.8.0)