# Follow-up survey for the binary black hole merger GW200224_222234 using
Subaru/HSC and GTC/OSIRIS

Takayuki Ohgami Division of Science, National Astronomical Observatory of
Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588, Japan  Josefa Becerra González
Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna, Tenerife, Spain
Universidad de La Laguna (ULL), Departamento de Astrofísica, E-38206 La
Laguna, Tenerife, Spain  Nozomu Tominaga Division of Science, National
Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo 181-8588,
Japan  Department of Astronomical Science, School of Physical Sciences, The
Graduate University of Advanced Studies, (SOKENDAI), 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba
277-8583, Japan  Department of Physics, Faculty of Science and Engineering,
Konan University, 8-9-1 Okamoto, Kobe, Hyogo 658-8501, Japan  Tomoki Morokuma
Planetary Exploration Research Center, Chiba Institute of Technology, 2-17-1
Tsudanuma, Narashino, Chiba 275-0016, Japan  Institute of Astronomy, Graduate
School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo
181-0015, Japan  Kavli Institute for the Physics and Mathematics of the
Universe (WPI), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba
277-8583, Japan  Yousuke Utsumi Kavli Institute for Particle Astrophysics and
Cosmology (KIPAC), SLAC National Accelerator Laboratory, Stanford University,
2575 Sand Hill Road, Menlo Park, CA 94025, USA  Yuu Niino Institute of
Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa,
Mitaka, Tokyo 181-0015, Japan  Research Center for the Early Universe,
Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Masaomi Tanaka Astronomical Institute, Tohoku
University, Sendai, Miyagi 980-8578, Japan  Smaranika Banerjee Astronomical
Institute, Tohoku University, Sendai, Miyagi 980-8578, Japan  Frédérick
Poidevin Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna,
Tenerife, Spain  Universidad de La Laguna (ULL), Departamento de Astrofísica,
E-38206 La Laguna, Tenerife, Spain  Jose Antonio Acosta-Pulido Instituto de
Astrofísica de Canarias (IAC), E-38200 La Laguna, Tenerife, Spain  Universidad
de La Laguna (ULL), Departamento de Astrofísica, E-38206 La Laguna, Tenerife,
Spain  Ismael Pérez-Fournon Instituto de Astrofísica de Canarias (IAC),
E-38200 La Laguna, Tenerife, Spain  Universidad de La Laguna (ULL),
Departamento de Astrofísica, E-38206 La Laguna, Tenerife, Spain  Teo Muñoz-
Darias Instituto de Astrofísica de Canarias (IAC), E-38200 La Laguna,
Tenerife, Spain  Universidad de La Laguna (ULL), Departamento de Astrofísica,
E-38206 La Laguna, Tenerife, Spain  Hiroshi Akitaya Planetary Exploration
Research Center, Chiba Institute of Technology, 2-17-1 Tsudanuma, Narashino,
Chiba 275-0016, Japan  Hiroshima Astrophysical Science Center, Hiroshima
University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan
Kenshi Yanagisawa Subaru Telescope, National Astronomical Observatory of
Japan, 650 North A’ohoku Place, Hilo, HI 96720, USA  Hiroshima Astrophysical
Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima,
Hiroshima 739-8526, Japan  Mahito Sasada Hiroshima Astrophysical Science
Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima
739-8526, Japan  Mizusawa VLBI Observatory, National Astronomical Observatory
of Japan, 2-12 Hoshigaoka, Mizusawa, Oshu, Iwate 023-0861, Japan  Michitoshi
Yoshida Subaru Telescope, National Astronomical Observatory of Japan, 650
North A’ohoku Place, Hilo, HI 96720, USA  Mirko Simunovic Subaru Telescope,
National Astronomical Observatory of Japan, 650 North A’ohoku Place, Hilo, HI
96720, USA  Ryou Ohsawa Institute of Astronomy, Graduate School of Science,
The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan  JASMINE
Project, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka,
Tokyo 181-8588, Japan  Ichi Tanaka  Subaru Telescope, National Astronomical
Observatory of Japan, 650 North A’ohoku Place, Hilo, HI 96720, USA  Tsuyoshi
Terai Subaru Telescope, National Astronomical Observatory of Japan, 650 North
A’ohoku Place, Hilo, HI 96720, USA  Yuhei Takagi Subaru Telescope, National
Astronomical Observatory of Japan, 650 North A’ohoku Place, Hilo, HI 96720,
USA  The J-GEM collaboration

###### Abstract

The LIGO/Virgo detected a gravitational wave (GW) event, named GW200224_222234
(a.k.a. S200224ca) and classified as a binary-black-hole coalescence, on
February 24, 2020. Given its relatively small localization skymap (71 deg2 for
a 90% credible region; revised to 50 deg2 in GWTC-3), we performed target-of-
opportunity observations using the Subaru/Hyper Suprime-Cam (HSC) in the
$r2$\- and $z$-bands. Observations were conducted on February 25 and 28 and
March 23, 2020, with the first epoch beginning 12.3 h after the GW detection.
The survey covered the highest probability sky area of 56.6 deg2,
corresponding to a 91% probability. This was the first deep follow-up
($m_{r}\gtrsim 24,~{}m_{z}\gtrsim 23$) for a binary-black-hole merger covering
$>$90% of the localization. By performing image subtraction and candidate
screening including light curve fitting with transient templates and examples,
we found 22 off-nucleus transients that were not ruled out as the counterparts
of GW200224_222234 with only our Subaru/HSC data. We also performed GTC/OSIRIS
spectroscopy of the probable host galaxies for five candidates; two are likely
to be located within the 3D skymap, whereas the others are not. In conclusion,
19 transients remain as possible optical counterparts of GW200224_222234;
however, we could not identify a unique promising counterpart. If there are no
counterparts in the remaining candidates, the upper limits of optical
luminosity are $\nu L_{\nu}<5.2^{+2.4}_{-1.9}\times 10^{41}$ erg s-1 and $\nu
L_{\nu}<1.8^{+0.8}_{-0.6}\times 10^{42}$ erg s-1 in the $r2$\- and $z$-bands,
respectively, at $\sim$12 h after GW detection. We also discuss improvements
in the strategies of optical follow-ups for future GW events.

Gravitational waves (678); Black holes (162); Surveys (1671)

††facilities: Subaru (HSC), GTC(OSIRIS)††software: hscPipe (Bosch et al.,
2018), astropy (Astropy Collaboration et al., 2013, 2018), Source Extractor
(Bertin & Arnouts, 1996)

##  1 Introduction

In general relativity, massive objects radiate energy via the distortion of
space time when their motion is accelerated. This energy radiation predicted
by Einstein is called a gravitational wave (GW; Einstein, 1916, 1918).
Astronomical objects or phenomena are expected to be sources of GW signals
with a large amplitude, which may be detected by current instruments. For
example, binary systems composed of compact objects such as black holes (BHs)
or neutron stars (NSs) emit strong GWs at their coalescence.

In the second observation run of the GW interferometers LIGO and Virgo, they
detected a GW signal from a binary-NS (BNS) coalescence using three detectors,
and the localization area was constrained to 28 deg2 for a 90% credible region
(GW170817; Abbott et al., 2017). The electromagnetic (EM) counterpart of GW
events was observed for the first time by multiple observatories across the EM
spectrum from radio to $\gamma$-rays (Andreoni et al., 2017; Arcavi et al.,
2017; Chornock et al., 2017; Coulter et al., 2017; Cowperthwaite et al., 2017;
Díaz et al., 2017; Drout et al., 2017; Evans et al., 2017; Kasliwal et al.,
2017; Lipunov et al., 2017; Nicholl et al., 2017; Pian et al., 2017; Smartt et
al., 2017; Soares-Santos et al., 2017; Troja et al., 2017; Tanvir et al.,
2017; Tominaga et al., 2018a; Utsumi et al., 2017; Valenti et al., 2017). It
was demonstrated that BNS mergers are accompanied by explosive EM emissions
called kilonova (e.g., Kasen et al., 2013; Shibata et al., 2017; Tanaka et
al., 2017; Kawaguchi et al., 2018; Perego et al., 2017; Rosswog et al., 2018;
Banerjee et al., 2020).

In contrast to NS mergers, binary-BH (BBH) mergers are not considered to be
accompanied by any EM emission. However, the Fermi Gamma-ray Burst Monitor
(GBM) reported the presence of a weak $\gamma$-ray transient after the
detection of GW150914 (Connaughton et al., 2016). Although the physical
association between GWs and $\gamma$-ray signals is ambiguous, and it is
unclear whether the $\gamma$-ray signal is truly astronomical because of the
low flux, various scenarios have been proposed, such as BBHs surrounded by
pre-existing material, for example, a circumbinary disc (Martin et al. 2018),
an accretion disc surrounding a galactic center BH (Bartos et al. 2017), and
remnants of gravitational collapse (Janiuk et al. 2017). Another case is the
detection of an optical counterpart candidate of GW190521 by the Zwicky
Transient Facility (ZTF19abanrhr). Graham et al. (2020) suggested that the EM
flare is consistent with the behavior expected from a kicked BBH merger in an
accretion disk of an active galactic nucleus (AGN, McKernan et al. 2019) and
ruled out other scenarios (for example, the intrinsic variability of AGN,
supernova, microlens, tidal disruption). Graham et al. (2022) comprehensively
searched for EM counterparts to BBH mergers and identified nine candidates. To
test various scenarios for EM radiation from a BBH merger, it is still
important to perform follow-ups of the BBH mergers.

The LIGO/Virgo collaboration began their third observation run (O3) in April
2019 and detected a GW event named GW200224_222234 (a.k.a. S200224ca) using
three detectors on February 24, 2020 at 22:22 UTC (The LIGO Scientific
Collaboration et al., 2021; LIGO Scientific Collaboration & Virgo
Collaboration, 2020). They released a preliminary localization skymap derived
using software called BAYESTAR (Singer & Price, 2016) on February 24, 2020 at
22:32 UTC. In this release, the GW event was classified as a BBH coalescence
with a $>99$% confidence level and a false alarm rate of $1.6\times 10^{-11}$
Hz (approximately one in 1975 years). The luminosity distance was $1583\pm
331$ Mpc, corresponding to a redshift of $0.30\pm 0.05$, and the 90%
localization sky area was as narrow as 71 deg2.

Upon receiving the alert for GW200224_222234, the Japanese collaboration for
Gravitational wave ElectroMagnetic follow-up (J-GEM; Morokuma et al., 2016;
Utsumi et al., 2017; Sasada et al., 2021) triggered a target-of-opportunity
(ToO) observation to search for its EM counterpart on February 25, 2020 at
10:43 UTC (Ohgami et al., 2020) using Hyper Suprime-Cam (HSC; Furusawa et al.,
2018; Kawanomoto et al., 2018; Komiyama et al., 2018; Miyazaki et al., 2012,
2018), which is a wide-field imager installed on the prime focus of the Subaru
Telescope. Its field of view (FoV) of 1.77 deg2 is the largest among current 8
m-class telescopes, which makes the Subaru/HSC the most efficient instrument
for the optical survey. The first exposure commenced at approximately 12.3 h
after GW detection, and the observation area reached 56.6 deg2. We also
conducted additional ToO observations using the Subaru/HSC on February 28 and
March 23, 2020, in the same fields as the first epoch observation. The LIGO
and Virgo collaboration published the GWTC-3 catalog (The LIGO Scientific
Collaboration et al., 2021) including the observations O1, O2, O3a, and O3b
and also released the reanalyzed the localization skymap. The 90% localization
sky area was updated from 71 deg2 to 50 deg2. Our observation area corresponds
to a cumulative probability of 91% in the updated skymap.

In this paper, we describe the details of the follow-ups of GW200224_222234
using the Subaru/HSC, the candidate selection algorithm, and a list of
candidates, including spectroscopic observations of the probable host galaxies
of the five candidates to measure their spectroscopic redshifts using the
Optical System for Imaging and low-Intermediate-Resolution Integrated
Spectroscopy111http://www.gtc.iac.es/instruments/osiris/ (OSIRIS), which is an
imager and spectrograph for the optical wavelength range, installed in the
10.4- m Gran Telescopio CANARIAS (GTC). All magnitudes are given in AB
magnitudes.

##  2 Observations with the Subaru/HSC and data analysis

###  2.1 ToO observation

We conducted optical imaging observations using the Subaru/HSC on February 25
(Day 1), 28 (Day 4), and March 23 (Day 28), 2020, in the $r2$\- and $z$-bands.
To trigger our ToO follow-up as rapidly as possible, we conducted the first
observations in the $r2$-band, which had already been set in the Subaru/HSC at
the time of the GW alert. The first exposure commenced on February 25, 2020 at
10:43 UTC, corresponding to 12.3 h after GW detection. We selected 60
observation pointings to cover the high probability area in the BAYSTAR skymap
for the HEALPix grid with a resolution of NSIDE $=64$, which corresponds to
0.84 deg2 pixel-2, enabling overlap of the FoVs. Figure 1 shows the survey
pointing map. Our survey area covered 56.6 deg2, corresponding to a cumulative
probability of 91% in the localization skymap refined in the GWTC-3 catalog
using the IMRPhenomXPHM model (Pratten et al., 2021), which models higher-
order spherical harmonics and spin precession. In this study, we used the
localization skymap created using the IMRPhenomXPHM model as the 3D
localization of GW200224_222234.

Figure 1: Observation pointings with the Subaru/HSC (red filled circles). The
left and right panels are the preliminary localization skymap of
GW200224_222234 (BAYESTAR; LIGO Scientific Collaboration & Virgo
Collaboration, 2020) and a refined version using the IMRPhenomXPHM model
(GWTC-3 catalog; The LIGO Scientific Collaboration et al., 2021),
respectively. The white contour lines indicate the 90% credible regions.

On February 25 and 28, 2020, we observed these pointings with 30 s of exposure
each in pointing ID order and revisited them in each band. The revisits were
conducted at least 1 h apart with a 1′ offset in each pointing to fill the CCD
gaps. Note that some areas were observed again at a time interval of less than
1 h because of the overlap between each pointing. On March 23, 2020, we
changed the exposure time and the order of the pointings. The pointings at
higher elevation were observed earlier. In the $z$-band, we took the shots
with a 35- s exposure each and revisited them. In the $r2$-band, we took the
shots with 50- s exposures for the first 32 pointings and 70- s exposures for
the remaining 28 pointings. The central coordinates and exposure times are
shown in Table 1 (This table is published in its entirety in the machine-
readable format).

Table 1: Central coordinates and exposure times of the observation pointings
with the Subaru/HSC.

| Central Coordinate | Exposure Time
---|---|---
Pointing ID | R.A. (J2000.0) | Decl. (J2000.0) | 2020-02-25 (Day 1) | 2020-02-28 (Day 4) | 2020-03-23 (Day 28)
| (HH:MM:SS.ss) | (DD:MM:SS.s) | $r2$ | $z$ | $r2$ | $z$ | $r2$ | $z$
00 | 11:40:18.75 |  $-$14:28:39.0 | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 50 s | 35 s$\times 2$
01 | 11:34:41.25 |  $-$12:01:28.9 | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 70 s | 35 s$\times 2$
02 | 11:31:52.50 |  $-$10:11:59.7 | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 50 s | 35 s$\times 2$
03 | 11:43:07.50 |  $-$15:05:41.2 | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 70 s | 35 s$\times 2$
04 | 11:37:30.00 |  $-$12:38:08.3 | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 30 s$\times 2$ | 70 s | 35 s$\times 2$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.

###  2.2 Data reduction and image subtraction

We reduced the observational data using hscPipe v4.0.5 (Bosch et al., 2018),
which is a standard analysis pipeline for the HSC. This pipeline provides full
packages for data analyses, which include image subtraction and source
detection. We evaluated the limiting magnitudes of the stacked images in each
epoch as follows: We assumed an aperture with a diameter of twice the full
width at half maximum (FWHM) of the point spread function (PSF) and
distributed it randomly, avoiding existing sources. Measuring the $5\sigma$
deviations of fluxes in this aperture, we obtained a map of the $5\sigma$
limiting magnitudes for each stack image. Table 2 lists the mode222The mode
values are the peak values of the histograms of the $5\sigma$ limiting
magnitudes. Here, we distributed the $5\sigma$ limiting magnitudes into ten
bins ranging between the maximum and minimum values., maximum, and minimum
values of the $5\sigma$ limiting magnitudes. The $5\sigma$ limiting magnitudes
scattered with a range of $\sim 1$ mag, and there were differences of
approximately 1 mag between those of the $r2$\- and $z$-bands. Note that the
limiting magnitudes in the $r2$-band on March 23, 2020, had a small positional
dependence owing to the different exposure time.

Table 2: Mode, maximum, and minimum values of the $5\sigma$ limiting magnitudes in the stacked images. Date | Filter | Limiting magnitude (AB)
---|---|---
(Days from GW detection) |  | mode | max | min
2020-02-25 | $r2$ | 25.29 | 25.71 | 24.51
(Day 1) | $z$ | 23.56 | 24.14 | 22.84
2020-02-28 | $r2$ | 24.61 | 25.13 | 23.98
(Day 4) | $z$ | 23.21 | 23.76 | 22.54
2020-03-23 | $r2$ | 25.53 | 26.25 | 24.64
(Day 28) | $z$ | 23.83 | 24.64 | 22.85

We performed image subtraction for the science images obtained in the first
(one day from GW detection; Day 1) and second (Day 4) epochs using the images
taken in the third (Day 28) epoch as reference images. The subtraction package
is included in hscPipe v4.0.5 and is based on an algorithm in which the FWHM
of the reference images is fitted to that in the science images via
convolution using kernels to make their PSFs equivalent, as proposed in Alard
& Lupton (1998) and Alard (1999). Table 3 lists the mode, maximum, and minimum
values of the $5\sigma$ limiting magnitudes evaluated from the difference
images. The positional dependence of the $5\sigma$ limiting magnitudes in the
$r2$-band is shown in Figure 2.

Table 3: Mode, maximum, and minimum values of the $5\sigma$ limiting magnitudes in the difference images. Difference Image | Filter | Limiting magnitude (AB)
---|---|---
|  | mode | max | min
Day 1$-$28 | $r2$ | 24.99 | 25.56 | 24.29
| $z$ | 23.28 | 23.89 | 22.54
Day 4$-$28 | $r2$ | 24.39 | 25.15 | 23.76
| $z$ | 22.97 | 23.57 | 22.23
Figure 2: Map (top) and histogram (bottom) of the $5\sigma$ limiting
magnitudes in the difference images within the first (Day 1 $-$ Day 28; left)
and second epochs (Day 4 $-$ Day 28; right) in the $r2$-band. The color
indicates larger/smaller (in red/blue) values than the mode value of each
epoch.

###  2.3 Source detection and screening

After image subtraction, we applied the following criteria to exclude bogus
detections (for example, caused by bad pixels and failure of image
subtraction) and select point sources from the difference images as in
Tominaga et al. (2018a) and Ohgami et al. (2021): (i) A signal-to-noise ratio
of PSF flux $(S/N)_{\rm PSF}>5$; (ii) $(b/a)/(b/a)_{\rm PSF}>0.65$, where $a$
and $b$ are the lengths of the major and minor axes, respectively, of the
shape of a source; (iii) $0.7<{\rm FWHM}/({\rm FWHM})_{\rm PSF}<1.3$; (iv)
PSF-subtracted residual with $<3\sigma$ standard deviation in the difference
image. Furthermore, we imposed the following criterion to exclude moving
objects such as minor planets: (v) Detection at least twice in the difference
images. As a result, we obtained 5213 variable point sources. To evaluate the
completeness of transient detection with image subtraction and our detection
criteria, we randomly injected artificial point sources with various
magnitudes into the observed images and detected them in the difference images
using the same detection criteria. Figure 3 shows the completeness of
transient detection in the difference images within the first (top: Day
1$-$28) and second epochs (bottom: Day 4$-$28). The vertical dotted lines
indicate the mode values of the $5\sigma$ limiting magnitudes, and these
magnitudes were approximately comparable to the magnitude for a completeness
of 23% (Day 1$-$28, $r$), 29% (Day 1$-$28, $z$), 39% (Day 4$-$28, $r$), and
32% (Day 4$-$28, $z$).

Figure 3: Completeness of transient detection in the difference images (top:
Day 1$-$28; bottom: Day 4$-$28). The dashed curves are the detection fractions
of artificial point sources as a function of magnitude. The solid curves are
the detection fractions of the remaining artificial point sources after
applying Criteria (i)$-$(iv). The vertical dotted lines indicate the mode
values of the $5\sigma$ limiting magnitudes.

Figure 4 shows a flowchart of the candidate screening and classification
process for the detected sources. We first checked whether the sources were
associated with known objects (for example, variable stars and AGNs) by
matching these sources with the Pan-STARRS1 (PS1) catalog (Flewelling et al.,
2020) within 1′′. Furthermore, to classify the PS1 objects, we used a flag,
objInfoFlag, that indicates if an object is extended or not. As a result, we
found that 2743 and 1883 sources were associated with stellar-like objects and
extended objects, respectively. The former are likely to have originated from
stellar variabilities, and thus we excluded them as candidates of an EM
counterpart of GW200224_222234. The latter could be variabilities of AGNs, an
EM signal from a BBH coalescence, or something else. The rate of BBH
coalescence may be enhanced in galactic nuclei (for example, McKernan et al.,
2019; Tagawa et al., 2020a, b). The light curves of AGNs are stochastic (for
example, Vanden Berk et al., 2004), and thus a clear difference between AGN
variability and a catastrophic phenomenon such as BBH coalescence is the long-
term variability. However, it was not possible to distinguish them from only
our three-nights’ observation over one month. Therefore, we also excluded the
1883 sources associated with extended objects in this study. Further
examination is left to future studies. After this selection, 587 sources
remained in the sample.

Figure 4: Flowchart of the candidate screening and classification process for
the detected sources after image subtraction. The number in each box
represents the number of remaining sources after each process.

As the next step, we performed a visual inspection to exclude bogus detections
and obtained 223 plausible candidates. Then, we classified these candidates by
whether there are extended objects in the PS1 catalog within an angular
separation, $\theta_{\rm sep}$. Here, we classified them into two groups:
“off-center” candidates ($1^{\prime\prime}\leq\theta_{\rm sep}\leq
15^{\prime\prime}$) and “hostless” candidates ($\theta_{\rm
sep}>15^{\prime\prime}$). We conservatively adopted the threshold
($\theta_{\rm sep}=15^{\prime\prime}$) corresponding to a separation of $\sim
60$ kpc at the distance of GW200224_222234 ($z\sim 0.3$). This separation is
larger than the typical size of galaxies (for example, Shen et al., 2003;
Kauffmann et al., 2003). Then, we investigated off-center candidates with
large separation individually. As a result, we obtained 118 off-center
candidates and 105 hostless candidates.

Furthermore, we calculated the probability $P_{\rm 3D}$ for 118 off-center
candidates. $P_{\rm 3D}$ represents how likely the extended objects were to be
located inside the 3D skymap of GW200224_222234. The definition and
calculation method are described in Tominaga et al. (2018a). We used the $r$\-
and/or $i$-band Kron-magnitudes (rMeanKronMag, iMeanKronMag in the PS1
catalog) of the extended objects to derive the absolute magnitudes of probable
host galaxy in the observer frame. We classified them with a threshold of
$P_{\rm 3D}=0.5$ and obtained 93 likely candidates within the $3\sigma$ error
of the arrival distance of GW200224_222234 ($P_{\rm 3D}\geq 0.5$) and 16
candidates that were likely outside ($P_{\rm 3D}<0.5$). The remaining nine
candidates had no information of the Kron-magnitudes in the $r$\- and
$i$-bands in the PS1 catalog (rMeanKronMag $=$ iMeanKronMag $=-999$) and were
classified as “No Info.” In the following, we investigate the nature of 93
candidates inside the 3D skymap, 16 candidates outside the 3D skymap, nine
candidates without host galaxy information, and 105 hostless candidates.

##  3 Light curve fitting

###  3.1 Method

We performed light curve fitting using transient templates and examples to
exclude common transients such as supernovae (SNe) from the 223 candidates. We
adopted a template set including the transient templates of Type Ia SNe (Hsiao
et al., 2007) and core-collapse SNe (CCSNe; Type Ibc, IIP, IIL, and IIn,
Nugent et al., 2002), and examples of rapid transients (RTs, Drout et al.,
2014) as in Tominaga et al. (2018b). RTs were included in the template set
despite their low rate (Drout et al., 2014) because their nature is still
under debate and they could have a similar timescale to a possible EM
counterpart of BBH coalescence (McKernan et al., 2019). The template set did
not include other types of transients, such as superluminous SNe (Inserra et
al., 2021) and SNe amplified by a gravitational lens (Oguri & Marshall, 2010),
because there is no indication that they are associated with BBH coalescence,
and their rates are sufficiently low that no detection was expected in our
observation (for example, Quimby et al., 2014; Moriya et al., 2019). Note that
the probability that the nature of detected candidates is that of such
transients cannot be excluded.

We derived the light curves of 223 candidates via forced PSF photometry of the
difference images at their location. The difference fluxes of the candidates
were measured using the difference images with the local background
subtracted. Note that the difference fluxes could not be the genuine fluxes of
candidates because the time interval between the science and reference images
was shorter than the typical timescales of SNe. Thus, forced PSF photometry
was also performed for the stacked images without the local background
subtracted. The measured fluxes included the fluxes of their host galaxies as
well as the genuine fluxes of candidates and thus were adopted as the upper
limits of the genuine fluxes of candidates. Then, we performed template
fitting not only to the difference fluxes but also to the upper limits.

We considered a variation in the explosion date and an intrinsic variation in
templates and examples. The light curve templates were derived with an
explosion date $t_{\rm exp}$, redshift $z$, and variations as conducted in
Tominaga et al. (2018b). To account for the variations in supernova
properties, we derived the peak absolute $B$-band magnitude $M_{B}$ using
Equation (4) in Barbary et al. (2012) by considering the stretch $s$, color
$c$, and intrinsic variation $I$ for the template of a Type Ia SN. For CCSNe
and RTs, we parameterized $M_{B}$ and the color excess of their host galaxy
$E_{B-V}$ characterizing the extinction of the host galaxy. We assumed that
the extinction curve of the host galaxy was the same as that in our Galaxy
(Pei, 1992).

To evaluate the difference between the observed light curves and the template,
we defined the following $\xi$ value:

$\displaystyle\xi:=$ $\displaystyle\sum_{i}^{N^{\rm obs}_{\rm
d}}\frac{\left(f^{\rm obs}_{{\rm d},i}-f^{\rm temp}_{{\rm
d},i}(p)\right)^{2}}{{\sigma^{\rm obs}_{{\rm d},i}}^{2}}$ (1)
$\displaystyle+\sum_{j}^{N^{\rm obs}_{\rm s}}\frac{\left(f^{\rm obs}_{{\rm
s},j}-f^{\rm temp}_{{\rm s},j}(p)\right)^{2}}{{\sigma^{\rm obs}_{{\rm
s},j}}^{2}},$

where the subscripts ${\rm d}$ and ${\rm s}$ indicate quantities derived from
the difference images (d) and stacked images (s), respectively, $N^{\rm obs}$
is the number of data points of the observation, $f^{\rm obs}$ and
$\sigma^{\rm obs}$ are the observed flux and its error, respectively, and
$f^{\rm temp}(p)$ is the template flux calculated using a template parameter
set $p$. When $f^{\rm obs}_{{\rm s},j}-f^{\rm temp}_{{\rm s},j}>0$, we set
$\sigma^{\rm obs}_{\rm s}$ to infinity because the flux in the stacked image
is an upper limit of the genuine flux of the candidate.

The $5\sigma$ limiting magnitudes in the difference images correspond to 0.37
$\mu{\rm Jy}$ (Day 1$-$28) and 0.64 $\mu{\rm Jy}$ (Day 4$-$28) in the
$r2$-band flux, and 1.77 $\mu{\rm Jy}$ (Day 1$-$28) and 2.36 $\mu{\rm Jy}$
(Day 4$-$28) in the $z$-band flux. Most of these candidates were detected in
the $r2$-band but not in the $z$-band. If the difference flux was lower than
the flux corresponding to the $5\sigma$ limiting magnitudes, we required the
template flux to be fainter than the flux corresponding to the $5\sigma$
limiting magnitudes.

We applied the Metropolis–Hastings (MH) algorithm, which is a type of Markov-
chain Monte Carlo (MCMC) method, for the template fitting. This method can
derive best-fit template parameter sets without being trapped by the local
minima of $\xi$ by probabilistically moving to a location where $\xi$ is
smaller in a parameter space. We applied the top-hat function as a prior
probability distribution with the same parameter ranges as adopted in Tominaga
et al. (2018b). The ranges were obtained from Barbary et al. (2012) for the
Type Ia SN, Dahlen et al. (2012) for CCSNe, and an assumption of $\pm 1$ mag
variation for RTs, as shown in Table 4. We derived a Markov chain using 50000
samples with each redshift from 0.00 to 0.50 every 0.05 using each template
model (Type Ia SN, 6 CCSNe, and 9 RTs). Each sample of the Markov chain was
proposed from a Gaussian distribution with the standard deviation shown in
Table 4, using the current value as the median. The initial values of each
parameter were set to the central value of each range. We adopted a large
number of burn-in steps to ensure a negligible dependence on the initial
values. To narrow down possible templates, we selected samples with a loose
threshold of $\xi<25$ from all MCMC samples.

Table 4: Template parameters applied in the MCMC method.

Parameter | Range† | SD of sampling‡
---|---|---
Common
Explosion Time $t_{\rm exp}$ (day) | $[-\infty,~{}\infty]$ | $0.1$
Type Ia SN
Color $c$ | $[-0.2,~{}0.8]$ | $0.01$
Stretch $s$ | $[0.6,~{}1.2]$ | $0.01$
Intrinsic Variation $I$ | $[-0.3,~{}0.3]$ | $0.01$
Type IIL normal
Peak Magnitude $M_{B}$ | $[-16.09,~{}-18.05]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Type IIL bright
Peak Magnitude $M_{B}$ | $[-17.92,~{}-19.96]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Type IIP
Peak Magnitude $M_{B}$ | $[-14.43,~{}-18.91]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Type IIn
Peak Magnitude $M_{B}$ | $[-16.98,~{}-20.66]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Type Ibc normal
Peak Magnitude $M_{B}$ | $[-16.09,~{}-18.05]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Type Ibc bright
Peak Magnitude $M_{B}$ | $[-18.46,~{}-20.30]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-10ah)
Peak Magnitude $M_{B}$ | $[-16.63,~{}-18.63]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-10bjp)
Peak Magnitude $M_{B}$ | $[-17.20,~{}-19.20]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-11bbq)
Peak Magnitude $M_{B}$ | $[-18.48,~{}-20.48]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-11qr)
Peak Magnitude $M_{B}$ | $[-18.03,~{}-20.03]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-12bb)
Peak Magnitude $M_{B}$ | $[-15.29,~{}-17.29]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-12bv)
Peak Magnitude $M_{B}$ | $[-18.44,~{}-20.44]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-12brf)
Peak Magnitude $M_{B}$ | $[-17.39,~{}-19.39]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-13ess)
Peak Magnitude $M_{B}$ | $[-17.49,~{}-19.49]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$
Rapid Transient (PS1-13duy)
Peak Magnitude $M_{B}$ | $[-17.77,~{}-19.77]$ | $0.01$
Color Excess $E_{B-V}$ | $[0,~{}1]$ | $0.01$

Note. — [${\dagger}$] Range of top hat distribution for prior distribution.
[${\ddagger}$] Standard deviation of the proposed distribution for sampling.

Additionally, for the candidates associated with the PS1 extended object (that
is, off-center candidates), we confirmed the consistency between the redshifts
of transient templates and the distance information of the host galaxies by
applying the following criteria:

Criterion 1-1

We eliminated templates whose redshifts were outside the redshift error range
($2\sigma$) of the PS1 extended object.

Criterion 1-2

If the PS1 objects exhibited a large angular separation, $\theta_{\rm
sep}>10^{\prime\prime}$, we did not apply Criterion 1-1 because the associated
PS1 objects were potentially not true host galaxies.

Criterion 2

When $\xi>25$ for $\theta_{\rm sep}<10^{\prime\prime}$ with Criterion 1-1, we
allowed the templates at a redshift larger than that of the PS1 extended
objects because a true faint host galaxy may exist behind the PS1 extended
object even if $\theta_{\rm sep}<10^{\prime\prime}$.

Here, we used the photometric redshifts (photo-$z$) from the Sloan Digital Sky
Survey (SDSS; York et al., 2000) catalog as the redshift of the PS1 extended
object if available. If not, we estimated a single-band redshift $z_{\rm
single}$ and its standard deviation $\sigma_{z}$ using the following formulae:

$\displaystyle z_{\rm single}$ $\displaystyle:=$
$\displaystyle\frac{\int^{\infty}_{0}\phi\,A\,z(D)~{}dD}{\int^{\infty}_{0}\phi~{}A~{}dD},$
(2) $\displaystyle\sigma_{z}^{2}$ $\displaystyle:=$
$\displaystyle\frac{\int^{\infty}_{0}\phi\,A\,(z(D)-z_{\rm
single})^{2}~{}dD}{\int^{\infty}_{0}\phi~{}A~{}dD},$ (3)

where $\phi=\phi(\lambda,~{}M)$ is the luminosity function of galaxies at a
rest wavelength $\lambda$ derived from the luminosity functions in the
$UBVRI$-band provided in Ilbert et al. (2005) and Planck cosmology (Planck
Collaboration et al., 2014), $A=A(D)$ is a surface area observed at a distance
$D$, $M=M(D;~{}m_{j})$ is an absolute magnitude derived with the $j$th-band
($r$ or $i$) apparent magnitude $m_{j}$ at $D$ in the observer frame, and
$\lambda=\lambda(D;~{}\lambda_{j})$ is the rest wavelength redshifted from an
observed wavelength $\lambda_{j}$ at $D$. We also estimated $z_{\rm single}$
and $\sigma_{z}$ for candidates with the SDSS photometric redshift (Table 5).
Their comparison illustrates their approximate consistency (Figure 5). Here,
we corrected the Galactic extinction (Schlafly & Finkbeiner,
2011)333http://irsa.ipac.caltech.edu/applications/DUST/ when we derived
absolute magnitudes from apparent magnitudes.

Figure 5: Comparison between the SDSS photometric redshifts and the redshift
estimated with Eqs. (2) and (3).

###  3.2 Results of fitting analysis

Among the 223 candidates for which we performed light curve fitting, 200
candidates were consistent with the templates of SNe, three candidates were
consistent only with the template of RTs, and the remaining 20 candidates were
not consistent with any templates or examples. Table 5 shows the coordinates,
$\theta_{\rm sep}$, $P_{\rm 3D}$, and probable transient templates of the 200
candidates consistent with the templates. A fitting result of JGEM20ewa with
$\xi=0.6$ is shown as an example of the candidates consistent with the
transient template (Figure 6). We concluded that the most probable origin of
this candidate is a Type Ia SN for a redshift of $z=0.25$. The best-fit light
curves are shown as solid lines in Figure 6. The top and bottom panels show
the PSF flux measured in the stacked and difference images, respectively.

Table 5: Information of the candidates that are consistent with the transient
templates or examples.

Name | Coordinate (J2000.0) | $\theta_{\rm sep}^{{\dagger}}$ | $z_{\rm single}$ | $\sigma_{z}$ | $z_{\rm SDSS}$ | $\sigma_{\rm SDSS}$ | $P_{\rm 3D}^{{\dagger}}$ | Probable transient template sets
---|---|---|---|---|---|---|---|---
| R.A. (HH:MM:SS.ss) | Decl. (DD:MM:SS.s) | (′′) |  |  |  |  |  |
JGEM20abe | 11:43:12.36 |  $-$15:39:25.6 | 8.4 | 0.37 | 0.13 | – | – | 77 | CCSNe (Ibc and IIL)
JGEM20abf | 11:43:11.52 |  $-$15:28:16.3 | – | – | – | – | – | – | CCSNe (Ibc, IIn and IIL)
JGEM20acf | 11:41:48.24 |  $-$14:55:29.3 | – | – | – | – | – | – | CCSNe (Ibc, IIn and IIL)
JGEM20adp | 11:50:40.82 |  $-$15:14:02.8 | – | – | – | – | – | – | Type Ia, CCSNe (Ibc, IIn and IIL)
JGEM20adq | 11:50:42.55 |  $-$15:13:57.4 | – | – | – | – | – | – | CCSNe (Ibc, IIn and IIL)
JGEM20ads | 11:50:52.06 |  $-$15:10:58.4 | 2.4 | 0.15 | 0.06 | – | – | 97 | Type Ia, CCSNe (Ibc, IIn, IIP and IIL)
JGEM20aej | 11:49:46.03 |  $-$15:01:31.8 | – | – | – | – | – | – | Type Ia, CCSNe (Ibc, IIn, IIP and IIL), RTs (13duy)
JGEM20afe | 11:48:22.75 |  $-$15:44:46.3 | 12.9 | 0.33 | 0.12 | – | – | 93 | CCSNe (Ibc and IIL)
JGEM20cvb | 11:40:58.56 |  $-$11:26:57.5 | 2.0 | – | – | – | – | – | Type Ia
JGEM20ewa | 11:37:13.22 |  $-$8:21:04.7 | 3.6 | 0.15 | 0.06 | – | – | 82 | Type Ia
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$

Note. — This table is published in its entirety in the machine-readable
format. A portion is shown here for guidance regarding its form and content.
[${\dagger}$] We indicate “–” if a candidate is not associated with PS1
extended objects.

Figure 6: Example of the fitting results. The light curve of JGEM20ewa (points
with error bars), and the best-fit light curve of the Type Ia SN template
(solid curves); $z=0.25$, $t_{\rm exp}=0.38$ d, $c=-0.15$, $s=0.60$, and
$I=-0.02$. The vertical axes of the top and bottom panels indicate the upper
limits measured in the stacked images and the difference fluxes measured in
the difference images, respectively.

Among the 200 candidates consistent with the templates of SNe, 105 candidates
were associated with the PS1 extended objects (that is, off-center
candidates), and the remaining 95 candidates had no PS1 extended objects
within $<15^{\prime\prime}$.

Criterion 1-1 was applied to 102 out of the 105 off-center candidates, and
these candidates were excluded from the final candidates. Two of the remaining
candidates (JGEM20fvn and JGEM20gqu) were associated with the PS1 extended
objects of $z=0.56\pm 0.05$ and $0.34\pm 0.04$, respectively. The positions of
the candidates and PS1 objects are shown in Figure 7; left (JGEM20fvn) and
middle (JGEM20gqu). Because JGEM20fvn and JGEM20gqu exhibited $\theta_{\rm
sep}=13^{\prime\prime}.8$ and $11^{\prime\prime}.0$ respectively, Criterion
1-2 was applied to them. We found that the templates of Type Ibc CCSNe with
$z=0.2$ and $0.15$ reproduced the light curves of JGEM20fvn and JGEM20gqu,
respectively. These results indicated that the PS1 objects associated with
JGEM20fvn and JGEM20gqu were unlikely to be true host galaxies and were likely
to be Type Ibc CCSNe with $z=0.2$ and $0.15$. Thus, we ruled out these two
objects from the final candidates. Criterion 2 was applied to the remaining
one candidate (JGEM20hgo). With Criterion 2, we found that a template of Type
Ibc CCSN with $z=0.1$ can reproduce the light curve of JGEM20hgo. The redshift
of the template was larger than the redshift of the PS1 object associated with
JGEM20hgo ($z=0.03\pm 0.01$). Therefore, we also ruled out JGEM20hgo from the
final candidates because it could be a Type Ibc CCSN behind the PS1 object.
The stacked and difference images of JGEM20hgo and the associated PS1 object
are shown in right panels of Figure 7. The parameters $\theta_{\rm sep}$ and
$P_{\rm 3D}$ of candidates consistent with the transient templates are shown
in Table 5.

Figure 7: Observed (top) and difference (bottom) $r2$-band images of
JGEM20fvn, JGEM20gqu, and JGEM20hgo on February 25, 2020. Each image displays
the $20^{\prime\prime}\times 20^{\prime\prime}$ region. Although the PS1
extended objects (plus-shaped markers) are associated with the candidates
(reverse L-shaped markers), these are unlikely to be the true host galaxies of
the candidates.

Among the 20 candidates inconsistent with the templates, one candidate
(JGEM20cvb) exhibited a large $\xi$ value of 2175. However, it apparently
matched the template light curve, as shown in Figure 8. This large $\xi$ could
be attributed to the underestimation of photometric errors or the
uncertainties on the template fluxes. For instance, $\xi$ was improved to
$\sim 2.4$ by assuming a template difference flux error of 0.1 mag. Thus, we
concluded that the origin of JGEM20cvb is a Type Ia SN and excluded it from
the final candidates.

We also performed template fitting without the single-band photometric
redshifts for the candidates associated with extended objects without SDSS
photometric redshifts. We found that the redshifts of the best-fit templates
were consistent for most of the candidates and confirmed that the consistent
templates did not change for candidates other than the candidates discussed
above, that is, JGEM20fvn, JGEM20gqu, and JGEM20hgo.

Figure 8: Light curve of JGEM20cvb (points with error bars), and the best-fit
light curve of the Type Ia SN template (solid curves); $z=0.05$, $t_{\rm
exp}=11.99$ d, ${\rm color}=-0.20$, $s=1.07$, and $I=0.05$.

As a result, 22 objects, that is, three objects consistent only with RTs and
19 objects inconsistent with all templates and examples, remained as the final
candidates of the optical counterpart of GW200224_222234. Their coordinates,
$\theta_{\rm sep}$, $P_{\rm 3D}$, and fitting results are shown in Table 6 and
Figure 9.

Table 6: Information on candidates inconsistent with the templates of SNe and
their associated PS1 extended objects.

Name | Coordinate (J2000.0) | $\theta_{\rm sep}^{{\dagger}}$ | $P_{\rm 3D}^{{\dagger}}$ | Best fitted transient templates
---|---|---|---|---
| R.A. (HH:MM:SS.ss) | Decl. (DD:MM:SS.s) | (′′) |  |
JGEM20acc | 11:41:52.70 |  $-$15:08:22.9 | – | – | CCSN Type Ibc ($\xi\sim 33$)
JGEM20aeh | 11:50:20.04 |  $-$15:12:54.0 | – | – | CCSN Type Ibc ($\xi\sim 55$)
JGEM20boq | 11:39:19.68 |  $-$12:44:10.0 | 8.9 | 75 | CCSN Type Ibc ($\xi\sim 334$)
JGEM20bsr | 11:37:26.83 |  $-$12:05:34.4 | – | – | RT 10bjp ($\xi\sim 16$)
JGEM20cpd | 11:35:05.23 |  $-$11:06:33.8 | 14.9 | 94 | CCSN Type Ibc ($\xi\sim 281$)
JGEM20dig | 11:41:48.72 |  $-$11:45:02.9 | 11.4 | 78 | CCSN Type Ibc ($\xi\sim 26$)
JGEM20dte | 11:27:56.54 |  $-$9:05:45.6 | – | – | CCSN Type Ibc ($\xi\sim 159$)
JGEM20ejv | 11:47:24.98 |  $-$9:40:32.5 | 6.3 | 75 | CCSN Type Ibc ($\xi\sim 62$)
JGEM20ekz | 11:45:32.81 |  $-$9:42:04.7 | – | – | CCSN Type Ibc ($\xi\sim 84$)
JGEM20env | 11:33:13.03 |  $-$8:10:35.0 | 14.3 | – | CCSN Type IIL ($\xi\sim 331$)
JGEM20enw | 11:32:50.98 |  $-$7:56:35.2 | 14.9 | 74 | RT 13duy ($\xi\sim 1$)
JGEM20eso | 11:28:42.14 |  $-$7:48:41.4 | – | – | RT 13duy ($\xi\sim 10$)
JGEM20fci | 11:44:18.89 |  $-$8:08:34.4 | – | – | Type Ia SN ($\xi\sim 1137$)
JGEM20fud | 11:31:23.86 |  $-$5:55:16.0 | 7.0 | 89 | Type Ia SN ($\xi\sim 1562$)
JGEM20fxo | 11:40:19.39 |  $-$6:36:34.9 | – | – | CCSN Type Ibc ($\xi\sim 77$)
JGEM20fyv | 11:40:00.36 |  $-$6:01:27.5 | 5.2 | 99 | CCSN Type Ibc ($\xi\sim 422$)
JGEM20gdm | 11:45:52.56 |  $-$6:13:49.4 | 1.6 | 91 | CCSN Type IIL ($\xi\sim 27$)
JGEM20hdq | 11:32:35.04 |  $-$3:56:15.7 | 3.0 | 60 | RT 10bjp ($\xi\sim 187$)
JGEM20hea | 11:32:20.45 |  $-$3:40:50.9 | 3.7 | 60 | CCSN Type Ibc ($\xi\sim 394$)
JGEM20hen | 11:32:21.91 |  $-$3:09:21.6 | 7.3 | 98 | Type Ia SN ($\xi\sim 348$)
JGEM20hfc | 11:31:30.55 |  $-$3:45:54.7 | 2.7 | 94 | CCSN Type Ibc ($\xi\sim 87$)
JGEM20hkw | 11:35:25.49 |  $-$4:11:36.6 | – | – | CCSN Type Ibc ($\xi\sim 37$)

Note. — [${\dagger}$] We indicate “–” if a candidate is not associated with
PS1 extended objects.

Figure 9: Light curves of 19 candidates inconsistent with the templates of SNe
and three candidates consistent only with RTs. The best-fit (but bad-fit for
the former 19 candidates) templates or examples are also shown for comparison.
The vertical axis of each panel indicates the upper limits measured in the
stacked images (top) and the difference fluxes measured in the difference
images (bottom).

##  4 Spectroscopic observations with the GTC/OSIRIS

Optical spectroscopic observations were performed using the OSIRIS instrument
mounted on the 10.4- m GTC telescope located at the Roque de los Muchachos
Observatory (Canary Islands, Spain). Five targets were selected for
spectroscopic follow-ups to determine the spectroscopic redshift (spec-$z$) of
the probable host galaxies of the candidates. The selected targets were the
PS1 extended objects associated with JGEM20fud, JGEM20fyv, JGEM20gdm,
JGEM20hdq, and JGEM20hen, which were sufficiently bright in the $z$-band for
short exposures. Owing to observation time constraints, we targeted only these
five objects.

The observations were performed on February 10, 2021. The instrumental setup
was the spectroscopic long-slit mode, with a R1000R grism and slit width of 1
arcsec. Each target observation was divided into three exposures of 400 s
(1200 s in total for JGEM20fud, JGEM20fyv, and JGEM20hen) and 600 s (1800 s in
total for JGEM20gdm and JGEM20hdq). In addition to the targets, the
spectrophotometric standard star G191-B2B was also observed using the same
instrumental setup and observation conditions. Standard calibration images for
bias, flat field, and the calibration lamp (HgAr+Xe+Ne) were also taken during
the same night. The data were reduced using standard procedures for bias
subtraction, flat-field correction, flux calibration, and atmospheric
extinction correction.

The flux calibrated spectra are shown in Figures 10 and 11. The spec-$z$ was
derived based on the identification of spectral lines. In addition to the
redshifted spectral features from the galaxies, we identified a possible
intervening system in the line of sight of JGEM20fud and JGEM20hen, which is
consistent with an Mg absorption line located at $z=0.33$. We show the
coordinate, the $z$-band magnitude, estimated redshift, and spec-$z$ of the
target galaxies in Table 7. If an SDSS photometric redshift is not available,
we show $z_{\rm single}$ and $\sigma_{z}$ calculated with Equations (2) and
(3). However, note that these have large uncertainties because they are
calculated from only a single band magnitude.

Figure 10: Optical spectra of candidates JGEM20fud, JGEM20fyv, JGEM20gdm,
JGEM20hdq, and JGEM20hen observed with GTC. The identified spectral features
are marked by vertical dotted lines. The Na interstellar absorption feature is
denoted by a solid vertical line. The emission lines identified as intervening
systems are marked by dashed lines.

Figure 11: Continuation of Figure 10.

Figure 12 shows the 2D location and distance of the five candidates. We
evaluated their associations with GW200224_222234 by comparing the 3D location
of their probable host galaxies with the 3D skymap of GW200224_222234. Here,
we define the probability of the association as follows:

$\displaystyle P_{\rm assoc}(x_{\rm gal}):=1-\int_{f_{x}>f_{x_{\rm
gal}}}f_{x}\,dV_{x},$ (4)

where $f_{x}$ is the probability distribution function (PDF) at the 3D
coordinates $x$, and $x_{\rm gal}$ is the location of the host galaxy. The
integration in Equation (4) represents the cumulative probability integrated
over the region with $f_{x}>f_{x_{\rm gal}}$. The PDF was obtained as
$f_{x}=\rho({\rm R.A.},~{}{\rm Decl.})\,N(D;~{}\mu,~{}\sigma^{2})$ using the
2D probability distribution $\rho=\rho({\rm R.A.},~{}{\rm Decl.})$ and normal
distribution $N(D;~{}\mu,~{}\sigma^{2})$ of the distance with the mean
$\mu=\mu({\rm R.A.},~{}{\rm Decl.})$ and standard deviation
$\sigma=\sigma({\rm R.A.},~{}{\rm Decl.})$ provided for each direction. To
perform the integration of Equation (4) numerically, we discretized $N(D)$
into 5000 bins from 0 to 50000 Mpc. The discretization of $\rho$ was performed
through pixelization using HEALPix. $P_{\rm assoc}$ is small when the host
galaxy is located outside the highly probable region. We set the threshold to
$P_{\rm assoc}=0.01$. Table 7 also shows $P_{\rm assoc}$ obtained for the five
targets. Although JGEM20gdm, JGEM20hdq, and JGEM20hen were located outside the
highly probable region ($P_{\rm assoc}<0.01$) and are unlikely to be related
to GW200224_222234, JGEM20fud and JGEM20fyv were close to the highly probable
region ($P_{\rm assoc}=0.76$ and $0.07$, respectively) and thus are possible
counterparts of GW200224_222234.

Figure 12: (Left) 2D location of the five candidates observed with the
GTC/OSIRIS and the localization skymap (IMRPhenomXPHM model) released in the
GWTC-3 catalog. The dots are the locations of the candidates. The gray contour
line indicates the 90% credible region. (Right) Distance of each candidate and
probability distribution at their location as a function of distance. The dots
are the distances of the candidates. The horizontal solid and dashed lines
indicate the mean and standard deviation of the probability distribution,
respectively. Table 7: Information of the target galaxies.

Name | Coordinate of the target galaxy (J2000.0) | ${m_{z}}^{\dagger}$ | Estimated | spec-$z$ | $P_{\rm assoc}$
---|---|---|---|---|---
| R.A. (HH:MM:SS.ss) | Decl. (DD:MM:SS.s) | (mag) | Redshift |  |
JGEM20fud | 11:31:24.32 |  $-$05:55:17.8 | 18.7 | $0.22\pm 0.03^{\ddagger}$ | 0.357 | $0.76$
JGEM20fyv | 11:40:00.36 |  $-$06:01:22.1 | 19.0 | $0.26\pm 0.07^{*}$ | 0.391 | $0.07$
JGEM20gdm | 11:45:52.63 |  $-$06:13:48.4 | 19.9 | $0.29\pm 0.09^{*}$ | 0.398 | $1\times 10^{-3}$
JGEM20hdq | 11:32:35.03 |  $-$03:56:18.6 | 20.1 | $0.56\pm 0.03^{\ddagger}$ | 0.632 | $1\times 10^{-8}$
JGEM20hen | 11:32:22.41 |  $-$03:09:22.9 | 18.4 | $0.36\pm 0.03^{\ddagger}$ | 0.438 | $4\times 10^{-3}$

Note. — [${\dagger}$] $z$-band magnitude of the target galaxy from the PS1
catalog. [${\ddagger}$] Photometric redshift $z_{\rm SDSS}$ from the SDSS
catalog. [∗] $z_{\rm single}$ and $\sigma_{z}$ calculated with Eqs. (2) and
(3).

##  5 Discussion

###  5.1 Contamination from supernovae

Next, we evaluated our detection criteria and screening process by comparing
our results with the expected number of SN detections. We estimated the
expected number of SN detections by summing up mock-SN samples weighted with
cosmological histories of SN rates using the observation depth, as in Niino et
al. (2014) and Ohgami et al. (2021). We assumed the SN rate of Okumura et al.
(2014) and Dahlen et al. (2012) for the Type Ia SN and CCSN, respectively. For
the Type Ia SN, the SN light curves were generated from the evolution of the
SN spectrum provided by Hsiao et al. (2007). For the CCSN, the light curves
were generated from the templates provided by Nugent et al. (2002)444The light
curve templates of Type Ia SNe and CCSNe are available at the website
https://c3.lbl.gov/nugent/nugent_templates.html.. We referred to the
luminosity distributions of SNe in Barbary et al. (2012) and Dahlen et al.
(2012).

Because our observation depth in the $z$-band was shallower than that in The
$r$-band, most of candidates were detected only in the $r$-band. Thus, we
sampled mock SNe whose brightness was decaying in the $r$-band, assuming the
reference images were taken 27 d after detection because the variation in the
SN magnitude during 3 d between Days 1 and 4 ($<0.1$ mag) was negligible
compared to the standard deviation of the $5\sigma$ limiting magnitude of our
observation. The number density as a function of limiting magnitude is shown
in Figure 13 as a dashed curve. The colored area represents the $\pm 50\%$
error from the $1\sigma$ error of the CCSN rate density. In the actual
observation, the detection completeness for fainter sources was lower (Figure
3). Therefore, we considered two detections based on the completeness of Day
1$-$28 and Day 4$-$28\. The solid curve in Figure 13 is the number density
derived by considering the completeness. The dot with error bars was estimated
with the number of transients consistent with SNe, $N_{\rm SN}=201$, and our
observed area of $S=56.6$ deg2. The vertical error bar was defined as
$\sqrt{N_{\rm SN}}/S$ by assuming that the number followed a Poisson
distribution with the expected value of $N_{\rm SN}$. The horizontal error bar
is the $1\sigma$ error of the $5\sigma$ limiting magnitude in the difference
images on Day 4 in the $r2$-band.

Our result is consistent with the expected number density of SN detection
within a $1\sigma$ error. This verified that our detection criteria and
screening process can be used to identify SNe with high completeness. Note
that our result was located at the $1\sigma$ lower end of the expected value.
This might be because we ignored nuclear transients in this search.

Figure 13: Expected number density of SN detection in the $r$-band as a
function of limiting magnitude. The colored area is the $\pm 50\%$ error
originating from the $1\sigma$ error of the core-collapse-SN rate density. The
dashed curve was obtained with 100% completeness, and the solid curve was
obtained with the actual completeness of our survey, shown in Figure 3, under
the assumption of twice detection. The dot with error bars was obtained using
the number density of the transients consistent with the templates of SNe. The
vertical error bar is defined as $\sqrt{N_{\rm SN}}/S$ by assuming that the
number follows a Poisson distribution with the expected value of $N_{\rm SN}$.
The horizontal error bar is the $1\sigma$ error of the $5\sigma$ limiting
magnitude in the difference images on Day 4 in the $r2$-band.

###  5.2 Implication in BBH merger models accompanied by EM emissions

We finally found 19 candidates of the optical counterpart of GW200224_222234
by performing light curve fitting and spectroscopic observations. Owing to the
limitation of observations, we must exclude sources located at the center of
galaxies, and thus these final candidates cannot be considered variabilities
of galactic nuclei, such as a kicked BBH merger in an accretion disk shown by
McKernan et al. (2019). Yamazaki et al. (2016) reported a BBH merger
potentially accompanied by $\gamma$-ray emissions caused by a relativistic
outflow. They showed that optical afterglow peaks appear $\sim 4.5\times
10^{4}$ s after GW detection (almost the same as the start of our first epoch
observation) if the ambient matter density is 0.01 cm-3; however, the peak
flux is expected to be lower than our observation limits (0.37 $\mu$Jy in
$r$-band) at a distance of 1710 Mpc. Therefore, such afterglow could not be
detected in our search.

Among the 19 candidates, 16 candidates inconsistent with all templates and
examples exhibited the increasing difference flux in one band and the
decreasing difference flux in another band (for example, JGEM20ejv), or a more
rapid decline in the $z$-band than in the $r$-band (for example, JGEM20boq).
If we perform further observations to obtain deeper images than the current
reference images after the candidates sufficiently decay, we can measure the
genuine fluxes of these candidates instead of the difference fluxes and hence
identify their natures.

The final 19 candidates are potentially unrelated to GW200224_222234. If there
are no counterparts of GW200224_222234, the upper limits of optical luminosity
are $\nu L_{\nu}<5.2^{+2.4}_{-1.9}\times 10^{41}$ erg s-1
($9.1^{+4.1}_{-3.3}\times 10^{41}$ erg s-1) and $\nu
L_{\nu}<1.8^{+0.8}_{-0.6}\times 10^{42}$ erg s-1 ($2.4^{+1.1}_{-0.9}\times
10^{42}$ erg s-1) on Day 1 (Day 4) in the $r2$-band and $z$-band,
respectively. Note that these upper limits are comparable with the luminosity
of a possible EM counterpart of the BBH event GW190521 ($10^{42}$ erg s-1,
McKernan et al. 2019).

###  5.3 Comparison with sources detected by other observations

For GW200224_222234, the Neil Gehrels Swift Observatory performed near-
UV/X-ray observations covering an area corresponding to 79.2%/62.4% of the GW
probability region and reported their detection of eight X-ray sources and
three near-UV sources (Klingler et al., 2021; Oates et al., 2021), whereas no
significant signal was detected in very high-energy $\gamma$-ray (Abdalla et
al., 2021). We searched our sources to find out whether these sources were
detected. The source JGEM20aoz was located at ${\rm R.A.}=11^{\rm h}41^{\rm
m}41^{\rm s}.81$, ${\rm Decl.}=-14^{\circ}07^{\prime}54^{\prime\prime}.78$
(J2000.0), 4.8 arcsec off center of the error circle of X-ray Source 7.
JGEM20aoz was detected four times (two epochs and two bands) and exhibited a
magnitude of 21.75 mag (21.59 mag) and 21.43 mag (21.25 mag) in the $r$\- and
$z$-bands, respectively, in the first (second) epoch difference images.
However, we classified it as a star-like object and ruled it out because the
coordinates matched that of a point source found in the PS1 catalog. In the
SIMBAD archival database, a BL Lac-type object (2FGL J1141.7-1404) is located
at a position consistent with JGEM20aoz. Therefore, we conclude that JGEM20aoz
is likely an AGN.

Our observation area also covered part of the area observed by the GW program
in the Dark Energy Survey (DESGW). They conducted observations of
GW200224_222234 using the Dark Energy Camera (DECam) in the $i$-band on
February 24, 25, 27, and March 5, 2020, with $10\sigma$ limiting magnitudes of
23.17, 23.19, 23.49, and 23.03, respectively, and reported eight transients
(Morgan et al., 2020). Comparing the transients found by the DESGW with our
candidates, JGEM20gfn matched the coordinates of AT2020ehw within
$1^{\prime\prime}$. JGEM20gfn was classified as an SN via light curve fitting.
We could not detect AT2020eho and AT2020eht because they were located at the
CCD gap of our difference images. Sources were detected at the location of
AT2020ehp and AT2020ehq in our difference images, but they were rejected
because of the large elongation, $(b/a)/(b/a)_{\rm PSF}=0.6$ and $0.5$,
respectively. Although no source was found at the position of AT2020ehy in
both the stacked and difference images, AT2020ehv and AT2020ehr were observed
in the stacked images but could not be detected in the difference images with
our criterion $(S/N)_{\rm PSF}>5$.

###  5.4 Comparison with a kilonova light curve

Although GW200224_222234 was classified as a BBH merger, we compared our
observation result with a kilonova model based on a radiative transfer
simulation with an ejecta mass $M_{\rm ej}=0.05\,M_{\odot}$ and an electron
fraction $Y_{\rm e}=0.30-0.40$ (Banerjee et al., 2020). This kilonova model
can explain the observed early multi-color light curve of GW170817/AT2017gfo.
The peak of bolometric luminosity approximately scales with $M_{\rm ej}$ to
the power of 0.35 (see Fernández & Metzger, 2016; Metzger, 2019; Tanaka,
2016). $Y_{\rm e}$ influences the light curves through the opacity. In the
low-$Y_{\rm e}$ case ($Y_{\rm e}<0.25$), the kilonova ejecta becomes
Lanthanide rich and the opacity becomes higher, and thus the light curves in
the early phase can be fainter by $2-3$ mag than that in the adopted model of
$Y_{\rm e}=0.30-0.40$.

Figure 14 shows a comparison between the kilonova light curve and the
magnitudes of two candidates whose host galaxies’ spectroscopic redshifts are
consistent with the 3D skymap of GW200224_222234 (JGEM20fud and JGEM20fyv).
All observed magnitudes were brighter and more slowly decaying than the
kilonova at the distance of the host galaxies. They did not exhibit rapid
decline in short-wavelength components, being a major feature of the kilonova
model; thus, they are unlikely to be kilonovae.

Figure 14: Comparison of light curves between the candidates whose host
galaxies’ spectroscopic redshifts are consistent with the 3D skymap of
GW200224_222234 and a kilonova model with an ejecta mass $M_{\rm
ej}=0.05\,M_{\odot}$ and an electron fraction $Y_{\rm e}=0.30-0.40$ (Banerjee
et al., 2020). We assume the kilonovae are located at the distance of each
candidate.

###  5.5 Future prospects

Here, we discuss the future prospects of follow-ups for kilonova events with
the Subaru/HSC in the era of next-generation GW interferometers, such as an
optimal upgrade of the LIGO facilities, known as “Voyager.” If a BNS merger
occurs at the distance of GW200224_222234 ($\sim$1710 Mpc), can we detect it?

Figure 15 shows the $r$-, $i$-, and $z$-band light curves under an assumption
that the kilonova is located at the estimated distance of GW200224_222234. The
colored area is the uncertainty caused by the 1$\sigma$ credible range of the
distance. The horizontal lines with the upward arrows are the 5$\sigma$
limiting magnitudes measured in the difference images in each epoch in each
band ($r$ and $z$). In the $r$-band, although the observation depth in the
first epoch was sufficiently deep to detect the rising phase of the light
curve, the model light curve decreased rapidly and was significantly fainter
than the limiting magnitude in the second epoch. The observation depths in the
$z$-band did not reach the brightness in either epoch. Thus, if a BNS merger
occurs at the distance of GW200224_222234, we might be able to detect the
light curve in the early phase in the $r$-band using the power of a wide-field
survey with an 8- m telescope.

Based on the follow-up of GW200224_222234, the number of sources detected only
on Day 1 and only in the $r2$-band (excluded using detection criterion (v))
reached 20137 555This number includes bogus detections because we did not
perform visual inspection.. It is difficult to determine their origin using
only single epoch and single band data.

Next, we consider what will improve the detection efficiency of future
surveys. (1) If we adopt a high-cadence (one day) and continuous (over three
days) observation, the kilonova at this distance will be detected on Days 1
and 2 and not detected on Day 3 in the $r$-band. This will illustrate the
rapidly evolving nature of the kilonova. (2) If we adopt the $i$-band instead
of the $z$-band, the $i$-band observation with Subaru/HSC will reach $\sim 24$
mag with 1 min of exposure and may detect the kilonova on Day 2. This will
constrain the color of the kilonova. These considerations illustrate the power
of a dedicated wide-field search with an 8m-class telescope, such as
Subaru/HSC and Rubin/Legacy Survey of Space and Time (LSST).

Figure 15: Comparison between a kilonova model located at the estimated
distance of the GW200224_222234 and the limiting magnitudes of our
observation. The solid curves are the light curves of the kilonova model, the
same as in Figure 14. The colored area represents the uncertainty caused by
the $1\sigma$ credible range of the distance.

##  6 Summary & Conclusion

The BBH coalescence GW200224_222234 was detected by the LIGO/Virgo
collaboration with their three detectors on February 24, 2020. We performed
ToO observations using the Subaru/HSC in the $r2$\- and $z$-bands during three
epochs: February 24, 28, and March 23, 2020. We selected the observation area
from the high-probability region in the preliminary localization skymap
covering 56.6 deg2. The integrated probability reaches 91% in the localization
skymap released in the GWTC-3 catalog.

We searched for the optical counterpart using the image subtraction technique.
We adopted the images taken on the third epoch as the reference images and
obtained the difference images of the first and second epochs. After screening
for the sources detected in the difference images via matching with the PS1
catalog and visual inspection, we found 223 candidates. We could not include
sources located at the galactic center owing to the limitations of
observation. Subsequently, we classified these candidates using their angular
separation from the nearby extended object and distance estimated from the
photometric data. Additionally, we investigated their nature using light curve
fitting with the transient template set. We adopted the templates of the Type
Ia SN and CCSN and examples of RTs and then found 201 candidates consistent
with the SN templates and likely not related to GW200224_222234.

To measure the spectroscopic redshifts of the probable host galaxies of the
final candidates, we also performed spectroscopic observations using the
GTC/OSIRIS for the extended PS1 objects associated with the five candidates.
We found that two targets (JGEM20fud and JGEM20fyv) are likely to be located
inside the highly probable region. The other targets are outside the highly
probable region and not related to the GW event.

As a result, we found 19 candidates as possible candidates of the optical
counterpart of GW200224_222234. The light curves of three candidates were
consistent only with those of RTs, and the light curves of the other 16
candidates were inconsistent with all transients and examples. These 19
candidates have a potential for being unrelated to GW200224_222234; however,
we could not establish their nature because of the lack of spectroscopic
observations for the candidates. If there is no counterpart of GW200224_222234
in the 19 final candidates, the upper limits of optical luminosity are
evaluated as $\nu L_{\nu}<5.2^{+2.4}_{-1.9}\times 10^{41}$ erg s-1
($9.1^{+4.1}_{-3.3}\times 10^{41}$ erg s-1) and $\nu
L_{\nu}<1.8^{+0.8}_{-0.6}\times 10^{42}$ erg s-1 ($2.4^{+1.1}_{-0.9}\times
10^{42}$ erg s-1) on Day 1 (Day 4) in the $r2$-band and $z$-band,
respectively, from the $5\sigma$ limiting magnitudes of our observation. These
upper limits are comparable with the brightness of a possible EM counterpart
of the BBH event GW190521 ($10^{42}$ erg s-1, McKernan et al. 2019).

We evaluated our detection criteria and screening process by comparing our
result with the expected number of SN detections, which was estimated from a
cosmic SN rate. The number density of the candidates consistent with the SN
templates was consistent with the expected number density within a $1\sigma$
error. This indicated that our method could identify SNe with high
completeness.

We also compared our sources with those found by the Swift and DESGW and found
that some of our sources were associated with their sources. However, these
identical sources are unlikely to be the optical counterpart of the GW.

Additionally, we discuss the implications of our result in several models of a
BBH merger accompanied by EM emissions. The kicked BBH merger in an accretion
disk reported by McKernan et al. (2019) cannot be a possible nature of the 16
candidates inconsistent with all transients because we excluded nucleus
transients. The relativistic outflow model reported in Yamazaki et al. (2016)
also cannot be a possible nature of the candidates because the expected flux
is fainter than the $5\sigma$ limiting magnitude of our observations. It is
important to perform further observations deeper than the observations on Day
28 to identify their nature.

We also compared the light curve of the two candidates with the kilonova light
curve by adopting a radiative transfer simulation that could explain the
multi-color light curve of GW170817/AT2017gfo. These candidates likely did not
originate from the same kilonova as AT2017gfo because these were brighter and
more slowly decaying than the kilonova. We also discuss future prospects by
comparing our observation depths with the kilonova light curve. If we perform
a high-cadence (one day) and continuous (over three days) observation, we can
reveal the rapidly evolving nature of transients. Moreover, our observation
can detect a kilonova or a possible EM counterpart of the BBH event GW190521
(McKernan et al., 2019) at the distance of GW200224_222234. This demonstrates
the power of 8m-class telescopes, such as Subaru/HSC and Rubin/LSST.

These collaborative observations were possible thanks to the leadership of the
National Astronomical Observatory of Japan (NAOJ) and the Instituto de
Astrofísica de Canarias (IAC). This research was based in part on data
collected at the Subaru Telescope, which is operated by the NAOJ, and in part
on observations made using the GTC telescope in the Spanish Observatorio del
Roque de los Muchachos of IAC under Director’s Discretionary Time. The
spectroscopic data were obtained using the instrument OSIRIS, built by a
Consortium led by the Instituto de Astrofísica de Canarias in collaboration
with the Instituto de Astronomía of the Universidad Autónoma de México. OSIRIS
was funded by GRANTECAN and the National Plan of Astronomy and Astrophysics of
the Spanish Government. We are grateful to the staff of the Subaru Telescope
for their help with the observations of this study. We are honored and
grateful for the opportunity to observe the Universe from Maunakea, which has
cultural, historical, and natural significance in Hawaii. This study was
supported by MEXT KAKENHI (JP17H06363) and JSPS KAKENHI (JP19H00694,
JP20H00158, JP20H00179, JP21H04997). JBG, FP, JAAP, IPF, and TMD acknowledge
financial support from the Spanish Ministry of Science and Innovation (MICINN)
through the Spanish State Research Agency, under Severo Ochoa Program
2020-2023 (CEX2019-000920-S) and the projects PID2019-107988GB-C22 and
PID2019-105552RB-C43. We would like to thank Editage (www.editage.com) for
English language editing.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abdalla et al. (2021) Abdalla, H., Aharonian, F., Ait Benkhali, F., et al. 2021, ApJ, 923, 109, doi: 10.3847/1538-4357/ac2e04
  * Alard (1999) Alard, C. 1999, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/9903111
  * Alard & Lupton (1998) Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325, doi: 10.1086/305984
  * Andreoni et al. (2017) Andreoni, I., Ackley, K., Cooke, J., et al. 2017, Publications of the Astronomical Society of Australia, 34, e069, doi: 10.1017/pasa.2017.65
  * Arcavi et al. (2017) Arcavi, I., McCully, C., Hosseinzadeh, G., et al. 2017, ApJ, 848, L33, doi: 10.3847/2041-8213/aa910f
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 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, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Banerjee et al. (2020) Banerjee, S., Tanaka, M., Kawaguchi, K., Kato, D., & Gaigalas, G. 2020, ApJ, 901, 29, doi: 10.3847/1538-4357/abae61
  * Barbary et al. (2012) Barbary, K., Aldering, G., Amanullah, R., et al. 2012, ApJ, 745, 31, doi: 10.1088/0004-637X/745/1/31
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5, doi: 10.1093/pasj/psx080
  * Chornock et al. (2017) Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19, doi: 10.3847/2041-8213/aa905c
  * Connaughton et al. (2016) Connaughton, V., Burns, E., Goldstein, A., et al. 2016, ApJ, 826, L6, doi: 10.3847/2041-8205/826/1/L6
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556, doi: 10.1126/science.aap9811
  * Cowperthwaite et al. (2017) Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJ, 848, L17, doi: 10.3847/2041-8213/aa8fc7
  * Dahlen et al. (2012) Dahlen, T., Strolger, L.-G., Riess, A. G., et al. 2012, ApJ, 757, 70, doi: 10.1088/0004-637X/757/1/70
  * Díaz et al. (2017) Díaz, M. C., Macri, L. M., Garcia Lambas, D., et al. 2017, ApJ, 848, L29, doi: 10.3847/2041-8213/aa9060
  * Drout et al. (2014) Drout, M. R., Chornock, R., Soderberg, A. M., et al. 2014, ApJ, 794, 23, doi: 10.1088/0004-637X/794/1/23
  * Drout et al. (2017) Drout, M. R., Piro, A. L., Shappee, B. J., et al. 2017, Science, 358, 1570, doi: 10.1126/science.aaq0049
  * Einstein (1916) Einstein, A. 1916, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 688
  * Einstein (1918) —. 1918, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), 154
  * Evans et al. (2017) Evans, P. A., Cenko, S. B., Kennea, J. A., et al. 2017, Science, 358, 1565, doi: 10.1126/science.aap9580
  * Fernández & Metzger (2016) Fernández, R., & Metzger, B. D. 2016, Annual Review of Nuclear and Particle Science, 66, 23, doi: 10.1146/annurev-nucl-102115-044819
  * Flewelling et al. (2020) Flewelling, H. A., Magnier, E. A., Chambers, K. C., et al. 2020, ApJS, 251, 7, doi: 10.3847/1538-4365/abb82d
  * Furusawa et al. (2018) Furusawa, H., Koike, M., Takata, T., et al. 2018, PASJ, 70, S3, doi: 10.1093/pasj/psx079
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Graham et al. (2022) Graham, M. J., McKernan, B., Ford, K. E. S., et al. 2022, arXiv e-prints, arXiv:2209.13004.  https://arxiv.org/abs/2209.13004
  * Hsiao et al. (2007) Hsiao, E. Y., Conley, A., Howell, D. A., et al. 2007, ApJ, 663, 1187, doi: 10.1086/518232
  * Ilbert et al. (2005) Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863, doi: 10.1051/0004-6361:20041961
  * Inserra et al. (2021) Inserra, C., Sullivan, M., Angus, C. R., et al. 2021, MNRAS, 504, 2535, doi: 10.1093/mnras/stab978
  * Janiuk et al. (2017) Janiuk, A., Bejger, M., Charzyński, S., & Sukova, P. 2017, New Astronomy, 51, 7, doi: 10.1016/j.newast.2016.08.002
  * Kasen et al. (2013) Kasen, D., Badnell, N. R., & Barnes, J. 2013, ApJ, 774, 25, doi: 10.1088/0004-637X/774/1/25
  * Kasliwal et al. (2017) Kasliwal, M. M., Nakar, E., Singer, L. P., et al. 2017, Science, 358, 1559, doi: 10.1126/science.aap9455
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54, doi: 10.1046/j.1365-8711.2003.06292.x
  * Kawaguchi et al. (2018) Kawaguchi, K., Shibata, M., & Tanaka, M. 2018, ApJ, 865, L21, doi: 10.3847/2041-8213/aade02
  * Kawanomoto et al. (2018) Kawanomoto, S., Uraguchi, F., Komiyama, Y., et al. 2018, PASJ, 70, 66, doi: 10.1093/pasj/psy056
  * Klingler et al. (2021) Klingler, N. J., Lien, A., Oates, S. R., et al. 2021, The Astrophysical Journal, 907, 97, doi: 10.3847/1538-4357/abd2c3
  * Komiyama et al. (2018) Komiyama, Y., Obuchi, Y., Nakaya, H., et al. 2018, PASJ, 70, S2, doi: 10.1093/pasj/psx069
  * LIGO Scientific Collaboration & Virgo Collaboration (2020) LIGO Scientific Collaboration, & Virgo Collaboration. 2020, GRB Coordinates Network, 27184, 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
  * Martin et al. (2018) Martin, R. G., Nixon, C., Xie, F.-G., & King, A. 2018, MNRAS, 480, 4732, doi: 10.1093/mnras/sty2178
  * McKernan et al. (2019) McKernan, B., Ford, K. E. S., Bartos, I., et al. 2019, ApJ, 884, L50, doi: 10.3847/2041-8213/ab4886
  * Metzger (2019) Metzger, B. D. 2019, Living Reviews in Relativity, 23, 1, doi: 10.1007/s41114-019-0024-0
  * Miyazaki et al. (2012) Miyazaki, S., Komiyama, Y., Nakaya, H., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84460Z, doi: 10.1117/12.926844
  * Miyazaki et al. (2018) Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1, doi: 10.1093/pasj/psx063
  * Morgan et al. (2020) Morgan, R., Palmese, A., Garcia, A., et al. 2020, GRB Coordinates Network, 27366, 1
  * Moriya et al. (2019) Moriya, T. J., Tanaka, M., Yasuda, N., et al. 2019, ApJS, 241, 16, doi: 10.3847/1538-4365/ab07c5
  * Morokuma et al. (2016) Morokuma, T., Tanaka, M., Asakura, Y., et al. 2016, PASJ, 68, L9, doi: 10.1093/pasj/psw061
  * Nicholl et al. (2017) Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L18, doi: 10.3847/2041-8213/aa9029
  * Niino et al. (2014) Niino, Y., Totani, T., & Okumura, J. E. 2014, PASJ, 66, L9, doi: 10.1093/pasj/psu115
  * Nugent et al. (2002) Nugent, P., Kim, A., & Perlmutter, S. 2002, PASP, 114, 803, doi: 10.1086/341707
  * Oates et al. (2021) Oates, S. R., Marshall, F. E., Breeveld, A. A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 507, 1296, doi: 10.1093/mnras/stab2189
  * Oguri & Marshall (2010) Oguri, M., & Marshall, P. J. 2010, MNRAS, 405, 2579, doi: 10.1111/j.1365-2966.2010.16639.x
  * Ohgami et al. (2020) Ohgami, T., Tominaga, N., Morokuma, T., et al. 2020, GRB Coordinates Network, 27205, 1
  * Ohgami et al. (2021) Ohgami, T., Tominaga, N., Utsumi, Y., et al. 2021, PASJ, 73, 350, doi: 10.1093/pasj/psab002
  * Okumura et al. (2014) Okumura, J. E., Ihara, Y., Doi, M., et al. 2014, PASJ, 66, 49, doi: 10.1093/pasj/psu024
  * Pei (1992) Pei, Y. C. 1992, ApJ, 395, 130, doi: 10.1086/171637
  * Perego et al. (2017) Perego, A., Radice, D., & Bernuzzi, S. 2017, ApJ, 850, L37, doi: 10.3847/2041-8213/aa9ab9
  * Pian et al. (2017) Pian, E., D’Avanzo, P., Benetti, S., et al. 2017, Nature, 551, 67, doi: 10.1038/nature24298
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591
  * Pratten et al. (2021) Pratten, G., García-Quirós, C., Colleoni, M., et al. 2021, Phys. Rev. D, 103, 104056, doi: 10.1103/PhysRevD.103.104056
  * Quimby et al. (2014) Quimby, R. M., Oguri, M., More, A., et al. 2014, Science, 344, 396, doi: 10.1126/science.1250903
  * Rosswog et al. (2018) Rosswog, S., Sollerman, J., Feindt, U., et al. 2018, A&A, 615, A132, doi: 10.1051/0004-6361/201732117
  * Sasada et al. (2021) Sasada, M., Utsumi, Y., Itoh, R., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A104, doi: 10.1093/ptep/ptab007
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Shibata et al. (2017) Shibata, M., Fujibayashi, S., Hotokezaka, K., et al. 2017, Phys. Rev. D, 96, 123012, doi: 10.1103/PhysRevD.96.123012
  * Singer & Price (2016) Singer, L. P., & Price, L. R. 2016, Phys. Rev. D, 93, 024013, doi: 10.1103/PhysRevD.93.024013
  * 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. (2017) Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, ApJ, 848, L16, doi: 10.3847/2041-8213/aa9059
  * Tagawa et al. (2020a) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020a, ApJ, 899, 26, doi: 10.3847/1538-4357/aba2cc
  * Tagawa et al. (2020b) Tagawa, H., Haiman, Z., & Kocsis, B. 2020b, ApJ, 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tanaka (2016) Tanaka, M. 2016, Advances in Astronomy, 2016, 634197, doi: 10.1155/2016/6341974
  * Tanaka et al. (2017) Tanaka, M., Utsumi, Y., Mazzali, P. A., et al. 2017, PASJ, 69, 102, doi: 10.1093/pasj/psx121
  * 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
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Tominaga et al. (2018a) Tominaga, N., Tanaka, M., Morokuma, T., et al. 2018a, PASJ, 70, 28, doi: 10.1093/pasj/psy007
  * Tominaga et al. (2018b) Tominaga, N., Niino, Y., Totani, T., et al. 2018b, PASJ, 70, 103, doi: 10.1093/pasj/psy101
  * Troja et al. (2017) Troja, E., Piro, L., van Eerten, H., et al. 2017, Nature, 551, 71, doi: 10.1038/nature24290
  * Utsumi et al. (2017) Utsumi, Y., Tominaga, N., Tanaka, M., et al. 2017, PASJ, 70, 1, doi: 10.1093/pasj/psx125
  * Valenti et al. (2017) Valenti, S., Sand, D. J., Yang, S., et al. 2017, ApJ, 848, L24, doi: 10.3847/2041-8213/aa8edf
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692, doi: 10.1086/380563
  * Yamazaki et al. (2016) Yamazaki, R., Asano, K., & Ohira, Y. 2016, Progress of Theoretical and Experimental Physics, 2016, 051E01, doi: 10.1093/ptep/ptw042
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513