Unusual Hard X-ray Flares Caught in NICER Monitoring of the Binary Supermassive Black Hole Candidate AT2019cuk/Tick Tock/SDSS J1430+2303

TL;DR

This study reports unique day-long hard X-ray flares in the changing-look AGN AT2019cuk, observed through high-cadence NICER monitoring, suggesting coronal variability as the likely cause and highlighting the importance of rapid X-ray observations.

Contribution

First detection of day-long hard X-ray flares in a changing-look AGN, providing new insights into AGN corona dynamics and variability mechanisms.

Findings

Repetitive hard X-ray flares with harder spectra observed

No evidence of periodic modulation in X-ray, UV, or optical bands

Variable corona/jet activity likely causes the flares

Abstract

The nuclear transient AT2019cuk/Tick Tock/SDSS J1430+2303 has been suggested to harbor a supermassive black hole (SMBH) binary near coalescence. We report results from high-cadence NICER X-ray monitoring with multiple visits per day from January-August 2022, as well as continued optical monitoring during the same time period. We find no evidence of periodic/quasi-periodic modulation in the X-ray, UV, or optical bands, however we do observe exotic hard X-ray variability that is unusual for a typical AGN. The most striking feature of the NICER light curve is repetitive hard (2-4 keV) X-ray flares that result in distinctly harder X-ray spectra compared to the non-flaring data. In its non-flaring state, AT2019cuk looks like a relatively standard AGN, but it presents the first case of day-long, hard X-ray flares in a changing-look AGN. We consider a few different models for the driving mechanism of these hard X-ray flares, including: (1) corona/jet variability driven by increased magnetic activity, (2) variable obscuration, and (3) self-lensing from the potential secondary SMBH. We prefer the variable corona model, as the obscuration model requires rather contrived timescales and the self-lensing model is difficult to reconcile with a lack of clear periodicity in the flares. These findings illustrate how important high-cadence X-ray monitoring is to our understanding of the rapid variability of the X-ray corona and necessitate further high-cadence, multi-wavelength monitoring of changing-look AGN like AT2019cuk to probe the corona-jet connection.

Tables

Table 1. Table 1: XMM-Newton/NuSTAR Spectral Fits

Model	Parameter (unit)	Value
zcut	$\Gamma$	${1.65}_{-0.01}^{+0.01}$
	$E_{\mathrm{cut}}$ (keV)	$>173$
zcut+zbb	$N_{H,\mathrm{neut}}$ (${10}^{22}$ cm-2)	${0.07}_{-0.01}^{+0.01}$
	$\Gamma$	${1.59}_{-0.03}^{+0.03}$
	$E_{\mathrm{cut}}$ (keV)	${75}_{-18}^{+33}$
	$k{T}_{\mathrm{bb}}$ (eV)	${81}_{-3}^{+3}$
(partcov*xstar)*zcut	$N_{H,\mathrm{ion}}$ (${10}^{22}$ cm-2)	${1.1}_{-0.2}^{+1.4}$
	$\log \left(\xi /\mathrm{erg}\mathrm{cm}{s}^{-1}\right)$	${1.91}_{-0.14}^{+0.11}$
	Covering Fraction	${0.88}_{-0.36}^{+0.12}$
	$\Gamma$	${1.73}_{-0.01}^{+0.01}$
	$E_{\mathrm{cut}}$ (keV)	$>281$

Full text

Active galactic nuclei (16)–High energy astrophysics (739)–Supermassive black holes (1663)–X-ray transient sources (1852)

††facilities: NICER, XMM-Newton, NuSTAR, Swift, ZTF, LOAO, MAO, CBNUO††software: XSPEC (v12.12.0; Arnaud 1996) HEASoft (v6.30.1; NASA HEASARC 2014) XMM-SAS (v20.0.0; Gabriel et al. 2004) HOTPANTS (Becker, 2015)

1 Introduction

In the theory of hierarchical structure growth, supermassive black hole (SMBH) binaries are thought to be a natural result of merging galaxies. Indeed, kiloparsec-separation dual SMBH systems have been identified in merging systems across the electromagnetic spectrum, including in the radio (e.g. Fu et al., 2011; Müller-Sánchez et al., 2015), optical (e.g. Liu et al., 2018; Silverman et al., 2020), and X-ray (e.g. Komossa et al., 2003; Koss et al., 2011; Foord et al., 2020). Some close-separation, parsec-scale binaries have also been resolved using very long baseline interferometry (VLBI; e.g. Rodriguez et al., 2006; Burke-Spolaor, 2011). However, imaging SMBH binaries at sub-parsec separations is difficult due to the extremely high angular resolution required, which is only possible with observatories like the Event Horizon Telescope or its upgrades (e.g. D’Orazio & Loeb, 2018; Kurczynski et al., 2022). Instead, spectroscopic signatures, like double-peaked or systematically shifted emission lines around the accreting black hole (e.g. Tsalmantza et al., 2011; Ju et al., 2013), or time-dependent photometric signatures, like periodic modulation of the light curve (e.g. Graham et al., 2015a, b; Liu et al., 2015; Charisi et al., 2016; Chen et al., 2022), have been argued to indicate the presence of a close-separation binary SMBH. While spectroscopic searches for binary candidates have been fruitful, this behavior is not unique to binaries, with other explanations including a recoiling SMBH (e.g. Blecha & Loeb, 2008; Komossa et al., 2008) or Keplerian motion in the disk (e.g. Eracleous & Halpern, 1994; Chornock et al., 2010; Gaskell, 2010). Likewise, finding periodic signals among stochastic AGN variability requires many periods of data (years to decades) to confirm a close-separation binary system (see e.g. Vaughan & Uttley, 2005; Vaughan et al., 2016). To date, our best close-separation binary SMBH candidates come from periodic light curves with extremely long temporal baselines (e.g. OJ 287; Sillanpaa et al., 1988; Lehto & Valtonen, 1996; Valtonen et al., 2008).

AT2019cuk (also known as Tick Tock, SDSS J1430, or ZTF18aarippg) is a nuclear transient discovered by the Zwicky Transient Facility (ZTF) on 2019 March 30 in the host galaxy SDSS J143016.05+230244.4 ($z=0.081$; Oh et al., 2015). It rose by approximately one magnitude in the optical $g$-band, reaching a peak magnitude of 17 and then declining with an apparently oscillatory pattern, reminiscent of what one would expect from a binary SMBH (see Figure 1 of this work and Figure 1 from Jiang et al., 2022). After including Swift X-ray and UV observations from late 2021, Jiang et al. (2022) found that the period of the oscillations rapidly decreased from $\sim$1 year to $\sim$1 month over the course of three years. The authors modeled this decreasing periodicity with a system of two merging SMBHs and suggested that they could merge within three years.

AT2019cuk may present a unique opportunity to potentially witness the coalescence of a binary SMBH for the first time. Hence, there has been extensive multi-wavelength follow-up of AT2019cuk, including new optical (e.g. Moiseev et al., 2022; Dotti et al., 2023), radio (e.g. An et al., 2022a, b; Bruni et al., 2022), and X-ray observations (e.g. Pasham et al., 2022; Dou et al., 2022). Continued optical observations have revealed an evolving spectrum, with the potential appearance of a He I emission line at 6678 Å  (Moiseev et al., 2022). Dotti et al. (2023) also present further optical photometry of AT2019cuk and argue against a binary SMBH model based on the decreasing amplitude of the oscillations. Instead, they suggest that the periodicity we see could be the result of Lens-Thirring precession of the accretion disk around only one SMBH. In the radio, milliarcsecond-resolution very long baseline interferometry revealed a compact ($<0.8$ pc) radio core with a flat radio spectrum producing $\sim$40% of the radio emission in AT2019cuk, which was suggested to be the result of either an optically-thick jet or corona (An et al., 2022b). However, so far, no radio outburst has been reported, as all reported flux levels are consistent with previous survey observations. Finally, X-ray observations offer unique insight into the inner accretion flow in this source, as, if AT2019cuk is indeed a binary SMBH, then the two putative SMBHs should be too close to host their own broad line regions, but could have their own mini accretion disks (e.g. Shen & Loeb, 2010; d’Ascoli et al., 2018). Extensive X-ray follow-up has been reported by Pasham et al. (2022) and Dou et al. (2022), which both present X-ray spectral analysis of initial NICER, XMM-Newton, NuSTAR, and Chandra observations of AT2019cuk, finding a potentially variable warm absorber and evidence for relativistic reflection.

In this work, we present high-cadence NICER X-ray observations and further optical monitoring of AT2019cuk. When Jiang et al. (2022) first announced AT2019cuk as a potential SMBH binary, we triggered a NICER guest observer program (ID: 4078, PI: Pasham) to perform high-cadence (two visits per day for several months) X-ray monitoring of AT2019cuk starting on January 20, 2022. This monitoring continued until the source reached solar exclusion in late August 2022. In this work, we present results up until September 8, 2022, but note that much of the NICER data is affected by the small sun angle beyond mid August 2022. We also carried out daily optical monitoring using multiple telescopes across the globe. In Section 2, we describe our data and their respective reduction methods. We present the spectral and timing studies in Section 3 and discuss the implications of our findings on the nature of the X-ray corona and with regards to the potential binary SMBH scenario in Section 4. Finally, we summarize our findings in Section 5. Throughout this paper, we assume a standard flat $\Lambda$CDM cosmology with $\Omega_{\Lambda}=0.73$ and $H_{0}$ = 70 km s*-1* Mpc*-1*. All uncertainties are quoted at the 90% confidence level, unless otherwise noted.

2 Observations and Data Reduction

2.1 NICER

NICER (Gendreau et al., 2016) has performed extensive high-cadence monitoring of AT2019cuk, with roughly 1-2 visits per day from January-August 2022. We reduced the NICER data using NICERDAS (version 9, as part of HEASoft version 6.30.1). We first ran nicerl2 to reprocess the data with standard NICERDAS tools and apply the default settings, with the following exceptions. We did not perform any undershoot or overshoot filtering (i.e. we use underonly_range=-, overonly_range=-, overonly_expr=NONE), and instead filtered on the background-subtracted products, using the 3C50 background detailed in Remillard et al. (2022). We utilized level 3 filtering, which implements a 0.2-0.3 keV background-subtracted count rate threshold of 2 cts s*-1* and a 13-15 keV background-subtracted count rate threshold of 0.05 cts s*-1*. Any good time intervals (GTIs) which exceed these values were removed from our analysis. Such filtering is designed to limit contamination of the energy range 0.3-12 keV, intended for investigations of the science target. In addition, we also filtered on the cut-off rigidity (COR) parameter in NICER data by using cor_range=1.5-* when running nicerl2 to ensure that any flares from either electron precipitation events in the geomagnetosphere or cosmic rays are removed from our data. We excluded data from detectors 14 and 34, which are known to be excessively noisy, as well as data from any “hot detectors,” defined as a detector whose 0-0.2 keV raw count rate exceeds the median 0-0.2 keV raw count rate for all detectors by $4\sigma$ using an iterative sigma-clipping procedure. Lastly, we looked at the distribution of the 15-18 keV raw count rates, which should all be background events, and removed any GTI which exceeds the median rate for all GTIs by $4\sigma$ using an iterative sigma-clipping procedure.

After all filtering is complete, we extracted a spectrum and estimated the background spectrum for each GTI using the nibackgen3C50 tool (Remillard et al., 2022). We utilized the nicerarf and nicerrmf tools to produce response files, with the appropriate weighting for any detectors we removed in filtering. We then grouped each spectrum using the optimal binning procedure (Kaastra & Bleeker, 2016) and to a minimum of 25 counts per bin to employ $\chi^{2}$ statistics when fitting. After this grouping, any spectrum with five or fewer channels was discarded as this leaves too few degrees of freedom to perform a fit to a spectral model including both a power law and a blackbody. We included data from January 20, 2022 when the NICER campaign began, up until September 08, 2022 at which point the source was in solar exclusion (ObsIDs 4202540102-108, 4578010101-299). However, much of the data from mid August onward was negatively affected by the small Sun angle. Hence, there is very little reliable data after 3C50 level 3 filtering in ObsIDs greater than 4578010280, although we include any GTI which passes all of the above described cuts. In total, this filtering resulted in 576 GTIs in our analysis for a total exposure time of 323.9 ks.

As Dou et al. (2022) first noted, there is a nearby AGN that is roughly 3.1 arcmin away from AT2019cuk, which is just within the field of view of NICER. Therefore, in Appendix B, we perform extensive checks to assess the impact of this nearby source. We utilize substantial contemporaneous observations with Swift to assess the impact of the nearby AGN on the NICER light curve, as Swift can separately image each of these sources. No other X-ray sources of comparable flux were detected near AT2019cuk in around 160 ks of stacked Swift data. The NICER light curve closely resembles the Swift light curve for AT2019cuk (see Figure 1 below and Figure 13 in the Appendix), and Swift also sees occasional hard X-ray spectra (see Figure 9 in the Appendix). Additionally, the nearby AGN’s contribution to the NICER measurements is further diluted by a factor of roughly 2-3 due to the reduced effective area of the NICER optics at the offset position (see Figure 11 in the Appendix). Thus, we conclude that the interloper in the NICER field of view contributes negligibly to NICER light curve and spectral analysis.

2.2 XMM-Newton

XMM-Newton (Jansen et al., 2001) performed two DDT observations of AT2019cuk in December 2021 and January 2022 (ObsIDs 0893810201 and 0893810401), which we use for comparison with the NICER data. We reduced the XMM-Newton EPIC-pn data from these two observations using the XMM-Newton Science Analysis System (XMM-SAS; v20.0.0). We followed standard data reduction processes, including running epproc to produce calibrated event files, screening for periods of high background flaring by removing time intervals with a count rate of $>$0.4 cts s*-1* in the 10-12 keV range, and producing response files using rmfgen and arfgen. All single and double events (PATTERN $\leq$ 4) were used when extracting spectra. Source spectra and light curves were extracted in a circular region with a radius of 35”, and background products were extracted from an off-source region on the same CCD chip with a radius of 40”. Spectra were grouped to a minimum of 25 counts per bin to employ $\chi^{2}$ statistics when fitting.

2.3 NuSTAR

In February 2022, following the initial announcement of AT2019cuk as a potential SMBH binary, we requested a 100 ks NuSTAR DDT observation (Harrison et al., 2013, ObsID 90801605002). All data was processed with NuSTARDAS (v2.1.2, as part of HEASoft v6.30.1) and with calibration files from NuSTAR CALDB v20220802. We reduced the NuSTAR data using nupipeline to produce calibrated event files and nuproducts to produce spectra and light curves for the FPMA and FPMB modules individually. Source products were extracted in a circular aperture with a radius of 60” and background products were extracted in a nearby off-source region with a radius of 80”. Spectra were grouped to a minimum of 25 counts per bin to employ $\chi^{2}$ statistics when fitting.

2.4 Swift

AT2019cuk has also been observed extensively with Swift (Gehrels et al., 2004), with a cadence of roughly one visit every two days from January-August 2022. We utilized the Swift XRT data to provide consistency checks for our NICER analysis and confirm that the NICER light curve is dominated by AT2019cuk, rather than a nearby AGN (see Appendix Section B). Using the online Swift XRT products generator111https://www.swift.ac.uk/user_objects/ (Evans et al., 2007, 2009), we created light curves on a per-GTI basis in the same energy bands as NICER (Total: 0.3-4 keV, Soft: 0.3-2 keV, Hard: 2-4 keV). We also extracted spectra on a per-ObsID basis and binned these spectra with the optimal binning procedure (Kaastra & Bleeker, 2016) and to a minimum of 1 count per bin for the use of $C$ statistic minimization when fitting. This binning was used due to significantly lower count rates from Swift compared to NICER.

During this monitoring, the Swift Ultra-Violet/Optical Telescope (UVOT; Roming et al. 2005) instrument performed observations with UVW1 as the primary filter, with only a few observations with all of the other Swift UVOT filters taken at the beginning of the observing period. The data was processed with the uvotsource command. A circular region of 5″, centered at the target position, was chosen for the source and another region of 40″ located at a nearby position free of point sources was used to estimate the background emission. UVW1 magnitudes were corrected for Galactic extinction using $E(B-V)=0.025$ (Schlafly & Finkbeiner, 2011).

2.5 Ground-Based Optical Photometry

In addition to X-ray and UV monitoring, ZTF has continued to monitor AT2019cuk. We performed point spread function (PSF) photometry on all publicly available ZTF data at the location of AT2019cuk using the ZTF forced-photometry service (Masci et al., 2019) in $g$- and $r$-band. Similar to UVOT, ZTF photometry was corrected for Galactic extinction.

We also collected more optical monitoring from observatories across the globe, including Maidanak Observatory (MAO, 1.5m, Ehgamberdiev 2018), Lemmonsan Optical Astronomy Observatory (LOAO, 1.0m), and Chungbuk National University Observatory (CBNUO, 0.6m), all of which are a part of the SomangNet network of telescopes (Im et al., 2021). Both $B$ and $R$ images from the SomangNet telescopes were obtained from February 23, 2022 to July 23, 2022 with roughly daily cadence. The data were reduced using our automatic image analysis pipeline, gpPy (G. S.-H. Paek et al., in preparation), which process includes bias and dark subtractions, flat-fielding, astrometry, and photometry calibration using point sources in the vicinity of AT2019cuk. We performed aperture photometry with a single aperture size to minimize the host galaxy contamination after convolving images from different epochs to have a common point-spread-function (PSF) using the HOTPANTS software (Becker 2015; see e.g. Kim et al. 2018, 2019). Image convolution is necessary due to different seeing values on different nights, which may sample different amounts of the host galaxy light (when done with a variable aperture size that is proportional to the PSF size) or the nuclear transient (when done with a fixed aperture size; see e.g. Kim et al. 2018). This process was done independently for each filter and each observatory and means that the optical light curve is contaminated by the host galaxy at different amounts. However, this contamination is negligible since the nuclear component dominates near the center where the aperture photometry was done. As with ZTF and Swift UVOT data, the optical data from the SomangNet network was corrected for Galactic extinction.

3 Analysis and Results

Figure 1 shows the long-term ZTF light curve spanning from 2019-2022, with a zoom-in on the observing period covered in this paper. The data shown include NICER and Swift XRT in the X-ray (top panel of the bottom zoom-in), Swift UVW1 in the UV (second panel of the bottom zoom-in), and ZTF, CBNUO, MAO, and LOAO in the optical (bottom two panels of the bottom zoom-in).

3.1 Identifying & Validating Hard X-ray Flares

One of the most striking features of the AT2019cuk NICER light curve is short, hard X-ray flares, seen most prominently in the 2-4 keV band. In Figure 2, we show the NICER light curves for AT2019cuk in three different bands (total: 0.3-4 keV, soft: 0.3-2 keV, hard: 2-4 keV), as well as a hardness ratio, defined as the 2-4 keV background-subtracted count rate divided by the 0.3-2 keV background-subtracted count rate. While the flares are evident by eye in the 2-4 keV light curve and hardness ratio, the GTIs are relatively short and at a low count rate. Thus, great care must be taken to ensure that the flares are indeed real, rather than an artifact of poor counting statistics or contamination from the NICER instrumental effects.

To assess the impact of counting statistics, we used the best fitting XMM-Newton/NuSTAR spectral model to simulate a NICER data set with the same distribution of exposure times and fluxes as the real data. The simulated data show essentially no hard X-ray flares and the real data show a tail of GTIs at significantly higher hardness ratio, indicating that the flares are not the result of poor counting statistics (see Figure 10 in the Appendix). Further details of these simulations are given in Section A.3 of the Appendix. Additionally, as the NICER instrument is most sensitive at soft energies ($\lesssim 2$ keV) and can have a significant background contribution, we are careful to ensure that these flares are astrophysical. We find no evidence for any instrumental or background effect that could lead to such flares, and we detail these extensive checks in Section A.1 of the Appendix. Similarly, we find consistency between our NICER findings and the contemporaneous Swift observations, including a few similarly hard Swift X-ray spectra (see Section A.2 of the Appendix). Thus, we find that neither counting statistics nor instrumental effects can account for the hard X-ray flares in AT2019cuk.

Identifying flares is best done using the hardness ratio, as it normalizes the long term variability using the soft band. We again utilize the simulated NICER data set to help identify flares, as it provides a less biased distribution with which to distinguish outliers. Therefore, we define a flare to be any GTI with a hardness ratio greater than $3\sigma$ above the mean hardness ratio of the simulated NICER data (which corresponds to a hardness ratio cut of 0.24). These flares are shown as orange stars in Figure 2 and throughout the rest of this paper. As evidenced by the third panel of Figure 2, this hardness ratio cut picks out the hard flares in the 2-4 keV band quite well. Additionally, we show a hardness-intensity diagram in Figure 3, which highlights the unique variability properties of the flares in AT2019cuk. For comparison, Figure 3 shows a blue arrow highlighting the expected coronal variability for AGN which follow the standard softer-when-brighter behavior (e.g. Sobolewska & Papadakis, 2009). The hard X-ray flares show variability that is opposite of what is usually seen in AGN.

Given the hard X-ray nature of these flares, we searched the Swift BAT 105-month catalog (Oh et al., 2018) and the Swift BAT transient monitoring catalog (Krimm et al., 2013) for AT2019cuk, but the source was not detected. Given the transient nature of this source and the relatively low X-ray count rate, it is not particularly surprising that it was not detected by Swift BAT, given the relatively high flux sensitivity limits of Swift BAT.

3.2 Timing Properties of the Hard X-ray Flares

The average duration of a flare is roughly 1 day, with measurements ranging from 0.9-1.6 days depending on how the flare times are measured and whether we include isolated flares (i.e. flares that last for only a single GTI). The average time between flares is roughly 5 days, as measured from the last GTI of one flare to the first GTI of the next flare. However, we note that the spread of the time between flares ranges from approximately 1 day to up to 15 days and is much larger than that of the flare duration. It is, however, worth noting that there is significant variability even within a given NICER GTI for AT2019cuk, as the longest GTIs (with more than 1000 seconds of exposure) tend to show a factor of two change in hardness on timescales as short as 200 seconds. Despite the repetitive nature of the flares, a Lomb-Scargle analysis on either the hard X-ray band (2-4 keV) or the hardness ratio over the entire observing period shows no significant periodicity or quasi-periodicity.

3.3 X-ray Spectral Modeling

All spectral fitting detailed in this section is performed in XSPEC (v12.12.0; Arnaud, 1996).

3.3.1 XMM-Newton/NuSTAR

We utilize the high quality XMM-Newton and NuSTAR spectra to perform spectral modeling that informs how we fit the short NICER observations. We simultaneously fit the two XMM-Newton spectra from 0.3-10 keV and the two spectra from NuSTAR FPMA and FPMB from 3-70 keV. Each redshifted model component is fixed at $z=0.081$, abundances are adopted from Wilms et al. (2000), and cross-sections are adopted from Verner et al. (1996). The column density of Galactic absorption (tbabs component in XSPEC) is fixed at the HI value of $N_{H}=2.26\times 10^{20}$ cm*-2* (HI4PI Collaboration et al., 2016). We allow for a constant offset between the XMM-Newton spectra and each NuSTAR spectrum for cross-calibration and to account for the lack of simultaneity of these observations. Otherwise, all parameters are tied across all four spectra when fitting.

A simple absorbed cutoff power law (i.e. tbabsztbabszcut in XSPEC notation) produces a poor spectral fit ($\chi^{2}_{\nu}=\chi^{2}/\mathrm{dof}=1.32$) with significant residuals in the soft part of the spectrum (see second panel of Figure 4). Adding an additional blackbody as a phenomenological model for the soft excess (i.e. tbabsztbabs(zcut+zbb) in XSPEC notation) improves the fit and provides a statistically good fit to the data ($\chi^{2}_{\nu}=1.09$). As a more physically-motivated approach, we also explored whether the soft features in the XMM-Newton spectra can be modeled with a warm absorber (similar to the modeling presented in Dou et al., 2022). To model this warm absorber, we produced a grid of photoionization models with variable column density and ionization parameter using xstar (Kallman & Bautista, 2001), assuming a power law ionizing spectrum with $\Gamma=2$ and a turbulent velocity of $v_{\mathrm{turb}}=300$ km s*-1*. This ionized absorption fit is statistically comparable to the zcut+zbb model with $\chi^{2}_{\nu}=1.11$, with no need for an additional neutral absorber or the blackbody component (i.e. tbabs*(partcov*xstar)*zcut model in XSPEC). In this model, we leave the column density, ionization parameter, covering fraction, and redshift of the ionized absorber free when fitting, but find that at 90% confidence, the absorber is consistent with being at the redshift of the source. The blueshift of the absorber could be determined with the higher energy resolution of the RGS instrument on XMM-Newton, but the low flux in these observations makes the RGS observations very close to the background level.

The fit results for each of these models are given in Table 1, and a ratio plot for each of these fits is shown in Figure 4. Overall, the XMM-Newton and NuSTAR data reveal that in early 2022 AT2019cuk had a fairly normal X-ray spectrum compared to other AGN. In particular, the photon index and cutoff energy are comparable to standard AGN (e.g. Ricci et al., 2017), and the ionized absorber is similar to the commonly-observed warm absorber in AGN.

3.3.2 NICER

To assess the spectral evolution during the hard X-ray flares in AT2019cuk during the NICER observing campaign, we perform spectral modeling of the NICER data on a per-GTI basis. We apply the same models that we fit to XMM-Newton/NuSTAR in the previous section, except that we utilize a power law rather than a cutoff power law when fitting NICER data as the NICER instrument has lower effective area at high energies and becomes background-dominated at much lower energy than we have access to with NuSTAR. All fits with NICER are performed in the 0.3-4 keV range, although we also fit each GTI in the 0.3-$E_{\mathrm{upper}}$ keV range, where $E_{\mathrm{upper}}$ denotes the energy at which the background begins dominating over the source emission, and found that there is no considerable difference to the fit results. However, at late times (MJD $>$ 59750) when the sun angle is decreasing and the overall flux has dropped, some of the GTIs are dominated by the background below 1 keV, which makes it very difficult to properly constrain fit parameters like $\Gamma$. Hence, we exclude the GTIs with $E_{\mathrm{upper}}<1$ keV in spectral fitting. Additionally, for each model, around 5% of GTIs are poorly fit with $\chi^{2}_{\nu}>2$, which also primarily occurred at late times when the soft part of the spectrum was potentially contaminated by the NICER noise peak due to low sung angle and the baseline model for AT2019cuk (set by the XMM-Newton/NuSTAR data from December 2021/January 2022) may have changed. Thus, we excluded these GTIs from our spectral analysis. When included, we bound the column density of the neutral absorber to be larger than $N_{H}=10^{19}$ cm*-2*, the column density of the ionized absorber to be within $N_{H}=10^{21}-10^{24}$ cm*-2*, the blackbody temperature to be within $kT_{\mathrm{bb}}=0.01-0.3$ keV, and the photon index to be within $\Gamma=0-3$. We also remove GTIs in which either the blackbody temperature, the photon index, or the column density of the ionized absorber was unconstrained between the two limits provided.

When modelled phenomenologically with a blackbody and power law (i.e. tbabsztbabs(zpo+zbb) in XSPEC notation), the largest difference between the flaring and non-flaring spectra is the distribution of photon indices. The weighted mean of the flaring photon indices is $\langle\Gamma_{\mathrm{flare}}\rangle=0.77\pm 0.17$ (where the error represents the $1\sigma$ scatter in the distribution), which is significantly lower than the non-flaring spectra ($\langle\Gamma_{\mathrm{non-flare}}\rangle=1.47\pm 0.25$) and the distribution that is seen in standard AGN ($\Gamma\approx 1.6-2.0$; e.g. Ricci et al., 2017). We also test a partial covering neutral absorber (i.e. tbabspcfabs(zpo+zbb)) with the photon index and blackbody temperature fixed at the values from the XMM-Newton/NuSTAR best-fit model to explore potential degeneracies between the photon index and neutral column density. This model yields comparable fit statistics with changing the photon index, but requires roughly a simultaneous increase in the column density and normalization of the power law, the implications of which we discuss further in Section 4.1.2.

We also explore whether this difference between flares and non-flares can be explained by changes in an ionized absorber by fitting each NICER GTI with the ionized absorption model (tbabs*(partcov*xstar)*zpo in XSPEC). Due to the limited spectral quality in individual NICER GTIs, we fix the covering fraction and ionization parameter at the XMM-Newton/NuSTAR values given in Table 1. We then proceed with two possible models—one in which we fix the photon index and fit for the column density of the ionized absorber, and vice versa. We find that these two models give similar fit statistics on average ($\langle\chi^{2}_{\nu,\,\mathrm{flare}}\rangle\approx 1.14-1.19$), comparable to the values from the phenomenological model previously discussed. With the column density free, both an increase in the intrinsic normalization of the power law and the column density of the ionized absorber are required, similar to the behavior seen in the partial covering neutral absorption model (see Section 4.1.2 for further discussion). This model produces similar fit statistics on average compared to the free photon index fits, but the free column density model struggles to produce the hardest X-ray spectra seen in the most extreme flares. With the photon index free, the flares have $\langle\Gamma_{\mathrm{flare}}\rangle=1.37\pm 0.19$, which is much more reasonable compared to the phenomenological model in which the blackbody may be leading to artificially low photon indices by trying to reproduce the features of the warm absorber. However, these photon indices are still harder than standard AGN and the non-flares, which have $\langle\Gamma_{\mathrm{non-flare}}\rangle=1.75\pm 0.14$. This intrinsic change to the spectral shape could result from increased magnetic activity in the corona, which we discuss further in Section 4.1.1. In Figure 5, we show an example of one non-flaring GTI (from ObsID 4578010170 on MJD 59672) and one flaring GTI (from ObsID 4578010169 on MJD 59671), along with various model fits. The results of the ionized absorber fits with a free photon index are shown in Figure 6, which highlights the clear difference between the photon indices in flares and non-flares.

In the case where we allow both the column density of the ionized absorber and the photon index to vary, there are degeneracies between a decreasing photon index and an increasing absorber column density, as both produce harder X-ray spectra, as seen in the flares. These degeneracies are difficult to decouple with the limited spectral quality of single NICER GTIs and limited energy range considered. We do find a statistically significant improvement in the average fit statistic ($\langle\chi^{2}_{\nu,\,\mathrm{flare}}\rangle=1.01$) when allowing both parameters to be free, but it is difficult to reconcile a simultaneous increase in the column density and hardening of the coronal emission physically. Similarly, allowing only the ionization parameter to be free, rather than the column density, requires a simultaneous decrease to the ionization parameter and an increase in the normalization (i.e. intrinsic luminosity) of the power law. Such changes are counter-intuitive based on the definition of the ionization parameter ($\xi=L_{\mathrm{ion}}/nr^{2}$), and thus indicate that either a significant increase in density or distance of the obscuring material would be required, neither of which are feasible on such short timescales.

3.4 Optical/UV Variability During NICER Coverage

The optical/UV light curves are shown in the bottom three panels of Figure 1, with the full ZTF light curve from 2019 onward shown in the top panel for reference. The most notable feature of the multi-wavelength data is the lack of variability on roughly 30 day timescales or shorter, as seen in and predicted by Jiang et al. (2022). There is an interesting spike in the X-rays around MJD 59740, followed by a significant decrease in X-ray flux, which is also seen in the UV and optical bands. However, the variability seen in the drop is similar to standard AGN behavior, whereby the strongest variability is in X-rays and decreases with increasing wavelength. This sort of behavior can be explained by the standard reprocessing scenario, whereby X-rays illuminate the disk and imprint their variability on short timescales, which gets dampened by reprocessing at longer wavelengths (e.g. Edelson et al., 1996; Cackett et al., 2007; Panagiotou et al., 2022).

4 Discussion

4.1 Possible Origins of the Hard X-ray Flares

X-ray spectra with photon indices as low as $\Gamma\sim 1$, as seen in the flares of AT2019cuk, are highly unusual for AGN. We propose three potential models for these exotic hard X-ray flares, detailed below and shown schematically in Figure 7.

4.1.1 A Variable Corona/Jet

While the photon indices reached in the flares of AT2019cuk are rather low compared to standard AGN, they do resemble some of the more extreme values seen in some jetted systems like blazars (e.g. Giommi et al., 2012; Tagliaferri et al., 2015; Bhatta et al., 2018) and jetted tidal disruption events (TDEs; e.g. Saxton et al., 2012). Similarly, rapid hard X-ray variability reminiscent of the flaring in AT2019cuk is seen in both blazars (e.g. Hayashida et al., 2015; Giommi et al., 2021), jetted TDEs (e.g. Burrows et al., 2011; Reis et al., 2012), and black hole binaries (e.g. Walton et al., 2017). In the jetted TDE Swift J1644+57, this hard X-ray variability has previously been attributed to the jet wobbling into and out of our line of sight (see Tchekhovskoy et al., 2014), on timescales on the order of hours to days, similar to the timescale of flares in AT2019cuk. However, AT2019cuk appears to be a radio-quiet AGN, with VLBI observations showing that the radio emission is rather weak (with radio-loudness parameter $R\approx$ 1) and contained within $\lesssim$100 pc (An et al., 2022b). Similarly, the optical data do not show rapid variability and the luminosity of AT2019cuk is too dim to be relativistically beamed, both of which are common properties of blazars.

Therefore, while we do not expect these flares to arise from a powerful jet pointed along our line of sight, we can draw on the similarity of the flares to the hard X-ray behavior of jetted systems. Magnetic reconnection has long been proposed as a potential heating mechanism for the X-ray corona in AGN and black hole binaries (e.g. Galeev et al., 1979) as well as for the acceleration of particles in black hole jets (e.g. Sironi et al., 2015; Sironi & Beloborodov, 2020). Simulations of particle acceleration in jets have shown that magnetic reconnection can produce relatively hard particle spectra ($p\sim 1.2$ in the limit of large magnetization, where the particle distribution is a power law with $n(\gamma)\propto\gamma^{-p}$). If these particles are emitting X-rays through a non-thermal radiative mechanism, then the resulting spectrum would be a power law with photon index $\Gamma=1+\frac{p-1}{2}\sim 1.1$ (e.g. Guo et al., 2014; Sironi & Spitkovsky, 2014), comparable to what we see in the flares of AT2019cuk. GRMHD simulations suggest that the base of the jet should have high magnetization (e.g. McKinney, 2006; Tchekhovskoy et al., 2011; Chatterjee et al., 2019), and is thus a natural place for magnetic reconnection to occur. Magnetic reconnection is therefore a viable method to produce hard X-ray spectra with non-thermal particle distributions, but such distributions are rather unusual in radio-quiet sources like AT2019cuk, which typically have coronae that produce a thermal Comptonization spectrum and are limited in spectral hardness by pair production.

A hardening of the observed coronal emission can also be explained by a dynamic and outflowing corona, as fewer disk photons reach and cool the corona, leading to extreme heating and photon-starvation of the corona (e.g. Beloborodov, 1999). Outflowing coronae have been claimed in a handful of AGN, but as evidenced from low reflection fractions rather than abnormally hard X-ray spectra (e.g. Markoff & Nowak, 2004; Markoff et al., 2005; Wilkins & Gallo, 2015; King et al., 2017; Wilkins et al., 2022). Unfortunately, we cannot simultaneously probe the reflection fraction and photon index with the existing NICER data due to the background-dominated nature of the spectrum in the Fe K band. However, the XMM-Newton/NuSTAR observations presented in this work show only a weak Fe K line, suggesting that a low reflection fraction is not unlikely in this system, even in the non-flaring spectrum.

We suspect that the hardening of the X-ray spectrum during the flares is likely related to increased magnetic activity in the corona driving repetitive reconnection events (see the left panel of 7). Despite being more extreme than standard corona variability, the timescales in which these flares occur are not atypical for coronal variability seen in typical AGN. One potential way to test this corona model for the X-ray flares is through high-cadence radio monitoring to probe the variability of the parsec-scale jet, as the jet and corona are likely closely related.

4.1.2 Variable Obscuration

When modeled with a simple phenomenological model such as a power law, highly obscured AGN can exhibit effective photon indices close to $\Gamma_{\mathrm{eff}}\sim 1$, as a result of not properly accounting for obscuration (e.g. Rovilos et al., 2014; Ricci et al., 2017; Wang et al., 2022). In general, the degeneracy between the obscuring column density and the photon index can only be broken by going to higher energies ($E>10$ keV), which are less affected by the obscuring material than softer energies. Thus, with its extended energy range and high effective area relative to other hard X-ray observatories, NuSTAR provides crucial constraints on the nature of highly obscured AGN (e.g. Arévalo et al., 2014; Baloković et al., 2014; Bauer et al., 2015; Brightman et al., 2015; Marchesi et al., 2018; LaMassa et al., 2019). The NuSTAR observation of AT2019cuk shows a rather typical X-ray spectrum without significant signs of heavy neutral obscuration. Likewise, both the NuSTAR and XMM-Newton observations lack a strong narrow Fe K$\alpha$ line, which is expected and commonly observed in Compton-thick AGN (e.g. Matt et al., 1996; Levenson et al., 2006). Thus, we don’t expect to see such heavy obscuration reminiscent of Compton-thick AGN in AT2019cuk, but lower levels of obscuration could still play a role in giving the flares their hard X-ray spectral shape.

Obscuration alone is expected to produce dips in the soft band while the hard band should remain relatively unchanged, since obscuration preferentially attenuates the soft X-ray band. Thus, to obtain the hard X-ray flares seen in AT2019cuk, something more complicated than a simple increase in the amount of obscuration must be going on. One way to increase the hard X-ray flux is to increase in the mass accretion rate, but this must be coupled to an increase in the level of obscuration that leads to a roughly steady soft X-ray flux. Such a scenario could occur if a sudden influx of material drives a wind that increases the amount of obscuring material (see middle panel of Figure 7). However, this scenario requires fine-tuning as the timescales for increasing the luminosity and obscuration must be comparable in every flare seen. The day-long timescales in the flares are also extremely fast for a black hole with a mass of roughly $10^{8}$ $M_{\odot}$ and are comparable to the dynamical timescale within roughly tens of gravitational radii. The rapid variability in the obscuration would require the wind to be driven from very close to the black hole, which is much closer than where either the neutral or mildly ionized absorbers are typically found in AGN (see e.g. Laha et al., 2021, and references therein).

Another possibility is that the flares could arise from the brief unveiling of a highly obscured sight line. In this scenario, the flares are the result of the column density of this sight line decreasing (down to $N_{H}\approx 10^{23}$ cm*-2*), temporarily revealing a hard X-ray spectrum. The increase in the observed luminosity could then be due to this additional light reaching us, rather than an intrinsic change to the mass accretion rate. To produce such rapid and repetitive spectral changes, the obscuring material must be in an optically-thick, clumpy torus, and to avoid a strong narrow Fe K$\alpha$ line that is commonly seen in highly obscured AGN, the obscuring material must either have a low covering fraction or be oriented close to edge-on such that the opposite side of the nucleus is obscured. However, this is difficult to reconcile with the fact that the system spends most of its time ($\sim 92\%$) in the non-flaring, typical unobscured AGN state.

4.1.3 Self-Lensing Events in a Binary System

Another possible model for the repeating flares in AT2019cuk is self-lensing in a binary SMBH, whereby emission from an accretion flow around one or both SMBHs in a binary can be gravitationally lensed once per orbit and generate repeating flares (D’Orazio & Di Stefano, 2018). One potential self-lensing flare was reported in the Kepler light curve of a binary SMBH candidate (Hu et al., 2020). For SMBHs with $M\approx 10^{6-9}~{}M_{\odot}$ and with orbital periods of days to years, the predicted lensing flares are bright ($10\%-1000\%$ magnification), last days to months in duration, and should occur in roughly 10% of accreting SMBH binaries (D’Orazio & Di Stefano, 2018).

The self-lensing model predicts that the only wavelength dependence in the flares should arise from a finite source size with respect to the Einstein radius of the lensing SMBH. The Einstein radius scales as $M_{\mathrm{BH}}^{1/2}$, but the emitting region at a fixed gravitational radius scales linearly with the black hole mass. Thus, finite source size effects are expected in higher mass systems, and in particular are possible for X-ray emission from within $R\approx 10R_{g}$ for a $M_{\mathrm{BH}}\approx 10^{8}\,M_{\odot}$ binary with a roughly equal mass ratio. Thus, these finite source effects could potentially produce the hardening of the X-ray spectrum, if the soft X-rays are produced at a larger radius than the hard X-rays. Additionally, the lack of flares in the UV and optical bands could be a result of most of that emission originating in a circumbinary disk, which is not lensed.

On the other hand, the lack of periodicity in the hard X-ray flares is inconsistent with the expectations of the self-lensing hypothesis, as the flares should arise on timescales set by the orbital period. One possible way around this requirement is with precession, but this would require the recurrence time of the flares to be correlated and no such behavior is seen in AT2019cuk. Hence, the irregular repetition of the flares in AT2019cuk poses a significant problem for the standard self-lensing hypothesis.

4.2 No Quasi-Periodic Behavior on Timescales of 30 days

None of the data collected since early 2022 shows a clear modulation on the order of 30 days or less, as would have been expected for the proposed binary SMBH system (Jiang et al., 2022). The NICER data are potentially contaminated by a nearby AGN (see Section B of the Appendix for more detail), whose variable nature makes it difficult to draw conclusions on the net overall flux from AT2019cuk. However, the Swift XRT and UVOT data, along with the ground-based optical data, also do not show any significant periodic variability on timescales of roughly 30 days. The lack of periodic behavior across the electromagnetic spectrum suggests that the initial behavior found in ZTF was due to stochastic red noise variability commonly observed in AGN. Despite the lack of periodic behavior, the large increase in the optical and infrared flux around 2018 and the corresponding change in the optical Balmer lines still make this a very compelling changing-look AGN.

4.3 Comparison to Other Nuclear Transients

4.3.1 Quasi-Periodic Eruptions (QPEs)

Flares on the order of hours-days are also seen in the recently discovered QPEs (Miniutti et al., 2019; Giustini et al., 2020; Arcodia et al., 2021; Chakraborty et al., 2021). However, there are numerous differences between the hard X-ray flares in AT2019cuk and QPEs, including that QPEs are primarily soft X-ray emitting sources, show much larger amplitude flares, and return to a stable low X-ray luminosity phase in between flares. Thus, the flares seen in AT2019cuk do not appear QPE-like.

4.3.2 Changing-Look AGN

The optical spectrum of AT2019cuk has shown clear changes since the SDSS spectrum taken in 2005 with a significant change to the broad H$\alpha$ line and the appearance of a broad component of the H$\beta$ complex (Jiang et al., 2022), which is the defining signature of a changing-look AGN. Similar to AT2019cuk, an accompanying rise in optical flux has been seen in numerous changing-look AGN with newly emerging broad lines (e.g. Shappee et al., 2014; MacLeod et al., 2016; Trakhtenbrot et al., 2019), which is commonly attributed to an increase in the mass accretion rate, although the mechanism driving this change is still debated (see Ricci & Trakhtenbrot, 2022, and references therein). In the X-ray band, many changing-look AGN exhibit spectral changes, but most are not outside of the typical spectral variations seen in standard AGN (e.g. Oknyansky et al., 2019; Parker et al., 2019; Wang et al., 2020; Guolo et al., 2021). Unfortunately, there are no X-ray observations that pre-date the optical outburst of AT2019cuk, except for ROSAT survey upper limits, so we cannot comment on the long-term X-ray variability. However, while AT2019cuk often looks like a relatively standard AGN in its non-flaring state, it shows the first evidence for day-long, hard X-ray flares in a changing-look AGN. The high-cadence NICER observations of AT2019cuk were crucial to the discovery of these hard X-ray flares and motivate more high-cadence X-ray monitoring of changing-look AGN to determine whether such X-ray behavior is common in these systems.

5 Summary

During a high-cadence monitoring program of the binary SMBH candidate and changing-look AGN AT2019cuk with NICER, we discovered hard X-ray flares unlike typical AGN flares. The flares appear most prominently in the 2-4 keV band, are roughly day-long, and are repetitive in nature, although they do not show signs of periodicity. We performed spectral modeling of all of the NICER data on a GTI-by-GTI basis, using both a blackbody and ionized absorber for the soft excess found in the XMM-Newton data. The flares have anomalously hard X-ray spectra with significantly lower photon indices than are typically seen in AGN and the non-flaring data. We present three possible models, namely corona/jet variability driven by increased magnetic activity, variable obscuration, and self-lensing in a binary SMBH, for these hard X-ray flares. We favor the corona variability model, given that the obscuration model requires significant fine-tuning of the timescales involved and the self-lensing model is difficult to reconcile with aperiodic flares. Finally, the multi-wavelength monitoring presented here shows no signs for 30 day periodicity, as was first reported in Jiang et al. (2022).

We thank the referee for feedback that helped improve this work. DRP was partially funded by the NASA grant 80NSSCK0962. DJD received funding from the European Union’s Horizon 2020 research and innovation programme under Marie Sklodowska-Curie grant agreement No. 101029157, and from the Danish Independent Research Fund through Sapere Aude Starting Grant No. 121587. ECF is supported by NASA under award number 80GSFC21M0002. MI, GSHP acknowledge the support from the National Research Foundation of Korea grants, No. 2020R1A2C3011091 & No. 2021M3F7A1084525, and the Korea Astronomy and Space Science Institute grant under the R&D program (Project No. 2022-1-860-03) supervised by the Ministry of Science and ICT.

We thank the staff of CBNUO and LOAO, Hyori Jeon, Haeun Kim, Jeongmuk Kim, Chaerin Kim, Jae-Hyuk Yoon, and In-Kyung Paek, for their help with the optical observations. This work made use of the data obtained at CBNUO and LOAO. LOAO is operated by the Korea Astronomy and Space Science Institute.

Appendix A Verifying the Hard X-ray Flares

A.1 NICER Data Quality Checks

Given the decreasing effective area of NICER at high energies and uncertainties in background modeling for NICER data, the hard X-ray flares seen in AT2019cuk require extreme scrutiny to ensure that they are neither instrumental nor a result of background contamination. In addition, periods of high background flaring may become more prevalent as we approach solar maximum and experience periods of increased solar activity, whereas NICER calibration data were primarily acquired early in the mission near solar minimum. Thus, we take numerous precautions in addition to the relatively restrictive data filtering presented in Section 2.1 to verify that these flares are not instrumental.

We first ensure that these flares are not a result of the ISS being too close to either the polar horn or the South Atlantic Anomaly (SAA) regions of Earth’s magnetosphere, both of which are known to produce significant background counts indistinguishable from cosmic X-rays. The flares are distributed throughout the orbital path of the satellite and are not preferentially located near either of these problematic regions. Next, we checked numerous other quantities associated with each GTI, including the exposure time, background rate, high-energy (15-18 keV) raw count rate, sunshine factor, and cutoff rigidity factor (COR), to compare the distribution of these quantities between flares and non-flares. We find that there is no significant difference between the distributions of any of these factors in flaring versus non-flaring GTIs. In addition, we check the undershoots and overshoots for each GTI and find that the flares and non-flares follow the same distributions and are reasonably distributed below the normal filtering levels. Additionally, the overshoot rates also fall below the COR-dependent value from the overonly_expr term in standard nicerl2 filtering. In Figure 8, we show the evolution of many of these parameters used in checking the validity of the hard flares, with a distribution of the non-flares and flares for each parameter to the right, each of which show consistency between the two populations.

Finally, another concern is that many of the flares occurred in only single, isolated GTIs. Although there is nothing in the data quality checks detailed above that hints at issues with these data points, we also checked for consistency among isolated flares (one GTI surrounded by non-flaring GTIs on either side) and non-isolated flares ($>2$ GTIs in a row marked as a flare). We find that there are no differences between the isolated and non-isolated flares, either in the data quality checks presented in Figure 8 or in the spectral fitting results. This suggests that the two are not independent distributions and the isolated flares can be treated as just as valid as the non-isolated flares. Likewise, we note that the number of isolated flares was greatest in the first 100 days or so of the NICER monitoring, in which the sampling rate was lower.

A.2 Simultaneous Swift Observations

Given the extensive Swift observations of AT2019cuk, we also check our NICER spectral modeling results by comparing to Swift spectral modeling. Swift has imaging capabilities that NICER lacks, which is beneficial because it allows us to ensure that the results are not driven by a nearby AGN (which is further discussed in Section B) nor the uncertainty in the empirical NICER background. Swift spectra and background spectra were created on a per-ObsID basis, using the online Swift XRT products generator (Evans et al., 2007, 2009). The spectra were then binned using the optimal binning procedure (Kaastra & Bleeker, 2016) and to a minimum of 1 count per bin. We fit each spectrum with $C$ statistic minimization, given the limited spectral quality and lower count rate of Swift relative to NICER. We fit the Swift spectra in the 0.3-10 keV range to the model tbabsztbabs(zpo+zbb), with the blackbody temperature fixed to the XMM-Newton/NuSTAR value of $kT_{\mathrm{bb}}=0.081$ keV and the column density of the neutral absorber fixed to the XMM-Newton/NuSTAR value of $N_{H}=6.9\times 10^{20}$ cm*-2*. We fixed these parameters to eliminate any degeneracies with the photon index due to the limited spectral quality of Swift XRT data. The resulting photon indices as a function of time are shown in Figure 9, with a comparison to all NICER fits as described in Section 3.3.

We find that the Swift and NICER data have very similar distribution of photon indices, albeit with large error bars on the Swift data. Some of the Swift data points have relatively low photon indices with $\Gamma\lesssim 1$, which is similar to what we find in the flaring NICER data with the same model. Thus, we find similar evidence in the Swift data for anomalously low photon indices and hard X-ray spectra, indicating that the flares seen in NICER are not a result of the contaminating source nor a result of poorly understood NICER background.

There is some evidence for a lower photon index preferentially found in the Swift data. This has two possible explanations–(1) NICER data are contaminated by the nearby AGN that has a softer spectrum with $\Gamma\approx 1.7$ (see Section B) or (2) Swift spectra are a combination of both flaring and non-flaring time periods, given the longer baseline in each data point due to the per-ObsID sampling. The first scenario is disfavored because the XMM-Newton/NuSTAR spectra, which also have the ability to clearly resolve the two sources, show a similar photon index to the average NICER value (when in a non-flaring state). Thus, we suspect that the lower photon indices seen in the Swift modeling are a result of coarser binning than we use with NICER, which could lead to some Swift data points including both flaring and non-flaring time periods. We require this rather coarse, ObsID-based binning rather than per-GTI binning with Swift in order to have sufficient signal to perform spectral fitting. However, when we bin the Swift light curve on a per-GTI basis (similar to NICER binning), we find that there are similar looking flares in the Swift light curve (see the bottom two panels of Figure 13).

A.3 Comparison to Simulated NICER Spectra

The NICER observations are also rather short and therefore contain a relatively low number of counts compared to the longer XMM-Newton and NuSTAR observations. Hence, to check whether the flares could arise from noise induced by counting statistics, we simulated a similar distribution of NICER observations from the best fit to the XMM-Newton/NuSTAR with the partial covering warm absorber. As each of the NICER GTIs has a different flux, exposure time, and background spectrum, we utilize the properties of each GTI to produce the simulated data, including rescaling the best fit XMM-Newton/NuSTAR spectrum to the flux of the NICER GTI. The simulated observations do not reveal any strong hard X-ray flares in the 2-4 keV band, as are seen in the real data, which can be seen qualitatively in the left panel of Figure 10 where we plot both the real and simulated hard X-ray light curves. In addition, the rightmost panel of Figure 10 shows that the hardness ratio seen in the real data is significantly harder than the simulated hardness ratio. Thus, the flares are not a result of noise due to counting statistics and cannot be attributed to standard AGN variability.

Appendix B Impact of a Nearby AGN on NICER Data

Upon investigating AT2019cuk with X-ray imaging telescopes including XMM-Newton, NuSTAR, and Swift, we noticed a nearby strong X-ray source located approximately 3.1’ away from AT2019cuk. This interloping source is a known AGN, located at ($\alpha$,$\delta$) = (14:30:08.65, +23:06:21.62) (Véron-Cetty & Véron, 2010). At 3.1’ away from AT2019cuk, this interloper is just within the field of view of NICER, but at a location such that the detector response has dropped by a factor of 2-3 in sensitivity, with the exact response depending on the exact pointing and which detectors are used. We compute the relative response of NICER’s X-ray Timing Instrument at the location of this interloper by using nicerarf with the RA and Dec set to the position of the interloper. We can then define the “effective ARF weight” as the mean of the weighting for each detector at the location of AT2019cuk divided by the mean of the weighting for each detector at the location of the interloper. A high effective ARF weighting means that AT2019cuk dominates the flux in the NICER pointings. In Figure 11, we show this effective ARF weighting over the entire observing period, together with the 2-4 keV light curve for easy comparison with the flares. From this it is clear that the flares are not simply an effect of pointing oscillations, as they appear in both times of high and low effective ARF weighting, with a similar distribution to the non-flaring points.

We also extracted a NuSTAR spectrum for both AT2019cuk and the interloper from the DDT observation in January 2022 to compare the spectral shapes between AT2019cuk and the interloper. Both spectra can be fit well with an absorbed power law. We leave the photon indices free between the two sources, and the best fit gives $\Gamma=1.74_{-0.07}^{+0.08}$ for AT2019cuk and $\Gamma=1.77_{-0.14}^{+0.17}$ for the interloper. Thus, as the spectral shape of the interloper is extremely similar to that of AT2019cuk and we see significant changes in the spectral shape during the hard X-ray flares, we do not expect the flares to be a result of the interloper.

Beyond the impact of the interloper on the NICER flares, we also need to assess the impact of the interloper on the overall variability of AT2019cuk. This is especially important for assessing any long-term periodicity of the light curve that could be suggestive of a binary SMBH. Given their ability to distinguish between the two sources, we used both NuSTAR and Swift to investigate the influence of the interloper relative to AT2019cuk. In the initial NuSTAR observation in January 2022, the interloper was a factor of three lower in flux than AT2019cuk. In Figure 12, we show the Swift light curve for these two sources during the NICER campaign, including a ratio of their count rates in the 0.3-4 keV band in the bottom panel. For the majority of the observing campaign, the interloper is brighter than AT2019cuk in the 0.3-4 keV band, although there are some time periods in which the interloper is brighter by a factor of no more than two. However, combined with the drop in the effective response at the location of the interloper of at least two, we find that the interloper never contributes more than 50% of the flux in the NICER band and usually contributes rather negligibly. Given the variability in the amount of flux from the interloper and the large gaps in Swift monitoring of this source during the NICER campaign, it is difficult to assess exactly how much impact this has on the overall variability of AT2019cuk with NICER.

In addition, we can utilize the Swift monitoring of AT2019cuk to assess the effect of the interloper on the NICER data. In Figure 13, we show both the NICER and the Swift light curves of AT2019cuk, where in each of the panels the NICER light curves have been scaled by a constant factor to account for the difference in the effective area between NICER and Swift. These factors were determined using the HEASARC WebPIMMS tool222https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl and assuming a $\Gamma=1.7$ power law for the conversion between the different count rates. In the 0.3-4 keV, 0.3-2 keV, and 2-4 keV bands respectively, these renormalization factors are 20, 25, and 10, where these factors represent the amount by which we divided the NICER data to accurately compare to the Swift data. From Figure 13, we can see that the NICER variability matches the Swift variability of AT2019cuk very well, including picking up key features of the light curve like the drop in flux around MJD 59750. Thus, this indicates that most of the variability in the NICER data is coming from AT2019cuk rather than the interloper.