1]Faculty of Mathematics and Physics, Institute of Science and Engineering,
Kanazawa University, Kakuma, Kanazawa, Ishikawa 920-1192, Japan

2]Graduate School of Science and Technology, Hirosaki University, Aomori
036-8561, Japan

3]Department of Physics and Astronomy, Louisiana State University, 202
Nicholson Hall, Baton Rouge, Louisiana 70803, USA

4]Department of Astronomy, Graduate School of Science, The University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

5]Department of Physics, Hiroshima University, 1-3-1 Kagamiyama, Higashi-
Hiroshima, Hiroshima 739-8526, Japan

6]Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5
Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

7]Department of Applied Physics, Research Institute of Explosive Stellar
Phenomena, Fukuoka University, 8-19-1, Jonan, Nanakuma, Fukuoka 814-0180,
Japan

8]Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan

9]Center for Gravitational Physics, Yukawa Institute for Theoretical Physics,
Kyoto University, 606-8502, Kyoto, Japan

10]Nishi-harima Astronomical Observatory, Center for Astronomy, University of
Hyogo, 407-2 Nishi-gaichi, Sayo, Hyogo 679-5313, Japan

11]Department of Physics, Osaka City University, Osaka, Osaka, 558-8585,
Japan, Nambu Yoichiro Institute of Theoretical and Experimental Physics, Osaka
City University, Osaka, Osaka, 558-8585, Japan

12]Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1
Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

13]Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama,
Meguro-ku, Tokyo 152-8551, Japan

14]Department of Physics, Rikkyo University, Toshima, Tokyo 171-8501, Japan

15]KEK and Sokendai, 1-1 Oho, Tsukuba 305-0801, Japan

16]Department of Physics, Okayama University, Okayama, Okayama 700-8530, Japan

17]Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, University of Tokyo,
Kashiwa, Chiba 277-8583, Japan

18]Department of Earth and Planetary Sciences, Faculty of Sciences, Kyushu
University, Nishi, Fukuoka 819-0395, Japan

19]High Energy Accelerator Research Organization (KEK), Oho 1-1, Tsukuba
305-0801, Japan

20]MAXI team, Institute of Physical and Chemical Research (RIKEN), 2-1,
Hirosawa, Wako, Saitama 351-0198, Japan

21]Institute of Astronomy, Graduate School of Science, The University of
Tokyo, 2-21-1, Osawa, Mitaka, Tokyo 181-0015, Japan

22]Faculty of Law, Ryukoku University, 67 Fukakusa Tsukamoto-cho, Fushimi-ku,
Kyoto 612-8577, Japan

23]Department of Physics, Nihon University, 1-8 Kanda-Surugadai, Chiyoda-ku,
Tokyo, 101-8308, Japan

24]Research Center for the Early Universe (RESCEU), School of Science, The
University of Tokyo, Tokyo 113-0033, Japan

25]Department of Physics, Faculty of Science and Engineering, Konan
University, 8-9-1 Okamoto, Kobe, Hyogo 658-8501, Japan

26]Astronomical Institute, Graduate School of Science, Tohoku University,
Aoba, Sendai 980-8578, Japan

27]College of Science and Engineering, Department of Physics and Mathematics,
Aoyama Gakuin University, 5-10-1 Fuchinobe, Chuo-ku, Sagamihara, Kanagawa
252-5258, Japan

28]Department of Physics, Faculty of Science, Toho University, 2-2-1 Miyama,
Funabashi, Chiba 274-8510, Japan

29]Department of Physics, Kobe University, Kobe 657-8501, Japan

30]National Institute of Technology, Numazu College, Ooka 3600, Numazu,
Shizuoka 410-8501, Japan

31]Research Center for Space Science, Advanced Research Laboratories, Tokyo
City University, 8-15-1, Todoroki, Setagaya, Tokyo 158-0082, Japan, and
Department of Information and Management Systems Engineering, Nagaoka
University of Technology, Niigata 940-2188, Japan

32]Division of Science, National Astronomical Observatory of Japan, 2-21-1,
Osawa, Mitaka, Tokyo 181-8588, Japan

33]Department of Earth Science and Astronomy, College of Arts and Sciences,
The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan

34]Kavli Institute for Particle Astrophysics and Cosmology, SLAC National
Accelerator Laboratory, Stanford University, 2575 Sand Hill Road, Menlo Park,
CA 94025, USA

35]Subaru Telescope, National Astronomical Observatory of Japan, 650 North
A’ohoku Place, Hilo, Hawaii 96720, USA

# Gravitational Wave Physics and Astronomy in the nascent era

Makoto Arimoto  Hideki Asada  Michael L. Cherry  Michiko S. Fujii  Yasushi
Fukazawa  Akira Harada  Kazuhiro Hayama  Takashi Hosokawa  Kunihito Ioka
Yoichi Itoh  Nobuyuki Kanda  Koji S. Kawabata  Kyohei Kawaguchi  Nobuyuki
Kawai  Tsutomu Kobayashi  Kazunori Kohri  Yusuke Koshio  Kei Kotake  Jun
Kumamoto  Masahiro N. Machida  Hideo Matsufuru  Tatehiro Mihara  Masaki Mori
Tomoki Morokuma  Shinji Mukohyama  Hiroyuki Nakano  Tatsuya Narikawa  Hitoshi
Negoro  Atsushi Nishizawa  Takayuki Ohgami  Kazuyuki Omukai  Takanori Sakamoto
Shigeyuki Sako  Mahito Sasada  Yuichiro Sekiguchi  Motoko Serino  Jiro Soda
Satoshi Sugita  Kohsuke Sumiyoshi  Hajime Susa  Teruaki Suyama  Hirotaka
Takahashi  Kazuya Takahashi  Tomoya Takiwaki  Takahiro Tanaka  Masaomi Tanaka
Ataru Tanikawa  Nozomu Tominaga  Nami Uchikata  Yousuke Utsumi  Mark R. Vagins
Kei Yamada  Michitoshi Yoshida  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [
[  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [

###### Abstract

The detections of gravitational waves (GW) by LIGO/Virgo collaborations
provide various possibilities to physics and astronomy. We are quite sure that
GW observations will develop a lot both in precision and in number owing to
the continuous works for the improvement of detectors, including the
expectation to the newly joined detector, KAGRA, and the planned detector,
LIGO-India. In this occasion, we review the fundamental outcomes and prospects
of gravitational wave physics and astronomy. We survey the development
focusing on representative sources of gravitational waves: binary black holes,
binary neutron stars, and supernovae. We also summarize the role of
gravitational wave observations as a probe of new physics.

xxxx, xxx

###### Contents

  1. 1 Introduction
    1. 1.1 Gravitational wave from binary black holes
    2. 1.2 Gravitational wave from binary neutron stars
    3. 1.3 Gravitational wave from Supernovae
    4. 1.4 Test of general relativity
    5. 1.5 Organization of this article
  2. 2 Binary Black Holes
    1. 2.1 Gravitational waves from binary black holes
    2. 2.2 Theoretical study on binary black hole formation
      1. 2.2.1 Formation of isolated stellar binaries
      2. 2.2.2 Binary black hole formation by dynamical interactions in dense star clusters
      3. 2.2.3 Outlook
  3. 3 Binary Neutron Stars
    1. 3.1 Exploring NS physics
      1. 3.1.1 BNS waveform modeling based on numerical-relativity simulations
      2. 3.1.2 Measuring tidal deformability from GW
    2. 3.2 Association of a short-duration gamma-ray burst
      1. 3.2.1 High-energy Observations
      2. 3.2.2 Theoretical Interpretation
    3. 3.3 Optical/Infrared Counterparts and Heavy Element Nucleosynthesis
      1. 3.3.1 J-GEM
      2. 3.3.2 GW170817
      3. 3.3.3 Advances of Theoretical Models
      4. 3.3.4 Future Prospects
  4. 4 Supernovae
    1. 4.1 GW signatures from core-collapse supernovae
      1. 4.1.1 The CCSN GW signatures
      2. 4.1.2 Circular Polarization of CCSN GWs
      3. 4.1.3 Rapid Rotation: GW and Neutrino signals from an exploding $27M_{\odot}$ star
      4. 4.1.4 Hydrodynamics and GW signals from a BH-forming massive star
      5. 4.1.5 Short Summary and Future Prospects
    2. 4.2 Understanding Supernovae via Neutrino Emissions
      1. 4.2.1 SK-Gd
      2. 4.2.2 Theoretical studies of supernova explosion and neutrino emission
      3. 4.2.3 Boltzmann-radiation-hydrodynamics simulation of the core-collapse supernovae
      4. 4.2.4 Light curve of supernova neutrinos
      5. 4.2.5 Information of supernova mechanism from supernova neutrino detectors
      6. 4.2.6 Neutrino bursts from black hole formation
      7. 4.2.7 Diffuse supernova neutrino background
      8. 4.2.8 Equation of state for supernovae and neutron star mergers
      9. 4.2.9 Supernovae simulations on GPUs
  5. 5 Probes of new physics
    1. 5.1 Test of gravity theories using gravitational wave data
      1. 5.1.1 Varieties of gravitational waveform within general relativity
      2. 5.1.2 Model independent test of waveform consistency with general relativity
      3. 5.1.3 Model independent test of extra polarizations
      4. 5.1.4 Gravitational waveform corresponding to various modifications of gravity
      5. 5.1.5 Dedicated test for particular models
    2. 5.2 Dark matter
      1. 5.2.1 Axion
      2. 5.2.2 Primordial black holes
  6. 6 Conclusion

##  1 Introduction

Table 1: Observing terms of LIGO/Virgo detectors. Detectors column: H1 at Hanford, Washington, and L1 at Livingston, Louisiana. V1 at Cascina, Italy.  | Observation Period | detectors | events reported
---|---|---|---
O1 | Sep 2015 to Jan 2016 | H1, L1 | 3 BBHs
O2 | Nov 2016 to Sep 2017 | H1, L1 | 3 BBHs
O2 | Aug 2017 | H1, L1, V1 | 4 BBHs, 1 BNS
O3a | Apr 2019 to Sep 2019 | H1, L1, V1 | many BBHs, 1 BNS
O3b | Nov 2019 to Apr 2020 | H1, L1, V1 | N/A

###  1.1 Gravitational wave from binary black holes

The first gravitational wave (GW) event from the binary black holes (BBHs),
GW150914, which was discovered by LIGO/Virgo collaboration brought us many
surprises GW150914 . Regarding binary neutron star (BNS) coalescences, it was
believed that there is not much uncertainty in estimating their event rate
compared with the case of BBHs. Even for BNS coalescences, the uncertainty in
the estimated event rate was expected to be about one order of magnitude. On
the other hand, although it is conceivable that BBHs should exist in our
universe, the estimation of the GW event rate had no observational support,
differently from the case of BNSs. Hence, BBHs were not considered to be a
promised GW source. It was somewhat surprising that such a system was the
first directly detected GW source. However, it was even more surprising that
those black holes (BHs) were heavier than many people expected, i.e., each
component BH has a mass about 30$M_{\odot}$. The mass of the BH candidates
found so far by electromagnetic observations is estimated to be at most about
20$M_{\odot}$ Farr:2010tu . Although there were some candidate objects that
did not contradict with even heavier mass within the observational error,
there is no object that can be said to have about 30$M_{\odot}$. The
observation of GWs has changed our understanding about the mass distribution
of BHs in the universe rosswog17 ; hotokezaka18 . Since then, the number of
BBH observations officially announced by the LIGO/Virgo collaborations has
been increasing LIGOScientific:2018mvr . Not all of them are massive BHs, but
the succeeding observations indicate that the first GW event, GW150914, was
not so special.

As a result of this discovery, how those BBHs detected by GWs were formed in
the history of the universe has become one of the great mysteries of
astrophysics. A variety of BBH formation scenarios have been proposed abbott16
. The most promising candidate for the formation channel would be the one
through the binary interaction between heavy star binaries formed from gases
with a low metallicity. Even for this leading scenario, there are various
uncertainties: the binary formation rate from the gas cloud, stellar evolution
including mass loss processes, binary interaction to form close binaries,
formation process of BHs, and so on. The binary evolution can be so complex,
changing its orbital radii and eccentricity. Therefore, it is a theoretical
challenge to predict the distributions of BH masses, orbital separation and
eccentricity at the formation of a BH binary, to make comparison with
observations. As alternative scenarios, we can think of the formation of BH
binaries in star clusters Rodriguez et al. (2015) and their formation in the
primordial universe Sasaki:2018dmp . Many variants of these scenarios can be
thought of. In the future, as the details of the theoretical model become
better understood, and as GW observations accumulate the information about the
event rate, the mass and distance distributions, the distribution of BH spins,
and so on, we would be able to identify the contributions from respective
formation channels of BBHs.

###  1.2 Gravitational wave from binary neutron stars

The first GW event from a BNS merger GW170817 observed on August 17, 2017,
marked a splendid opening of GW astronomy GW170817 . BNS coalescence was
counted as a guaranteed GW source. The discovery was followed by simultaneous
observation of gamma-ray bursts, identification of the electromagnetic
counterpart by optical and infrared telescopes abbott17MMA , long-term follow-
up observations by radio waves and X-rays, detection of apparent superluminal
motion Mooley:2018qfh ; Ghirlanda:2018uyx . How BNS coalescence proceed was
revealed to be something close to what was theoretically predicted. In the
study of BNS coalescence, weak gravity approximation is not appropriate any
more. Hence, numerical relativity is an indispensable research tool with the
knowledge of equation of state for the high density nuclear matter. At the
same time various radiation fields, photons and neutrinos, originating from
high density matter must be solved and the magnetic fields also can play an
important role in the simulations. Therefore combining various techniques for
numerical simulation is required. In the case of GW170817, from its optical
and infrared light curves, it is thought that the ejecta mass was about
$0.05M_{\odot}$ and a part of them possessed high neutron excess, resulting in
production of significant amount of $r$-process elements. The amount of the
produced $r$-process elements may be sufficient to explain almost all of them
in the universe. In addition, numerical simulations show that the merged
binary becomes a massive neutron star (NS) supported by rotation, and the
material stripped around it forms a disk. The simultaneous (1.7 second after
the coalescence time) observation of a short gamma-ray burst (SGRB) indicates
the connection between SGRB and BNS coalescence. Although the observed SGRB
was extraordinary weak, this can be understood because the SGRB jet was not
pointing at the line of sight toward us Ioka:2017nzl ; Ioka:2019jljb . This
interpretation gained a strong support because apparent superluminal motion
was observed by radio waves Mooley:2018qfh ; Ghirlanda:2018uyx . The detection
of apparent superluminal motion is also a strong evidence for that a
relativistic jet has been launched. In addition, the GW observations at this
event also provided restrictions on the equation of state of dense nuclear
matter, mostly from the constraint on the tidal deformability of NSs just
before the coalescence GW170817 . Here again, numerical relativity technique
played an important role. It is also an important theoretical challenge to
derive the equation of state that can be compared with phenomena from the
first principle of QCD. High-energy gamma-rays, neutrinos, and cosmic rays
were not detected in GW170817 Abdalla et al. (2017); Albert et al. (2017);
Ajello et al. (2018); Abe et al. (2018), although BNS mergers and SGRBs are
considered to be sources of these particles Murase et al. (2018); Kimura et
al. (2018); Biehl et al. (2018); Kyutoku & Ioka (2016). We do not further
discuss these multi-messenger signals in this review.

###  1.3 Gravitational wave from Supernovae

Although GWs from supernovae have not been detected yet, they are important
targets of gravitational wave observations. After more than a half-century of
continuing efforts ever since the seminal work by Colgate and White (1966)
colgate , theory and neutrino radiation-hydrodynamics simulations of core-
collapse supernovae (CCSNe) are now converging to a point that multi-
dimensional (multi-D) hydrodynamics instabilities play a key role in the
neutrino mechanism bethe , the most favoured scenario to trigger explosions.
In fact, self-consistent simulations in three spatial dimensions (3D) now
report shock revival of the stalled bounce shock into explosion by the multi-D
neutrino mechanism (see bernhard16 ; janka16 ; burrows13 ; Kotake12  for
reviews).

From observational point of view, mounting evidence from multi-wavelength
electromagnetic-wave (EM) observations (such as blast/ejecta morphologies,
line profiles, polarizations) are pointing toward CCSNe being generally
aspherical (e.g., casA ; tanaka17  and references therein). Nevertheless,
these EM signals are secondary as a probe into the supernova inner-workings,
providing information only after a shock break-out from the massive stars. On
the other hand, the GWs from CCSNe are direct observables, which imprint a
live information of the central engine. Furthermore coincident detection of
CCSN neutrinos with GWs is important not only for significantly enhancing the
detectability of the GWs (e.g., naka016 ), but also for unraveling the
hydrodynamics evolution from the onset of core-collapse toward the forming
compact objects (neutron stars and black holes, see mirizzi  for a review).
Current estimates of CCSN rates in our Milky way predict one event every $\sim
40\pm 10$ year horiuchi_kneller . While rare, they provide us unique
opportunity to study CCSNe using a full set of multi-messenger observables
including GWs, neutrinos, and nuclear gamma-rays gehrels .

###  1.4 Test of general relativity

The direct detection of GWs has opened a new window for tests of general
relativity (GR). In order to detect GW signals, we need theoretical prediction
for the GW waveform in advance. In more complex theory of gravitation that
extends GR, theoretically predicting the waveform is more difficult. Our
current knowledge about the waveform in modified gravity is limited, but some
possible deviations in the waveform from GR have been predicted. GW
observations so far have not suggested any deviation from GR
LIGOScientific:2019fpa . However, the lack of deviation itself can have
important physical meaning. Besides constraints on the deviation in the
waveform, an important constraint was also derived from the speed of
propagation of GWs, i.e., from the detection of short gamma-ray burst
associated with GW170817 event at a distance of about 40 Mpc. The fractional
deviation of the propagation speed of GWs from that of light is less than
$\approx 10^{-15}$ Monitor:2017mdv , which excludes various possible
extensions of GR Creminelli:2017sry ; Sakstein:2017xjx ; Baker:2017hug ;
Langlois:2017dyl ; Crisostomi:2017lbg ; Dima:2017pwp .

GWs from binaries can be used as the standard siren Schutz:1986gp , providing
an alternative method to measure the cosmic expansion history. GW170817
already gives us an estimate of the Hubble constant Abbott:2017xzu ;
Hotokezaka:2018dfi , but we do not further discuss the aspects of testing
cosmological models using GWs in this review.

###  1.5 Organization of this article

In this review paper we overview the vast area of newly developing
gravitational wave physics and astronomy mentioned above. This paper is
organized as follows in Sec. 2, 3 and 4, we focus on three major sources of
gravitational waves, i.e., BBHs, BNSs, and SNe. We discuss various topics
emanating from the detection of gravitational waves, including the future
prospects. In Sec. 5 we focus more on the new physics examined by using GWs as
a probe. We discuss the test of gravity and implication to the dark matter
physics. Sec. 6 is devoted to a short conclusion.

##  2 Binary Black Holes

###  2.1 Gravitational waves from binary black holes

LIGO/Virgo has reported GWs from BBH mergers as a result of O1 and O2
observing run LIGOScientific:2018mvr  that is called “Gravitational-Wave
Transient Catalog” (GWTC)-1. The best fit parameters for these binaries are
listed in Table 2. We have large errors in the parameter estimation,
especially for the individual masses and the effective spin parameter, because
the signal-to-noise ratios (SNRs) are not so high. There are results reported
by alternative searches. One of them reports other candidates of GW events
Venumadhav:2019lyq .

Table 2: BBH merger events observed during LIGO/Virgo O1 and O2. We show only the median values for the individual masses $m_{1}$ and $m_{2}$, chirp mass ${\cal M}$, effective spin $\chi_{\rm eff}$, final BH mass $M_{{\rm fi}}$, and its angular momentum parameter $a_{{\rm fi}}$. These numbers are of source frame. The estimated distance in Mpc and red-shift $z$ are also shown. The 90% credible intervals have been reported in Table III of Ref. LIGOScientific:2018mvr . Event | $m_{1}/M_{\odot}$ | $m_{2}/M_{\odot}$ | ${\cal M}/M_{\odot}$ | $\chi_{\rm eff}$ | $M_{{\rm fi}}/M_{\odot}$ | $a_{{\rm fi}}$ | distance (Mpc,$z$)
---|---|---|---|---|---|---|---
GW150914 | 35.6 | 30.6 | 28.6 | -0.01 | 63.1 | 0.69 | 430 Mpc | 0.09
GW151012 | 23.2 | 13.6 | 15.2 | 0.05 | 35.7 | 0.67 | 1060 Mpc | 0.20
GW151226 | 13.7 | 7.7 | 8.9 | 0.18 | 20.5 | 0.74 | 440 Mpc | 0.09
GW170104 | 30.8 | 20.0 | 21.4 | -0.04 | 49.1 | 0.66 | 960 Mpc | 0.19
GW170608 | 11.0 | 7.6 | 7.9 | 0.03 | 17.8 | 0.69 | 320 Mpc | 0.07
GW170729 | 50.2 | 34.0 | 35.4 | 0.37 | 80.3 | 0.81 | 2750 Mpc | 0.48
GW170809 | 35.0 | 23.8 | 24.9 | 0.08 | 56.4 | 0.70 | 990 Mpc | 0.20
GW170814 | 30.6 | 25.2 | 24.1 | 0.07 | 53.4 | 0.72 | 580 Mpc | 0.12
GW170818 | 35.4 | 26.7 | 26.5 | -0.09 | 59.8 | 0.67 | 1020 Mpc | 0.20
GW170823 | 39.5 | 29.0 | 29.2 | 0.09 | 65.6 | 0.71 | 1850 Mpc | 0.34

The obtained mass and spin distributions are summarized in Fig. 1. The chirp
mass is defined by

${\cal M}=\mu^{3/5}M^{2/5}\,,$ (1)

with $M=m_{1}+m_{2}$ and $\mu=m_{1}m_{2}/M$, where $m_{1}$ and $m_{2}$ are the
component masses. The effective spin parameter is defined by

$\chi_{\rm
eff}=\frac{1}{M}\left(\frac{\mbox{\boldmath{$S$}}_{1}}{m_{1}}+\frac{\mbox{\boldmath{$S$}}_{2}}{m_{2}}\right)\cdot\hat{\mbox{\boldmath$L$}}\,,$
(2)

where $\mbox{\boldmath{$S$}}_{1}$ and $\mbox{\boldmath{$S$}}_{2}$ are spin
angular momentum vectors of the respective bodies and
$\hat{\mbox{\boldmath$L$}}$ is the unit vector pointing at the direction of
the orbital angular momentum of the binary.

Figure 1: Distribution of the median values of the chirp mass ${\cal M}$, mass
ratio $q=m_{2}/m_{1}$, effective spin $\chi_{\rm eff}$, and individual masses
(which are redundant though) of BBHs observed during LIGO/Virgo O1 and O2. It
should be noted that there are large errors in the current observations. Even
for the chirp mass which is the most sensitive parameter in the GW waveform of
BBH mergers, the 90% credible interval gives $\lesssim$ 20% difference in the
estimation.

The reported event rate by the LIGO/Virgo Collaborations is ranging from $9.7$
to $101$ $\mbox{Gpc}^{-3}\mbox{yr}^{-1}$ for BBH mergers. The masses of the
observed BHs are relatively large compared with the masses of stellar mass BH
candidates found by X-ray observations (see, e.g., Ref. Farr:2010tu  and
references therein). The origin of BBHs is an important new issue of debate
raised by the GW observations. A promising formation path is the one through
the stellar evolution of low metallicity stars, which is the main topic of
Section 2.2. There are various formation channels and the main focus is
whether or not the stellar dynamics in the star cluster plays a crucial role.
An alternative is the formation of BH binaries in the early universe, which is
the so-called primordial BH scenario. The formation process of primordial BH
binaries is discussed in Section 5.2.

###  2.2 Theoretical study on binary black hole formation

Mergers of BBHs are the dominant population of gravitational wave sources.
Currently, two evolutionary pathways are envisaged as origins of BBHs. First
is the so-called in-situ BBH formation scenario, where massive star binaries
are formed and each member star becomes a black hole (BH) after stellar
evolution. Second is the so-called ex-situ scenario, where single stars formed
in a dense star cluster become BHs after stellar evolution and those BHs
become binaries as a result of dynamical interaction. Here, we describe our
current understanding of those scenarios in this order.

####  2.2.1 Formation of isolated stellar binaries

##### A

Binary star formation at low metallicities

Figure 2:  An example of the latest 3D radiation-hydrodynamics simulations of
the primordial star formation, in which the metallicity is exactly zero (from
Sugimura et al. submitted (Sugimura et al., 2020)), showing a snapshot at the
epoch when massive protostars in a binary system accrete the gas under the
radiative feedback. Presented are the volume-rendered 3D distribution of the
mass density, together with ionization fronts (yellow surfaces) and protostars
(black dots).

Whereas massive stellar binaries are possible astrophysical progenitors of BH-
BH binaries that finally merge, their formation process is still uncertain. In
particular, major focus is put on formation of such systems at low
metallicities $Z\lesssim 0.1~{}Z_{\odot}$, which are required for yielding
massive ($\simeq$ a few $\times 10~{}M_{\odot}$) BHs avoiding significant mass
loss during the stellar evolution. Observations suggest that a large fraction
of massive stars are born in binaries in the solar neighborhood. Some of them
are in very tight systems with separations of only $\lesssim 0.1$ AU, and are
possible candidates that finally evolve into BH binaries that merge within the
Hubble time. On the other hand, since the delay time until the BBH merger is
potentially as long as the Hubble timescale, depending on the initial
separation, a part of the progenitor stars may have been born in the early
universe. Therefore, our goal is to elucidate the formation process of massive
and tight binary stars in a variety of low-metallicity environments, including
the local and early universe.

A straightforward method which ultimately leads us to the goal is directly
following the evolution with numerical simulations. One of the key processes
is the gas gravitational fragmentation that potentially leads to bearing
multiple stars, some of which are to be in binaries. Magnetic fields, if any,
should be important because it transfers the angular momentum of the gas,
affecting the orbital evolution of newly forming binaries. The stellar
radiative feedback is also critical to determine how massive binaries finally
appear, halting the mass accretion onto the stars (McKee & Tan, 2008; Hosokawa
et al., 2011, 2016). We consequently need to solve the interplay between gas
dynamics, magnetic fields, and radiative fields in full 3D spatially resolving
a very wide dynamic range. Ultimately integrating the evolution over $\sim$
Myr is required to cover the whole evolution of the star formation process.
Although this is still challenging, recent studies are advancing to
incorporate the above processes as much as possible with limited computational
resources. For instance, Fig. 2 shows a snapshot in our recent 3D radiation
hydrodynamic simulations using the adaptive-mesh-refinement code SFUMATO
(Sugimura et al., 2020). In this particular case, we follow the long-term
($\sim 0.1$ Myr) evolution in the protostellar accretion stage of the
primordial star formation. We see a multiple stellar system including a
massive binary system that has appeared after the disk fragmentation in the
earlier stage. The double bipolar H II regions are growing from the two
massive stars that are still accreting the surrounding gas under the radiative
feedback. Although this provides a glimpse of the binary formation in the
early universe, the newborn binary has a large separation: the two massive
stars exciting the bipolar H II regions are separated by more than $10^{3}$
AU. The formation process of the massive and tight binaries, some of which
should evolve to the gravitational wave sources, is yet to be revealed.
Obviously further studies and simulations are necessary for explaining how
such systems are formed. We describe our recent studies and future prospects
in what follows.

##### B

Role of magnetic fields in close binary formation

Observations of nearby star-forming regions indicate that a large fraction of
stars are members of binary systems (Duchêne & Kraus, 2013) and the binary
frequency is even higher for high-mass stars than for low-mass ones (Sana et
al., 2012). Theoretically, numerical simulations have shown that a binary
stellar system forms as a result of fragmentation of a molecular cloud core
during the gravitational collapse (Matsumoto & Hanawa, 2003). In rapidly
rotating and gravitationally unstable cores, fragmentation occurs more easily
(Miyama et al., 1984; Tsuribe & Inutsuka, 1999) and the accretion rate is
higher after the protostars are formed. The rotation velocity increases and
the cloud becomes disk-like during the collapse owing to the angular momentum
conservation. This promotes fragmentation and subsequent binary formation.

According to numerical simulations, the binary separation continues to
increase with time as a parcel of gas with larger angular momentum falls onto
the central region or proto-binary system at later times (Matsumoto & Hanawa,
2003; Machida et al., 2008) and the outcome is a wide binary system with a
separation typically greater than $\gg 100$ AU. In contrast, observationally
high-mass stars tend to be in close-binary systems with the separation 1-10 AU
(Duchêne & Kraus, 2013). This discrepancy between theory and observations can
be circumvented if we take non-ideal magnetohydrodynamics (MHD) effect into
account.

Recent non-ideal MHD simulations for present-day star formation showed that
the excess angular momentum is transported from the proto-binary system by
magnetic braking and by outflow launching and, as a consequence, the close
binary system can be formed in a highly unstable magnetized cloud. Figures 3
is a snapshot of a proto-binary system in an example of such simulations
(Saiki & Machida, 2020). In the figure, we can see that powerful jet is driven
by each of the protostars. The separation of the binary system is about 10 AU
at this epoch. On the equatorial plane, both the circumstellar and
circumbinary disks can be confirmed (the left panel of Fig. 3). Two cavity-
like structures created by the jets are seen in the middle panel of Fig. 3.
The velocity of the jets driven by the protostars exceeds $>100$ km s-1. Those
high-velocity components are enclosed by the low-velocity outflow of $\sim 10$
km s-1 driven by the circumbinary disk. The angular momentum of the binary
system is extracted from the system both by the large-scale wide-angle
outflows and small-scale collimated jets. The binary jet/outflow structure is
also reported from recent ALMA observations (Tobin et al., 2019).

Figure 3:  Proto-binary formation in the MHD simulation. Shown are density
(color) and velocity (arrows; middle and right panels) distributions on the
$z=0$ (left) and $y=0$ (middle and right) plane at three different spatial
scales. The elapsed time after the onset of prestellar cloud $t$ and that
after protostar formation $t_{\rm ps}$ are described in the left panel.

##### C

Stellar system formed from metal-free gas

Next we see the nature of first stellar systems formed from the primordial gas
and discuss possible formation of first star binaries. The collapse of a
slightly gravitationally unstable gas cloud ($n_{\rm H}=1.4\times 10^{4}{\rm
cm^{-3}}$ and $T=200$K) is simulated by the SPH method. Armed with elaborate
sink particle creation scheme and usage of the barotropic equation of state
pre-calculated by a one-zone model, we succeeded in extending the time
integration up to 25000 yr in the accretion phase, $\sim 17$ times longer than
the previous calculation (Susa, 2019).

Figure 4 shows the number of fragments, i.e., sinks, as functions of the
scaled time normalized by the free-fall time at the density $n_{\rm
th}=10^{15}{\rm cm}^{-3}$, above which the temperature rises adiabatically.
The results can be scaled using the scaled time according to the scale-free
nature of gravity (Susa, 2019). Note that the integration time of the present
simulation is also the longest so far in terms of the scaled time. We can see
a simple trend that the number of fragments is proportional to the elapsed
time after the central density exceeds $n_{\rm th}$. This trend has already
been known previously (Susa, 2019), but now we confirm that the relation holds
even after 10 times longer integration time, although the simulation is done
for one particular realization. This result is also generally consistent with
results of other simulations in the literature.

Figure 4:  Number of fragments in a first-star-forming cloud as function of
the scaled time. We assume $\rho_{\rm ad}$ is $2\times
10^{-5}{\rm(g/cm^{-3})}$, above which the gas behaves adiabatic. Red line with
small triangles denotes the present result, and the thin lines correspond to
our previous result(Susa, 2019). Results in the literature(Stacy et al., 2012;
Smith et al., 2011; Clark et al., 2011; Greif et al., 2011, 2012; Vorobyov et
al., 2013; Susa, 2013; Machida & Nakamura, 2015; Hartwig et al., 2015a; Stacy
et al., 2016; Hirano & Bromm, 2017) are also shown. Figure 5:  The separations
of binaries identified in the present simulation as functions of time.
Different colors denote different binaries. Thick red solid curve shows the
mean square average of the distance between the stars and the center of mass
of the stellar system whereas the green thin curve shows the distance in
proportion to the elapsed time since the first sink formation.

Several systems of binaries are identified in this cluster. Figure 5 shows the
time evolution of their separations. Different colors correspond to different
pair of binaries. The members of the binaries easily change as a result of the
close encounter/merger with other stars. The thick red curve at the top of the
panel shows the typical size of this cluster, which grows in proportion to the
elapsed time after the first protostar formation, as is predicted by the
conservation of angular momentum and the gas distribution at the end of the
collapse phase. On the other hand, the binary separations also grow when the
binary systems acquire the angular momentum as a result of the gas accretion,
which is, however, counterbalanced by the effect of encounters with other
stars that makes the binary tighter. Hence, the separations are kept rather
constant during the computation. The separations observed in this particular
run are $\sim$ a few 10 AU, which is too wide to end up with BBHs that
coalesce within a Hubble time. In the present simulation, however, we adopt a
lower threshold density than the true adiabatic density for computational
reasons and the actual binary separation should be scaled from our nominal
value. After the scaling, the time should be 100 times shorter and the
separation $\sim 70$ times less. Thus, our binary separation corresponds to a
few times $0.1$ AU after scaling, although the integrated time is only 250 yr.
If this separation is kept constant for another $\sim 5000$ yr when the
surrounding gas is swept away by the radiative feedback, these binaries could
be progenitors of merging BBHs observed in gravitational waves.

To draw firm conclusions about the nature of formed stellar systems, we need
to extend the numerical simulation in time until the radiative feedback from
the protostar terminates the accretion to the system at several thousand years
after the formation of the first protostar. We should also consider the
effects of turbulence and magnetic fields. In fact, the presence of a coherent
magnetic field in primordial environment induces a strong magnetic breaking
effect (Machida & Doi, 2013), leading to less fragmentation in the accretion
disk. Also recent simulations predict more fragmentation in case with strong
turbulence. Those are the issues to be addressed in the future.

####  2.2.2 Binary black hole formation by dynamical interactions in dense
star clusters

##### A

Binary black holes formed in star clusters

Dense star clusters are one of promising formation sites of BBHs. Especially,
globular clusters (GCs) have been expected as a source of BBHs merging within
the Hubble time. In these days, most studies on BBH formation in GCs are being
done with the Monte-Carlo method (Rodriguez et al., 2015; Askar et al., 2017).
Those works suggested a local merger rate density of $\sim 5$ Gpc-3 yr-1
(Askar et al., 2017; Rodriguez et al., 2016). They also predicted BBHs formed
in GCs have several distinct features from those formed under isolated
environments. The GC-originated BBHs have more isotropic spin-orbit
misalignment angles (Rodriguez et al., 2016). Some of those BBHs may contain
massive ($\sim 100M_{\odot}$) BHs locating in the so-called “the pair-
instability mass gap” (Rodriguez et al., 2019), or may have finite
eccentricities at a GW frequency of $\sim 10$ Hz (Samsing et al., 2018;
Rodriguez et al., 2018). Such Monte-Carlo simulations, however, need careful
calibrations of some parameters from comparison with the $N$-body results
(Giersz et al., 2013; Rodriguez et al., 2016). On the other hand, direct
$N$-body simulations with a realistic number of particles (i.e., $N\sim
10^{6}$) are still numerically too expensive with cost of $N^{2}$ and with
very small binary orbital period compared with the GC lifetime. As a result,
currently, there is only a single example of million-body runs performed using
a GPU cluster (Wang et al., 2016), although the situation may change in near
future by adoption of tree-direct hybrid methods (Oshino et al., 2011; Iwasawa
et al., 2015) and novel schemes for binary integration (Wang et al., 2020).

On the other hand, open clusters (OCs) have not been expected to produce such
tight BBHs because they are less massive and less dense compared with GCs.
From the mass function of star clusters, however, there are hundred times more
OCs than GCs. Therefore, we thought it might be interesting to examine the
formation rate of BBHs in OCs seriously.

We performed a series of $N$-body simulations of OCs using a direct $N$-body
code, NBODY6++GPU, which includes both stellar and binary evolution models
Wang et al. (2015). As an initial condition, we adopted a Plummer model
Plummer (1911) with a total mass of $2.5\times 10^{3}M_{\odot}$ with 0.1 solar
metallicity. We calculated 360 realizations per a given set of the parameters
and obtained $\sim 1$ BBH per cluster on average. The merger time of the BBHs
due to gravitational wave emission is given as a function of the semi-major
axis ($a$) and eccentricity ($e$) of BBHs Peters&Mathews (1963). In Fig. 6,
the distributions of $a$ and $e$ of BBHs formed in the simulations are shown.
We can see that most of BBHs distribute in red-colored region (Region 3),
which correspond to those dynamically formed and the characteristic value of
their semi-major axis is set by the condition that the potential depth of
binaries is of the order of that of the cluster. Owing to shallow potential
depth of OCs, their separation is too large for the BBHs to merge within the
Hubble time. On the other hand, we found some BBHs that can merge within the
Hubble time in blue-colored region (Region 1), in which BBHs have very low
eccentricity and small semi-major axis. Those are formed via binary stellar
evolution such as stable and unstable mass transfer. Some of them experience a
dynamical encounter after the binary stellar evolution and have large
eccentricity (yellow-colored, Region 2). Three formation channels above are
summarized in Fig. 6.

In a GC, more BBHs are dynamically formed and, owing to the deeper potential
well, they tend to have smaller value of semi-major axis and some of them
merge within the Hubble time. On the other hand, in open clusters, the major
channel for BBH mergers is different: a binary consisting of two massive main-
sequence (MS) stars are first formed and becomes tight as a result of the
binary stellar evolution. If this separation is kept until BBHs are formed
after the stellar lifetime, they can merge within the Hubble time.

From our results, we can estimate a merger rate of BBHs formed in OCs at 10
Gyr ago to be 0.3 yr-1 Gpc-3 Kumamoto et al. (2019). This value corresponds to
20–50 % of those formed in GCs citeRodriguez16. This means BBHs formed in OCs
contribute to GW events much more than previously expected.

Figure 6:  Semi-major axis ($a$) and eccentricity ($e$) distribution of BBHs
obtained from our simulations (from Fig. 1 of (Kumamoto et al., 2019)). Star
symbols show BBHs with a merger time less than 10 Gyr. Colors of the symbols
indicate the primary mass ($M_{1}$) of the BBHs. Colored background indicates
classification depending on their formation mechanisms illustrated in 7.
Colored curves show the relation between $a$ and $e$ to give $t_{\rm GW}=10$
Gyr for BBHs with masses of 15, 25, 35, and 45$M_{\odot}$ (red, yellow, green,
and blue).  Figure 7:  Schematic figure of the formation processes of three
different types of BBHs. (From Fig. 3. of (Kumamoto et al., 2019))

We also repeat similar simulations but with different metallicities ranging
from zero to one solar metallicity. Using the BBH formation rate and
distribution of their orbital parameters, and assuming cosmic metallicity
evolution with some dispersion Chruslinka&Melemans (2019) and a cosmic star
formation rate density (Madau&Fragos, 2017), we can estimate the merger rate
of BBHs formed in OCs at the present universe by summing up the contributions
from all OCs. The estimated event rate is $\sim 70~{}{\rm yr}^{-1}{\rm
Gpc}^{-3}$ Kumamoto et al. (2020), comparable to that estimated from
LIGO/Virgo O2 and O3, $64.0_{-33.0}^{+73.5}~{}{\rm yr}^{-1}{\rm Gpc}^{-3}$
Abbott et al. (2019).

Our simulations also give a prediction on BH-MS binary distribution observable
with Gaia (Gaia Collaboration, 2016). In the next data release of Gaia (Gaia
DR3), binaries will be included. Using the results of our simulations, we
estimated that $\sim 10$ BH-MS binaries would be detectable with Gaia
Shikauchi et al. (2020), even if we assume relatively stringent observational
constraints Yamaguchi et al. (2018). In contrast to BH-MS binary formed in
galactic fields, those formed in OCs would have smaller mass of the MS
companion ($<5M_{\odot}$), longer orbital period ($>1$ yr), and higher
eccentricity ($>0.1$).

##### B

Binary population synthesis model of extremely low-metallicity star

Binary population synthesis (BPS) is a powerful tool to predict merger rates
of BBHs formed in galactic fields (e.g., Belczynski et al. (2016)). Coupled
with star cluster simulations, BPS can also predict the merger rate of BBHs
formed in dense stellar clusters (e.g., Rodriguez et al. (2016)). So far, most
studies have focused only on BBHs formed under Pop I and II environments
($0.01\lesssim Z/Z_{\odot}\lesssim 1$). On the other hand, it has been
suggested that Pop III-originated BBHs have distinct features from Pop I/II-
originated BBHs (Kinugawa et al., 2014). Thus, it is important to bridge the
metallicity gap (i.e. $0<Z/Z_{\odot}\lesssim 0.01$) in order to elucidate
origins of merging BBHs. Hereafter, we call stars in this metallicity gap the
extreme metal-poor (EMP) stars. In order to bridge the metallicity gap, we
have made a table for evolutionary tracks of stars with
$-8\leq\log(Z/Z_{\odot})\leq-2$. Note that there has been no evolution track
computation for EMP stars before. We have compiled those evolutionary tracks
as a set of fitting formulae, and incorporated them into BSE and NBODY6
(Hurley et al. (2002); Aarseth (2003), respectively). The BSE code is based on
several BPS calculations (e.g., Belczynski et al. (2016); Kinugawa et al.
(2014)), and the NBODY6 code and its derivatives are widely used for obtaining
BBH merger rates in dense star clusters (e.g., Banerjee et al. (2010);
Tanikawa (2013); Bae et al. (2014); Fujii et al. (2017); Park et al. (2017);
Hong et al. (2018); Kumamoto et al. (2019)). Here, we describe these fitting
formulae briefly (see Tanikawa et al. (2019), in detail).

We first make stellar evolution models as a reference for fitting formulae by
means of 1 dimensional (1D) simulation. For the 1D simulation, we use the
HOSHI code (Takahashi et al., 2016, 2018, 2019; Yoshida et al., 2019). We
follow the time evolution of stars with $8\leq M/M_{\odot}\leq 160$ for
$-8\leq\log(Z/Z_{\odot})\leq-2$. We start the calculation from the zero-age
main-sequence (ZAMS) and stop the calculation at the carbon-ignition time. We
do not include stellar wind mass loss. We take into account the post carbon-
ignition evolution (e.g. supernova) and stellar wind mass loss as post
processing. This is the same as in the BSE and NBODY6 codes incorporating the
evolution of Pop I/II stars. The evolution of EMP stars differs from that of
Pop I/II in the following respects: neither Hertzsprung-gap (HG) nor blue-loop
phase exists and the envelope of massive ($M\sim 20M_{\odot}$) stars remains
radiative until the end of their lives. This has significant impact on binary
evolution, such as the mass accretion rate in the Roche-lobe overflow phase,
and occurrence of post-MS mergers and common envelope evolution.

Figure 8: HR diagram for comparison between stellar evolution models (black)
and fitting formulae (red). The stellar evolution in this figure ends at the
carbon-ignition time, and stellar wind mass loss is not taken into account.

We use the above stellar evolution models as reference, and make fitting
formulae of EMP stars. Figure 8 shows stellar evolution on HR diagram in the
cases of $\log(Z/Z_{\odot})=-5$ and $-8$ for comparison with the results
obtained from 1D calculations and the fitting formulae. Note that the stellar
evolution in this figure is terminated at the carbon-ignition time, and
stellar wind mass loss is not taken into account. The tracks by fitting
formulae show good agreement with those by 1D calculations, with the maximum
deviation less than $20$%. This error is satisfactory for BPS and star cluster
simulations coupled with BPS, considering larger uncertainties coming from
uncertainties in tidal interaction, common envelope evolution, stellar wind
mass loss, natal kicks of neutron star and black holes, supernova mass loss
including pair instability, pulsational pair instability supernovae, and so
on.

####  2.2.3 Outlook

To distinguish the possible formation channels described above, we need
further statistical constraints with increasing the number of BBHs identified
by GW detections. For instance, binary population synthesis models predict
that the Pop III progenitor stars provide characteristic features in resulting
chirp-mass distributions. A most prominent feature is a peak around
$30M_{\odot}$ (Kinugawa et al., 2014), which may be still consistent with the
latest LIGO/Virgo O3a results. Future GW observations will impose more
constraints on other properties of BBHs, such as mass ratio, separation, spin,
and eccentricity, etc. Maturing both in-situ and ex-situ formation channels is
indispensable for understanding the origins of BBHs, making full use of the
observational constraints.

##  3 Binary Neutron Stars

###  3.1 Exploring NS physics

####  3.1.1 BNS waveform modeling based on numerical-relativity simulations

GWs from binary neutron stars (BNSs) contain rich information of the NSs. In
particular, the tidal deformability of NSs has been proposed as one of the
most promising quantities related to the equation of state that can be
extracted from the GW observation Lai:1993pa ; Mora:2003wt ; Flanagan:2007ix ;
Read:2009yp ; Damour:2009wj ; Hinderer:2009ca ; Vines:2011ud ; Damour:2012yf ;
Bini:2012gu ; Favata:2013rwa ; Yagi:2013baa ; Read:2013zra ; Bini:2014zxa ;
Bernuzzi:2014owa ; Wade:2014vqa ; Lackey:2014fwa . The tidal deformability
$\lambda$ is defined through $Q_{ij}=-\lambda E_{ij}$, where $Q_{ij}$ and
$E_{ij}$ are the induced quadrupole moment and the external tidal field,
respectively. The tidal deformability can also be written in terms of the
radius of the neutron star, $R$, and its Love number, $k_{2}$, as
$\lambda=2k_{2}R^{5}$. Therefore, the behavior of $\lambda$ depends on the
inner structure of a NS, its EOS. For example, for the soft (stiff) EOSs, the
radius of a NS is small (large), therefore, the NS is less (largely) deformed,
i.e. $\lambda$ takes small (large) value. The simultaneous measurement of the
mass and the tidal deformability of the NSs provides a substantial constraint
on the equation of state of nuclear matter which is yet poorly understood
Lattimer:2012nd ; Lattimer:2015nhk .

To extract the tidal deformability of NSs from the observed signals, an
accurate theoretical waveform template is crucial. For this purpose, many
efforts have been made to model the gravitational waveforms of a NS merger in
the past few decades. The waveforms including the linear-order tidal effects
are derived by post-Newtonian (PN) calculation for the early inspiral stage
Flanagan:2007ix ; Vines:2011ud . Furthermore, the waveforms employing the
effective-one-body (EOB) formalism are derived to incorporate higher-order PN
effects Damour:2009wj ; Bini:2012gu ; Damour:2012yf ; Bini:2014zxa ;
Bernuzzi:2014owa ; Hinderer:2016eia ; Steinhoff:2016rfi ; Dietrich:2017feu ;
Nagar:2018zoe ; Messina:2019uby . However, these waveform model could not be
accurate enough for the estimation of the tidal deformability due to a
significant systematic error of the unknown higher-order PN terms
Favata:2013rwa ; Yagi:2013baa ; Wade:2014vqa ; Lackey:2014fwa . Furthermore,
the tidal correction to the waveforms is derived based on the PN calculation
which is only valid in the regime where the orbital separation is sufficiently
large. Since the tidal effect is most significant just before the merger, it
is necessary to examine the accuracy of these waveform models.

Numerical-relativity (NR) simulation is the unique method to predict the tidal
effects in a regime where the non-linear effect of hydrodynamics should be
taken into account in the framework of general relativity Thierfelder:2011yi ;
Baiotti:2011am ; Bernuzzi:2012ci ; Radice:2013hxh ; Hotokezaka:2015xka ;
Haas:2016cop ; Hotokezaka:2016bzh ; Dietrich:2017feu ; Dietrich:2017aum ;
Kiuchi:2017pte ; Dietrich:2019kaq . The recent progress in numerical
techniques and computational resources enabled us to perform the high-
precision NR simulations for BNS mergers and to generate the waveforms with
systematic errors sufficiently small for the GW data analysis.

In our recent work Kiuchi:2017pte ; Kawaguchi:2018gvj ; Kiuchi:2019kzt , we
perform NR simulations for NS-NS mergers systematically varying the mass and
the tidal deformability of NSs: 6 and 4 models that cover the symmetric mass
ratio, $\eta$, in the range of $0.247$–$0.25$, the binary with the total mass
of $\approx 2.7$ and $2.5\,M_{\odot}$, and 5 different EOSs which cover the
binary tidal deformability in the range of $\approx 300$–$1900$. Here, the
binary tidal deformability, ${\tilde{\Lambda}}$, is defined by

$\displaystyle{\tilde{\Lambda}}$
$\displaystyle=\frac{8}{13}\left[\left(1+7\eta-31\eta^{2}\right)\left(\Lambda_{1}+\Lambda_{2}\right)\right.$
$\displaystyle-\left.\sqrt{1-4\eta}\left(1+9\eta-11\eta^{2}\right)\left(\Lambda_{1}-\Lambda_{2}\right)\right],$
(3)

where $\Lambda_{i}=\lambda_{i}/m_{i}^{5}\,~{}(i=1,2)$ are the tidal
deformability of each individual star (Note that $m_{1}<m_{2}$). By this work
the waveforms from NS-NS mergers for more than $15$ inspiral orbits are
obtained for a wide range of binary parameters.

We carefully estimate the error budgets of the waveforms obtained by the NR
simulations. In our NR simulation code, the field equations are solved under
discretization. We found that the error due to the finite grid resolution is
the largest error budget among the error sources of the numerical simulation.
To estimate the error due to the finite grid resolution, we perform the
simulations with various gird resolutions and checked the convergence property
of the obtained waveform data. As the results, we find that the phase error of
the waveforms is a sub-radian order Kiuchi:2017pte ; Kawaguchi:2018gvj ;
Kiuchi:2019kzt .

Figure 9: Comparison between NR and SEOBNRv2T waveforms for the models with
the total mass of $\approx 2.5\,M_{\odot}$ and the symmetric mass ratio of
$0.247$. The top and bottom panels show the amplitude of NR and SEOBNRv2T
waveforms and the phase difference between them, respectively. Here, $D$ and
$m_{0}$ denote the distance to the source and the total mass of the binary.
The vertical dashed lines are the peak amplitude time when the GW amplitude
reaches a peak in each model.

Employing the obtained high-precision NR waveforms, we examined the validity
of the analytical waveform models for NS binary. We show that the latest
tidal-EOB (SEOBNRv2T Hinderer:2016eia ; Steinhoff:2016rfi ) waveforms can be
accurate even up to $\approx 3\,{\rm ms}$ before the onset of merger
Kiuchi:2017pte  (see Fig. 9 for an example). However, the phase difference
between the tidal-EOB waveforms and the NR results is still larger than
$\approx 1\,{\rm rad}$ after two NSs come into contact for the case that the
NS radii are larger than $\approx 13\,{\rm km}$. This finding indicated that
further improvement of the waveform model would be needed to suppress the
systematic error in the measurement of the tidal deformability.

For this purpose, we develop a phenomenological model for GWs from inspiraling
NSs based on our high-precision NR simulations. Our phenomenological waveform
model is derived in the frequency domain as in the Phenom-series for binary
black holes Khan:2015jqa  for convenience in data analysis. We decompose a
frequency-domain gravitational waveform, ${\tilde{h}}\left(f\right)$, into the
frequency-domain amplitude, $A\left(f\right)$, and phase,
$\Psi\left(f\right)$, (with an ambiguity in the origin of the phase) by

$\displaystyle{\tilde{h}}\left(f\right)=A\left(f\right){\rm
e}^{-i\Psi\left(f\right)},$ (4)

and we define the corrections of the NS tidal deformation to GW amplitude and
phase by

$\displaystyle A_{\rm tidal}\left(f\right)=A\left(f\right)-A_{\rm
BBH}\left(f\right)$ (5)

and

$\displaystyle\Psi_{\rm tidal}\left(f\right)=\Psi\left(f\right)-\Psi_{\rm
BBH}\left(f\right),$ (6)

respectively. Here, $A_{\rm BBH}\left(f\right)$ and $\Psi_{\rm
BBH}\left(f\right)$ are the GW amplitude and phase of a binary black hole with
the same mass as the BNS, respectively (hereafter referred to as the point-
particle parts). In this work, we employ the SEOBNRv2 waveforms
Taracchini:2013rva  as the fiducial point-particle part of GWs.

For the tidal part phase and amplitude of our phenomenological waveform model,
we employ the following function forms motivated by the 2.5 PN order post-
Newtonian formula Refs. Damour:2012yf  and Hotokezaka:2016bzh  for the tidal-
part amplitude and phase correction:

$\displaystyle\Psi^{\rm
tidal}=\frac{3}{128\eta}\left[-\frac{39}{2}{\tilde{\Lambda}}\left(1+a\,{\tilde{\Lambda}}^{2/3}x^{p}\right)\right]x^{5/2}$
$\displaystyle\times\left(1+\frac{3115}{1248}x-\pi
x^{3/2}+\frac{28024205}{3302208}x^{2}-\frac{4283}{1092}\pi x^{5/2}\right)$ (7)

for the phase correction and

$\displaystyle A^{\rm tidal}$
$\displaystyle=\sqrt{\frac{5\pi\eta}{24}}\frac{m_{0}^{2}}{D_{\rm
eff}}{\tilde{\Lambda}}x^{-7/4}$
$\displaystyle\times\left(-\frac{27}{16}x^{5}-\frac{449}{64}x^{6}-b\,x^{q}\right)$
(8)

for the amplitude correction. Here, $a$, $p$, $b$, and $q$ are the free
parameters of the models.

Since our NR waveforms only contain the waveforms of which frequency is larger
than $\approx 400\,{\rm Hz}$, they are too short to calibrate the
phenomenological waveform model only by themselves. Thus, for this purpose, we
construct hybrid waveforms employing our NR waveforms and the effective-one-
body waveforms of Refs. Hinderer:2016eia ; Steinhoff:2016rfi  (SEOBNRv2T) as
the high and low-frequency parts, respectively, and use them to calibrate the
phenomenological waveform model. The hybridization of the waveforms is
performed in the time-domain by the procedure described in Ref.
Hotokezaka:2016bzh ; Kawaguchi:2018gvj .

Figure 10: (Left) Difference in the tidal-part phase between the hybrid
waveforms and our analytical model for the various binary configurations (see
Kiuchi:2017pte  and Kawaguchi:2018gvj  for the details of the models). Phase
differences are plotted after the alignment in the frequency range of 10–1000
Hz. (Right) Relative difference of amplitude between the hybrid waveforms and
our analytical model.

To focus on the inspiral-phase waveform and to avoid the contamination from
the post-merger waveforms, which would have large uncertainties, we restrict
the GW frequency range in $10$–$1000\,{\rm Hz}$. The fitting parameters are
determined by employing the hybrid waveforms with the largest value of binary
tidal deformability in the models studied in this work. By performing the
least square fit with respect to the phase difference and relative difference
of the amplitude, we obtained $a=12.55$, $p=4.240$, $b=4251$, and $q=7.890$.
Employing these parameters, we find that the phase and amplitude corrections
to the frequency domain waveforms are reproduced within $0.1\,{\rm rad}$ and
$25\%$ for all the binary configurations studied in this work (see Fig. 10).

To validate our phenomenological waveform model more quantitatively from the
view point of the data analysis, we calculate the mismatch between the our
waveform models and the hybrid waveforms, ${\bar{F}}$, defined by

$\displaystyle{\bar{F}}=1-\max_{\phi_{0},t_{0}}\frac{\left({\tilde{h}}_{1}\middle|{\tilde{h}}_{2}{\rm
e}^{2\pi
ift_{0}+i\phi_{0}}\right)}{||{\tilde{h}}_{1}||\,||{\tilde{h}}_{2}||},$ (9)

where $(\cdot|\cdot)$ and $||\cdot||$ are defined by

$\displaystyle\left({\tilde{h}}_{1}\middle|{\tilde{h}}_{2}\right)=4{\rm
Re}\left[\int_{f_{\rm min}}^{f_{\rm
max}}\frac{{\tilde{h}}_{1}\left(f\right){\tilde{h}}^{*}_{2}\left(f\right)}{S_{\rm
n}\left(f\right)}df\right],$ (10)

and

$\displaystyle\rho=||{\tilde{h}}||=\sqrt{\left({\tilde{h}}\middle|{\tilde{h}}\right)}.$
(11)

Here $S_{\rm n}$ denotes the one-sided noise spectrum density of the detector,
and we employ the noise spectrum density of the ZERO_DETUNED_HIGH_POWER
configuration of advanced LIGO aLIGOnoise  for it. We found that the mismatch
between our phenomenological waveform model and the hybrid waveforms is always
within $10^{-5}$ Kawaguchi:2018gvj ; Kiuchi:2019kzt .

Figure 11: Tidal phase in the frequency domain normalized by the leading,
Newtonian (relative 5PN-order) tidal phase formula. Our analytical model
(referred to as Kyoto tidal) assuming ${\tilde{\Lambda}}=$ 1000 (dot-dashed,
blue), 400 (dashed, blue), and 100 (dotted, blue), the NRTidal model (solid,
red Dietrich:2017aum ) and the NRTidalv2 model (solid, cyan Dietrich:2019kaq )
are shown in the plot. Note that latter two models are independent of
${\tilde{\Lambda}}$ when normalized by the leading tidal phase.

Gravitational waveform models for BNSs based on NR waveforms are also derived
in Dietrich:2017aum  and Dietrich:2019kaq  in a similar manner. The main
difference between our and their work is the difference of the NR waveforms
and the tidal-EOB waveforms used for the model calibration. Figure 11 compares
our tidal phase correction model and those of Dietrich:2017aum  and
Dietrich:2019kaq  normalized by the leading-order PN term (see also
Narikawa:2019xng  and Dietrich:2019kaq ). Figure 11 shows that those models
agree with each other within $10$–$20\%$ for 10–1000 Hz. Indeed, those models
give consistent results for the current detector sensitivity as shown in the
next subsection (see also Dietrich:2019kaq ).

####  3.1.2 Measuring tidal deformability from GW

As mentioned above, NSs provide unique laboratories to study the properties of
ultra-dense matter. To constrain the equation of state of NSs effectively, we
need accurate waveform models for the data analysis. We assume that the
waveform of the GWs can be decomposed into the point-particle, spin and tidal
parts. The point-particle and the spin parts are the same as those of the
binary black hole case, on the other hand, tidal part is that of the BNS case
because of the matter effects of NSs. The standard waveform, so-called
TaylorF2 (hereafter TF2), is derived by the post-Newtonian (PN) approximation
Sathyaprakash:1991mt . The point-particle Buonanno:2009zt ; Blanchet:2013haa ,
and spin parts Bohe:2013cla ; Arun:2008kb ; Mikoczi:2005dn  are calculated up
to the 3.5PN-order, while the effect of the tidal part (let us call it
PNTidal) enters from the 5PN-order, calculated up to the 7.5 PN order
Damour:2012yf . Since the tidal effects become important at the late inspiral
stage where the PN approximation becomes invalid, we need better waveform
models to estimate tidal deformability accurately.

Dietrich et al. Dietrich:2017aum ; Dietrich:2019kaq  (NRTidal) and Kawaguchi
et al. Kawaguchi:2018gvj  (KyotoTidal) independently construct new waveform
models in which the tidal part is improved by the calibration using NR
waveforms. In KyotoTidal, nonlinear effects of $\tilde{\Lambda}$ is taken into
account. The lack of the point-particle and spin parts beyond the 4PN-order
also bias the measurement of binary quantities. Kawaguchi et al. have also
constructed higher PN-order correction terms of the point-particle part by
fitting the SEOBNRv2 waveform model Purrer:2015tud ; Taracchini:2013rva , so-
called TF2+. Here, TF2+_KyotoTidal is calibrated in the frequency range of
10-1000 Hz, because they focus on the inspiral phase to avoid uncertainties in
the post-merger phase. Several studies have also derived constraints on the NS
EOS via measuring tidal deformability from GW170817 Margalit:2017dij ;
Bauswein:2017vtn ; Ruiz:2017due ; Annala:2017llu ; Zhou:2017pha ;
Fattoyev:2017jql ; Paschalidis:2017qmb ; Nandi:2017rhy ; Most:2018hfd ;
Raithel:2018ncd ; Landry:2018prl .

Narikawa et al. separately analyze the GW data of a binary-neutron-star
merger, GW170817, by each of the Advanced LIGO twin detectors Narikawa:2018yzt
. They found that the posterior probability distributions of the binary tidal
deformability for the Hanford and Livingston detectors are distinctively
different as shown in Fig. 12. They have suggested that the difference between
the detectors cannot be understood by statistical error only, since only the
posterior with the Livingston detector does not change smoothly as the upper
cutoff frequency increases. Their findings suggest that further research of
the noise properties in the high-frequency region of the Livingston data is
needed to extract information on the tidal deformability of GW170817. Narikawa
et al. also found that there is a difference in the estimates of
$\tilde{\Lambda}$ for GW170817 between NR calibrated waveform models (the
KyotoTidal and NRTidalv2 models) as shown in Fig. 13 Narikawa:2019xng . Here,
NRTidalv2 model is an upgrade of the NRTidal model Dietrich:2019kaq . The
order of peak values of $\tilde{\Lambda}$ for the different waveform models in
Fig. 13 can be understood by the difference in the magnitude of the tidal
phase at around $\tilde{\Lambda}\approx 600-900$ shown in Fig. 11. However,
the difference is smaller than the statistical error.

The second BNS merger event GW190425 was observed during O3a Abbott:2020uma .
The total mass of the system was about 3.4$M_{\odot}$, which is larger than
the total mass of any other known BNS system. Since the LIGO Hanford detector
was offline at the time of the event, parameter estimation was performed by
using the data from the LIGO Livingston and Virgo detectors. The EOS with
large tidal deformability $\tilde{\Lambda}>1100$ was disfavored, which is
consistent with the several analyses done for GW170817. For the low-spin prior
(component spin is enforced as $\chi_{1,2}\sim U[-0.05,~{}0.05]$), the 90%
upper limit of GW190425 becomes much smaller than that of GW170817,
$\tilde{\Lambda}\leq 600$. Since the total mass of the system was large, the
possibility of BH-NS system is discussed in Abbott:2020uma .

Figure 12: Marginalized posterior probability distribution of binary tidal
deformability, $\tilde{\Lambda}$ for GW170817, derived by data of different
detectors with $f_{\rm max}$ = 2048 Hz. (Reprinted with permission from
Narikawa:2018yzt . ©(2017) by the American Physical Society.)  Figure 13:
Marginalized posterior PDFs of binary tidal deformability, $\tilde{\Lambda}$
for GW170817, estimated by different waveform models for $f_{\rm
max}=1000~{}{\rm Hz}$.

###  3.2 Association of a short-duration gamma-ray burst

####  3.2.1 High-energy Observations

The observations of GW170817 by high-energy instruments, Fermi, MAXI, Swift,
CALET and Chandra, are reviewed in this subsection.

##### Fermi

On August 17, 2017, at 12:41:06.5 (UTC), hard X-ray emission from GRB 170817A
triggered Fermi Gamma-ray Burst Monitor (GBM with a coverage of 8keV–40MeV;
Meegan, et al. (2009)) near the detection limit (Goldstein:2017mmi, ). The
emission with a short hard pulse and a soft tail lasted for $\sim$2 s, as
shown in Fig. 14. From the obtained time duration, this GRB is classified as a
short GRB. Also the SPectrometer on board INTEGRAL Anti-Coincidence Shield
(SPI-ACS) detected a hard pulse completely coincident with the Fermi/GBM one
Savchenko:2017ffs .

Amazingly, the binary coalescence detected by the Laser Interferometer
Gravitational-wave Observatory (LIGO) occurred $\sim$1.7s before the Fermi/GBM
trigger (GW170817, ) and the sky position determined by the LIGO GW
observation is consistent with the Fermi/GBM one. The gamma-ray spectrum for
the initial hard pulse is well fitted by a power-law function with a photon
index of $\alpha$ = -1.6$\pm$0.4 and an exponential cutoff at $E_{\rm peak}$ =
185 $\pm$ 62keV. For the soft tail emission, a blackbody with $k_{\rm B}T$ =
10.3 $\pm$ 1.5keV well represents the spectrum. The fluence over the total
duration is $\sim$3 $\times$ 10-7erg cm-2, which falls in the $\sim$50th
percentile of the fluence distribution obtained by Fermi/GBM
(Goldstein:2017mmi, ). Assuming a distance to the host galaxy NGC 4993 as
43Mpc, the isotropic equivalent energy and peak luminosity in the 1keV – 10MeV
band are $\sim$3 $\times$ 1046erg and $\sim$2 $\times$ 1047erg s-1,
respectively, which are by 3–6 orders of magnitude smaller than those of other
short GRBs detected by Fermi/GBM. The observed features for GRB 170817A would
suggest that the hard X-ray emission came from a misaligned jet to the
observer (i.e., off-axis jet).

The higher-energy gamma-ray observation of GRB 170817A by Fermi Large Area
Telescope (LAT; Atwood, et al. (2009)) covering the 0.1–300GeV band was
unfortunately not available due to entry to the South Atlantic Anomaly at the
time of occurrence of the merger time Ajello et al. (2018). Fermi/LAT resumed
scientific operation $\sim$103s after the merger time and put a flux upper
limit of 4.5 $\times$ 10-10erg cm-2 s-1 (95% confidence level) in the 0.1–1GeV
band using the initial 1000s observation after the operation resumption.

Motivated by study of GRB 170817A, investigation of previous short GRBs with a
similar signature has been performed (Burns, et al., 2018; Troja, et al.,
2018; von Kienlin, et al., 2019; Matsumoto & Piran, 2020). In particular, the
nearby short GRB 150101B detected by Fermi/GBM (z = 0.134 Levan, et al.
(2015)), which could accompany the kilonova and the off-axis jet emission
suggested by the optical and X-ray observations Troja, et al. (2018), has a
hard spike and a soft tail. Thus, the two-component sigunature could be a
common feature for short GRBs Burns, et al. (2018). The isotropic equivalent
energy and luminosity in the 1keV – 10MeV band are $\sim$2 $\times$ 1049erg
and $\sim$4 $\times$ 1050erg s-1, respectively, which are in a fiducial range
of short GRBs. One of the possible interpretation for these phenomena suggests
that GRB 150101B could originate from a more on-axis jet than GRB 170817A.

Figure 14:  Light curve of GRB 170817A detected by GBM in the 8–50keV and
50–300keV bands. Emission consists of a hard pulse and a weak soft tail. The
vertical solid line represents the time of the binary neutron star merger. The
horizontal dashed line represents the averaged background level.

##### MAXI

Monitor of All-sky X-ray Image (MAXI; 2009PASJ…61..999M ) is an instrument for
monitoring X-ray sky, which was launched in 2009 and has been observing from
the International Space Station (ISS). Because MAXI does not have a moving
mechanism and is mounted to the ISS, we cannot actively control the pointing
direction. The orbital period of the ISS is 92 minutes and thus an X-ray
source is usually observed once in the same period. Since MAXI covers about
85% of the sky every orbit, the source activity can trace back to the time
before the event.

In studying X-ray counterparts of GW events, we use the Gas Slit Camera (GSC;
2011PASJ…63S.623M ) GSC consists of 12 position sensitive proportional
counters. Six of them covers 3 degree $\times$ 160 degree field of view (FOV)
toward the horizontal direction and the others compose zenithal FOV with the
same dimension. These FOVs move as the ISS rotates. A scan observation
(transit) for a point source lasts 40–150 seconds depending on the source
position in the FOV 2011PASJ…63S.635S .

The effective area of the GSC for a source changes with time linearly like a
triangle. For a source with constant flux (or that can be assumed as constant
in a scan transit), we calculate an average flux by dividing the source count
by the effective exposure, that is a sum of the effective area along the time.

For GW170817 MAXI did not detect the X-ray emission. The upper limits are
summarized in 2018PASJ…70…81S . Fig. 15 shows the luminosity and upper limits
for early ($<10$ days) X-ray afterglows swift:phil ; 2019ApJ…886L..17H  in the
case of GW170817. The position of the EM counterpart lay at the gap between
the horizontal and the zenithal cameras for the first three orbits and the
first scan started after 4.7 hours. Moreover the first three observations were
incomplete and the upper limits were 1–2 order of magnitude higher than the
usual. However, it is important to mention that this was the earliest
observation of the X-ray counterpart of the GW event.

Although the observation start time was late in this case, if MAXI’s
observation starts $\sim$ 100 seconds after the trigger, MAXI could observe an
X-ray afterglow equivalent to the level observed by Swift at several hours
after the trigger.

The best lesson of the observation of GW170817 was that it was important to
continuously observe the sky as long as possible. In O3 we have increased the
observation efficiency to reduce the gap in the observation.

Figure 15:  The observed luminosity and upper limits of the X-ray counterpart
of GW170817 in early time. MAXI upper limits are based on the result of
2018PASJ…70…81S . The first six points are for each scan transit and the last
two are for one day and 10-day observations. The data of Swift XRT and Chandra
are converted to the luminosity at 2-10keV assuming the spectral parameters in
their original papers swift:phil ; 2019ApJ…886L..17H . The dotted, dashed, and
dash-dotted lines are models for luminosity $L$ as a function of the time $t$,
$L\propto t^{\alpha}$, and $\alpha$ are -1, -1.5, and -2, respectively. The
horizontal gray line is a typical GSC upper limit for one scan transit, scaled
to the distance of GW170817 (40Mpc).

##### Swift

The Neil Gehrels Swift Observatory (Swift; (swift:neil, )) consists of three
scientific instruments. The Burst Alert Telescope (BAT; (swift:bat, )) is a
wide field of view telescope monitoring about 1/6 of the sky in the 15–150keV
band. The X-ray Telescope (XRT; (swift:xrt, )) is an X-ray telescope covering
the energy range in the 0.3–10keV band in 24′ field of view. The UV-Optical
Telescope (UVOT; (swift:uvot, )) observes the wavelength range of 170–600nm in
17′ field of view.

At the merger time of GW170817, the error region of GW170817 was occulted by
the Earth and not visible by BAT (swift:bat_gw170817, ; swift:phil, ).
Therefore, no useful information can be derived from the BAT data. The XRT and
UVOT started their observations about an hour after the merger time around the
center of the Fermi GBM error region ($\sim 1^{\circ}$ radius). After the
localization information became available from the LIGO and Virgo at 4.8 hours
after the merger time, Swift started a series of 120s observations at the
positions of known galaxies inside the 33.6 deg2 error region of the
gravitational-wave detectors. No new X-ray sources ($\geq$10-12erg
cm${}^{-2}\,$s-1) were found during those observations (swift:phil, ).

Swift started the observations at the position of the optical transient of
GW170817 about 14.4 hours after the merger time. Although the XRT found no
X-ray from the position of the transient, the UVOT detected a fast decaying UV
emission which is interpreted as the blue kilonova associated with GW170817.
The non-detection of X-ray at the early phase by XRT places a crucial
difference to a typical short GRB afterglow which is dominated by an on-axis
emission (swift:phil, ).

##### CALET

The Calorimetric Electron Telescope (CALET; (calet:torii, ; calet:asaoka, ))
is the scientific instrument attached on the International Space Station
(ISS). The CALET Gamma-ray Burst Monitor (CGBM;(calet:yamaoka, )) is the
gamma-ray burst monitor onboard CALET and is capable to observe emissions from
7keV to 1MeV (Hard X-ray Monitor; HXM, FoV of $\sim$$60^{\circ}$ from the
boresight) and from 40keV to 20MeV (Soft Gamma-ray Monitor; SGM, FoV of
$\sim$$110^{\circ}$ from the boresight) utilizing two different detectors. No
significant emission was observed by CGBM at the time of the merger of
GW170817. The 90% upper limit is $1.3\times 10^{-7}$erg cm${}^{-2}\,$s-1 in
the 10-1000keV band assuming no shielding by the structure of ISS
(Monitor:2017mdv, ).

In addition, the main CALET instrument, the high-energy Calorimeter (CAL),
observes gamma-rays from $\sim$1GeV up to 10TeV with a field of view of nearly
2sr. The field of view of CAL did not include the location of GW170817 at the
time of the merger. A search for delayed emission over the 2-month period
following the event resulted in 90% C. L. upper limits on the energy flux of
$1.2\times 10^{-10}$erg cm${}^{-2}\,$s-1 for gamma-rays above 1GeV and
$4.0\times 10^{-10}$erg cm${}^{-2}\,$s-1 above 10GeV (calet:cal, ).

##### Chandra

The Chandra X-ray Observatory, which provides the best sensitivity in X-ray
for a point source, first targeted to the optical transient of GW170817 at 2.2
days after the merger time. Although no X-ray emission was detected in this
first visit, Chandra detected the significant X-ray emission on the 2nd visit
at 8.9 days after the merger in the isotropic luminosity of $9\times
10^{38}$erg s-1 at a distance of 40Mpc chan:nora1 .

According to the year-long monitoring observations by Chandra, the X-ray flux
continues to increase up to $\sim$160 days after the merger time by a power-
law index of 0.9 vanEerten:2018vgj ; chan:hajela . After the peak, the flux
shows the decline by a power-law index of $\sim-2$ (Fig. 16). There is a hint
of an X-ray flaring feature around 155 days after the merger time chan:piro .
Chandra has been collecting the data up to 2.5 years after the merger and
continues monitoring. This unique X-ray light curve behavior along with the
radio data will be discussed in the following section.

Figure 16:  The X-ray light curve of the counterpart of GW170817 observed by
Chandra. The source count-rates are converted to the unabsorbed flux in the
0.3-8keV band using the Galactic absorption column $N_{H}=7.6\times
10^{20}$cm-2 and the photon index of 1.59 vanEerten:2018vgj .

####  3.2.2 Theoretical Interpretation

The GW event GW170817 (GW170817, ) and the associated short GRB 170817A
(Monitor:2017mdv, ; Goldstein:2017mmi, ; Savchenko:2017ffs, ) have provided,
for the first time, clear evidence for a relativistic jet from a binary
neutron star merger observed from off-axis (as schematically shown in Fig.
17). The observations also show that the jet is structured, excluding a top-
hat distribution of energy and Lorentz factor of the jet. We also discuss the
current issues about the exact structure of the jet, the origin of the gamma-
ray emission, and future prospects.

Figure 17:  Electromagnetic counterparts associated with the off-axis jet from
short GRB 170817A (from Fig. 1 of Ref. (Ioka:2017nzl, )).

##### A

Evidence of an off-axis jet

The afterglow observations of GW170817A verify the existence of a relativistic
jet: the superluminal motion of the radio image indicates a relativistic
collimated source (Mooley:2018qfh, ; Ghirlanda:2018uyx, ) and the closure
relation between the spectral index and the light curve slope after the
luminosity peak is consistent with the standard afterglow model of a
relativistic jet (vanEerten:2018vgj, ; Mooley:2018clx, ; Lamb:2018qfn, ),
excluding a failed-jet scenario that the jet is choked by the ejecta from the
neutron star merger (Kasliwal:2017ngb, ; Gottlieb:2017pju, ; Nakar:2018cbe, ).

The jet should be strong enough for successfully breaking out the merger
ejecta in this event (Nagakura:2014hza, ; Murguia-Berthier:2014pta, ). The
breakout time since the jet launch ($t=t_{0}$) is approximately estimated as

$\displaystyle t_{b}-t_{0}\sim 0.17\,{\rm s}\left(\frac{t_{0}-t_{m}}{0.1\,{\rm
s}}\right)^{1/2}\left(\frac{L_{{\rm iso},0}}{10^{51}\,{\rm erg}\,{\rm
s}^{-1}}\right)^{-1/2}+0.078\,{\rm s}\left(\frac{L_{{\rm
iso},0}}{10^{51}\,{\rm erg}\,{\rm s}^{-1}}\right)^{-1},$ (12)

where $L_{{\rm iso},0}$ is the isotropic luminosity at the base of the jet and
$t_{0}-t_{m}$ is the time delay of the jet launch since the merger (modifying
Eq. (32) of Ref. (Hamidani:2019qyx, )). Here the second term is necessary for
this event because the ejecta expansion velocity is comparable to the jet head
velocity, in contrast to collapsars with a static envelope. In order to have
the short GRB 170817A $\sim 1.7$s after the merger, the jet luminosity should
be large enough $L_{{\rm iso},0}\gtrsim 3\times 10^{49}$erg s-1 and the delay
time should be short enough $t_{0}-t_{m}\lesssim 1.3$ s (Hamidani:2019qyx, ).

The observed isotropic energy of the gamma-ray emission is very small, several
orders of magnitude smaller than typical (Monitor:2017mdv, ;
Goldstein:2017mmi, ; Savchenko:2017ffs, ). Therefore the jet should be off-
axis for a reasonable radiative efficiency (Ioka:2017nzl, ; Ioka:2019jljb, ).
An off-axis observer receives photons emitted outside the relativistic beaming
cone, and the apparent energy of the off-axis jet becomes faint (see also
below).

##### B

Diverse jet structures

The afterglow observations, in particular the slowly rising light curves from
radio to X-ray, are not consistent with a top-hat jet (Mooley:2017enz, ), but
indicate a structured jet (DAvanzo:2018zyz, ; Ghirlanda:2018uyx, ;
Lazzati:2017zsj, ; Lyman:2018qjg, ; Margutti:2018xqd, ; Ruan:2017bha, ;
Troja:2018ruz, ; vanEerten:2018vgj, ; Lamb:2018qfn, ). These are the first
confirmation of previous suspicion that realistic GRB jets are structured
(Meszaros:1997je, ; Zhang:2001qt, ).

However, there remains confusion about the jet structure as different authors
give different structures, such as Gaussian (Lyman:2018qjg, ;
vanEerten:2018vgj, ; Lamb:2018qfn, ) and power-law profiles (DAvanzo:2018zyz,
; Ghirlanda:2018uyx, ) (see Fig. 1 of Ref. (Ioka:2017nzl, )). These structures
are obtained by assuming a functional form of the jet structure with model
parameters, and adjusting the parameters through fitting the theoretical light
curves to the observations.

We instead invent a novel method to determine the functional form of the jet
structure itself (Takahashi:2019otc, ). This method solves an inverse problem,
uniquely reconstructing a jet structure from a given off-axis GRB afterglow by
integrating an ordinary differerential equation, which is formulated based on
the standard theory of GRB afterglows. Applying the inversion method to GRB
170817A, we clarify that the current observational errors still allow various
jet structures, discovering by-product that a hollow-cone jet is also
consistent with the observations as shown in Fig. 18.

Figure 18:  Left: Examples of the jet structure that are consistent with the
afterglow observations of GRB 170817A: a hollow-cone (magenta line) and
Gaussian jet (dashed line). The green region is not in principle determined by
the current data. Right: The corresponding light curves for a hollow-cone
(thick dashed lines) and Gaussian jet (thin dashed line) with the data. The
blue line is used for the inverse reconstruction (modifying Fig. 8 of Ref.
(Takahashi:2019otc, )).

An important caveat is that the current afterglow observations only constrain
the jet core (magenta in Fig. 18), not in principle the outer jet structure,
(green; which is crucial for the gamma-ray emission as discussed below). Early
afterglow observations are required to obtain the outer jet structure
(Takahashi:2019otc, ).

The origin of the jet structure is still unclear. First, the jet may have
intrinsic structure from the beginning of the launch. Second, the jet may be
structured during the propagation if the outer part loads baryon from the
ejecta. Third, the cocoon component (the collided jet and ejecta) may
contribute to the outer structure since a part of the cocoon is relativistic
due to partial mixing of the baryon. Early macronova/kilonova observations are
important for probing the cocoon, which has crucial information on the jet.

##### C

Origin of the gamma-ray emission

The prompt gamma-ray emission of GRB 170817A generally comes from an off-
center jet, neither the jet core nor the line of sight but the middle, as in
Fig. 19 (Ioka:2019jljb, ), given the jet structure with a rapidly decaying
tail as in Fig. 18. This is because the emission from the jet core is
suppressed by the relativistic de-beaming, while the line-of-sight emission
becomes too faint to dominate the off-axis emission from the inner jet. The
off-center emission solves seemingly puzzles that the observed prompt spectrum
is inconsistent with the spectral Amati relation (Amati:2002ny, ;
Yonetoku:2003gi, ) (because the off-center is different from the usually-
observed jet core) and causes the compactness problem (Matsumoto:2018xnd, ;
Matsumoto:2019sor, ) (because the off-centre jet is much less energetic and
much closer to the line of sight than the jet core).

Figure 19:  The surface brightness distribution of the prompt emission from a
Gaussian jet (from Fig. 3 of Ref. (Ioka:2019jljb, ))

The off-center structure still has large uncertainties because it is not
constrained by the observed afterglows (as shown by the green region in Fig.
18). Whether it is a jet or cocoon is also unknown. The observed gamma-rays
might be the jet emission scattered by the cocoon (Kisaka:2017mgz, ). Future
observations of off-axis events will reveal the jet structure, where roughly
$\sim 10\%$ events are expected to be brighter at smaller viewing angles than
GRB 170817A.

###  3.3 Optical/Infrared Counterparts and Heavy Element Nucleosynthesis

Figure 20:  Optical, infrared, and radio telescopes in the J-GEM network.

The origin of the elements in the Universe is a subject of interest in
astronomy and astrophysics. In particular, the origin of the elements
synthesized by rapid neutron capture, or so-called $r$-process, is one of the
long-standing, unsolved problems. Theoretically, NS mergers are one of the
promising candidate sites of $r$-process nucleosynthesis lattimer74 ;
lattimer76 ; eichler89 ; freiburghaus99 ; roberts11 ; goriely11 ; korobkin12 ;
bauswein13 ; wanajo14 . However, there has been no observational evidence of
$r$-process by NS mergers.

When NS mergers occur, a small fraction of mass of NSs is ejected into
interstellar space rosswog99 ; rosswog05 ; goriely11 ; rosswog13a ;
hotokezaka13 ; bauswein13 . In the ejected material, $r$-process
nucleosynthesis is expected to take place. Then, we expect optical and
infrared emission powered by radioactive decays of freshly synthesized nuclei
li98 ; kulkarni05 ; metzger10 ; kasen13 ; barnes13 ; tanaka13 . This
phenomenon is called “kilonova” or “macronova” (hereafter we call it
kilonova). In other words, if we can detect kilonova after detection of GWs
from a NS merger, we can test $r$-process nucleosynthesis by the NS merger.

The event rate of NS mergers can be measured from GW observations. The
ejection of $r$-process elements per event can be estimated from
electromagnetic (optical/infrared) observations. Therefore, by these multi-
messenger observations, we can study the production rate of $r$-process
elements by NS mergers and test if NS mergers can be the origin of $r$-process
elements in the Universe.

To achieve observations of kilonova, we coordinated the Japanese collaboration
for Gravitational wave ElectroMagnetic follow-up (J-GEM, Section 3.3.1).
Thanks to this observational network, we succeeded in observations of the
electromagnetic counterpart of GW170817 (Section 3.3.2). By combining our
state-of-art numerical simulations (Section 3.3.3), we advance our knowledge
about the physics of NS mergers and the origin of $r$-process elements in the
Universe.

####  3.3.1 J-GEM

J-GEM is a network of optical, infrared, and radio telescopes for EM follow-up
observations of GW sources. Figure 20 shows a summary of telescopes in the
J-GEM. The roles of the optical telescopes can be divided into two categories:
one is the wide-field surveys (Kiso, MOA-II, and Subaru Hyper Suprime-Cam
(HSC)) and galaxy targeted surveys (other telescopes). For more details of
each telescope, see Morokuma et al. (2016) morokuma16 . J-GEM signed
Memorandum of Understanding with LIGO/Virgo collaboration in 2014, and started
to receive the alert of GW detection from the first observing run of Advanced
LIGO (O1).

We demonstrated coordinated observations using this observing network for the
first GW source GW150914 (BBH, GW150914 ). Within 4 days from the GW detection
(2 days from the initial alert), about 24deg2 area has been observed with Kiso
Schmidt telescopes. Also, within 6-7 days from the GW detection, 18 nearby
galaxies were observed with the B&C 61-cm telescope. Although no EM
counterpart has been identified, this was the first attempt of multi-messenger
observations of GW sources morokuma16 .

More intensive observing campaign has been performed for the second GW source
GW151226 (BBH, GW151226PRL ). For this event, 238 nearby galaxies were
observed by galaxy targeted surveys yoshida17 , which is great improvement
compared with the observations of GW150914. Also, Kiso Schmidt telescope and
MOA-II telescopes covered 778deg2 and 145deg2 area, respectively yoshida17 .
Furthermore, Subaru/HSC covered 63.5deg2 area with unprecedented sensitivity
(limiting magnitudes of 24.6 and 23.8 for the $i$ and $z$ bands, respectively,
utsumi18 ). The total area covered by the wide-field surveys was 986.5deg2,
which corresponds to $\sim$29% of the probability map of GW151226. It is worth
emphasizing that the sensitivity of Subaru/HSC observations was deep enough to
detect expected kilonova emission even at about 200Mpc, which is about the
design sensitivity of Advanced LIGO, Advanced Virgo, and KAGRA.

Subaru/HSC is one of the most powerful wide-field imager: its field-of-view
(1.8deg2) is the largest among 8-10m class telescopes. Table 3 summarized our
follow-up observations using Subaru/HSC. Details of follow-up observations of
S190510g are presented by Ohgami et al. (2021) ohgami21 .

Table 3: Subaru/HSC observations for GW sources
1 Detected with low-letency pipeline but not recovered with the offline
analysis.
2 Initially reported as BNS with a high probability but can be terrestrial
noise.
3 Initially reported as Mas Gap object but later classified as BBH.  Target | Type | Disc. date | Obs. date | Filter | Area | Probability | Lim. mag
---|---|---|---|---|---|---|---
|  | (UUT) | (UT) |  | (deg2) | (%) |
GW151226 | BBH | 2015-12-26 | 2016-01-07 | HSC-i | 63.1 | 12 | 24.4
|  |  | 2016-01-07 | HSC-z | 63.3 | 12 | 24.0
|  |  | 2016-01-13 | HSC-i2 | 63.1 | 12 | 24.4
|  |  | 2016-01-13 | HSC-z | 63.2 | 12 | 24.0
|  |  | 2016-02-06 | HSC-i2 | 65.0 | 12 | 24.5
|  |  | 2016-02-06 | HSC-z | 65.0 | 12 | 24.0
G2754041 | BH-NS | 2017-02-25 | 2017-02-27 | HSC-z | 26.0 | 0.5 | 23.2
GW170814 | BBH | 2017-08-14 | 2017-08-16 | HSC-Y | 32.5 | 35 | 22.2
|  |  | 2017-08-17 | HSC-Y | 32.6 | 35 | 22.2
GW170817 | BNS | 2017-08-17 | 2017-08-18 | HSC-z | 26.2 | 61 | 21.7
|  |  | 2017-08-19 | HSC-z | 27.3 | 61 | 21.6
S190510g | BNS2 | 2019-05-10 | 2019-05-10 | HSC-Y | 118.8 | 1.1 | 22.9
S191216ap | BBH3 | 2019-12-16 | 2019-12-20 | HSC-z | 2.3 | 0.035 | 25.2
S200224ca | BBH | 2020-02-24 | 2020-02-25 | HSC-r2 | 57.8 | 80.8 | 25.3
|  |  |  | HSC-z | 57.8 | 80.8 | 23.6
|  |  | 2020-02-28 | HSC-r2 | 57.8 | 80.8 | 24.6
|  |  |  | HSC-z | 57.8 | 80.8 | 23.2
|  |  | 2020-03-23 | HSC-r2 | 57.8 | 80.8 | 25.5
|  |  |  | HSC-z | 57.8 | 80.8 | 23.8
Figure 21:  Left: Localization map of GW170817 and Subaru/HSC pointing
tominaga18 . Right: three-color image of the counterpart of GW170817 utsumi17
.

####  3.3.2 GW170817

The first GW detection from a NS merger has been achieved for GW170817
GW170817 . The detection triggered EM follow-up observations covering entire
wavelength range abbott17MMA . Thanks to the observations with three GW
detectors, the position of the GW source is localized to $\sim 30$deg2. In
optical/infrared wavelengths, a counterpart AT2017gfo was identified by
several groups andreoni17 ; arcavi17 ; chornock17 ; coulter17 ;
cowperthwaite17 ; diaz17 ; drout17 ; evans17 ; kasliwal17 ; kilpatrick17 ;
lipunov17 ; muccully17 ; nicholl17 ; pian17 ; shappee17 ; siebert17 ; smartt17
; soares-santos17 ; tanvir17 ; tominaga18 ; troja17 ; utsumi17 ; valenti17 .

We have performed wide-field survey using Subaru/HSC, covering 23.6deg2
tominaga18 . The observations recovered the detection of AT2017gfo.
Furthermore, thanks to the wide coverage and good sensitivity, we
statistically ruled out that other transients in the localization area are
associated with GW170817, or in other words, we concluded that AT2017gfo is
the most likely counterpart.

Intensive optical and infrared observations of AT2017gfo have been performed
in the framework of J-GEM network utsumi17 . Observations with Subaru/HSC show
that the $z$-band brightness quickly declines. On the other hand, the
observations with IRSF/SIRIUS shows that near-infrared light curves evolve
more slowly. This is fully consistent with the expected behavior of kilonova.

Using radiative transfer simulations, we find that the observed light curves
can be explained by $\sim 0.03M_{\odot}$ of ejecta including lanthanide
elements, which have a high opacity tanaka17a . Also, the blue optical
component, revealed by observations with Subaru/HSC, B&C/Tripol5, and MOA-
II/MOA-cam3, requires an ejecta component with a smaller fraction of
lanthanide. Since the estimated total ejecta mass is higher than the expected
ejection by the first dynamical ejection from the NS merger, we concluded that
the observed signals are mainly produced by post-merger mass ejection.

|
---|---
Figure 22:  Density distribution and electron fraction distribution in the
post-merger ejecta, obtained by viscous-radiation hydrodynamics simulations
fujibayashi18 .  Figure 23:  Multi-color light curves of GW170817/AT2017gfo
(circles) compared with the results of two-dimensional radiative transfer
simulations (solid lines) kawaguchi18 ; kawaguchi20

####  3.3.3 Advances of Theoretical Models

The observed properties of AT2017gfo are consistent with the expectation from
numerical relativity simulations shibata17 . In particular, to have an ejecta
component with a moderate lanthanide fraction, at least temporal presence of
hypermassive neutron star is suggested. If the merger results in the prompt
collapse to a BH, there is no source of neutrino radiation, and thus, almost
entire ejecta are expected to be lanthanide-rich, which is not consistent with
the observations. In this way, numerical relativity simulations play crucial
roles to connect the initial conditions estimated from GWs (mass and mass
ratio) with the final outcome of kilonova.

Mass ejection from NS mergers can be roughly divided into two phases: (1)
dynamical mass ejection mainly by the tidal force and shock interaction in the
dynamical time scale and (2) subsequent post-merger mass ejection from the
accretion torus. Observations of GW170817/AT2017gfo highlighted the importance
of post-merger mass ejection from NS mergers. However, the detailed
simulations have been difficult. We have performed long-term general
relativistic neutrino radiation hydrodynamics simulations, by taking into
account the effects of viscosity fujibayashi18 . We find that post-merger
ejecta ($>10^{-2}M_{\odot}$) dominate over the dynamical ejecta
($<10^{-2}M_{\odot}$). Thanks to the inclusion of neutrino transfer, we could
also provide the distribution of electron fraction in the ejecta. We showed
that the ejecta have high electron fraction ($>0.25$), which means that post-
merger ejecta is lanthanide-poor (Fig. 22).

Depending on the mass and mass ratios of the NS merger, various ejecta
components with different properties (mass, velocity, and electron fractions)
are expected. Therefore, it is important to study the observable outcome using
the results of numerical relativity. We have performed multi-dimensional
radiative transfer simulations by consistently taking the interplay of
multiple ejecta components into account kawaguchi18 . We showed that AT2017gfo
can be naturally explained by the geometry predicted by numerical relativity
simulations (Figure 22).

We then extended our radiative transfer simulations for various initial
conditions. Using the updated atomic data tanaka20 , Kawaguchi et al. (2020)
kawaguchi20  showed that the properties of kilonova should have diversity
depending mainly on the total mass of the system. Since this is multi-
dimensional simulations, we can also predict the variety of the luminosity for
different viewing angles. It is shown that optical emission is greatly
suppressed when NS merger is observed from the equatorial plane while infrared
emission is almost unchanged. These simulations will be useful to facilitate
future EM follow-up observations.

Figure 24:  Expected number of NS merger as a function of distance assuming
the event rate estimated from the O2 run LIGOScientific:2018mvr . The dashed
lines show the range for 90% confidence level.

####  3.3.4 Future Prospects

Our observational and theoretical results indicate that NS mergers synthesize
a wide range of $r$-process elements. With the current estimate of the event
rate and mass ejection, NS merger can explain the total amount of $r$-process
elements in our Galaxy rosswog17 ; hotokezaka18 . However, the event rate is
still largely uncertain. Also, we still do not fully understand that the mass
ejection is universal or not. Varieties in the mass ejection and
nucleosynthesis are expected by numerical simulations. Therefore, it is
important to observe more NS mergers to accurately derive the event rate and
to obtain the general picture of the mass ejection and nucleosynthesis.

The sensitivity of GW detectors are still improving. In the third observing
run (O3), a typical sensitivity corresponds to the range of the NS merger up
to 130Mpc (LIGO Livingstone), 110Mpc (LIGO Hanford), and 45Mpc (Virgo). Thanks
to these high sensitivity, two reliable GW events including at least one NS
have been reported. However, due to the poor positional localization or large
distances, any EM counterpart has not been identified. We refer the reader to
Sasada et al. (2021) sasada21  and Ohgami et al. (2021) ohgami21  for our
observing effort in O3. The situation will be improved when the sensitivity of
Virgo becomes comparable. Furthermore, when KAGRA joins the network with a
comparable sensitivity, positional localization will be greatly improved (down
to a few deg2 for the signals such as GW150914, gaebel17 ).

Figure 24 shows the expected event number per year as a function of the
distance. Once all the sensitivity of GW detectors reach $\sim 200$Mpc, we
expect $>10$ events every year. Kilonova associated with such events is
expected to be faint, $\sim 22$ mag in optical within 2-3 days after the
merger. Therefore, observations with Subaru/HSC will play important roles.
Also, deep near-infrared imaging and spectroscopic observations with large-
aperture telescopes, such as Subaru, Keck, Gemini, and upcoming TAO 6.5m
telescope doi18 , will be important for firm identification of kilonova. Such
multi-wavelength observations will provide us with unique information of
$r$-process nucleosynthesis.

##  4 Supernovae

###  4.1 GW signatures from core-collapse supernovae

####  4.1.1 The CCSN GW signatures

Unlike the GW signals from compact binary coalescence where a template-based
search is best suited (e.g., rasio99 ), gravitational waveforms from CCSNe are
essentially of stochastic nature. This is because the waveforms are affected
by turbulence in the postbounce core, which is governed by the multi-D, non-
linear hydrodynamics. In order to clarify the GW emission mechanisms,
extensive numerical simulations have been done so far in different contexts
(e.g., Dimmelmeier08 ; Scheidegger10 ; CerdaDuran13 ; Ott13 ; Yakunin15 ;
KurodaT14 ; pan18 ; viktoriya18  and Kotake13  for a review). For canonical
supernova progenitors Heger05 , core rotation is generally too slow to affect
the dynamics (e.g., takiwaki16 ; summa18 ). For such progenitors, the GW
emission takes place in the postbounce phase, which is characterized by prompt
convection, neutrino-driven convection, proto-neutron star (PNS) convection,
the standing accretion shock instability (SASI, thierry15 ), and the
$g$(/$f$)-mode oscillation of the PNS (e.g., EMuller97 ; EMuller04 ; Murphy09
; Kotake09 ; BMuller13 ; viktoriya18 ).

The most generic GW emission feature seen in recent self-consistent 3D CCSN
models is the one from the PNS oscillation kuroda16 ; andresen17 ; radice19 ;
alex19 . The characteristic GW frequency increases almost monotonically with
time due to an accumulating accretion to the PNS, which ranges approximately
from $\sim 100$ to $1000$Hz. On the other hand, the typical frequency of the
SASI-induced GW signals is concentrated in the lower frequency range of $\sim
100$ to $260$Hz and persists when the SASI dominates over neutrino-driven
convection kuroda16 ; andresen17 ; andresen18 ; radice19 . The detection of
these distinct GW features is considered as the key to infer which one is
working more predominantly in the preexplosion supernova core, namely
neutrino-driven convection or the SASI andresen17 .

Figure 25: An example gravitational waveform of the $+$ (blue line) and
$\times$ (orange line) mode from a 3D-GR CCSN model of $15M_{\odot}$ star
(model SFHx in kuroda16 ). The time is measured after core bounce (:$t=0$). A
source distance of $10$kpc is assumed. The waveform was extracted via a
standard quadrupole formula with GR corrections KurodaT14 . Since the waveform
is weakly dependent on the source orientation, an unbiased direction is chosen
(e.g., from the north pole ($(\theta,\phi)=(0,0)$). This figure is taken from
kawahara18 .

Shown Fig. 25 is an example waveform taken from a 3D CCSN simulation in
kuroda16 . For the model, the hydrodynamics evolution is self-consistently
followed in full general relativity (GR), starting from the onset of core-
collapse of a $15M_{\odot}$ star WW95 , through core bounce, up to $\sim$
350ms after bounce. As consistent with the outcomes from recent 3D models
(e.g., andresen17 ; radice19 ), the hydrodynamic evolution and the associated
GW waveform is characterized by the prompt convection phase shortly after
bounce ($T_{\rm pb}\lesssim 20$ms with $T_{\rm pb}$ the postbounce time, shown
simply as ”time” in Fig. 25), then the linear (or quiescent) phase, which is
followed by the non-linear phase ($T_{\rm pb}\gtrsim 140$ms) when the vigorous
SASI activity (as well as the growing GW amplitudes) was observed for the
model (see kuroda16  for more details). The dominance of the SASI over
neutrino-driven convection persists over $140\mbox{ms}\lesssim T_{\rm
pb}\lesssim 300$ms, after which neutrino-driven convection dominates over the
SASI.

In order to discuss the detectability of the CCSN GW signals (like of Fig.
25), previous studies have traditionally relied on a GW spectrogram analysis,
aiming to specify the GW feature by the excess power in the time-frequency
domain (by taking the square norm of the short-time Fourier transform).
However, the spectrogram analysis has a trade-off relation between the time
and frequency resolution, leading to a limited time-frequency localization.
Because of this drawback, some features could have been potentially overlooked
in the previous spectrogram analysis in the context of CCSN GWs. Recent
developments in the time-frequency analysis (TFA) have yielded various time-
frequency representation alternatives to the spectrogram. In particular,
quadratic time-frequency representations, such as the Wigner–Ville
distribution and its modified forms, enable us to perform high-resolution TFA
(e.g., cohen95 ; stankovic2013time ; Boashash2015 ).

Recently, the S-method utilizing the Wigner–Ville distribution was applied in
kawahara18 , which successfully extracted several key (CCSN) GW features in
the time-frequency domain (left panel of Fig. 26, where the waveform (Fig. 25)
was used). In comparison with the conventional TFA based on the Fourier
transform (e.g., Fig. 1 in kuroda16 ), the S-method significantly improves the
sharpness of the modes. Furthermore by defining the instantaneous frequency
(IF) of the excess in the time-frequency domain (corresponding to the bright-
color regions in the left panel of Fig. 26, see kawahara18  for more details),
five modes (A B, C, $\mathrm{C}^{\\#}$, and D, the right panel of Fig. 26)
were identified.

Mode A in the right panel of Fig. 26 corresponds to the PNS $g$(/$f$)-mode
oscillation, which is excited in the vicinity of the PNS surface by the
downflows and/or directly originates from the deceleration of infalling
convective plums (BMuller13, ). Modes B and C are quasi-static modes, in the
sense that these two modes barely change with time in the time-frequency
domain. Mode B ($f_{B}\sim$ 130Hz) is originated from SASI, the frequency of
which is twice of the SASI frequency ($f_{\rm SASI}\sim$ 65Hz) because of the
quadrupole nature of the GW emission (see also andresen18 ). Mode D is the
overtone of the mode B ($f_{D}\approx 2f_{B}$), which is most likely to come
from the PNS core oscillation (see Fig. 8 in kawahara18 ). Mode C intersects
mode A at $T_{\rm pb}\sim 180$ms, after which this mode is denoted as
$\mathrm{C}^{\\#}$. Although yet to be investigated elaborately, modes C and
$\mathrm{C}^{\\#}$ are likely to arise in the innermost region ($\sim 10$km)
alex18 . Note that the typical frequency of the SASI-induced GWs is in the
range of $\sim$100 – 260Hz (e.g., modes B and D in the left panel of Fig. 26),
which is in the best sensitivity range of the currently running
interferometers. For a Galactic event, the signal-to-noise ratio (SNR) of the
GW signals predicted in the most recent 3D CCSN models was estimated in the
range of $\sim$ 4 – 10 for the advanced LIGO andresen17 . With the third-
generation detectors on line (e.g., Cosmic Explorer and Einstein Telescope),
these signals would be never missed for the CCSN events throughout the Milky
way.

Figure 26: The left panel is the time-frequency representation of the sample
waveform (Fig. 25), which is extracted by the S-method. The right panel shows
the GW mode identification of the spectrogram based on the instantaneous
frequency identification method via the S-method (see kawahara18  for more
details). In the right panel, the white dashed line represents the theoretical
prediction of the peak GW frequency of the PNS $g$-mode oscillation BMuller13
. Note that the deviation of mode A from the white dashed line (especially in
the later postbounce phase) was previously observed in the literature
andresen17 ; viktoriya18 . These figures are taken from kawahara18 .

####  4.1.2 Circular Polarization of CCSN GWs

Recently, yet another GW feature has been reported, which is circular
polarization of the CCSN GWs. The importance of detecting the GW circular
polarization was first pointed out by Hayama16  in the context of rapidly
rotating core-collapse. More recently, Hayama et al. (2018) Hayama18
presented analysis of the GW circular polarization using results from core-
collapse of a non-rotating $15$$M_{\odot}$ star kuroda16 , the waveform of
which is already shown in Fig. 25.

Figure 27: Shown are the trajectory of the GW circular polarization on the
$h_{+}-h_{\times}$ plane of models SFHx (left panel) and TM1 (right panel) of
3D-GR simulations using different equations of state (see text). The two
characteristic epoch with the right-handed ($140\mbox{ms}\lesssim T_{\rm
pb}\lesssim 200$ms) and left-handed ($220\mbox{ms}\lesssim T_{\rm pb}\lesssim
320$ms) polarization are highlighted with the red and blue color. The gray
line denotes the whole trajectory during the simulation time. The source
distance is assumed at 10kpc. These figures are taken from Hayama18 .

Figure 27 visualizes the time evolution of the GW circular polarization
(hereafter, CP) of models SFHx (left panel) and TM1 (right panel) from
Hayama18 . Note that SFHx SFH  and TM1 denotes the name of the nuclear
equation of state (EOS) employed in the 3D-GR simulation kuroda16 . The key
difference is that the softer EOS (SFHx) makes the PNS radius and the shock
radius at the shock-stall more compact than those of TM1. This leads to much
stronger activity of the SASI for SFHx compared to TM1.

In fact, one can see a clearer polarization signature for SFHx (left panel of
Fig. 27) characterized by the bigger GW amplitude with the right-handed (red
line) and left-handed mode (blue line) than those for TM1 (right panel). Note
in each panel that the two SASI-dominant phases of SFHx are colored by red or
blue (see also bottom panels of Fig. 1 in Hayama18 ). The non-axisymmetric
flows associated with spiral SASI give rise to the circular polarization.
Before $T_{\rm pb}\sim 140$ms, the CP is small because the SASI activity is
still weak. But, after $T_{\rm pb}\sim 140$ms when the (spiral) SASI activity
begins to be vigorous (e.g., kuroda16 ), the right-handed CP firstly emerged,
which is followed by the left-handed CP (the blue line) until $T_{\rm pb}\sim
320$ms, after which neutrino-driven convection dominates over the SASI.

To explore the detectability of the CP, Monte Carlo simulations of the
coherent network analysis hayama15  were performed for the signal
reconstruction, where the network of LIGO Hanford, LIGO Livingston, VIRGO and
KAGRA was considered. They reported the enhancement of the signal-to-noise
ratio of the GW circular polarization (${\rm SNR}_{\rm CP}$) relative to that
of the GW signal itself (${\rm SNR}_{\rm TF}$). Although this is likely
because the Gaussian noise has little component of the circular polarization,
more detailed analysis is required to draw a robust conclusion whether or not
the GW CP could provide a new probe to clarifying the CCSN inner-workings such
as the SASI and the PNS oscillations (see Hayama18  for more detail).

####  4.1.3 Rapid Rotation: GW and Neutrino signals from an exploding
$27M_{\odot}$ star

Rapid rotation in the iron core (with the initial rotation frequency typically
greater than $0.5$rad/s) leads to significant rotational flattening of the
collapsing and bouncing core, which produces a time-dependent quadrupole (or
higher) GW emission (e.g., ott_rev  for a review). For the bounce signals
having a strong and characteristic signature, the iron core must rotate enough
rapidly111Although rapid rotation and the strong magnetic fields in the core
are attracting great attention as the key to solve the dynamics of collapsars
and magnetars, one should keep in mind that recent stellar evolution
calculations predict that such an extreme condition can be realized only in a
special case woos06  ($\lesssim$ 1% of massive star population).. The GW
frequency associated with the rapidly rotating collapse and bounce is in the
range of $\sim 600-1000$Hz Dimmelmeier08 . A current estimate based on a
coherent network analysis using predictions from a set of 3D models shows that
these GW signals could be detectable up to about $\sim 20$kpc for the most
rapidly rotating model (hayama15 , see also powell19 ; gossan16 ).

Figure 28 shows the neutrino and GW signals takikotake  obtained in 3D core-
collapse supernova simulation of a rapidly rotating 27$M_{\odot}$ star that is
exploding due to the growth of the so-called low-T/$|W|$ instability
takiwaki16 . The time modulation seen in the left panel corresponds to the
neutrino light-house effect where the spinning of strong neutrino emission
regions around the rotational axis leads to quasi-periodic modulation in the
neutrino signal. Depending on the observer’s viewing angle, the time
modulation will be clearly detectable in IceCube and the future Hyper-
Kamiokande. The GW emission is also anisotropic where the GW signal is
emitted, as previously identified (see Kotake13  for a review), most strongly
toward the equator at rotating core-collapse and bounce, and the non-
axisymmetric instabilities in the postbounce phase lead to stronger GW
emission toward the spin axis. The right panel in Fig. 28 shows that these GW
signals can be a target of LIGO-class detectors for a Galactic event. The
origin of the postbounce GW emission (e.g., the bar-mode ($m=2$) deformation
of the PNS due to the low T/$|W|$ instability), naturally explains the reason
that the peak GW frequency is about twice of the neutrino modulation
frequency. These results demonstrate that the simultaneous detection of the
rotation-induced neutrino and GW signatures could provide a smoking-gun
signature of a rapidly rotating PNS at the birth (see also walk18 ; walk19 ).

Figure 28: The left panel shows detection rates of $\bar{\nu}_{e}$ at 10kpc
for a rapidly rotating $27M_{\odot}$ model as a function of time after bounce
takikotake . The red and green lines correspond to the event rates (per 1ms
bin) for an observer along the equator and the pole (e.g., parallel to the
rotation axis), respectively. A quasi-periodic modulation (red line)
corresponds to $\sim 120$Hz (red lines) seen for the observer along the
equator. The right panel shows the characteristic GW spectra relative to the
sensitivity curves of advanced LIGO, advanced Virgo (advv ) and KAGRA (aso13
). The peak GW frequency (the vertical red band) is about twice of the
neutrino modulation frequency.

####  4.1.4 Hydrodynamics and GW signals from a BH-forming massive star

Successful observations of GWs from binary black holes (BBHs) by the LIGO-
VIRGO collaborations are now posing a new challenge; what is the origin of the
massive BHs ? One of the most plausible scenarios is a binary stellar
evolution in a low-metallicity environment (see abbott16  for a review). It
has been proposed that two massive stars in the approximate range of $40$ to
$100M_{\odot}$ lead to the formation of a massive helium core after
experiencing the Roche-lobe overflow and common envelope phase (e.g., bel10 ;
Langer12 ; Kinugawa16  for collective references therein). The gravitational
collapse of the massive core ($\sim 30M_{\odot}$) could account for some of
the relevant BH mass ranges (at least in the high-mass end) in the GW events,
although the formation path to the massive core and further to the BH is still
very uncertain due to the complexity of the binary evolution and the fallback
dynamics (e.g., Fryer99 ).

In order to clarify the formation process of the BH, one requires GR neutrino
radiation-hydrodynamics core-collapse simulations of such massive stars in 3D.
Due to the high numerical cost, most of the previous studies with BH formation
have been done assuming spherical symmetry (1D) (e.g., janka16 ; Kotake12  and
collective references therein). In the context of multi-D simulations with
multi-energy neutrino transport, Ref. Pan17  reported 1D and two-dimensional
(2D) core-collapse simulations of a solar-metallicity 40$M_{\odot}$ star using
two-flavor IDSA scheme (idsa ) and a post-Newtonian gravity to include GR
effects. More recently, a BH-forming 3D-GR simulation of a zero-metallicity
40$M_{\odot}$ star with an approximate neutrino transport (FMT) scheme was
reported in chan18 . It is only recently that the first 3D-GR simulation
(using a 70$M_{\odot}$ star Takahashi14 ) with detailed neutrino transport
shibata11  following the dynamics up to BH formation was published Kuroda18 .
In the simulation, the evolution equations of metric, hydrodynamics are solved
based on the Baumgarte-Shibata-Shapiro-Nakamura formalism Shibata95 ;
Baumgarte99  and the evolution of neutrino radiation field based on the multi-
energy M1 scheme Kuroda18 . Casting the helium core of $\sim 31$$M_{\odot}$,
the zero-metallicity 70$M_{\odot}$ star was chosen as one of the possible
progenitor candidates of the BBHs.

Figure 29: The left panel shows time evolution of the maximum rest-mass
density (solid line) and minimum lapse function (dotted line) for a
70$M_{\odot}$ star (denoted as Z70, red line) and a 40$M_{\odot}$ star (as
S40, black line), respectively. The right panel shows angle-averaged profiles
of the (rest-mass) density (black line) and the ratio $r_{\rm s}/R$ (red
lines, see text) for Z70 as a function of the enclosed baryon mass for three
representative time slices near the final simulation time. ”-0.321” and
”-2.359”ms denote the time before the final simulation time ($T_{\rm
fin}=293$ms after bounce). The green line corresponds to $r_{\rm s}/R$ for
S40. These figures are from Kuroda18 .

The left panel of Fig. 29 shows temporal evolution of the maximum density
$\rho_{\rm max}$ (solid lines) and the minimum lapse $\alpha_{\rm min}$
(dotted lines) for the 70$M_{\odot}$ star (denoted as Z70, red line) and the
40$M_{\odot}$ star (as S40, black line), respectively. Note that S40 is shown
for comparing with previous results. The compactness parameter of S40 is much
smaller ($\xi_{2.5}\sim 0.26$) than that of Z70. The maximum density at bounce
($\rho_{\rm max}\sim 4\times 10^{14}$g cm-3) is quite similar between Z70 and
S40. After bounce, the increase of the maximum density of Z70 (red solid line)
is significantly faster than that of S40 (black solid line). In both Z70 and
S40, the minimum lapse (dotted lines) shows a gradual decrease after bounce.
At around $T_{\rm fin}=293$ms after bounce, it shows a drastic drop to
$\alpha_{\rm min}=0.0645$ for Z70, indicating that the PNS core starts to
collapse rapidly toward a BH formation.

The right panel of Fig. 29 explains how $T_{\rm fin}$ is related to the BH
formation time. Shown in the panel is the profiles of (angle-averaged) density
(black lines) and a diagnostics to measure the BH formation as a function of
the enclosed baryon mass $M_{\rm b}$ at some representative snapshots near
$T_{\rm fin}$ for Z70 (red lines) and at $T_{\rm pb}=400$ms for S40 (green
line). As a BH-formation diagnostics, the ratio of $r_{\rm s}/R$ is shown,
where $r_{\rm s}$ and $R$ denotes the Schwarzschild radius and the radial
coordinate, respectively. One can see that the maximum $r_{\rm s}/R$ is $\sim
0.3$ at 2.359ms before $T_{\rm fin}$ (the thin red line labelled by ”-2.359
[ms]”) and rapidly increases with time, approaching to unity (precisely,
$0.932$) at $T_{\rm fin}$ (thickest red line), which was judged as the epoch
of the BH formation in Kuroda18 . It should be noted that for the unambiguous
definition of the BH formation, one requires the implementation of the so-
called apparent horizon finder (e.g., thorn04 ) in numerical relativity
simulation, which is still left as a future work.

At the (fiducial) BH formation time, the mass and the radius is $M_{\rm
b(g),BH}\sim 2.60(2.51)$$M_{\odot}$ and $R_{\rm iso}\sim 4$km, respectively.
By contrast, S40 shows significantly less compact structure (green line) at
the final simulation time ($T_{\rm pb}=400$ms). The BH formation should occur
much later, possibly when the mass shell at $R(M_{\rm b}=2.6M_{\odot})\sim
10^{9}$cm accretes to the stalled shock. Using the same EOS (LS220 LS ), this
expectation is in line with Pan17 , which reported the BH formation at $T_{\rm
pb}\sim 700$ms and with chan18  at $T_{\rm pb}\sim 1$ s.

The time evolution of the (angle-averaged) shock radii $R_{\rm s}$, the gain
radius $R_{\rm g}$, the ratio of the advection timescale to the neutrino-
heating timescale in the gain region $\tau_{\rm adv}/\tau_{\rm heat}$, and the
mass in the gain region $M_{\rm gain}$ are presented in Fig. 2 of Kuroda18 ,
respectively. For model Z70, the shock revival was obtained after $T_{\rm
pb}\,\gtrsim\,260$ms. At this time, the maximum temperature in the core
becomes as high as $T\sim 100$ MeV at a slightly off-center region at $R_{\rm
iso}\sim 10$km (equivalently at $M_{\rm b}\sim 1.0M_{\odot}$). Subsequently
the high temperature region propagates outward in the mass coordinate,
although spatially inward, due to the continuous mass accretion. The maximum
temperature reaches $\sim 170$ MeV at $R_{\rm iso}\sim 1$km ($M_{\rm b}\sim
1.4M_{\odot}$) just before $T_{\rm fin}$. In this second collapse phase to the
forming BH, the high neutrino emission makes the heating timescale shorter
than the competing advection timescale in the gain region. Aided by strong
convection behind the shock, the stalled shock is revived at $T_{\rm
pb}\,\gtrsim\,260$ms ($\tau_{\rm adv}/\tau_{\rm heat}\geq 1$). This also
results in the increase in the gain mass (see the blue line) due to the shock
expansion.

Figure 30: A snapshot of the entropy distribution ($k_{\rm B}$ baryon-1) for
Z70 at $T_{\rm pb}\sim 294$ms just before the BH formation (Kuroda18 ). The
sheet represents the lapse function ($\alpha$) on the $z=0$ plane. This figure
is taken from Kuroda18 .

Figure 30 visualizes a 3D hydrodynamics feature near at the BH formation.
During the first $\sim 160$ms after bounce, the neutrino heating is still weak
and high entropy bubbles are yet to appear. Only after $T_{\rm pb}\gtrsim
230$ms, the formation of high entropy plumes ($s\gtrsim 15k_{\rm B}$ baryon-1)
was seen due to the intense neutrino heating from the hot PNS. At this time,
the mass in the gain region $M_{\rm gain}$ also starts to increase. The
expansion of the (merging) high entropy plumes was observed, leading to the
shock revival. The lapse function shows a steepest drop in the center (see the
cusp in the plane of Fig. 30), which corresponds to the BH formation. By
expanding the shock radius into the spherical harmonics, we find that the
deviation of the shock from spherical symmetry (in the low-modes $\ell=1,2$)
is less than $\sim 2\%$. This clearly indicates that neutrino-driven
convection dominates over the SASI in this case.

Finally, Fig. 31 shows the GW prediction for the Z70 star. The waveform is
extracted along positive $z$-axis via the quadrupole formula. The GW amplitude
(multiplied by the distance to the source) stays at small value $\lesssim
10$cm before the shock expansion occurs ($T_{\rm pb}\sim 260$ms). The strong
GW emission thereafter mainly originates from strong convection motion behind
the shock. Just before the final simulation time, the GW amplitude reaches
$\sim 100$cm, though only for a short duration. Although more quantitative
discussion is apparently needed, a rough estimate shows that the GW signal
could be targeted by third-generation GW detectors such as the Einstein
Telescope (ET) and the Cosmic Explorer (CE) (ET ; CE ) for the source distance
less than $\sim$ Mpc. Regarding the neutrino signals, both electron type
(anti-)neutrinos show decreasing and plateau trend for $T_{\rm pb}\gtrsim
260$ms, whereas heavy-lepton neutrinos show a rapid increase both in the
luminosity and energy. These features are consistent with those previously
identified in 1D full-GR simulations with Boltzmann neutrino transport
(matthias04 ), which is due to rapid contraction of the PNS to the forming BH
(see also, Sumiyoshi07 ; tobias09  ). The detection of the short-live ($\sim
300$ms after bounce) neutrino signals are basically limited to Galactic events
(see mirizzi  for a review). However, further study would be needed to clarify
the contribution of these BH forming massive stars to the prediction of
diffuse neutrino supernova background (e.g., Lunardini09 ; horiuchi18 ).

Figure 31: A gravitational waveform for the BH-forming core-collapse of a
70$M_{\odot}$ star. Note that $h_{+}$ and $D$ denote the GW amplitude of the
$+$ polarization and the distance to the source, respectively.

####  4.1.5 Short Summary and Future Prospects

To summarize, each phase of a CCSN has a range of characteristic GW signatures
and the GW polarization that can provide diagnostic constraints on the
evolution and physical parameters of a CCSN and on the explosion
hydrodynamics. The GW signatures described in this section commonly appear in
the frequency range of $100$ \- $1000$Hz in recent simulations of the CCSN.
Among them, the PNS oscillatory modes are currently recognized as the model-
independent GW signature. The characteristic frequency ($f_{p}\propto
M/R^{2}$) of the fundamental ($f/g$) modes is predominantly dependent on the
PNS mass ($M$) and the radius ($R$) bernhard13 . In order to break the
degeneracy, the detection of other eigen-modes ($p,w$ modes) of the PNS
oscillations is mandatory. In fact, the GW asteroseismology of the PNS has
just started forne19 ; morozova18 ; sotani17 , the outcomes of which should
reveal the requirement of the next-generation detectors to detect the whole of
the eigen-modes (presumably extending up to several kHz) for the future CCSN
event.

Regarding the SASI-induced GWs for a Galactic supernova source, they are
predicted to be detectable by the advanced LIGO with the singal-to-noise ratio
(SNR) in the range of $\sim 4$–$10$ by the most recent 3D CCSN models
andresen17 . With the third-generation detectors (e.g. Cosmic Explorer and
Einstein Telescope) online, these signals would never be missed for CCSN
events throughout the Milky Way. A current estimate of SNR for GWs from
rapidly rotating collapse and bounce based on predictions from a set of 3D
models using a coherent network analysis shows that these GW signals could be
detectable up to about $\sim 20$kpc for the most rapidly rotating model
(hayama15 ; see also Refs. powell19 ; gossan16 ). See a recent review by Ref.
ernazar  for a more comprehensive review regarding the detectability of the
CCSN GW signals by the future detectors.

After decades of progress, the CCSN theory is finally converging, but not
without detailed numerical effort and much theoretical sweat. But at last not
in the least, we have to draw a caution that the current numerical results and
the associated GW predictions that we have reported in this Section should
depend on the next-generation calculations by which more sophistication can be
made not only in the neutrino transport schemes (such as Boltzmann transport
sumi12 ; harada19 ; nagakura19 ), but also in the treatment of GR with the BH
formation. Therefore we provide here only a snapshot of the moving (long-run)
documentary film that records our endeavours for making our dream of ”GW
astronomy of the next Galactic CCSN” come true.

###  4.2 Understanding Supernovae via Neutrino Emissions

One of the primary sources of near-field GWs, core collapse supernova
explosions are among the most dramatic and important events to take place in
the Universe. Agents of great destruction – anything within tens of light-
years is utterly annihilated by the dying star – they are nevertheless
responsible for the existence of life itself, for they (and the collisions of
neutron stars they produce) are the sole source of all elements heavier than
helium. Therefore, understanding these complex explosive events is necessary
if we are to understand why we are here, and why the Universe looks the way it
does today.

Supernova neutrinos, famously observed from SN1987A, provide a unique and
vital probe into the inner dynamics of these events. Released together with
GWs during the initial stellar collapse, neutrinos and GWs are both certain to
travel through any obscuring dust or gas and remain undiminished upon their
arrival at the Earth. Neutrinos also carry information regarding the end state
of the star: for explosions within our galaxy, collapses into neutron stars or
black holes, the eventual sources of far-field GWs, can be differentiated via
observations of neutrino emissions.

Our goal is to make theory and experiment work together so that we will be
ready to make the best possible observations of the next galactic explosion
and maximize our extraction of information on the explosion mechanism,
progenitor, and nuclear physics ingredients.

In order to both predict and make sense of the neutrino signals from the next
explosion in the Milky Way galaxy, the world’s most advanced supernova
numerical simulations are needed. Using nuclear data such as equation of state
and neutrino reactions, and by solving the 6D Boltzmann equation in 2D/3D, we
can provide detailed predictions of the time distribution and energy spectra
of supernova neutrinos.

To verify these complex computational models, we will enrich the famous Super-
Kamiokande (SK) detector with gadolinium [Gd] salt. Doing so will turn it into
the world’s most advanced supernova neutrino detector, capable of real-time
tagging and identification of individual supernova neutrino interactions with
nanosecond-scale time resolution. This will make the diffuse supernova
neutrino flux from all past supernova explosions visible for the first time.
Starting around the middle of 2020, a gadolinium-loaded SK will collect a
steady stream of supernova neutrino data, the first such new data in over 30
years.

The new experimental measurements made possible by gadolinium will then be fed
into the theoretical models, testing their predictions using real supernova
neutrino data and allowing us to better prepare for the next nearby supernova
explosion. The models will be refined as needed and as indicated by the past
supernova data. The gadolinium-loading of SK will also greatly improve the
detector’s response to a nearby supernova. Therefore, both experimentally and
theoretically, we will be well-prepared and ready for it.

####  4.2.1 SK-Gd

Water Cherenkov detectors such as Kamiokande, IMB, and SK, have been used for
decades as effective detectors for neutrino interactions and nucleon decay
searches. While many important measurements have been made with these very
large detectors, a major drawback has been their inability to efficiently
detect the presence of thermal neutrons.

If a water Cherenkov detector could be improved to observe the neutrons
produced by the inverse beta process ($\bar{\nu}_{e}p\rightarrow e^{+}n$),