Follow-up Survey for the Binary Black Hole Merger GW200224_222234 Using Subaru/HSC and GTC/OSIRIS

We conducted optical imaging observations using the Subaru/HSC on 2020 February 25 (day 1), 28 (day 4), and March 23 (day 28) 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 2020 February 25 at 10:43 UTC, corresponding to 12.3 hr 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 deg² pixel⁻², enabling overlap of the FoVs. Figure 1 shows the survey pointing map. Our survey area covered 56.6 deg², 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.

Zoom In Zoom Out Reset image size

On 2020 February 25 and 28, 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 hr 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 hr because of the overlap between each pointing. On 2020 March 23, 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 machine-readable format).

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σ deviations of fluxes in this aperture, we obtained a map of the 5σ limiting magnitudes for each stack image. Table 2 lists the mode, ¹⁷  maximum, and minimum values of the 5σ limiting magnitudes. The 5σ limiting magnitudes scattered with a range of ∼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 2020 March 23 had a small positional dependence owing to the different exposure time.

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σ limiting magnitudes evaluated from the difference images. The positional dependence of the 5σ limiting magnitudes in the r2 band is shown in Figure 2.

Zoom In Zoom Out Reset image size

After image subtraction, we applied the following criteria to exclude bogus detections (e.g., 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 the PSF flux (S/N)_PSF > 5; (ii) (b/a)/(b/a)_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 < FWHM/(FWHM)_PSF < 1.3; and (iv) PSF-subtracted residual with <3σ 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σ 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).

Zoom In Zoom Out Reset image size

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 (e.g., 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 whether an object is extended. As a result, we found that 2743 and 1883 sources were associated with star-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 (e.g., McKernan et al. 2019; Tagawa et al. 2020a, 2020b). The light curves of AGNs are stochastic (e.g., 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 our three-night observation over one month alone. 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.

Zoom In Zoom Out Reset image size

As the next step, we performed a visual inspection to exclude bogus detections and obtained 223 plausible candidates. Then, we classified these candidates according to whether the PS1 catalog included extended objects within an angular separation, θ_sep, from them. Here, we classified them into two groups: "off-center" candidates (1'' ≤ θ_sep ≤ 15'') and "hostless" candidates (θ_sep > 15''). We conservatively adopted a threshold (θ_sep = 15'') corresponding to a separation of ∼60 kpc at the distance of GW200224_222234 (z ∼ 0.3). This separation is larger than the typical size of galaxies (e.g., Kauffmann et al. 2003; Shen 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_3D for 118 off-center candidates. P_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_3D = 0.5 and obtained 93 likely candidates within the 3σ error of the arrival distance of GW200224_222234 (P_3D ≥ 0.5) and 16 candidates that were likely outside (P_3D < 0.5). The remaining 9 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, 9 candidates without host galaxy information, and 105 hostless candidates.

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 their timescale might be similar to that of 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 (e.g., Quimby et al. 2014; Moriya et al. 2019). Note that the probability that the nature of detected candidates is that of these 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 on the difference fluxes, but also on 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}_{\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 ξ value:

$$\begin{eqnarray}\begin{array}{rcl}\xi & := & \displaystyle \sum _{i}^{{N}_{{\rm{d}}}^{\mathrm{obs}}}\displaystyle \frac{{\left({f}_{{\rm{d}},i}^{\mathrm{obs}}-{f}_{{\rm{d}},i}^{\mathrm{temp}}(p)\right)}^{2}}{{{\sigma }_{{\rm{d}},i}^{\mathrm{obs}}}^{2}}\\ & & +\,\displaystyle \sum _{j}^{{N}_{{\rm{s}}}^{\mathrm{obs}}}\displaystyle \frac{{\left({f}_{{\rm{s}},j}^{\mathrm{obs}}-{f}_{{\rm{s}},j}^{\mathrm{temp}}(p)\right)}^{2}}{{{\sigma }_{{\rm{s}},j}^{\mathrm{obs}}}^{2}},\end{array}\end{eqnarray} \tag{ 1 }$$

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

The 5σ limiting magnitudes in the difference images correspond to 0.37 μJy (day 1−28) and 0.64 μJy (day 4−28) in the r2-band flux, and 1.77 μJy (day 1−28) and 2.36 μ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σ limiting magnitudes, we required the template flux to be fainter than the flux corresponding to the 5σ 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 ξ by probabilistically moving to a location where ξ 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, from Dahlen et al. (2012) for CCSNe, and an assumption of ±1 mag variation for RTs, as shown in Table 4. We derived a Markov chain using 50,000 samples with each redshift from 0.00 to 0.50 every 0.05 using each template model (Type Ia SN, six CCSNe, and nine 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 ξ < 25 from all MCMC samples.

Table 4. Template Parameters Applied in the MCMC Method

Parameter	Range†	SD of Sampling‡
Common
Explosion Time ${t}_{\exp }$ (day)	[−∞ , ∞]	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 MB	[−16.09, −18.05]	0.01
Color Excess EB−V	[0, 1]	0.01
Type IIL bright
Peak Magnitude MB	[−17.92, −19.96]	0.01
Color Excess EB−V	[0, 1]	0.01
Type IIP
Peak Magnitude MB	[−14.43, −18.91]	0.01
Color Excess EB−V	[0, 1]	0.01
Type IIn
Peak Magnitude MB	[−16.98, −20.66]	0.01
Color Excess EB−V	[0, 1]	0.01
Type Ibc normal
Peak Magnitude MB	[−16.09, −18.05]	0.01
Color Excess EB−V	[0, 1]	0.01
Type Ibc bright
Peak Magnitude MB	[−18.46, −20.30]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-10ah)
Peak Magnitude MB	[−16.63, −18.63]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-10bjp)
Peak Magnitude MB	[−17.20, −19.20]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-11bbq)
Peak Magnitude MB	[−18.48, −20.48]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-11qr)
Peak Magnitude MB	[−18.03, −20.03]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-12bb)
Peak Magnitude MB	[−15.29, −17.29]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-12bv)
Peak Magnitude MB	[−18.44, −20.44]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-12brf)
Peak Magnitude MB	[−17.39, −19.39]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-13ess)
Peak Magnitude MB	[−17.49, −19.49]	0.01
Color Excess EB−V	[0, 1]	0.01
Rapid Transient (PS1-13duy)
Peak Magnitude MB	[−17.77, −19.77]	0.01
Color Excess EB−V	[0, 1]	0.01

Note. [†] Range of top-hat distributions for prior distribution. [‡] Standard deviation of the proposed distribution for sampling.

Download table as:  ASCIITypeset images: 1 2

Additionally, for the candidates associated with the PS1 extended object (i.e., 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σ) of the PS1 extended object.
• Criterion 1-2: If the PS1 objects exhibited a large angular separation, θ_sep > 10'', we did not apply criterion 1–1 because the associated PS1 objects were potentially not true host galaxies.
• Criterion 2: When ξ > 25 for θ_sep < 10'' with criterion 1–1, we allowed the templates at a redshift higher than that of the PS1 extended objects because a true faint host galaxy may exist behind the PS1 extended object even if θ_sep < 10''.

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 when available. If not, we estimated a single-band redshift z_single and its standard deviation σ_z using the following formulae:

$$\begin{eqnarray}&&{z}_{\mathrm{single}}:= \displaystyle \frac{{\int }_{0}^{\infty }\phi \,A\,z(D)\,{dD}}{{\int }_{0}^{\infty }\phi \,A\,{dD}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\sigma }_{z}^{2}:= \displaystyle \frac{{\int }_{0}^{\infty }\phi \,A{(z(D)-{z}_{\mathrm{single}})}^{2}\,{dD}}{{\int }_{0}^{\infty }\phi \,A\,{dD}},\end{eqnarray} \tag{ 3 }$$

where ϕ = ϕ(λ, M) is the luminosity function of galaxies at a rest wavelength λ derived from the luminosity functions in the U, B, V, R, and I bands provided in Ilbert et al. (2005) and Planck cosmology (Planck Collaboration et al. 2014), A = A(D) is the surface area observed at a distance D, M = M(D; m_j ) is the absolute magnitude derived with the jth band (r or i) apparent magnitude m_j at D in the observer frame, and λ = λ(D; λ_j ) is the rest wavelength redshifted from an observed wavelength λ_j at D. We also estimated z_single and σ_z for candidates with an SDSS photometric redshift (Table 5). Their comparison illustrates their approximate consistency (Figure 5). Here, we corrected for the Galactic extinction (Schlafly & Finkbeiner 2011) ¹⁸  when we derived absolute magnitudes from apparent magnitudes.

Zoom In Zoom Out Reset image size

Table 5. Information of the Candidates that Are Consistent with the Transient Templates or Examples

Name	Coordinate (J2000.0)		${\theta }_{\mathrm{sep}}^{\dagger }$	zsingle	σz	zSDSS	σSDSS	${P}_{3{\rm{D}}}^{\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
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮

Note. [†] We indicate with a ellipsis when a candidate is not associated with PS1 extended objects.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Of 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, θ_sep, P_3D, and probable transient templates of the 200 candidates consistent with the templates. A fitting result of JGEM20ewa with ξ = 0.6 is shown as an example of the candidates that are 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.

Zoom In Zoom Out Reset image size

Of the 200 candidates consistent with the templates of SNe, 105 candidates were associated with the PS1 extended objects (i.e., off-center candidates), and the remaining 95 candidates had no PS1 extended objects within <15''.

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 ± 0.05 and 0.34 ± 0.04, respectively. The positions of the candidates and PS1 objects are shown in Figure 7 in the left (JGEM20fvn) and middle panels (JGEM20gqu). Because JGEM20fvn and JGEM20gqu exhibited θ_sep = 13''8 and 11''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 one remaining 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 higher than the redshift of the PS1 object associated with JGEM20hgo (z = 0.03 ± 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 θ_sep and P_3D of candidates consistent with the transient templates are shown in Table 5.

Zoom In Zoom Out Reset image size

Of the 20 candidates that were inconsistent with the templates and examples, one candidate (JGEM20cvb) exhibited a high ξ value of 2175. However, it apparently matched the template light curve, as shown in Figure 8. This high ξ could be attributed to the underestimation of photometric errors or the uncertainties on the template fluxes. For instance, ξ was improved to ∼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.

Zoom In Zoom Out Reset image size

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.

As a result, 22 objects, that is, 3 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, θ_sep, P_3D, and fitting results are shown in Table 6 and Figure 9.

Zoom In Zoom Out Reset image size

Table 6. Information on Candidates that Are Inconsistent with the Templates of SNe and Their Associated PS1 Extended Objects

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

Note. [†] We indicate with an ellipsis when a candidate is not associated with PS1 extended objects.

Download table as:  ASCIITypeset image

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). ¹⁹  We refer 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 r2 band, most of the candidates were detected only in the r2 band. Thus, we sampled mock SNe whose brightness was decaying in the r band, assuming the reference images were taken 27 days after detection because the variation in the SN magnitude during 3 days between days 1 and 4 (<0.1 mag) was negligible compared to the standard deviation of the 5σ 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 ±50% error from the 1σ 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}_{\mathrm{SN}}=201$, and our observed area of S = 56.6 deg². The vertical error bar was defined as $\sqrt{{N}_{\mathrm{SN}}}/S$ by assuming that the number followed a Poisson distribution with the expected value of ${N}_{\mathrm{SN}}$. The horizontal error bar is the 1σ error of the 5σ limiting magnitude in the difference images on day 4 in the r2 band.

Zoom In Zoom Out Reset image size

Our result is consistent with the expected number density of the SN detection within a 1σ 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σ lower end of the expected value. This might be because we ignored nuclear transients in this search.

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 γ-ray emissions caused by a relativistic outflow. They showed that optical afterglow peaks appear ∼4.5 × 10⁴ s after GW detection (almost the same as the start of our first epoch of observations) if the ambient matter density is 0.01 cm⁻³; but the peak flux is expected to be lower than our observation limits (0.37 μJy in r band) at a distance of 1710 Mpc. Therefore, an afterglow like this could not be detected in our search.

Of the 19 candidates, 16 candidates inconsistent with all templates and examples exhibited an increasing difference flux in one band and a decreasing difference flux in the other band (e.g., JGEM20ejv), or a more rapid decline in the z band than in the r2 band (e.g., 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 }\lt {5.2}_{-1.9}^{+2.4}\times {10}^{41}$ erg s⁻¹ (${9.1}_{-3.3}^{+4.1}\times {10}^{41}$ erg s⁻¹) and $\nu {L}_{\nu }\lt {1.8}_{-0.6}^{+0.8}\times {10}^{42}$ erg s⁻¹ (${2.4}_{-0.9}^{+1.1}\times {10}^{42}$ erg s⁻¹) 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⁴² erg s⁻¹, McKernan et al. 2019).

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 γ-ray (Abdalla et al. 2021). We searched our sources to find out whether these sources were detected. The source JGEM20aoz was located at R. A. = 11^h41^m4181, $\mathrm{decl}.\,=\,-14^\circ 07^{\prime} 54^{\prime\prime} 78$ (J2000.0), 48 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 r2 and z bands, respectively, in the first-epoch (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 likely is 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 2020 March 5, with 10σ 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''. 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)_PSF = 0.6 and 0.5, respectively. Although no source was found at the position of AT2020ehy in either the stacked or the 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)_PSF > 5.

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_ej = 0.05 M_⊙ and an electron fraction Y_e = 0.30–0.40 (Banerjee et al. 2020). This kilonova model can explain the observed early multicolor light curve of GW170817/AT2017gfo. The peak of the bolometric luminosity approximately scales with M_ej to the power of 0.35 (see Tanaka 2016; Fernández & Metzger 2016; Metzger 2019). Y_e influences the light curves through the opacity. In the low-Y_e case (Y_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_e = 0.30–0.40.

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

Zoom In Zoom Out Reset image size

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 (∼1710 Mpc), can we detect it?

Figure 15 shows the r-, i-, and z-band light curves under the assumption that the kilonova is located at the estimated distance of GW200224_222234. The colored area is the uncertainty caused by the 1σ credible range of the distance. The horizontal lines with the upward arrows are the 5σ limiting magnitudes measured in the difference images in each epoch in each band (r2 and z). In the r2 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.

Zoom In Zoom Out Reset image size

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. ²⁰  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 ∼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 the Rubin/Legacy Survey of Space and Time (LSST).