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

Megan Masterson MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Erin Kara MIT
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA  Dheeraj R. Pasham MIT Kavli Institute
for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Daniel J. D’Orazio Niels Bohr International Academy,
Niels Bohr Institute, Blegdamsvej 17, DK-2100 Copenhagen, Denmark  Dominic J.
Walton Institute of Astronomy, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK  Centre for Astrophysics Research, University of
Hertfordshire, College Lane, Hatfield AL10 9AB, UK  Andrew C. Fabian Institute
of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
Matteo Lucchini MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Ronald A.
Remillard MIT Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Zaven
Arzoumanian  X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center,
Code 662, Greenbelt, Maryland 20771, USA  Otabek Burkhonov Ulugh Beg
Astronomical Institute, Astronomy Street 33, Tashkent 100052, Uzbekistan
Hyeonho Choi SNU Astronomy Research Center, Astronomy Program, Department of
Physics and Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu,
Seoul, 08826, Republic of Korea  Shuhrat A. Ehgamberdiev Ulugh Beg
Astronomical Institute, Astronomy Street 33, Tashkent 100052, Uzbekistan
National University of Uzbekistan, Tashkent 100174, Uzbekistan  Elizabeth C.
Ferrara Department of Astronomy, University of Maryland, College Park, MD,
20742, USA  Center for Exploration and Space Studies (CRESST), NASA/GSFC,
Greenbelt, MD 20771, USA  NASA Goddard Space Flight Center, Greenbelt, MD
20771, USA  Muryel Guolo Department of Physics and Astronomy, Johns Hopkins
University, 400 N. Charles St., Baltimore, MD 21218, USA  Myungshin Im SNU
Astronomy Research Center, Astronomy Program, Department of Physics and
Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul, 08826,
Republic of Korea  Yonggi Kim Department of Astronomy and Space Science,
Chungbuk National University Observatory, Chungdae-ro 1, Seowon-Gu, Cheongju,
Chungbuk 28644, Republic of Korea  Davron O. Mirzaqulov Ulugh Beg Astronomical
Institute, Astronomy Street 33, Tashkent 100052, Uzbekistan  Gregory S. H.
Paek SNU Astronomy Research Center, Astronomy Program, Department of Physics
and Astronomy, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul,
08826, Republic of Korea  Hyun-Il Sung  Korea Astronomy and Space Science
Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
Joh-Na Yoon  Department of Astronomy and Space Science, Chungbuk National
University Observatory, Chungdae-ro 1, Seowon-Gu, Cheongju, Chungbuk 28644,
Republic of Korea Megan Masterson<EMAIL_ADDRESS>
###### 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.

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

Figure 1: Multi-wavelength light curve for AT2019cuk, with zoom in on the
duration of the high-cadence NICER monitoring campaign. Top panel: Full ZTF
light curve from January 2019 to September 2022. Red diamonds show the ZTF
$r$-band data, offset by 0.8 mag, and yellow squares show the ZTF $g$-band
data. The orange shaded region shows the prediction of the merger time from
Jiang et al. (2022), using post-Newtonian modeling of both optical and X-ray
data (see Figure 3/Extended Figure 9 of Jiang et al., 2022). The grey box
represents the zoom in region shown in the lower four panels. Second panel:
NICER (blue circles) and Swift XRT (green X’s) background-subtracted light
curves in the 0.3-4 keV band. Swift data are scaled by a factor of 20 to
account for the difference in effective area between the two telescopes. Third
panel: Swift UVW1 light curve. Data have been corrected for Galactic
extinction. Bottom two panels: Optical light curves, with the fourth panel
showing the ZTF $g$-band (yellow squares) and LOAO, MAO, and CBNUO $B$-band
(light blue pentagons, triangles, and plus signs respectively) data and the
fifth panel showing ZTF $r$-band (red diamonds) LOAO, MAO, and CBNUO $R$-band
(light pink pentagons, triangles, and plus signs respectively). There is a
$\approx$ 9-10% drop in the $g$\- and $B$-band data in the fourth panel around
MJD 59750 and a $\approx$ 3-5% drop in the $r$\- and $R$-band data at the same
time, which corresponds to the time at which a large drop in the X-ray flux
($\approx$ 50%) is also seen.

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

Figure 2: Top 3 panels: NICER light curves in three different energy bands
(from top to bottom: 0.3-4 keV, 0.3-2 keV, and 2-4 keV) for AT2019cuk. Each
light curve is background-subtracted, rescaled to the rate for 52 detectors,
and shown in counts per second. Bottom panel: NICER hardness ratio, defined as
the count rate in 2-4 keV divided by the count rate in 0.3-2 keV. The orange
dashed line shows the boundary between our definition of flaring vs. non-
flaring GTIs, set at a hardness ratio of 0.2. In all four panels, all dark
blue circles correspond to non-flaring GTIs and orange stars correspond to
flaring GTIs.

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.

Figure 3: Hardness-intensity diagram for the NICER data. Flaring GTIs are
shown as orange stars (any point with a hardness ratio greater than 0.24), and
non-flaring GTIs are shown as blue circles. The blue arrow shows a
representative track for the standard softer-when-brighter behavior seen in
AGN (e.g. Sobolewska & Papadakis, 2009), which is similar to what is seen in
the non-flaring points for AT2019cuk. The orange arrow shows the harder-when-
brighter behavior seen in the hard X-ray flares of AT2019cuk, which are
starkly different from the standard AGN variability that is expected.

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. tbabs*ztbabs*zcut 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. tbabs*ztbabs*(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.

Figure 4: Spectral fits to the XMM-Newton observations from December 2021
(Obs. 1, shown in light green) and January 2022 (Obs. 2, shown in dark blue)
and the NuSTAR observation from January 2022 (FPMA/FPMB shown in dark/light
orange, respectively). Top panel: Unfolded spectrum fit to the
tbabs*ztbabs*(zcut+zbb) model. Second panel: Ratio plot for the
tbabs*ztbabs*zcut model. Clear residuals are seen in the soft part of the
spectrum. Third panel: Ratio plot for the tbabs*ztbabs*(zcut+zbb) model.
Bottom panel: Ratio plot for the tbabs*(partcov*xstar)*zcut model. Both the
partial covering ionized absorber and the blackbody give statistically
comparable fits and significant improvement over the cutoff power law alone.

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.

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}$
|  $kT_{\mathrm{bb}}$ (eV) | $81_{-3}^{+3}$
(partcov*xstar)*zcut |  $N_{H,\mathrm{ion}}$ ($10^{22}$ cm-2) | $1.1_{-0.2}^{+1.4}$
| $\log(\xi/\mathrm{erg}\,\mathrm{cm}\,\mathrm{s}^{-1})$ | $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$

When modelled phenomenologically with a blackbody and power law (i.e.
tbabs*ztbabs*(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. tbabs*pcfabs*(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.

Figure 5: Spectral fits to example non-flaring and flaring NICER GTIs. The
non-flaring spectrum is shown in blue, and the flaring spectrum is shown in
orange. Top panel: Unfolded spectrum fit to the tbabs*(partcov*xstar)*zpo
model with the photon index free between the two GTIs. Second panel: Ratio
plot for the tbabs*ztbabs*(zpo+zbb) model. Third panel: Ratio plot for the
tbabs*(partcov*xstar)*zpo) model with the column density of the ionized
absorber free. The photon index and the ionization parameter and covering
fraction of the ionized absorber are fixed at the best fit values from the
XMM-Newton/NuSTAR fit. Fourth panel: Ratio plot for the
tbabs*(partcov*xstar)*zpo) model with the photon index of the power law free.
The column density, the ionization parameter, and covering fraction of the
ionized absorber are fixed at the best fit values from the XMM-Newton/NuSTAR
fit. Bottom panel: Ratio plot for the tbabs*(partcov*xstar)*zpo model, but
with both the column density of the ionized absorber and the photon index
free.

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.

Figure 6: NICER spectral fits to the ionized absorption model with a variable
photon index (i.e. tbabs*(partcov*xstar)*zpo). In each panel, the dark blue
circles are the non-flaring data and orange stars the flaring data (as defined
by the hardness ratio cut detailed in Section 3.1). The left shows all of the
NICER data, with each point corresponding to a single NICER GTI, whereas the
right shows a histogram of the flaring and non-flaring data over the entire
observing period. Top Panel: Hard X-ray luminosity in the 2-4 keV band. Middle
Panel: Photon index of the power-law component. The horizontal light green
dashed line shows the best-fit value for $\Gamma$ from the XMM-Newton/NuSTAR
fit, with the shaded region showing the errors. There is generally good
consistency between the non-flares and the XMM-Newton/NuSTAR value, but the
flaring data have significantly lower photon indices on average. Bottom Panel:
Reduced $\chi^{2}$ statistic. Both the flares and the non-flares are similarly
well-fit by this model.

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.

Figure 7: Schematics for potential flare scenarios. Panel (A): Flares driven
by intrinsic changes in the X-ray corona. Increased magnetic activity in the
X-ray corona can lead to magnetic reconnection events causing the observed
spectrum to harden. Panel (B): Flares driven by changes in the line-of-sight
absorption. The flares correspond to an increase in mass accretion rate, which
drives a simultaneous increase in the observed luminosity and a wind that
increases the observed column density of the absorber. Panel (C): Flares are
the result of gravitational self-lensing in a binary SMBH.

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

Figure 8: NICER instrument parameters for the flares and non-flares. The rows
correspond to: (1) Hard count rate (2-4 keV), (2) Exposure time for the GTI,
(3) Raw high energy count rate (15-18 keV), (4) Satellite latitude, (5) NICER
sunshine factor, (6) Undershoots, (7) Overshoots, (8) Cut-off rigidity. In all
8 panels, the blue circles correspond to non-flaring data points and orange
stars correspond to flares, as defined in Section 3.1. In the lower five
panels, corresponding to NICER instrumental parameters, the error bars
indicate the minimum and maximum values of each parameter reached during a
given GTI. The right-most panels are normalized distributions for the non-
flares in blue and the flares in orange for each of the 8 parameters
previously described. No clear difference in seen between flares and non-
flares, indicating that none of these NICER instrument parameters are
responsible for the flaring behavior.

###  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 tbabs*ztbabs*(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).

Figure 9: Left: Distribution of photon indices as a function of time for both
NICER (blue) and Swift XRT (green) observations. Both sets of data are fit
with an absorbed blackbody plus power law model, but in the Swift fitting, the
column density of the absorber and the blackbody temperature are both fixed at
the value of the XMM-Newton/NuSTAR fits due to limited spectral constraints.
Any upward/downward facing triangles indicate lower/upper limits,
respectively, where the upper limit of $\Gamma$ was set at 3 and the lower
limit was set at 0. Right: Histogram of photon indices across all observations
with both NICER (blue) and Swift (green). The distributions are similar, with
both extending down to quite low photon indices compared to standard AGN. In
both panels, the orange dashed line is the XMM-Newton/NuSTAR value for the
same model, with the shaded area showing the uncertainty on that value.

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

Figure 10: Left: GTI-based hard X-ray light curves (2-4 keV) for both the real
NICER observations (top) and the simulated NICER observations (bottom) with
the same distribution exposure times as the real observations. The simulated
data are based on the joint XMM-Newton/NuSTAR fit with the partial covering
warm absorber. The real data show clear, strong flares that are not present in
the simulated data, indicating that the flares are not due to noise induced
from counting statistics. Right: Distribution of hardness ratios for both the
real (dark blue line) and simulated (light green line) data. There is a clear
skew in the real data toward harder spectra.

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

Figure 11: Top panel: NICER 2-4 keV light curves for the AT2019cuk field of
view. Bottom panel: Effective ARF weight, defined 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, both of which
are computed using nicerarf. In both panels the blue circles represent non-
flaring GTIs and the orange stars represent flaring GTIs, defined by the
hardness ratio cut detailed in Section 3.1. The right side shows histograms
for both the 2-4 keV rate and the effective ARF weight for the two different
populations (flaring vs. non-flaring) in their respective colors. As the
flares and non-flares have similar effective ARF weights, it is evident that
the flares are not a result of pointing differences or jitter.

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.

Figure 12: Top 3 panels: Swift XRT light curves in three different energy
bands (from top to bottom: 0.3-4 keV, 0.3-2 keV, and 2-4 keV). Each light
curve is grouped by GTI, is background-subtracted, and shown in counts per
second. The blue circles and green squares correspond to AT2019cuk and the
interloper, respectively. Fourth panel: Swift hardness ratio, defined as the
count rate in 2-4 keV divided by the count rate in 0.3-2 keV. The two sources
have relatively comparable hardness ratios, meaning that the hard X-ray flares
seen in NICER are not likely due to picking up more flux from the interloper.
Bottom panel: Ratio of the 0.3-4 keV count rates of AT2019cuk and the
interloper. The black dashed line shows an equal 0.3-4 keV rate between the
two sources, and any data points above the line indicate that AT2019cuk was
brighter than the interloper. In combination with the difference in the
response at the location of the interloper (see Figure 11), the bottom panel
shows decisively that the interloper contributes at worst roughly half of the
observed NICER rate for AT2019cuk.

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.

Figure 13: Top 3 panels: NICER (blue circles) and Swift XRT (green X’s) light
curves in three different energy bands (from top to bottom: 0.3-4 keV, 0.3-2
keV, and 2-4 keV) for AT2019cuk. Both the NICER and Swift data are binned by
GTI. Each light curve is background-subtracted and show in counts per second.
Bottom panel: NICER and Swift hardness ratios, defined as the count rate in
2-4 keV divided by the count rate in 0.3-2 keV. The orange dashed line shows
the corresponding hardness ratio cut-off for the flares, defined in Section
3.1. In each panel, the NICER data have been rescaled by a constant factor to
account for the effective area difference between the two instruments. There
is generally good agreement between the Swift and NICER data after this
rescaling, which indicates that AT2019cuk is contributing most of the flux
that NICER sees, rather than the interloper. Similarly, there are some
relatively hard Swift data points coincident with some of the NICER flares.

## References

  * An et al. (2022a) An, T., Zhang, X., Wang, A., et al. 2022a, The Astronomer’s Telegram, 15267, 1
  * An et al. (2022b) An, T., Zhang, Y., Wang, A., et al. 2022b, A&A, 663, A139, doi: 10.1051/0004-6361/202243821
  * Arcodia et al. (2021) Arcodia, R., Merloni, A., Nandra, K., et al. 2021, Nature, 592, 704, doi: 10.1038/s41586-021-03394-6
  * Arévalo et al. (2014) Arévalo, P., Bauer, F. E., Puccetti, S., et al. 2014, ApJ, 791, 81, doi: 10.1088/0004-637X/791/2/81
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Baloković et al. (2014) Baloković, M., Comastri, A., Harrison, F. A., et al. 2014, ApJ, 794, 111, doi: 10.1088/0004-637X/794/2/111
  * Bauer et al. (2015) Bauer, F. E., Arévalo, P., Walton, D. J., et al. 2015, ApJ, 812, 116, doi: 10.1088/0004-637X/812/2/116
  * Becker (2015) Becker, A. 2015, HOTPANTS: High Order Transform of PSF ANd Template Subtraction, Astrophysics Source Code Library, record ascl:1504.004.  http://ascl.net/1504.004
  * Beloborodov (1999) Beloborodov, A. M. 1999, ApJ, 510, L123, doi: 10.1086/311810
  * Bhatta et al. (2018) Bhatta, G., Mohorian, M., & Bilinsky, I. 2018, A&A, 619, A93, doi: 10.1051/0004-6361/201833628
  * Blecha & Loeb (2008) Blecha, L., & Loeb, A. 2008, MNRAS, 390, 1311, doi: 10.1111/j.1365-2966.2008.13790.x
  * Brightman et al. (2015) Brightman, M., Baloković, M., Stern, D., et al. 2015, ApJ, 805, 41, doi: 10.1088/0004-637X/805/1/41
  * Bruni et al. (2022) Bruni, G., Rodi, J., Panessa, F., et al. 2022, The Astronomer’s Telegram, 15306, 1
  * Burke-Spolaor (2011) Burke-Spolaor, S. 2011, MNRAS, 410, 2113, doi: 10.1111/j.1365-2966.2010.17586.x
  * Burrows et al. (2011) Burrows, D. N., Kennea, J. A., Ghisellini, G., et al. 2011, Nature, 476, 421, doi: 10.1038/nature10374
  * Cackett et al. (2007) Cackett, E. M., Horne, K., & Winkler, H. 2007, MNRAS, 380, 669, doi: 10.1111/j.1365-2966.2007.12098.x
  * Chakraborty et al. (2021) Chakraborty, J., Kara, E., Masterson, M., et al. 2021, ApJ, 921, L40, doi: 10.3847/2041-8213/ac313b
  * Charisi et al. (2016) Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MNRAS, 463, 2145, doi: 10.1093/mnras/stw1838
  * Chatterjee et al. (2019) Chatterjee, K., Liska, M., Tchekhovskoy, A., & Markoff, S. B. 2019, MNRAS, 490, 2200, doi: 10.1093/mnras/stz2626
  * Chen et al. (2022) Chen, Y.-J., Zhai, S., Liu, J.-R., et al. 2022, arXiv e-prints, arXiv:2206.11497.  https://arxiv.org/abs/2206.11497
  * Chornock et al. (2010) Chornock, R., Bloom, J. S., Cenko, S. B., et al. 2010, ApJ, 709, L39, doi: 10.1088/2041-8205/709/1/L39
  * d’Ascoli et al. (2018) d’Ascoli, S., Noble, S. C., Bowen, D. B., et al. 2018, ApJ, 865, 140, doi: 10.3847/1538-4357/aad8b4
  * D’Orazio & Di Stefano (2018) D’Orazio, D. J., & Di Stefano, R. 2018, MNRAS, 474, 2975, doi: 10.1093/mnras/stx2936
  * D’Orazio & Loeb (2018) D’Orazio, D. J., & Loeb, A. 2018, ApJ, 863, 185, doi: 10.3847/1538-4357/aad413
  * Dotti et al. (2023) Dotti, M., Bonetti, M., Rigamonti, F., et al. 2023, MNRAS, 518, 4172, doi: 10.1093/mnras/stac3344
  * Dou et al. (2022) Dou, L., Jiang, N., Wang, T., et al. 2022, A&A, 665, L3, doi: 10.1051/0004-6361/202244450
  * Edelson et al. (1996) Edelson, R. A., Alexander, T., Crenshaw, D. M., et al. 1996, ApJ, 470, 364, doi: 10.1086/177872
  * Ehgamberdiev (2018) Ehgamberdiev, S. 2018, Nature Astronomy, 2, 349, doi: 10.1038/s41550-018-0459-3
  * Eracleous & Halpern (1994) Eracleous, M., & Halpern, J. P. 1994, ApJS, 90, 1, doi: 10.1086/191856
  * Evans et al. (2007) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379, doi: 10.1051/0004-6361:20077530
  * Evans et al. (2009) —. 2009, MNRAS, 397, 1177, doi: 10.1111/j.1365-2966.2009.14913.x
  * Foord et al. (2020) Foord, A., Gültekin, K., Nevin, R., et al. 2020, ApJ, 892, 29, doi: 10.3847/1538-4357/ab72fa
  * Fu et al. (2011) Fu, H., Zhang, Z.-Y., Assef, R. J., et al. 2011, ApJ, 740, L44, doi: 10.1088/2041-8205/740/2/L44
  * Gabriel et al. (2004) Gabriel, C., Denby, M., Fyfe, D. J., et al. 2004, in Astronomical Society of the Pacific Conference Series, Vol. 314, Astronomical Data Analysis Software and Systems (ADASS) XIII, ed. F. Ochsenbein, M. G. Allen, & D. Egret, 759
  * Galeev et al. (1979) Galeev, A. A., Rosner, R., & Vaiana, G. S. 1979, ApJ, 229, 318, doi: 10.1086/156957
  * Gaskell (2010) Gaskell, C. M. 2010, Nature, 463, E1, doi: 10.1038/nature08665
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005, doi: 10.1086/422091
  * Gendreau et al. (2016) Gendreau, K. C., Arzoumanian, Z., Adkins, P. W., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, T. Takahashi, & M. Bautz, 99051H, doi: 10.1117/12.2231304
  * Giommi et al. (2012) Giommi, P., Polenta, G., Lähteenmäki, A., et al. 2012, A&A, 541, A160, doi: 10.1051/0004-6361/201117825
  * Giommi et al. (2021) Giommi, P., Perri, M., Capalbi, M., et al. 2021, MNRAS, 507, 5690, doi: 10.1093/mnras/stab2425
  * Giustini et al. (2020) Giustini, M., Miniutti, G., & Saxton, R. D. 2020, A&A, 636, L2, doi: 10.1051/0004-6361/202037610
  * Graham et al. (2015a) Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015a, Nature, 518, 74, doi: 10.1038/nature14143
  * Graham et al. (2015b) —. 2015b, MNRAS, 453, 1562, doi: 10.1093/mnras/stv1726
  * Guo et al. (2014) Guo, F., Li, H., Daughton, W., & Liu, Y.-H. 2014, Phys. Rev. Lett., 113, 155005, doi: 10.1103/PhysRevLett.113.155005
  * Guolo et al. (2021) Guolo, M., Ruschel-Dutra, D., Grupe, D., et al. 2021, MNRAS, 508, 144, doi: 10.1093/mnras/stab2550
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103, doi: 10.1088/0004-637X/770/2/103
  * Hayashida et al. (2015) Hayashida, M., Nalewajko, K., Madejski, G. M., et al. 2015, ApJ, 807, 79, doi: 10.1088/0004-637X/807/1/79
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116, doi: 10.1051/0004-6361/201629178
  * Hu et al. (2020) Hu, B. X., D’Orazio, D. J., Haiman, Z., et al. 2020, MNRAS, 495, 4061, doi: 10.1093/mnras/staa1312
  * Im et al. (2021) Im, M., Kim, Y., Lee, C.-U., et al. 2021, Journal of Korean Astronomical Society, 54, 89
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1, doi: 10.1051/0004-6361:20000036
  * Jiang et al. (2022) Jiang, N., Yang, H., Wang, T., et al. 2022, arXiv e-prints, arXiv:2201.11633.  https://arxiv.org/abs/2201.11633
  * Ju et al. (2013) Ju, W., Greene, J. E., Rafikov, R. R., Bickerton, S. J., & Badenes, C. 2013, ApJ, 777, 44, doi: 10.1088/0004-637X/777/1/44
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Kallman & Bautista (2001) Kallman, T., & Bautista, M. 2001, ApJS, 133, 221, doi: 10.1086/319184
  * Kim et al. (2019) Kim, J., Im, M., Choi, C., & Hwang, S. 2019, ApJ, 884, 103, doi: 10.3847/1538-4357/ab40cd
  * Kim et al. (2018) Kim, J., Karouzos, M., Im, M., et al. 2018, Journal of Korean Astronomical Society, 51, 89, doi: 10.5303/JKAS.2018.51.4.89
  * King et al. (2017) King, A. L., Lohfink, A., & Kara, E. 2017, ApJ, 835, 226, doi: 10.3847/1538-4357/835/2/226
  * Komossa et al. (2003) Komossa, S., Burwitz, V., Hasinger, G., et al. 2003, ApJ, 582, L15, doi: 10.1086/346145
  * Komossa et al. (2008) Komossa, S., Zhou, H., & Lu, H. 2008, ApJ, 678, L81, doi: 10.1086/588656
  * Koss et al. (2011) Koss, M., Mushotzky, R., Treister, E., et al. 2011, ApJ, 735, L42, doi: 10.1088/2041-8205/735/2/L42
  * Krimm et al. (2013) Krimm, H. A., Holland, S. T., Corbet, R. H. D., et al. 2013, ApJS, 209, 14, doi: 10.1088/0067-0049/209/1/14
  * Kurczynski et al. (2022) Kurczynski, P., Johnson, M. D., Doeleman, S. S., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12180, Space Telescopes and Instrumentation 2022: Optical, Infrared, and Millimeter Wave, ed. L. E. Coyle, S. Matsuura, & M. D. Perrin, 121800M, doi: 10.1117/12.2630313
  * Laha et al. (2021) Laha, S., Reynolds, C. S., Reeves, J., et al. 2021, Nature Astronomy, 5, 13, doi: 10.1038/s41550-020-01255-2
  * LaMassa et al. (2019) LaMassa, S. M., Yaqoob, T., Boorman, P. G., et al. 2019, ApJ, 887, 173, doi: 10.3847/1538-4357/ab552c
  * Lehto & Valtonen (1996) Lehto, H. J., & Valtonen, M. J. 1996, ApJ, 460, 207, doi: 10.1086/176962
  * Levenson et al. (2006) Levenson, N. A., Heckman, T. M., Krolik, J. H., Weaver, K. A., & Życki, P. T. 2006, ApJ, 648, 111, doi: 10.1086/505735
  * Liu et al. (2015) Liu, T., Gezari, S., Heinis, S., et al. 2015, ApJ, 803, L16, doi: 10.1088/2041-8205/803/2/L16
  * Liu et al. (2018) Liu, X., Guo, H., Shen, Y., Greene, J. E., & Strauss, M. A. 2018, ApJ, 862, 29, doi: 10.3847/1538-4357/aac9cb
  * MacLeod et al. (2016) MacLeod, C. L., Ross, N. P., Lawrence, A., et al. 2016, MNRAS, 457, 389, doi: 10.1093/mnras/stv2997
  * Marchesi et al. (2018) Marchesi, S., Ajello, M., Marcotulli, L., et al. 2018, ApJ, 854, 49, doi: 10.3847/1538-4357/aaa410
  * Markoff & Nowak (2004) Markoff, S., & Nowak, M. A. 2004, ApJ, 609, 972, doi: 10.1086/421099
  * Markoff et al. (2005) Markoff, S., Nowak, M. A., & Wilms, J. 2005, ApJ, 635, 1203, doi: 10.1086/497628
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Matt et al. (1996) Matt, G., Brandt, W. N., & Fabian, A. C. 1996, MNRAS, 280, 823, doi: 10.1093/mnras/280.3.823
  * McKinney (2006) McKinney, J. C. 2006, MNRAS, 368, 1561, doi: 10.1111/j.1365-2966.2006.10256.x
  * Miniutti et al. (2019) Miniutti, G., Saxton, R. D., Giustini, M., et al. 2019, Nature, 573, 381, doi: 10.1038/s41586-019-1556-x
  * Moiseev et al. (2022) Moiseev, A. V., Kozlova, D. V., & Kotov, S. S. 2022, The Astronomer’s Telegram, 15319, 1
  * Müller-Sánchez et al. (2015) Müller-Sánchez, F., Comerford, J. M., Nevin, R., et al. 2015, ApJ, 813, 103, doi: 10.1088/0004-637X/813/2/103
  * NASA HEASARC (2014) NASA HEASARC. 2014, HEAsoft: Unified Release of FTOOLS and XANADU, Astrophysics Source Code Library, record ascl:1408.004.  http://ascl.net/1408.004
  * Oh et al. (2015) Oh, K., Yi, S. K., Schawinski, K., et al. 2015, ApJS, 219, 1, doi: 10.1088/0067-0049/219/1/1
  * Oh et al. (2018) Oh, K., Koss, M., Markwardt, C. B., et al. 2018, ApJS, 235, 4, doi: 10.3847/1538-4365/aaa7fd
  * Oknyansky et al. (2019) Oknyansky, V. L., Winkler, H., Tsygankov, S. S., et al. 2019, MNRAS, 483, 558, doi: 10.1093/mnras/sty3133
  * Panagiotou et al. (2022) Panagiotou, C., Papadakis, I., Kara, E., Kammoun, E., & Dovčiak, M. 2022, ApJ, 935, 93, doi: 10.3847/1538-4357/ac7e4d
  * Parker et al. (2019) Parker, M. L., Schartel, N., Grupe, D., et al. 2019, MNRAS, 483, L88, doi: 10.1093/mnrasl/sly224
  * Pasham et al. (2022) Pasham, D., Fabian, A., Walton, D., et al. 2022, The Astronomer’s Telegram, 15225, 1
  * Reis et al. (2012) Reis, R. C., Miller, J. M., Reynolds, M. T., et al. 2012, Science, 337, 949, doi: 10.1126/science.1223940
  * Remillard et al. (2022) Remillard, R. A., Loewenstein, M., Steiner, J. F., et al. 2022, AJ, 163, 130, doi: 10.3847/1538-3881/ac4ae6
  * Ricci & Trakhtenbrot (2022) Ricci, C., & Trakhtenbrot, B. 2022, arXiv e-prints, arXiv:2211.05132.  https://arxiv.org/abs/2211.05132
  * Ricci et al. (2017) Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, ApJS, 233, 17, doi: 10.3847/1538-4365/aa96ad
  * Rodriguez et al. (2006) Rodriguez, C., Taylor, G. B., Zavala, R. T., et al. 2006, ApJ, 646, 49, doi: 10.1086/504825
  * Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95, doi: 10.1007/s11214-005-5095-4
  * Rovilos et al. (2014) Rovilos, E., Georgantopoulos, I., Akylas, A., et al. 2014, MNRAS, 438, 494, doi: 10.1093/mnras/stt2228
  * Saxton et al. (2012) Saxton, C. J., Soria, R., Wu, K., & Kuin, N. P. M. 2012, MNRAS, 422, 1625, doi: 10.1111/j.1365-2966.2012.20739.x
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Shen & Loeb (2010) Shen, Y., & Loeb, A. 2010, ApJ, 725, 249, doi: 10.1088/0004-637X/725/1/249
  * Sillanpaa et al. (1988) Sillanpaa, A., Haarala, S., Valtonen, M. J., Sundelius, B., & Byrd, G. G. 1988, ApJ, 325, 628, doi: 10.1086/166033
  * Silverman et al. (2020) Silverman, J. D., Tang, S., Lee, K.-G., et al. 2020, ApJ, 899, 154, doi: 10.3847/1538-4357/aba4a3
  * Sironi & Beloborodov (2020) Sironi, L., & Beloborodov, A. M. 2020, ApJ, 899, 52, doi: 10.3847/1538-4357/aba622
  * Sironi et al. (2015) Sironi, L., Petropoulou, M., & Giannios, D. 2015, MNRAS, 450, 183, doi: 10.1093/mnras/stv641
  * Sironi & Spitkovsky (2014) Sironi, L., & Spitkovsky, A. 2014, ApJ, 783, L21, doi: 10.1088/2041-8205/783/1/L21
  * Sobolewska & Papadakis (2009) Sobolewska, M. A., & Papadakis, I. E. 2009, MNRAS, 399, 1597, doi: 10.1111/j.1365-2966.2009.15382.x
  * Tagliaferri et al. (2015) Tagliaferri, G., Ghisellini, G., Perri, M., et al. 2015, ApJ, 807, 167, doi: 10.1088/0004-637X/807/2/167
  * Tchekhovskoy et al. (2014) Tchekhovskoy, A., Metzger, B. D., Giannios, D., & Kelley, L. Z. 2014, MNRAS, 437, 2744, doi: 10.1093/mnras/stt2085
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * Trakhtenbrot et al. (2019) Trakhtenbrot, B., Arcavi, I., MacLeod, C. L., et al. 2019, ApJ, 883, 94, doi: 10.3847/1538-4357/ab39e4
  * Tsalmantza et al. (2011) Tsalmantza, P., Decarli, R., Dotti, M., & Hogg, D. W. 2011, ApJ, 738, 20, doi: 10.1088/0004-637X/738/1/20
  * Valtonen et al. (2008) Valtonen, M. J., Lehto, H. J., Nilsson, K., et al. 2008, Nature, 452, 851, doi: 10.1038/nature06896
  * Vaughan & Uttley (2005) Vaughan, S., & Uttley, P. 2005, MNRAS, 362, 235, doi: 10.1111/j.1365-2966.2005.09296.x
  * Vaughan et al. (2016) Vaughan, S., Uttley, P., Markowitz, A. G., et al. 2016, MNRAS, 461, 3145, doi: 10.1093/mnras/stw1412
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435
  * Véron-Cetty & Véron (2010) Véron-Cetty, M. P., & Véron, P. 2010, A&A, 518, A10, doi: 10.1051/0004-6361/201014188
  * Walton et al. (2017) Walton, D. J., Mooley, K., King, A. L., et al. 2017, ApJ, 839, 110, doi: 10.3847/1538-4357/aa67e8
  * Wang et al. (2022) Wang, C., Luo, B., Brandt, W. N., et al. 2022, ApJ, 936, 95, doi: 10.3847/1538-4357/ac886e
  * Wang et al. (2020) Wang, J., Xu, D. W., & Wei, J. Y. 2020, ApJ, 901, 1, doi: 10.3847/1538-4357/abaa48
  * Wilkins & Gallo (2015) Wilkins, D. R., & Gallo, L. C. 2015, MNRAS, 449, 129, doi: 10.1093/mnras/stv162
  * Wilkins et al. (2022) Wilkins, D. R., Gallo, L. C., Costantini, E., Brandt, W. N., & Blandford, R. D. 2022, MNRAS, 512, 761, doi: 10.1093/mnras/stac416
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016