# Searching for electromagnetic emission in an AGN from the gravitational wave
binary black hole merger candidate S230922g111based on observations made with
the Southern African Large Telescope (SALT)222Some of the data presented
herein were obtained at Keck Observatory, which is a private 501(c)3 non-
profit organization operated as a scientific partnership among the California
Institute of Technology, the University of California, and the National
Aeronautics and Space Administration. The Observatory was made possible by the
generous financial support of the W. M. Keck Foundation.

Tomás Cabrera McWilliams Center for Cosmology and Astrophysics, Department of
Physics, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213
Tomás Cabrera<EMAIL_ADDRESS>Antonella Palmese McWilliams Center for
Cosmology and Astrophysics, Department of Physics, Carnegie Mellon University,
5000 Forbes Avenue, Pittsburgh, PA 15213  Lei Hu McWilliams Center for
Cosmology and Astrophysics, Department of Physics, Carnegie Mellon University,
5000 Forbes Avenue, Pittsburgh, PA 15213  Brendan O’Connor McWilliams Fellow
McWilliams Center for Cosmology and Astrophysics, Department of Physics,
Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213
K.E.Saavik Ford Center for Computational Astrophysics, Flatiron Institute, 162
5th Ave, New York, NY 10010, USA  Department of Astrophysics, American Museum
of Natural History, New York, NY 10024, USA  Department of Science, BMCC, City
University of New York, New York, NY 10007, USA  Barry McKernan Center for
Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY
10010, USA  Department of Astrophysics, American Museum of Natural History,
New York, NY 10024, USA  Department of Science, BMCC, City University of New
York, New York, NY 10007, USA  Igor Andreoni Joint Space-Science Institute,
University of Maryland, College Park, MD 20742, USA  Department of Astronomy,
University of Maryland, College Park, MD 20742, USA  Astrophysics Science
Division, NASA Goddard Space Flight Center, Mail Code 661, Greenbelt, MD
20771, USA  Tomás Ahumada Division of Physics, Mathematics, and Astronomy,
California Institute of Technology, Pasadena, CA 91125, USA  Ariel Amsellem
McWilliams Center for Cosmology and Astrophysics, Department of Physics,
Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Malte
Busmann University Observatory, Faculty of Physics, Ludwig-Maximilians-
Universität München, Scheinerstr. 1, 81679 Munich, Germany  Peter Clark
Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth,
PO1 3FX, UK  Michael W. Coughlin School of Physics and Astronomy, University
of Minnesota, Minneapolis, Minnesota 55455, USA  Ekaterine Dadiani McWilliams
Center for Cosmology and Astrophysics, Department of Physics, Carnegie Mellon
University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Veronica Diaz McWilliams
Center for Cosmology and Astrophysics, Department of Physics, Carnegie Mellon
University, 5000 Forbes Avenue, Pittsburgh, PA 15213  Matthew J. Graham
Division of Physics, Mathematics, and Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA  Daniel Gruen University Observatory,
Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1,
81679 Munich, Germany  Excellence Cluster ORIGINS, Boltzmannstr. 2, 85748
Garching, Germany  Keerthi Kunnumkai McWilliams Center for Cosmology and
Astrophysics, Department of Physics, Carnegie Mellon University, 5000 Forbes
Avenue, Pittsburgh, PA 15213  Jake Postiglione Department of Science, BMCC,
City University of New York, New York, NY 10007, USA  Julian S. Sommer
University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität
München, Scheinerstr. 1, 81679 Munich, Germany  Francisco Valdes NSF National
Optical-Infrared Research Laboratory, 950 N. Cherry Ave., Tucson, AZ 85719,
USA

###### Abstract

We carried out long-term monitoring of the LIGO/Virgo/KAGRA binary black hole
(BBH) merger candidate S230922g in search of electromagnetic emission from the
interaction of the merger remnant with an embedding active galactic nuclei
(AGN) accretion disk. Using a dataset primarily composed of wide-field imaging
from the Dark Energy Camera (DECam) and supplemented by additional photometric
and spectroscopic resources, we searched $\sim 70\%$ of the sky area
probability for transient phenomena, and discovered 6 counterpart candidates.
One especially promising candidate - AT 2023aagj - exhibited temporally
varying asymmetric components in spectral broad line regions, a feature
potentially indicative of an off-center event such as a BBH merger. This
represents the first _live_ search and multiwavelength, photometric, and
spectroscopic monitoring of a GW BBH optical counterpart candidate in the disk
of an AGN.

††facilities: Blanco (DECam), Keck:I (LRIS), Gemini:South (GMOS), SALT (RSS),
Hale (DBSP), WO:2m (3KK)††software:  astropy (Astropy Collaboration et al.,
2013, 2018, 2022), DRAGONS (Labrie et al., 2019), dustmaps (Green, 2018),
gwemopt (Coughlin et al., 2018), healpy (Górski et al., 2005; Zonca et al.,
2019), ppxf (Cappellari, 2023), Fritz SkyPortal (Coughlin et al., 2023),
ligo.skymap (Singer et al., 2016), LPIPE (Perley, 2019), matplotlib (Hunter,
2007), numpy (Harris et al., 2020), pandas (McKinney, 2010), ParSNIP (Boone,
2021), pyraf-dbsp (Bellm & Sesar, 2016), Redrock
(https://github.com/desihub/redrock), SExtractor (Bertin & Arnouts, 1996),
SFFT (Hu et al., 2022), SWarp (Bertin et al., 2002)

##  1 Introduction

As multimessenger astronomy continues to develop, the branch of this field
dedicated to studying compact binary coalescences (CBCs) stands as an
important forerunner as the first to include direct observations of
gravitational waves (GWs) (Abbott et al., 2017a). These events manifest
phenomena at scales presently unobtainable in man-made environments, and hence
are important natural laboratories in which to study extreme physics such as
heavy element nucleosynthesis (Symbalisty & Schramm, 1982) and measure key
physical quantities such as the expansion rate of the universe (Schutz, 1986;
Holz & Hughes, 2005).

This potential was realized following the first successful multimessenger
detection of the binary neutron star merger GW170817/GRB 170817A (Abbott et
al., 2017b), wherein a global electromagnetic astronomy campaign succeeded in
locating a kilonova (KN) (Coulter et al., 2017; Arcavi et al., 2017; Andreoni
et al., 2017; Drout et al., 2017; Nicholl et al., 2017; Evans et al., 2017;
Lipunov et al., 2017; Pian et al., 2017; Smartt et al., 2017; Utsumi et al.,
2017; Valenti et al., 2017; Kasliwal et al., 2017; Troja et al., 2017;
Margutti et al., 2017; Shappee et al., 2017; Soares-Santos et al., 2017;
Chornock et al., 2017; Tanvir et al., 2017) incident with the initially
identified GW and gamma ray burst (GRB) events (Abbott et al., 2017a;
Savchenko et al., 2017; Goldstein et al., 2017). This single event by itself
was enough to result in hundreds of analyses, including constraints on neutron
star (NS) physics (e.g. Abbott et al., 2018) and the Hubble constant $H_{\rm
0}$ (Abbott et al., 2017c). In the case of $H_{\rm 0}$, it is predicted that
$\mathcal{O}(100)$ binary neutron star (BNS) events with electromagnetic (EM)
counterparts are required to make a $\sim$few% measurement of the parameter
(Chen et al., 2018), although joint multimessenger constraints on the binary
viewing angle can further improve precision by a factor of a few (e.g.
Hotokezaka et al. 2019; Palmese et al. 2024).

While BNS mergers are the only class of GW event with a confirmed EM
counterpart, it is also predicted that binary black hole (BBH) mergers can
produce counterparts in certain circumstances. The majority of proposed
counterpart mechanisms involve the interaction of the binary or merger remnant
with a coexistent medium, usually the gaseous accretion disk of an active
galactic nucleus (AGN), whether through accretion (Bartos et al., 2017b), ram-
pressure stripping of gas about the remnant or jetted Bondi accretion
(McKernan et al., 2019), or breakout emission from a post-merger jet (Tagawa
et al., 2023, 2024). A cartoon of the AGN-associated EM counterpart mechanism
is show in Figure 1.

Figure 1:  Multi-panel schematic showing the mechanism believed to underpin
luminous EM counterparts to BBH mergers in AGN disks. In panel 1, the pre-
merger BBH accretes from mini-disks within its Hill sphere in the AGN disk
midplane and blows a cocoon within the disk via feedback. In panel 2, the
merger happens, forming a highly spinning BH (dimensionless spin parameter
$a\sim 0.7$ typically). A jet is presumed to form at this stage (although it
has yet to be established whether such a jet can persist for long, or whether
it is choked off by high mass accretion). Mass and spin asymmetries in the
progenitor black holes lead to a kick at merger (depicted by the arrow in
panel 2). Panels 3 and 4 show the development of BHL accretion as the newly
merged BH exits its original Hill sphere into the rest of the AGN disk,
powering a luminous transient. In panel 5 the BH emerges from the AGN disk,
dragging disk gas with it. In panel 6 the EM emission fades as the disk
material is consumed and the BH continues on an inclined orbit around the
SMBH, and will re-enter the AGN disk on half the orbital timescale.

Counterpart candidates for BBH mergers from previous LIGO/Virgo/KAGRA (LVK)
observing runs have been proposed (Graham et al., 2020, 2023), although
confirmation of a counterpart has remained challenging. There are several
reasons for this difficulty. First, even if we know that a given merger
happens in an AGN, in principle there is only a $\mathcal{O}(1/4)$ chance that
we could detect an EM counterpart: our view of around half of AGN (the Type
2s) is obscured, and an unfortunate kick-direction out the opposite side of a
Type 1 AGN could still obscure any emission (Graham et al., 2023). Second,
even if a flare emerges from an AGN, it must be discerned in the presence of
other AGN variability. This makes our task of searching for counterparts more
challenging in brighter AGN, and our search is biased against less luminous
counterparts. Third, possible flare parameters are weakly constrained, both
because AGN disk properties are uncertain to orders of magnitude, and because
the properties of the post-merger flare depend on these uncertain parameters.

One possible means to confirm a counterpart is through optical spectroscopy:
because the BBH merger occurs off-center in the AGN disk, the resulting flare
could unevenly illuminate the broad-line region (BLR) of the AGN, resulting in
asymmetric broad-line features in the optical spectrum (McKernan et al.,
2019). Synchronized evolution of a transient with such spectral features can
be a key piece of evidence in favor of the classification of a transient as a
EM counterpart to a BBH merger. At the very least, an evolving asymmetric line
profile suggests an off-center flaring event in the AGN, ruling out most
sources of AGN variability near the SMBH. Flare energetics and lightcurve
profiles can allow us to rule out e.g. embedded supernovae, leaving few
candidates for a sufficiently energetic off-center AGN flare, including a BBH
merger.

Certainly, the search for EM counterparts to BBH mergers has unique challenges
versus the search for NS merger counterparts (intrinsic AGN variability, delay
time uncertainties, etc.); however, as the historical rate of detection of BBH
mergers has been $\sim\mathcal{O}(100)\times$ greater than that for NS
mergers, there are many more opportunities to search for a counterpart to the
former kind of event. Because multimessenger observations of BBH mergers are
also useful in making cosmological measurements (Palmese et al., 2021; Bom &
Palmese, 2023; Alves et al., 2024), the pursuit of EM counterparts to BBH
mergers has the potential to significantly contribute to the multimessenger
observations required to make the first 2% measurement of the Hubble constant
with standard sirens, which the level of precision needed to help us
understand the Hubble tension (Chen et al., 2018).

Currently, GW follow-up of distant events is best enabled with technology
capable of meeting the colloquial requirements of “wide, fast, deep”, as to
ensure the rapid and thorough coverage of event volumes. The Dark Energy
Camera (DECam) (Flaugher et al., 2015) on the 4m Victor M. Blanco Telescope st
the Cerro Tololo Inter-American Observatory (CTIO) is one of the premiere
instruments that can address these requirements. DECam already has a
respectable history in GW follow-up, having been used for this purpose from
the first GW event (Soares-Santos et al., 2016; Annis et al., 2016) and many
events since then (Goldstein et al., 2019; Andreoni et al., 2019; Herner et
al., 2020; Andreoni et al., 2020; Morgan et al., 2020; Garcia et al., 2020),
and importantly being one of the first instruments to detect the counterpart
KN to GW170817 (Soares-Santos et al., 2017). The new survey program
Gravitational Wave MultiMessenger Astronomy DECam Survey (GW-MMADS) (PIs:
Andreoni & Palmese) is designed to find EM counterparts to BNS, NS BH, and BBH
mergers via rapid DECam follow-up of GW events during the fourth gravitational
wave observing run (O4).

In this work, we present results from our follow up of the GW event S230922g.
S230922g is a LVK BBH merger candidate detected at 2023-09-22 02:03:44.886 UTC
(Ligo Scientific Collaboration et al., 2023a, b). The 90% credible region of
324 deg2 is one of the smallest areas of O4 a (the first part of O4, completed
in January 2024); this, combined with the localization being visible from
Chile in the months following the trigger, made an effective follow-up with
DECam possible in pursuit of potential AGN-linked counterparts like those
mentioned previously. In this paper we present our long-term monitoring of
S230922g; the paper is organized as follows: §2 presents our program’s generic
strategy for GW follow-up and the data collected, highlighting components
especially relevant for S230922g; §3 describes the methods we used to distill
our population of detected transients into a shortlist of counterpart
candidates; §4 describes our most favored candidate, and to a lesser extent
additional candidates of interest; in §5 we estimate parameters for our most
favored counterpart candidate, finding it well within the confines of existing
BBH counterpart theory; in §6 we summarize our findings and highlight
considerations relevant for future efforts. Throughout this work we use a flat
$\Lambda$CDM model with $H_{0}=70$ km/s/Mpc and matter density of
$\Omega_{m}=0.3$,

##  2 Data

###  2.1 GW data

S230922g is an event of interest because it was detected with high
significance (False Alarm Rate - FAR - 1 per $1.6\times 10^{-16}$ years from
gstlal; Messick et al. 2017; Sachdev et al. 2019) by both LIGO Livingston and
LIGO Hanford, it is well localized compared to the O4a population, and it has
$\sim$100% probability of being a BBH. It was also identified with high
significance by the Burst CWB (Klimenko et al., 2016) search pipeline,
potentially indicating a loud, short burst, as one may expect for a massive
BBH spending a short fraction of the late inspiral phase in the LVK band. The
luminosity distance of the event, marginalized over the entire sky, is
$d_{L}=1491\pm 443$ Mpc. These quantities are reported with the skymap from
the Bilby (Ashton et al., 2019) reduction (Ligo Scientific Collaboration et
al., 2023c); we use this skymap for our work.

Higher mass ($M_{\rm tot}\gtrsim 50~{}M_{\odot}$) BBH mergers are more likely
to have originated from the disks of AGNs (e.g. Gayathri et al. 2021), and so
we also inform our decision to trigger BBH follow up by estimating the total
mass of the binary. We follow a similar calculation to that in Graham et al.
(2020), and consider that $A_{90}\propto~{}{\rm SNR}^{-2}$ (Berry et al.,
2015), as well as ${\rm SNR}\propto\mathcal{M}_{c}^{5/6}/d_{L}$ (Finn &
Chernoff, 1993), where $\mathcal{M}_{c}$ is the chirp mass of the binary. For
this event, we assume a 1.4-1.4 $M_{\odot}$ binary neutron star merger
detection horizon of 150 and 152 Mpc for Hanford and Livingston, respectively,
an equal mass system, and the luminosity distance marginalized over the sky;
we derive a total rest frame mass of $M_{\rm tot}\sim 90~{}M_{\odot}$ for
S230922g. Note that this estimate is highly uncertain with an error bar of at
least a factor of 2.

###  2.2 DECam observations

Our team was notified of S230922g through a General Coordinates Network (GCN)
listener333Adapted from https://github.com/scimma/slackbot. that sent a digest
of the event to our Slack workspace.444Our team also uses the Fritz science
data platform (Coughlin et al., 2023) for notifications of GW events via phone
call, but we reserve this kind of notification for exceptionally time-
sensitive events such as BNS and NS BH mergers. Our first response was to
generate an initial observing strategy with gwemopt (Coughlin et al., 2018) to
asses the probability coverage possible with a single night of DECam
observations. gwemopt selects exposure pointings from a preset sky tiling.
Because our image subtraction pipeline uses archival DECam exposures as
templates, we implemented a custom tiling based on images calibrated by the
NSF National Optical-Infrared Research Laboratory (NOIRLab) DECam Community
Pipeline (CP) (Valdes et al., 2014) and made available through the NOIRLab
Astro Data Archive555https://astroarchive.noirlab.edu/ (McManus & Olsen, 2021)
to ensure our observing plans include pointings with usable templates. During
each observing night, the CP calibrates the response and astrometry of each
exposure shortly after it is taken which is then available to the image
subtraction pipeline.

We determined that 70% of the GW localization area could be covered with 60 s
and 80 s DECam exposures in $g$\- and $i$-bands, respectively, with $\sim$4.7
hours of telescope time. Our DECam tiling is shown in Figure 2, overlaid on
the LVK Bibly skymap (Ligo Scientific Collaboration et al., 2023b) for
S230922g. When this plan was finalized, we communicated our pointings to the
larger astronomical community though the GW Treasure
Map666https://treasuremap.space/ (Wyatt et al., 2020), and we updated
information appropriately as we executed our plan.

Figure 2:  Our observation plan for follow up of S230922g. The LVK skymap for
the event is plotted on the lower layers of the figure, and the 50% (90%)
credible interval regions are outlined in the solid (dashed) line. The
candidates composing our shortlist are shown as black X’s, with our favored
candidate indicated with a green star.

Leading BBH counterpart models predict delays of 1-3 weeks after the merger
before the manifestation of any EM counterpart, due to the need for the signal
to first emerge from the depths of the AGN disk. Even so, it is valuable to
initiate EM follow up for these kinds of events within a few days of the GW
trigger, as high-cadence archival data are generally not available in the
localization area of a given GW event and a baseline of sources in the
localization area can be useful when vetting candidates. Accordingly, we
decided to trigger our Target of Opportunity (ToO) program for the night
following the event (the evening of September 22, 2023) in order to establish
a baseline of activity at the GW merger time and to facilitate the
identification of novel phenomena in subsequent observing epochs.

On the first night we were only able to complete $\sim$1 hour of observations
before inclement weather began, but we were able to observe our full plan in
its entirety the following night (September 23). We publicly reported
transients in the area through GCN within 24 hours of these observations
(Cabrera et al., 2023). Further observations were conducted 11, 33, 41, and 70
days after the LVK trigger (October 3, October 25, November 2, and December
1). As we processed the data throughout the campaign, we took additional
observations of the most interesting candidates, adding up to eight epochs
altogether. We submitted all transient sources discovered through this
campaign to the Fritz SkyPortal (van der Walt et al., 2019; Coughlin et al.,
2023). After composing our final candidate shortlist (as detailed in §3), we
reported the respective sources to the Transient Name Server (TNS).

###  2.3 Spectroscopic observations

BBH counterpart models have some degeneracy with those for stochastic AGN
variability and other transient events like tidal distruption events (TDEs),
but spectroscopic data has been predicted to serve as a key discriminator
among these phenomena: specifically, the off-center location of a BBH merger
in the accretion disk of the AGN as opposed to the central location of TDEs or
accretion-based phenomena is expected to induce asymmetry in the broad lines
of AGN spectra as the transient event “washes” over the BLR in an asymmetric
manner (McKernan et al., 2019). The detection of a BLR asymmetry that evolves
in concert with the light curve of a transient is a smoking gun that locates
the transient in an off-center position in the accretion disk expected for BBH
mergers, and the coincidence of such an event with a GW merger is a strong
evidence standard astronomical techniques can provide linking BBH mergers and
AGN flares.

In pursuit of such evidence, we triggered several spectroscopic resources for
additional follow-up of transients whose light curves demonstrated proposed
counterpart features (see §3.1). A total of five spectra were collected during
our follow-up campaign: two with Gemini GMOS, and one each with Keck LRIS,
SALT RSS, and P200 DBSP. We discuss the methodology for each spectrum in the
context of the respective candidate in §4.2.

##  3 Methods

###  3.1 Automated vetting

Detection-based cuts
---
Initial photometry detections | 25,392,140 | 1.00
Data quality masking | 14,053,259 | 0.553
Real/bogus score $\geq$ 0.7 | 1,501,800 | 0.059
Source-based cuts
Initial transients | 233,313 | 1.00
Remove LS variable stars | 140,230 | 0.601
Remove MPC objects | 18,528 | 0.079
$\geq$1 real/bogus score $\geq$0.9 | 17,515 | 0.075
$\geq$2 detections with $\Delta t\geq 30$ min | 3558 | 0.015
LS galaxy sep $\leq$ 1” | 2388 | 0.010
Additional vetting | 6 | 2.6e-5

We analyze our data with our difference photometry pipeline described in Hu et
al. (2024). We perform image subtraction with Saccadic Fast Fourier Transform
(SFFT) (Hu et al., 2022), a scalable image subtraction algorithm and GPU-
enabled implementation of the same that produces difference images up to an
order of magnitude faster than widely used tools with comparable accuracy.
Aperture photometry is conducted on the resulting difference images using
SExtractor (Bertin & Arnouts, 1996) with fluxes calibrated to the DESI Legacy
Survey (LS) source catalog and corrected for extinction using the dustmaps
package (Green, 2018) and the $E(B-V)_{\rm SFD}$ coefficients from Schlafly &
Finkbeiner 2011. Our image differencing pipeline detected over 25 million
transient features in our campaign data, which we distilled into a tractable
list of astrophysical transients through a series of cuts. Table
LABEL:tab:vetting summarizes the vetting steps we applied, which are described
in this section.

We first remove any detections on difference images contaminated by bad pixels
recorded in the DECam CP data quality mask products. Surviving features were
then scored using a Convolutional Neural Network (CNN) trained on archival
DECam data products from our pipeline to perform real/bogus classification.
Possible CNN scores range from 0 to 1, with higher scores associated with more
realistic astrophysical objects; at this step we kept any detections with
scores greater than 0.7. Roughly 6% of the initial list of detections passed
these steps to form our list of real astrophysical detections.

These detections were then cross-matched among themselves to identify features
present in data from multiple epochs and filters; the resulting list of
transients consisted of 233,313 objects. This list was subsequently
crossmatched with the LS DR10.1 catalog (Dey et al., 2019) and the Minor
Planet Center (MPC) service, and any transients respectively matched to star-
like sources (defined as a source assigned a LS morphological type of “PSF”)
or minor planets were removed from our search. The LS crossmatch was also used
to sort the sources into three categories based on proximity to the nearest
galaxy-type777We consider all LS sources with PSF types other than “PSF”
(stellar) and “DUP” (extended source components) as “galaxy-type” sources. LS
source: sources within 0.3” of a galaxy-type source were labeled as “A”
sources, sources greater than 1.0” away from the nearest such source were
labeled as “T” sources, sources falling in the middle ground between these two
categories were labeled as “C” sources.

Two final cuts were then applied to limit our search to the most realistic
persistent sources, requiring sources to have least two distinct observations
separated by $>$30 minutes, and at least one high-fidelity detection
(real/bogus CNN score $\geq$0.9). Of the remaining 3558 sources, 2388 of them
were within 1” of a galaxy-type LS object (“A” and “C” sources); this final
group composed our automated list of candidates based on our photometric DECam
data. See Table LABEL:tab:vetting for a summary of cuts applied and resulting
numbers of candidates at each cut.

###  3.2 Additional vetting

Our 2388 automatically-identified sources were then examined by eye to
identify candidates for further study. Sources flagged for further
investigation included those that had brightened since the time of the GW
event and did not exhibit early reddening (this latter criteria was motivated
by the expectation that the luminous counterpart emerges from the depths of
the accretion disk, see Figure 1, and we do not expect it to redden as the
optical depth decreases). Other considerations included whether there was a
perceptible delay time from the GW event time (reflective of an EM counterpart
needing time to escape the embedding accretion disk before being observed) and
whether sources appeared similar to different known transients such as
supernovae. Image stamps were also examined for each source to exclude any
artifacts and the like that survived the CNN cut.

The resulting list was trimmed to remove any transients over 2$\sigma$ away
from the GW distance posterior, as such information was available. To
determine distances to each transient, the LS galaxy-type object matched to
each transient was cross-matched with several galaxy catalogs (with a search
radius of 1”) in search of a redshift measurement. Galaxy catalogs used
included the NASA/IPAC Extragalactic Database (NED) Local Volume Sample (Cook
et al., 2023) 888https://ned.ipac.caltech.edu/, the Dark Energy Spectroscopic
Instrument (DESI) galaxy catalog (DESI Collaboration et al., 2023), and Quaia
(Storey-Fisher et al., 2024); LS DR10.1 photometric redshifts were also
included. If the host galaxy matched in multiple catalogs, then the best
redshift measurement was used, preferring spectroscopic redshifts over Quaia
“spectrophotometric” redshifts (Storey-Fisher et al., 2024), and the latter
over photometric redshifts. Comparing the available redshift measurements to
the distance posteriors (taken from the skymap
HEALPix999https://healpix.sourceforge.io/ tile containing the respective
transient) informed the elimination of sources beyond the 2$\sigma$ threshold
cut, using the larger of the two uncertainties between the GW and galaxy
catalog measurements. Note that transients whose host galaxy lacked any
redshift measurement were not subject to this cut.

Candidates in our shortlist were photometrically classified with
Parameterization of SuperNova Intrinsic Properties (ParSNIP) (Boone, 2021),
assuming a model trained on the PLAsTiCC (Kessler et al., 2019) simulations;
because ParSNIP uses redshift as an input parameter, this classification was
limited to only those sources with a redshift measurement. When classifying a
light curve, ParSNIP assigns a probability to each of a set list of transient
classes so that the total probability sums to 1, with each probability
reflecting the relative likeness of the light curve to each class. We note
that ParSNIP does not contain a catch-all class such as “Other” for use when
classifying transients of unfamiliar phenomenology: that is, if a transient is
strongly dissimilar to all but one of the classes, it could receive a high
classification probability for that one class, even if it not a strong match
in and of itself. We consider the “TDE” ParSNIP class as the one most similar
to the AGN flares of interest (all other classes concern some kind of SN,
except for the KN class, whose typical timescales are considerably shorter
than those we are interested in). Accordingly, in our use case we broadly
interpret transients classified as “TDE” as those that are not identifiable as
SN-like events, and we refer to the “TDE” class as “Non-SN” to better reflect
this perspective. We find that most of our classified transients receive this
classification.

The final shortlist of 23 transients surviving these cuts is shown in Table
LABEL:tab:candidates. We separate several subsets from this list to identify
transients that we exclude from our search via asynchronous analysis after our
follow-up campaign; these subsets appear in labeled sections of the table. 6
of the transients in our shortlist showed significant brightening, but did not
peak during the time they were observed (so potentially consistent with
longer-term AGN variability), and so we are unable to consider their full
nature with our present dataset. 9 transients, while initially interesting for
further monitoring, exhibited reddening in later epochs, and so are excluded
through disagreement with our assumed counterpart model that predicts a
signature that becomes more blue with time. Finally, 2 transients were
excluded as counterpart candidates through spectroscopic classification (see
§4.2.2 and §4.2.3). The remaining 6 transients are those that cannot be
excluded as counterparts to S230922g, and are listed at the top of the table.
Table LABEL:tab:candidates includes host redshift information (where
available), along with GW skymap localization information: the 1$\sigma$
distance posterior expressed as redshift under the assumed cosmology, and the
2D and 3D confidence interval (CI) as calculated with the
ligo.skymap.postprocess.crossmatch.crossmatch
routine101010https://lscsoft.docs.ligo.org/ligo.skymap/postprocess/crossmatch.html#ligo.skymap.postprocess.crossmatch.crossmatch
(Singer et al., 2016). The transients are ranked by ascending 2D CI, s.t. the
objects in the highest probability regions are listed first. The highest-
probability ParSNIP class and the associated probability are listed in the
last two columns of the table, where available.

AT 2023adwb | - | - | $1360\pm 430$ | 0.136 | - | - | -
---|---|---|---|---|---|---|---
AT 2023uho | - | - | $1250\pm 552$ | 0.484 | - | - | -
AT 2023adwp | - | - | $1274\pm 433$ | 0.652 | - | - | -
AT 2023aagj | 0.184 | Specz (this work) | $1478\pm 425$ | 0.754 | 0.829 | Non-SN | 0.944
AT 2023adwt | - | - | $1429\pm 413$ | 0.756 | - | - | -
AT 2023uea | $0.195\pm 0.084$ | Quaia SPz | $1459\pm 426$ | 0.790 | 0.814 | Non-SN | 0.991
Transients without peak
AT 2023adio | $0.212\pm 0.097$ | Quaia SPz | $1410\pm 446$ | 0.242 | 0.363 | SNII | 0.901
AT 2023adwd | - | - | $1328\pm 451$ | 0.285 | - | - | -
AT 2023adwe | $1.528\pm 0.671$ | LS photz | $1293\pm 463$ | 0.285 | - | Non-SN | 0.927
AT 2023adwn | $0.349\pm 0.078$ | LS photz | $875\pm 486$ | 0.617 | 0.810 | Non-SN | 0.676
AT 2023adwq | $0.391\pm 0.112$ | Quaia SPz | $1545\pm 454$ | 0.683 | 0.784 | Non-SN | 0.968
AT 2023adfo | - | - | $1471\pm 422$ | 0.729 | - | - | -
Reddening transients
AT 2023adwc | 0.216 | NED specz | $1228\pm 445$ | 0.143 | 0.080 | Non-SN | 0.996
AT 2023adwf | 0.180 | NED specz | $1148\pm 417$ | 0.300 | 0.155 | Non-SN | 0.985
AT 2023adwj | - | - | $1345\pm 457$ | 0.429 | - | - | -
AT 2023uec | - | - | $1127\pm 469$ | 0.459 | - | - | -
AT 2023uos | - | - | $1448\pm 458$ | 0.475 | - | - | -
AT 2023adwl | - | - | $1051\pm 441$ | 0.583 | - | - | -
AT 2023unj | - | - | $1127\pm 439$ | 0.642 | - | - | -
AT 2023unl | 0.248 | Specz (this work) | $1535\pm 444$ | 0.677 | 0.652 | - | -
AT 2023adws | 0.178 | NED specz | $1451\pm 419$ | 0.750 | 0.827 | Non-SN | 0.928
Excluded via spectra
AT 2023aden | - | - | $1154\pm 451$ | 0.658 | - | - | -
AT 2023uab | 0.128 | Specz (this work) | $919\pm 479$ | 0.785 | 0.554 | SNIa | 0.814

##  4 Candidate counterparts

In this section we describe candidate counterparts that received particular
scrutiny during our campaign, including one we identify as the most likely to
be a counterpart to S230922g. While we generally refer to candidates by their
TNS name, we include the internal name for each of our transients as well, as
these names were used for initial reports. Internal names are assembled from
the date of discovery and RA/dec of the source; for example, the internal name
of C202309242206400m275139 refers to a transient first detected on September
24, 2023 with an RA/dec of 22h 06m 40.0s, -27∘ 51m 39s (the letter between the
RA and dec can be either “p” or “m”, denoting a positive or negative
declination, respectively). The leading character in the object names
indicates the proximity of the transient to the nearest galaxy-type LS source,
in the same convention explained in §3.1 (“A” for nuclear sources, “T” for
non-nuclear sources, and “C” for marginally nuclear sources). The internal
names for all transients in our shortlist can be found in the subfigure titles
of Figure 5.

###  4.1 AT 2023aagj (C2023092422006400m275139)

We consider AT 2023aagj to be the transient from our sample most likely to be
an optical counterpart to S230922g. This transient is distinguished by a
$\sim$1 magnitude brightening over a period of a month, remaining blue in
color throughout its evolution. It is within 0.5” of the host galaxy centroid,
and as such may be associated with activity in the host nucleus. Difference
photometry, spectra, and sample stamps for this candidate are visible in
Figure 3. ParSNIP classifies AT 2023aagj as a “TDE”/“Non-SN”-type object with
a probability of 94.4%, which in this context distinguishes the transient from
SNe and the like.

Figure 3:  Observational data for AT 2023aagj. Upper left: Difference
photometry for the transient. $g(i)$-band data are shown in green (yellow).
The colored vertical lines in the light curve plot indicate the times
additional data (spectra and X-ray) were taken. Upper right: Stamps from the
template, science, and difference images for a sample epoch; the epoch from
which the sample stamps are taken is outlined in the photometry plot with a
gray square. Center left: WISE photometry for the host of AT 2023aagj. The
span of the difference photometry plot is shaded in gray. Bottom: Spectra
taken for the transient. Telluric regions are shaded in green. Several known
emission lines (shifted according to the measured redshift of the host) are
labeled and shown as dashed gray lines.

Two spectra of AT 2023aagj were taken before the transient faded. The initial
Keck Low-Resolution Imaging Spectrograph (LRIS) spectrum, taken on December 7
(PI: Kasliwal, PID: C360), and reduced with the standard LPIPE routine
(Perley, 2019), revealed ionized gas emission lines and broad line features
for the H$\alpha$ and MgII regions, indicating a possible AGN host. The
redshift of the host was calculated to be 0.184, placing the transient within
the 83% credible volume for the S230922g skymap and consistent at the
1.4$\sigma$ level with the luminosity distance posterior of the GW event. An
additional spectrum of the object was taken with Gemini North Gemini Multi-
Object Spectrograph (GMOS) on December 25, 2023 (PI: Cabrera, PID:
GN-2023B-DD-109), using a B480 grating with 2x900s dithered exposures and a
central wavelength of 6250 Å. The data were reduced with Data Reduction for
Astronomy from Gemini Observatory North and South (DRAGONS) (Labrie et al.,
2019), and flexure corrections were applied by re-running the reduction
software without sky subtraction, and calculating a scalar offset by comparing
to a sky line catalog (Rousselot et al., 2000). An interesting feature
observed in the Keck spectrum is some degree of asymmetry (skewed towards the
redder wavelengths) in both the H$\alpha$ and MgII broad lines, as seen in
Figure 4; however, the H$\alpha$ asymmetry does not appear in the later Gemini
North-GMOS spectrum. The latter spectrum does not cover the MgII feature, and
so no comparison between the two times is possible in this regime. Given that
there is no similarly-resolved archival spectrum for this source, we are not
able to distinguish transient spectral features from persistent host features
at this time.

Figure 4:  pPXF model and data for the two AT 2023aagj spectra around the two
broad lines of interest. The bottom panels show the residuals between the
model and data. An asymmetry redward of the line appears to be present in the
Keck spectrum, taken closer to peak, but disappears in the later Gemini
spectrum.

We fit the Gemini spectrum using the Penalized Pixel-Fitting software (ppxf;
Cappellari 2023) in order to study the possible presence of an AGN. We fit the
spectrum with stellar templates from the E-MILES SSP library (Vazdekis et al.,
2016) along with a number of narrow and broad Gaussian peaks corresponding to
various emission and absorption lines. Using these measured gas line
strengths, we place this object in the Baldwin–Phillips–Terlevich (BPT)
diagram, finding that it lands in the “composite” region (according to the
delineations of Kewley et al. 2001; Kauffmann et al. 2003; Kewley et al.
2006). This indicates that it is likely that the observed lines are due to a
combination of star formation and AGN activity. We further measure the mass of
the AGN using the fitted H$\alpha$ velocity dispersion, which we find to be
$2282\pm 41$ km/s. The resulting BH mass according to the method presented in
(Greene & Ho, 2005) is $M_{\rm SMBH}\sim 2\times 10^{7}~{}M_{\odot}$. Note
that this measurement may be inaccurate if the transient is significantly
contributing to the BLR emission.

Wide-field Infrared Survey Explorer (WISE) (Wright et al., 2010a) infrared
data for the host of this event is shown in Figure 3. Following Clark et al.
(2024), we retrieved this data from the NASA/IPAC infrared science archive
(IRSA)111111https://irsa.ipac.caltech.edu/, utilizing data from both the
ALLWISE (Wright et al., 2010b) and NEOWISE Reactivation Releases (NEOWISE-R;
Mainzer et al., 2011, 2014). Individual observations are combined following
quality control cuts to obtain one mean magnitude per filter for each
observation period, providing a $\sim$ 6 month cadence. Notably, there are
variations in the host system of similar amplitude to our detected transient,
but it must be recognized that such variations occur on much longer ($\sim
10\times$) timescales than the timescale of AT 2023aagj. We compare this
behaviour with publicly available PanSTARRS (Flewelling et al., 2020) and
DECam images at this location, which, although sparse, do not show any
significant change in flux. We observe that a significant jump of $\sim 0.3$
mag in $g$ band 5” aperture magnitude occurred according to the SkyMapper data
release (Onken et al., 2024) between MJD $\sim$57650 and $\sim$57850\. It is
possible that an associated optical signature existed in that time period, but
was missed by the relatively sparse nature of the optical observations, and
the infrared (IR) component we observed was delayed from the optical via dust
reprocessing. We do not observe a rebrightening of the source in the WISE
bands over the timescale of the optical transient, but there appears to be a
$\sim 0.1$ increase in magnitude in December 2023, potentially indicating the
start of a reprocessed, delayed emission from the event.

We also produce a WISE color-color plot for W1-W2 versus W2-W3, although this
is not shown as W3 is only available in one epoch, which may not be
representative of the entire color evolution of the object over the past
years. By comparison with populations of extragalactic objects (as defined in
Wright et al. 2010a) we find that while the host most clearly lies in the WISE
population regions associated with spiral and luminous infrared galaxies, it
is near the boundary into the Seyfert and quasar parameter spaces, and has
been observed to temporarily visit these regions.

We carried out a ToO observation of AT 2023aagj with the Neil Gehrels Swift
Observatory (Gehrels et al., 2004) X-ray Telescope (XRT; Burrows et al.,
2005). Our observation began on 2023-12-25 at 06:02 UT with a total exposure
of 3265 s in Photon Counting (PC) mode. We used the Swift-XRT automated tools
(Evans et al., 2009, 2023) to analyze our data and all previous data covering
the source position from May 2006 and April 2021 (totaling 7.8 ks exposure in
PC mode, including the latest data). We measure a source position of RA, DEC
(J2000) = $22^{h}06^{m}40^{s}.26$, $-27^{\circ}51\arcmin 41.8\arcsec$ with an
uncertainty radius of $4.7\arcsec$ (90% confidence). In our latest observation
the source displays a soft spectrum with all photon counts below 5 keV,
similar to the previous observations. The time-averaged spectra were fit with
an absorbed powerlaw model using XSPEC v12.14.0 (Arnaud, 1996) within HEASoft
v6.32. The inferred photon index is relatively unconstrained but favors
$\Gamma\approx 2.5$ with hydrogen column density $N_{H}\approx 2\times
10^{20}$ cm-2. The unabsorbed flux ($0.3-10$ keV) is
$F_{X}=(3.2^{+8.8}_{-1.0})\times 10^{-13}$ erg cm-2 s-1 in our latest
observation. We find a factor of $\sim$ $4$ increase in flux between the 2006
and 2021 observations, but no significant change (a factor of $\sim$ $1.2\pm
0.5$) between 2021 and 2023 (due also to the larger errors on the flux). We
therefore conclude that there is no ongoing X-ray outburst from AT 2023aagj,
although we note that these observations were taken $\sim 2$ months after the
optical peak at a time when any additional flux from the transient may have
been too faint to be detected. It is interesting to note how the WISE flare
peaked around 2018, notably coincident with the significant change in X-ray
flux. We speculate this may be related to a changing look AGN event, although
this is hard to confirm due to the lack of an older spectrum.

###  4.2 Other candidates of scrutiny

We took spectra for three additional candidates beyond AT 2023aagj during our
follow-up campaign. The spectrum for one transient (AT 2023unl) provides
inconclusive evidence for the nature of the source; however, the source was
observed to redden in later epochs, and so is disfavored by our assumed
counterpart model. The spectra for the other two transients (AT 2023aden and
AT 2023uab) contain SN-like features; both of these transients are
subsequently excluded as possible counterparts to S230922g. The three spectra
are shown in Figure 9 (Appendix A).

####  4.2.1 AT 2023unl (C202310042207549m253435)

The spectrum for AT 2023unl was taken by the Robert Stobie Spectrograph (RSS)
(Burgh et al., 2003; Kobulnicky et al., 2003) on the South African Large
Telescope (SALT) (Buckley et al., 2006). The RSSMOSPipeline (Hilton et al.,
2018)121212https://github.com/mattyowl/RSSMOSPipeline was used for the
reduction of SALT-RSS data. The pipeline automatically identifies and extracts
one dimensional spectra using the data products delivered by the SALT
observation team. Using the known redshift of the observed object and the
reduced spectra, PYQSOFit (Guo et al.,
2018)131313https://github.com/legolason/PyQSOFit was then used to fit and
identify possible emission lines in the data. The signal-to-noise ratio (SNR)
for the spectrum is relatively low, with few strong features. The present
H$\beta$ line does not demonstrate any unusual characteristics such as
asymmetry. A possible broad line feature is present around an observed
wavelength of 5525 Å, albeit a chip gap appears to have ended up at the peak
of the feature, and so it is uncertain how much of the observed variation is
due to imperfect data reduction; that said, there is no similar feature around
the second chip gap near 6575 Å. We conclude there is no conclusive evidence
for this transient being in an AGN.

####  4.2.2 AT 2023aden (A202310262246341m291842)

A spectrum was taken for this transient using Gemini North GMOS on November
16, 2023 (PI: Cabrera, PID: GN-2023B-DD-103), and was acquired and reduced in
a similar manner as the Gemini spectrum for AT 2023aagj, with a central
wavelength of 7600 Å. The spectrum is blue and has strong absorption features
characteristic of SNe; subsequently, we do not consider this transient as a
likely counterpart to S230922g.

####  4.2.3 AT 2023uab (C202309242248405m134956)

We used the Double Spectrograph (DBSP) mounted at the 200-inch Hale Telescope
at Palomar Observatory. For this source, we used a 1.5 arcsec slitmasks, a D55
dichroic, a blue grating of 600/4000 and red grating of 316/7500. The data
were reduced using a custom PyRAF DBSP reduction pipeline (Bellm & Sesar,
2016). The spectrum exhibits clear H$\alpha$ and H$\beta$ features which
enable a redshift measurement of 0.128 (computed with
Redrock141414https://github.com/desihub/redrock; Bailey et al., in prep.). The
continuum of the spectrum peaks around $\sim$5000 Å, with several broad
absorption features present around the peak. This, in concert with ParSNIP
classifying the object as a type Ia SN, leads us to exclude this object as a
possible counterpart.

##  5 Discussion

Taking our favored counterpart candidate of AT 2023aagj, we work under the
assumption that the transient is related to S230922g and try to infer whether
the the properties of the flare may be consistent with expectations for the
counterpart signal, mechanism, and host.

What can be noted is that the persistent blue color is unlike typical
supernovae, and structure function (SF) arguments as in Palmese et al. (2021)
show that this transient is distinguished from typical AGN activity at the
5.9$\sigma$ level (corresponding to a one-sided probability $p_{\rm
flare}\sim\mathcal{O}(10^{-16})$). Taking the 90% CI volume of the event as
$3.73\times 10^{8}$ Mpc3 (as calculated with ligo.skymap) and a fiducial AGN
number density of $10^{-4.75}$ Mpc-3 (Bartos et al., 2017a; Greene & Ho,
2007), we estimate a count of 6600 AGN in the volume. Alternatively,
crossmatching the Milliquas (Flesch, 2023) quasar catalog with the skymap
(also with ligo.skymap) yields a 3D match of 826 objects. Given the discovery
window of our follow-up campaign $t_{\rm window}=92$ days (calculated as the
difference in time between the second epoch and the last), we calculate the
probability of a chance discovery of a similar flare as

$p_{cc}=n_{\rm AGN}\frac{t_{\rm window}}{\Delta t_{\rm rise}}p_{\rm
flare}\approx\left(4.48\times
10^{-16}\frac{\text{flares}}{\text{AGN}}\right)n_{\rm AGN},$ (1)

which predicts $p_{cc}\sim\mathcal{O}(10^{-12})$ and
$p_{cc}\sim\mathcal{O}(10^{-13})$ for the number density and Milliquas
estimates, respectively.

We compare our candidate flare with the sample of Zwicky Transient Facility
(ZTF) transients from Graham et al. (2023). Assuming a roughly constant
spectral flux density across the $g$\- and $B$-band frequency domain, we apply
the $B$ bandpass bolometric correction of 5.15 from Duras et al. (2020) to
derive a total flare energy of

$\begin{split}E_{\rm tot}&\approx 4\pi d_{L}^{2}\left(5.15\int\mathrm{d}\nu
T_{B}(\nu)\right)\int_{t_{d}}^{t_{e}}\mathrm{d}t10^{-0.1m_{g,\rm
AB}(t)+48.6}\\\ &\approx 1.2\times 10^{50}\text{ erg},\end{split}$ (2)

where $T_{B}(\nu)$ is the transmission function for the $B$ bandpass (Bessell
1990), and $t_{d}$ and $t_{e}$ are the times of the first and last
observations of the flare. Fitting a Gaussian rise-exponential decay model to
the $g$-band light curve yields a rise time parameter (the standard deviation
of the Gaussian) of $t_{\rm rise}\approx 5.25\text{ days}$. This places this
transient in an intermediate regime between the SN and TDE populations from
Graham et al. 2023 (c.f. Figure 2 in that paper) among fast TDEs and high-
energy SNe. Notably, the Graham et al. 2023 proposed flare counterparts are
all slower and more energetic than AT 2023aagj; this may be associated with
the measurement of our flare parameters being based on difference photometry,
while the ZTF light curves used in the compared study are assembled via direct
photometry that includes the flux of the host AGN.

We estimate physical parameters of the system following the methodology of
Graham et al. (2020), which tests the hypothesis that the flare is generated
from a Bondi accreting BH as the remnant emerges from the accretion disk due
to a post-merger kick. With a flare onset time delay of $\sim 2$ days from the
GW trigger and total mass estimate of $M_{\rm tot}\sim 90M_{\odot}$, under the
assumption of association we can conclude one or more of the following: either
(i) the kick velocity $v_{\rm kick}$ is large, or (ii) the disk is
geometrically thin (aspect ratio $h\leq\mathcal{O}(10^{-3})$) or (iii) the
disk is not optically thick away from a geometrically thin mid-plane, or (iv)
there is a cavity in the disk due to feedback from the pre-merger BBH.

We continue by assuming the kick velocity is not specially large or small. If
we assume a kick velocity $v_{\rm kick}\sim\mathcal{O}(200){\rm km/s}$ (which
is consistent with the peak of the prior in Varma et al. 2022) then we can
constrain the approximate disk height from (Graham et al., 2023)

$t_{\rm exit}\sim\frac{H\sqrt{2{\rm ln}\tau_{\rm mp}}}{v_{\rm kick}}$ (3)

where $H$ is the disk height and $\tau_{\rm mp}$ is the mid-plane optical
depth. Since $t_{\rm exit}\sim 2\text{ days}$, $\tau_{\rm mp}\sim 10^{[3,6]}$
yields a factor $4-5$ from the square root, so $H\sim 0.3r_{g}(M_{\rm
SMBH}/2\times 10^{7}M_{\odot})$. This implies a thin disk, possibly similar to
the thin regions of a Thompson et al. (2005) model, which also has a lower
$\tau_{\rm mp}$ than other models such as Sirko & Goodman (2003). The thinner
regions of AGN disks are where gas damping is most efficient and consequently
where mergers are more likely to occur (McKernan & Ford, 2024).

If we assume from Figure 3 that the flare begins 2 days post-S230922g (denoted
by a vertical line), and the flare lasts through the spectra denoted by the
blue and red lines, then the overall flare duration is $t_{\rm
flare}\sim\mathcal{O}(10^{2})$days. Given $E_{\rm tot}\sim 10^{50}{\rm erg}$,
the average luminosity of the flare is ${\overline{L}_{\rm flare}}\sim
10^{43}{\rm erg/s}$, which can be parameterized as $\overline{L}_{\rm
flare}\sim 10^{3}L_{\rm Edd}(M_{\rm BBH}/90M_{\odot})$ where $L_{\rm Edd}$ is
the Eddington luminosity. If we parameterize the Bondi accretion rate onto the
merged BH ($\dot{M}_{\rm BHL}$) as:

$\dot{M}_{\rm BHL}\sim\frac{0.03M_{\odot}}{{\rm yr}}\left(\frac{M_{\rm
BBH}}{90M_{\odot}}\right)^{2}\\\ \times\left(\frac{v_{\rm rel}}{200{\rm
km/s}}\right)^{-3}\left(\frac{\rho_{\rm disk}}{10^{-11}{\rm g/cm^{3}}}\right)$
(4)

where $\rho_{\rm disk}$ is the disk gas density and we assume the relative
velocity ($v_{\rm rel}$) between the BH and the AGN disk gas is $v_{\rm
rel}\sim v_{\rm kick}$, i.e. the sound speed in the AGN disk gas is less than
$v_{\rm kick}$. We have chosen $\rho_{\rm disk}$ to be comparable to the
typical density in a Thompson et al. (2005) model where the disk is near its
thinnest ( aspect ratio $h\sim 10^{-3}$). Then $L_{\rm flare}\sim
10^{-2}L_{\rm BHL}$, where $L_{\rm BHL}=\eta\dot{M}_{\rm BHL}c^{2}$.
Interestingly, this accretion rate is within an order of magnitude of the
inferred accretion rate onto embedded BBH in AGN from recent GRMHD simulations
(Dittmann et al., 2024).

Associating the line asymmetry with an off-center flare implies that the
signature persists in the BLR for $\sim 70\text{ days}$ post event, but is
gone $\sim 90\text{ days}$ post event. A light-travel time across the full BLR
of $\mathcal{O}(90\text{ days})$ ($\sim 80$ in rest frame), corresponds to a
distance scale of $0.07{\rm pc}\sim 7\times 10^{4}r_{g}(M_{\rm SMBH}/2\times
10^{7}M_{\odot})$, implying that the merger occurred around $4\times
10^{4}r_{g}(M_{\rm SMBH}/2\times 10^{7}M_{\odot})$. This is located further
out than the thinnest regions of a Sirko & Goodman (2003) model disk, so we
can rule out an origin in a disk similar to that model. This location could be
consistent with the thinnest ($h\sim 10^{-3}$) regions of a Thompson et al.
(2005) model disk, which extend $\mathcal{O}[0.01,0.1]{\rm pc}$.

An alternative candidate for an off-center luminous flare in an AGN is a
$\mu$-TDE, where a star is tidally disrupted by a stellar mass BH embedded in
an AGN disk. Perna et al. (2021) show that in the outer regions ($\sim 7\times
10^{4}r_{g}(M_{\rm SMBH}/2\times 10^{7}M_{\odot}$) of a Thompson et al. (2005)
disk model, a $\mu$-TDE will peak in afterglow very quickly ($\sim$hours),
which is inconsistent with our observations. However, if $M_{\rm SMBH}$ is
actually a factor of a few more massive, then the diffusion time becomes
longer, potentially months, which could be consistent with our observations.

We also note that the lightcurve of AT 2023aagj resembles those of afterglows
in AGN disks as found in Wang et al. (2022). The rise of an afterglow from a
BBH merger is not implausible, assuming that these objects are able to launch
jets, as explored by Tagawa et al. (2023). Future work could investigate the
possibility of fitting those afterglow models in high-density environments.

Under the assumption that AT 2023aagj is the counterpart to S230922g, our
requirements are that the AGN disk must be relatively thin, similar to a
Thompson et al. (2005) model, and that the accretion rate onto the kicked BH
is significantly super-Eddington, but well below the Bondi rate (and requires
an associated jet for the radiation to emerge). One of our most likely false
positive for this event is a $\mu$-TDE in the disk if our estimate of $M_{\rm
SMBH}$ is too small by a factor of a few.

##  6 Conclusions

In this work we present the results of the GW-MMADS follow-up of S230922g. We
discuss in detail our most likely candidate AT 2023aagj, noting it is mostly
constant blue color, AGN-hosted nature, and especially the presence of
asymmetry in BLR spectral features as evidence in favor of recognizing it as a
counterpart to S230922g. However, we do not find the current dataset
sufficient to clearly identify the transient as such.

It is worth noting that several aspects lead to an inconclusive outcome
regarding the significance of this association. The variable nature of the AGN
in question, showing a significant flare in the infrared and a change in X-ray
flux over the past decade, exemplifies that there is some inherent activity in
the host, and the discovery of a GW counterpart is only possible if a means of
distinguishing the transient signal from any host variability. One possible
discriminator is the observation of asymmetric components in the broad lines
of the AGN spectrum, indicative of an off-center phenomenon; however, other
sources such as disk winds and tilted dust obscuration may cause similar
profiles. A smoking gun for the BBH association would be the evolution of the
asymmetric component in sync with the transient light curve, with the
asymmetry shifting towards bluer wavelengths and eventually disappearing as
the transient fades. For AT 2023aagj, due to the lack of wavelength coverage
for the Gemini spectrum and subsequent lack of multiple observations of MgII,
a chip gap on H$\beta$ in the Keck spectrum, the telluric region blueward of
H$\alpha$, and the target setting in December, the extraction of such an
evolution is not feasible with the present data. A spectrum of this object in
2024 as the target rises again may be informative about the nature of the
asymmetry. A persistent line asymmetry long after the merger event may point
to similarly persistent asymmetric illumination in this source, such as a
warped disk, and could rule out the association between AT 2023aagj and
S230922g. Finally, given the blue color and the ParSNIP prediction, we cannot
exclude that this transient may have originated from a TDE.

We consider the uncertainty of our result as another indicator of the
challenging nature of the search for BBH EM counterparts and successively the
need for the continuation and growth of future efforts for this purpose. More
so than for NS merger counterparts, predictions of BBH counterparts share a
parameter space with existing transient families (e.g. that of TDEs) and AGN
variability, and distinguishing between the different classes of events at the
present time may be challenging. Spectroscopic data are expected to be the
most helpful towards this end, and so it is important to combine the
appropriate resources with photometric follow-up while transients are still
active, for as long as it takes to develop models mature enough to perform
well without requiring such expensive resources. This is especially relevant
for the remainder of O4 follow-up, as all efforts made now will be the last
uptime for GW follow-up until the fifth gravitational wave observing run (O5)
begins in 2027. With Virgo having joined the GW detector network for O4 b, we
expect future O4 searches to yield better localized events, a trend that is
expected to continue with O5; these lower search volumes, in concert with
forthcoming powerful resources such as Legacy Survey of Space and Time (LSST)
(Ivezić et al., 2019), will help future counterpart searches be more efficient
and effective, and will help deepen and broaden the impact of GW and
multimessenger astronomy in the years to come.

## Acknowledgements

TC, AP, and LH acknowledge that this material is based upon work supported by
NSF Grant No. 2308193. BO gratefully acknowledges support from the McWilliams
Postdoctoral Fellowship at Carnegie Mellon University. BM & KESF are supported
by NSF AST-2206096 and NSF AST-1831415 and Simons Foundation Grant 533845 as
well as Simons Foundation sabbatical support. The Flatiron Institute is
supported by the Simons Foundation. AP thanks Rosalba Perna, Armin Rest, and
Stephen Smartt for useful discussion.

This research used resources of the National Energy Research Scientific
Computing Center, a DOE Office of Science User Facility supported by the
Office of Science of the U.S. Department of Energy under Contract No. DE-
AC02-05CH11231 using NERSC award HEP-ERCAP0029208 and HEP-ERCAP0022871. This
work used resources on the Vera Cluster at the Pittsburgh Supercomputing
Center (PSC). We thank T.J. Olesky and the PSC staff for help with setting up
our software on the Vera Cluster. MWC acknowledges support from the National
Science Foundation with grant numbers PHY-2308862 and PHY-2117997.

This project used data obtained with the Dark Energy Camera (DECam), which was
constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES
Projects has been provided by the US Department of Energy, the US National
Science Foundation, the Ministry of Science and Education of Spain, the
Science and Technology Facilities Council of the United Kingdom, the Higher
Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli
Institute for Cosmological Physics at the University of Chicago, Center for
Cosmology and Astro-Particle Physics at the Ohio State University, the
Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M
University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho
de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de
Desenvolvimento Científico e Tecnológico and the Ministério da Ciência,
Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the
Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University
of California at Santa Cruz, the University of Cambridge, Centro de
Investigaciones Enèrgeticas, Medioambientales y Tecnològicas–Madrid, the
University of Chicago, University College London, the DES-Brazil Consortium,
the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH)
Zürich, Fermi National Accelerator Laboratory, the University of Illinois at
Urbana-Champaign, the Institut de Ciéncies de l’Espai (IEEC/CSIC), the
Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory,
the Ludwig-Maximilians Universit at M unchen and the associated Excellence
Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, the Ohio State University, the OzDES Membership Consortium, the
University of Pennsylvania, the University of Portsmouth, SLAC National
Accelerator Laboratory, Stanford University, the University of Sussex, and
Texas A&M University.

Based on observations at Cerro Tololo Inter-American Observatory, NSF’s
NOIRLab (NOIRLab Prop. ID 2022B-715089; PI: Palmese; 2023B-851374 , PI:
Andreoni & Palmese; 2023B-735801, PI: Palmese & Wang), which is managed by the
Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation. We thank Kathy
Vivas, Alfredo Zenteno, and CTIO staff for their support with DECam
observations.

Based on observations obtained at the international Gemini Observatory (Prop.
ID GN-2023B-DD-103, PI: Cabrera; GN-2023B-109, PI: Cabrera), a program of NSF
NOIRLab, which is managed by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the U.S. National Science
Foundation on behalf of the Gemini Observatory partnership: the U.S. National
Science Foundation (United States), National Research Council (Canada),
Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia,
Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia,
Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science
Institute (Republic of Korea). Gemini data were processed using DRAGONS (Data
Reduction for Astronomy from Gemini Observatory North and South).

This work was enabled by observations made from the Gemini North telescope,
located within the Maunakea Science Reserve and adjacent to the summit of
Maunakea. We are grateful for the privilege of observing the Universe from a
place that is unique in both its astronomical quality and its cultural
significance.

Some of the observations reported in this paper were obtained with the
Southern African Large Telescope (SALT).

The Legacy Surveys consist of three individual and complementary projects: the
Dark Energy Camera Legacy Survey (DECaLS; Proposal ID 2014B-0404; PIs: David
Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID
2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey
(MzLS; Prop. ID 2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together
include data obtained, respectively, at the Blanco telescope, Cerro Tololo
Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward
Observatory, University of Arizona; and the Mayall telescope, Kitt Peak
National Observatory, NOIRLab. Pipeline processing and analyses of the data
were supported by NOIRLab and the Lawrence Berkeley National Laboratory
(LBNL). The Legacy Surveys project is honored to be permitted to conduct
astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular
significance to the Tohono O’odham Nation. LBNL is managed by the Regents of
the University of California under contract to the U.S. Department of Energy.

BASS is a key project of the Telescope Access Program (TAP), which has been
funded by the National Astronomical Observatories of China, the Chinese
Academy of Sciences (the Strategic Priority Research Program “The Emergence of
Cosmological Structures” Grant XDB09000000), and the Special Fund for
Astronomy from the Ministry of Finance. The BASS is also supported by the
External Cooperation Program of Chinese Academy of Sciences (Grant
114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant
12120101003, 11433005).

The Legacy Survey team makes use of data products from the Near-Earth Object
Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet
Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by
the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director,
Office of Science, Office of High Energy Physics of the U.S. Department of
Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract; and by the U.S. National Science Foundation, Division of
Astronomical Sciences under Contract No. AST-0950945 to NOAO.

This research has made use of the NASA/IPAC Extragalactic Database (NED),
which is funded by the National Aeronautics and Space Administration and
operated by the California Institute of Technology.

##  Appendix A Additional data

This section contains light curves (Figures 5-8) and spectra (Figure 9) for
all candidates other than AT 2023aagj.

For some transients we took additional data with the Three Channel Camera
(3KK) (Lang-Bardl et al., 2016) at the 2.1m-Fraunhofer Wendelstein Telescope
(Hopp et al., 2014) located in the German alps. The instrument is capable of
taking three images at the same time in two optical passbands and one near-
infrareds (NIRs) band. For the monitoring of the transients we used the $g$,
$i$ and $J$ bands. The data were reduced with a pipeline developed at the
University Observatory Munich. For the detrending of the images the pipeline
makes use of the tools from Gössl & Riffeser (2002). Then the images are
calibrated and coadded using SExtractor (Bertin & Arnouts, 1996), SCAMP
(Bertin, 2006) and SWarp (Bertin et al., 2002), difference photometry is
conducted with SFFT (Hu et al., 2022), and extinction corrections are applied
as described in §3.1. During this process, the astrometric solution is
computed against the Gaia EDR3 catalog (Gaia Collaboration et al., 2016, 2021)
and the zero points for the optical bands are calibrated with the PanSTARRS 1s
(PS1s) catalog (Magnier et al., 2013). The NIRs zero points are matched to the
Two Micron All Sky Surveys (2MASSs) catalog (Skrutskie et al., 2006) and
converted to AB magnitudes (Blanton & Roweis, 2007). Wendelstein data are
plotted as outlined diamonds in Figure 5, with $g$, $i$, and $J$ band data
plotted in green, yellow, and red, respectively.

Figure 5:  Light curves for our remaining 22 candidates (continued in
following figures). The dashed line indicates the S230922g event time. The
sample stamps for each transient are taken from the exposure with the highest
SNR, indicated with a gray square. Data taken with Wendelstein appear as small
diamonds, where relevant (the red point for AT 2023uab is $J$-band).  Figure
6: Light curves for our remaining 22 candidates (cont.). Figure 7: Light
curves for our remaining 22 candidates (cont.). Figure 8: Light curves for our
remaining 22 candidates (cont.). Figure 9:  Additional spectra taken as a part
of our follow-up campaign. A selection of spectral lines and telluric features
are marked as in Figure 3.

## References

  * Abbott et al. (2017a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, The Astrophysical Journal, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2017b) —. 2017b, The Astrophysical Journal, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017c) —. 2017c, Nature, 551, 85, doi: 10.1038/nature24471
  * Abbott et al. (2018) —. 2018, Physical Review Letters, 121, 161101, doi: 10.1103/PhysRevLett.121.161101
  * Alves et al. (2024) Alves, L. M. B., Sullivan, A. G., Yang, Y., et al. 2024, Monthly Notices of the Royal Astronomical Society, 531, 3679–3683, doi: 10.1093/mnras/stae1360
  * Andreoni et al. (2017) Andreoni, I., Ackley, K., Cooke, J., et al. 2017, PASA, 34, e069, doi: 10.1017/pasa.2017.65
  * Andreoni et al. (2019) Andreoni, I., Goldstein, D. A., Anand, S., et al. 2019, The Astrophysical Journal, 881, L16, doi: 10.3847/2041-8213/ab3399
  * Andreoni et al. (2020) Andreoni, I., Goldstein, D. A., Kasliwal, M. M., et al. 2020, The Astrophysical Journal, 890, 131, doi: 10.3847/1538-4357/ab6a1b
  * Annis et al. (2016) Annis, J., Soares-Santos, M., Berger, E., et al. 2016, The Astrophysical Journal Letters, 823, L34, doi: 10.3847/2041-8205/823/2/L34
  * Arcavi et al. (2017) Arcavi, I., Hosseinzadeh, G., Howell, D. A., et al. 2017, Nature, 551, 64, doi: 10.1038/nature24291
  * 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
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, The Astrophysical Journal Supplement Series, 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, Astronomy and Astrophysics, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, The Astronomical Journal, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, The Astrophysical Journal, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bartos et al. (2017a) Bartos, I., Haiman, Z., Marka, Z., et al. 2017a, Nature Communications, 8, 831, doi: 10.1038/s41467-017-00851-7
  * Bartos et al. (2017b) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017b, The Astrophysical Journal, 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bellm & Sesar (2016) Bellm, E. C., & Sesar, B. 2016, pyraf-dbsp: Reduction pipeline for the Palomar Double Beam Spectrograph.  http://ascl.net/1602.002
  * Berry et al. (2015) Berry, C. P. L., Mandel, I., Middleton, H., et al. 2015, The Astrophysical Journal, 804, 114, doi: 10.1088/0004-637x/804/2/114
  * Bertin (2006) Bertin, E. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 112
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astronomy and Astrophysics Supplement Series, 117, 393, doi: 10.1051/aas:1996164
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bertin et al. (2002) Bertin, E., Mellier, Y., Radovich, M., et al. 2002, Astronomical Data Analysis Software and Systems XI, 281, 228
  * Bertin et al. (2002) Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 281, Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228
  * Bessell (1990) Bessell, M. S. 1990, Publications of the Astronomical Society of the Pacific, 102, 1181, doi: 10.1086/132749
  * Blanton & Roweis (2007) Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734, doi: 10.1086/510127
  * Bom & Palmese (2023) Bom, C. R., & Palmese, A. 2023, Standard Siren Cosmology with Gravitational Waves from Binary Black Hole Mergers in Active Galaxy Nuclei, arXiv, doi: 10.48550/arXiv.2307.01330
  * Boone (2021) Boone, K. 2021, The Astronomical Journal, 162, 275, doi: 10.3847/1538-3881/ac2a2d
  * Buckley et al. (2006) Buckley, D. A. H., Swart, G. P., & Meiring, J. G. 2006, Ground-based and Airborne Telescopes, 6267, 62670Z, doi: 10.1117/12.673750
  * Burgh et al. (2003) Burgh, E. B., Nordsieck, K. H., Kobulnicky, H. A., et al. 2003, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, 4841, 1463, doi: 10.1117/12.460312
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Science Reviews, 120, 165, doi: 10.1007/s11214-005-5097-2
  * Cabrera et al. (2023) Cabrera, T., Hu, L., Andreoni, I., et al. 2023, GRB Coordinates Network, 34763, 1
  * Cappellari (2023) Cappellari, M. 2023, MNRAS, 526, 3273, doi: 10.1093/mnras/stad2597
  * Chen et al. (2018) Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545, doi: 10.1038/s41586-018-0606-0
  * Chornock et al. (2017) Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19, doi: 10.3847/2041-8213/aa905c
  * Clark et al. (2024) Clark, P., Graur, O., Callow, J., et al. 2024, Monthly Notices of the Royal Astronomical Society, 528, 7076, doi: 10.1093/mnras/stae460
  * Cook et al. (2023) Cook, D. O., Mazzarella, J. M., Helou, G., et al. 2023, Completeness of the NASA/IPAC Extragalactic Database (NED) – Local Volume Sample.  https://arxiv.org/abs/2306.06271
  * Coughlin et al. (2018) Coughlin, M. W., Tao, D., Chan, M. L., et al. 2018, Monthly Notices of the Royal Astronomical Society, 478, 692, doi: 10.1093/mnras/sty1066
  * Coughlin et al. (2023) Coughlin, M. W., Bloom, J. S., Nir, G., et al. 2023, A Data Science Platform to Enable Time-Domain Astronomy, arXiv, doi: 10.48550/arXiv.2305.00108
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556, doi: 10.1126/science.aap9811
  * DESI Collaboration et al. (2023) DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2023, arXiv e-prints, arXiv:2306.06308, doi: 10.48550/arXiv.2306.06308
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168, doi: 10.3847/1538-3881/ab089d
  * Dittmann et al. (2024) Dittmann, A. J., Dempsey, A. M., & Li, H. 2024, ApJ, 964, 61, doi: 10.3847/1538-4357/ad23ce
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 1570, doi: 10.1126/science.aaq0049
  * Duras et al. (2020) Duras, F., Bongiorno, A., Ricci, F., et al. 2020, Astronomy and Astrophysics, 636, A73, doi: 10.1051/0004-6361/201936817
  * Evans et al. (2023) Evans, P. A., Page, K. L., Beardmore, A. P., et al. 2023, MNRAS, 518, 174, doi: 10.1093/mnras/stac2937
  * Evans et al. (2009) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177, doi: 10.1111/j.1365-2966.2009.14913.x
  * Evans et al. (2017) Evans, P. A., Cenko, S. B., Kennea, J. A., et al. 2017, Science, 358, 1565, doi: 10.1126/science.aap9580
  * Finn & Chernoff (1993) Finn, L. S., & Chernoff, D. F. 1993, Physical Review D, 47, 2198–2219, doi: 10.1103/physrevd.47.2198
  * Flaugher et al. (2015) Flaugher, B., Diehl, H. T., Honscheid, K., et al. 2015, The Astronomical Journal, 150, 150, doi: 10.1088/0004-6256/150/5/150
  * Flesch (2023) Flesch, E. W. 2023, The Million Quasars (Milliquas) Catalogue, V8, arXiv, doi: 10.48550/arXiv.2308.01505
  * Flewelling et al. (2020) Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2020, The Astrophysical Journal Supplement Series, 251, 7, doi: 10.3847/1538-4365/abb82d
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, Astronomy and Astrophysics, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Garcia et al. (2020) Garcia, A., Morgan, R., Herner, K., et al. 2020, The Astrophysical Journal, 903, 75, doi: 10.3847/1538-4357/abb823
  * Gayathri et al. (2021) Gayathri, V., Yang, Y., Tagawa, H., Haiman, Z., & Bartos, I. 2021, The Astrophysical Journal Letters, 920, L42, doi: 10.3847/2041-8213/ac2cc1
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, The Astrophysical Journal, 611, 1005, doi: 10.1086/422091
  * Goldstein et al. (2017) Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJ, 848, L14, doi: 10.3847/2041-8213/aa8f41
  * Goldstein et al. (2019) Goldstein, D. A., Andreoni, I., Nugent, P. E., et al. 2019, The Astrophysical Journal, 881, L7, doi: 10.3847/2041-8213/ab3046
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, The Astrophysical Journal, 622, 759, doi: 10.1086/427976
  * Gössl & Riffeser (2002) Gössl, C. A., & Riffeser, A. 2002, A&A, 381, 1095, doi: 10.1051/0004-6361:20011522
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Physical Review Letters, 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Graham et al. (2023) Graham, M. J., McKernan, B., Ford, K. E. S., et al. 2023, The Astrophysical Journal, 942, 99, doi: 10.3847/1538-4357/aca480
  * Green (2018) Green, G. M. 2018, Journal of Open Source Software, 3, 695, doi: 10.21105/joss.00695
  * Greene & Ho (2005) Greene, J. E., & Ho, L. C. 2005, ApJ, 630, 122, doi: 10.1086/431897
  * Greene & Ho (2007) Greene, J. E., & Ho, L. C. 2007, Astrophys. J., 667, 131, doi: 10.1088/0004-637X/704/2/1743
  * Guo et al. (2018) Guo, H., Shen, Y., & Wang, S. 2018, PyQSOFit: Python code to fit the spectrum of quasars, Astrophysics Source Code Library.  http://ascl.net/1809.008
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Herner et al. (2020) Herner, K., Annis, J., Brout, D., et al. 2020, Astronomy and Computing, 33, 100425, doi: 10.1016/j.ascom.2020.100425
  * Hilton et al. (2018) Hilton, M., Hasselfield, M., Sifón, C., et al. 2018, ApJS, 235, 20, doi: 10.3847/1538-4365/aaa6cb
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, The Astrophysical Journal, 629, 15, doi: 10.1086/431341
  * Hopp et al. (2014) Hopp, U., Bender, R., Grupp, F., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9145, Ground-based and Airborne Telescopes V, ed. L. M. Stepp, R. Gilmozzi, & H. J. Hall, 91452D, doi: 10.1117/12.2054498
  * Hotokezaka et al. (2019) Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, Nature Astronomy, 3, 940, doi: 10.1038/s41550-019-0820-1
  * Hu et al. (2024) Hu, L., Cabrera, T., Palmese, A., et al. 2024, in prep.
  * Hu et al. (2022) Hu, L., Wang, L., Chen, X., & Yang, J. 2022, The Astrophysical Journal, 936, 157, doi: 10.3847/1538-4357/ac7394
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, The Astrophysical Journal, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Kasliwal et al. (2017) Kasliwal, M., Nakar, E., Singer, L., et al. 2017, Science, 358, 1559
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055, doi: 10.1111/j.1365-2966.2003.07154.x
  * Kessler et al. (2019) Kessler, R., Narayan, G., Avelino, A., et al. 2019, Publications of the Astronomical Society of the Pacific, 131, 094501, doi: 10.1088/1538-3873/ab26f1
  * Kewley et al. (2001) Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121, doi: 10.1086/321545
  * Kewley et al. (2006) Kewley, L. J., Groves, B., Kauffmann, G., & Heckman, T. 2006, MNRAS, 372, 961, doi: 10.1111/j.1365-2966.2006.10859.x
  * Klimenko et al. (2016) Klimenko, S., Vedovato, G., Drago, M., et al. 2016, Physical Review D, 93, doi: 10.1103/physrevd.93.042004
  * Kobulnicky et al. (2003) Kobulnicky, H. A., Nordsieck, K. H., Burgh, E. B., et al. 2003, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, 4841, 1634, doi: 10.1117/12.460315
  * Labrie et al. (2019) Labrie, K., Anderson, K., Cárdenes, R., Simpson, C., & Turner, J. E. H. 2019, ASP Conference Series, 523, 321
  * Lang-Bardl et al. (2016) Lang-Bardl, F., Bender, R., Goessl, C., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 990844, doi: 10.1117/12.2232039
  * Ligo Scientific Collaboration et al. (2023a) Ligo Scientific Collaboration, VIRGO Collaboration, & Kagra Collaboration. 2023a, GRB Coordinates Network, 34757, 1
  * Ligo Scientific Collaboration et al. (2023b) —. 2023b, GRB Coordinates Network, 34758, 1
  * Ligo Scientific Collaboration et al. (2023c) —. 2023c, GRB Coordinates Network, 34758, 1
  * Lipunov et al. (2017) Lipunov, V. M., Gorbovskoy, E., Kornilov, V. G., et al. 2017, ApJ, 850, L1, doi: 10.3847/2041-8213/aa92c0
  * Magnier et al. (2013) Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20, doi: 10.1088/0067-0049/205/2/20
  * Mainzer et al. (2011) Mainzer, A., Grav, T., Bauer, J., et al. 2011, ApJ, 743, 156, doi: 10.1088/0004-637X/743/2/156
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30, doi: 10.1088/0004-637X/792/1/30
  * Margutti et al. (2017) Margutti, R., Berger, E., Fong, W., et al. 2017, ApJ, 848, L20, doi: 10.3847/2041-8213/aa9057
  * McKernan & Ford (2024) McKernan, B., & Ford, K. E. S. 2024, MNRAS, 531, 3479, doi: 10.1093/mnras/stae1351
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, The Astrophysical Journal, 884, L50, doi: 10.3847/2041-8213/ab4886
  * McKinney (2010) McKinney, W. 2010, Proceedings of the 9th Python in Science Conference, 56, doi: 10.25080/Majora-92bf1922-00a
  * McManus & Olsen (2021) McManus, S., & Olsen, K. 2021, The NOIRLab Mirror, 2, 33
  * Messick et al. (2017) Messick, C., Blackburn, K., Brady, P., et al. 2017, Physical Review D, 95, doi: 10.1103/physrevd.95.042001
  * Morgan et al. (2020) Morgan, R., Soares-Santos, M., Annis, J., et al. 2020, The Astrophysical Journal, 901, 83, doi: 10.3847/1538-4357/abafaa
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L18, doi: 10.3847/2041-8213/aa9029
  * Onken et al. (2024) Onken, C. A., Wolf, C., Bessell, M. S., et al. 2024, SkyMapper Southern Survey: Data Release 4.  https://arxiv.org/abs/2402.02015
  * Palmese et al. (2021) Palmese, A., Fishbach, M., Burke, C. J., Annis, J., & Liu, X. 2021, The Astrophysical Journal, 914, L34, doi: 10.3847/2041-8213/ac0883
  * Palmese et al. (2024) Palmese, A., Kaur, R., Hajela, A., et al. 2024, Physical Review D, 109, 063508
  * Perley (2019) Perley, D. A. 2019, PASP, 131, 084503, doi: 10.1088/1538-3873/ab215d
  * Perna et al. (2021) Perna, R., Lazzati, D., & Cantiello, M. 2021, ApJ, 906, L7, doi: 10.3847/2041-8213/abd319
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67, doi: 10.1038/nature24298
  * Rousselot et al. (2000) Rousselot, P., Lidman, C., Cuby, J. G., Moreels, G., & Monnet, G. 2000, Astronomy and Astrophysics, 354, 1134
  * Sachdev et al. (2019) Sachdev, S., Caudill, S., Fong, H., et al. 2019, The GstLAL Search Analysis Methods for Compact Binary Mergers in Advanced LIGO’s Second and Advanced Virgo’s First Observing Runs.  https://arxiv.org/abs/1901.08580
  * Savchenko et al. (2017) Savchenko, V., Ferrigno, C., Kuulkers, E., et al. 2017, ApJ, 848, L15, doi: 10.3847/2041-8213/aa8f94
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, The Astrophysical Journal, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310, doi: 10.1038/323310a0
  * Shappee et al. (2017) Shappee, B. J., Simon, J. D., Drout, M. R., et al. 2017, Science, 358, 1574, doi: 10.1126/science.aaq0186
  * Singer et al. (2016) Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016, The Astrophysical Journal, 829, L15, doi: 10.3847/2041-8205/829/1/L15
  * Sirko & Goodman (2003) Sirko, E., & Goodman, J. 2003, MNRAS, 341, 501, doi: 10.1046/j.1365-8711.2003.06431.x
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Smartt et al. (2017) Smartt, S. J., Chen, T. W., Jerkstrand, A., et al. 2017, Nature, 551, 75, doi: 10.1038/nature24303
  * Soares-Santos et al. (2016) Soares-Santos, M., Kessler, R., Berger, E., et al. 2016, The Astrophysical Journal Letters, 823, L33, doi: 10.3847/2041-8205/823/2/L33
  * Soares-Santos et al. (2017) Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, The Astrophysical Journal, 848, L16, doi: 10.3847/2041-8213/aa9059
  * Storey-Fisher et al. (2024) Storey-Fisher, K., Hogg, D. W., Rix, H.-W., et al. 2024, The Astrophysical Journal, 964, 69, doi: 10.3847/1538-4357/ad1328
  * Symbalisty & Schramm (1982) Symbalisty, E., & Schramm, D. N. 1982, Astrophysical Letters, 22, 143
  * Tagawa et al. (2023) Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2023, ApJ, 950, 13, doi: 10.3847/1538-4357/acc4bb
  * Tagawa et al. (2024) Tagawa, H., Kimura, S. S., Haiman, Z., Perna, R., & Bartos, I. 2024, Astrophys. J., 966, 21, doi: 10.3847/1538-4357/ad2e0b
  * Tanvir et al. (2017) Tanvir, N. R., Levan, A. J., González-Fernández, C., et al. 2017, ApJ, 848, L27, doi: 10.3847/2041-8213/aa90b6
  * Thompson et al. (2005) Thompson, T. A., Quataert, E., & Murray, N. 2005, ApJ, 630, 167, doi: 10.1086/431923
  * Troja et al. (2017) Troja, E., Piro, L., Van Eerten, H., et al. 2017, Nature, 551, 71
  * Utsumi et al. (2017) Utsumi, Y., Tanaka, M., Tominaga, N., et al. 2017, PASJ, 69, 101, doi: 10.1093/pasj/psx118
  * Valdes et al. (2014) Valdes, F., Gruendl, R., & DES Project. 2014, in Astronomical Society of the Pacific Conference Series, Vol. 485, Astronomical Society of the Pacific Conference Series, ed. N. Manset & P. Forshay, 379
  * Valenti et al. (2017) Valenti, S., Sand, D. J., Yang, S., et al. 2017, ApJ, 848, L24, doi: 10.3847/2041-8213/aa8edf
  * van der Walt et al. (2019) van der Walt, S. J., Crellin-Quick, A., & Bloom, J. S. 2019, Journal of Open Source Software, 4, doi: 10.21105/joss.01247
  * Varma et al. (2022) Varma, V., Biscoveanu, S., Islam, T., et al. 2022, Phys. Rev. Lett., 128, 191102, doi: 10.1103/PhysRevLett.128.191102
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409, doi: 10.1093/mnras/stw2231
  * Wang et al. (2022) Wang, Y.-H., Lazzati, D., & Perna, R. 2022, Monthly Notices of the Royal Astronomical Society, 516, 5935–5944, doi: 10.1093/mnras/stac1968
  * Wright et al. (2010a) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010a, The Astronomical Journal, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Wright et al. (2010b) —. 2010b, The Astronomical Journal, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Wyatt et al. (2020) Wyatt, S. D., Tohuvavohu, A., Arcavi, I., et al. 2020, The Astrophysical Journal, 894, 127, doi: 10.3847/1538-4357/ab855e
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, The Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298