Gravitational Wave Physics and Astronomy in the nascent era

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.


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 [1]. 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 30M . The mass of the BH candidates found so far by electromagnetic observations is estimated to be at most about 20M [2]. 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 30M . The observation of GWs has changed our understanding about the mass distribution of BHs in the universe [3,4]. Since then, the number of BBH observations officially announced by the LIGO/Virgo collaborations has been increasing [5]. 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 [6]. 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 [7] and their formation in the primordial universe [8]. Many variants of these scenarios can be 5/89 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.

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 [9]. 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 [10], longterm follow-up observations by radio waves and X-rays, detection of apparent superluminal motion [11,12]. 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 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 [13,14]. This interpretation gained a strong support because apparent superluminal motion was observed by radio waves [11,12]. 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 [9]. 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 [15][16][17][18], although BNS mergers and SGRBs are considered to be sources of these particles [19][20][21][22]. We do not further discuss these multi-messenger signals in this review. point that multi-dimensional (multi-D) hydrodynamics instabilities play a key role in the neutrino mechanism [24], the most favoured scenario to trigger explosions. In fact, selfconsistent simulations in three spatial dimensions (3D) now report shock revival of the stalled bounce shock into explosion by the multi-D neutrino mechanism (see [25][26][27][28] 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., [29,30] 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., [31]), but also for unraveling the hydrodynamics evolution from the onset of core-collapse toward the forming compact objects (neutron stars and black holes, see [32] for a review). Current estimates of CCSN rates in our Milky way predict one event every ∼ 40 ± 10 year [33]. 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 [34].

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 [35]. 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 ≈ 10 −15 [36], which excludes various possible extensions of GR [37,38,[40][41][42][43].
GWs from binaries can be used as the standard siren [44], providing an alternative method to measure the cosmic expansion history. GW170817 already gives us an estimate of the Hubble constant [45,46], but we do not further discuss the aspects of testing cosmological models using GWs in this review.

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. 7/89 We discuss the test of gravity and implication to the dark matter physics. Sec. 6 is devoted to a short conclusion. 8/89 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 [5] 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-tonoise ratios (SNRs) are not so high. There are results reported by alternative searches. One of them reports other candidates of GW events [47]. 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 M, effective spin χ eff , final BH mass M fi , and its angular momentum parameter a 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. [5]. The obtained mass and spin distributions are summarized in Fig. 1. The chirp mass is defined by with M = m 1 + m 2 and µ = m 1 m 2 /M , where m 1 and m 2 are the component masses. The effective spin parameter is defined by where S 1 and S 2 are spin angular momentum vectors of the respective bodies andL is the unit vector pointing at the direction of the orbital angular momentum of the binary. The reported event rate by the LIGO/Virgo Collaborations is ranging from 9.7 to 101 Gpc −3 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. [2] 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 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 20% difference in the estimation. 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.

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.  Fig. 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 [48]), 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).

A.
Binary star formation at low metallicities 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 0.1 Z , which are required for yielding massive ( a few ×10 M ) 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 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 [49][50][51]. 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 ∼ Myr is required to cover the whole evolution of 11/89 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 [48]. In this particular case, we follow the long-term (∼ 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 [52] and the binary frequency is even higher for high-mass stars than for low-mass ones [53]. 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 [54]. In rapidly rotating and gravitationally unstable cores, fragmentation occurs more easily [55,56] 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 [54,57] and the outcome is a wide binary system with a separation typically greater than 100 AU. In contrast, observationally high-mass stars tend to be in close-binary systems with the separation 1-10 AU [52]. 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 [58]. 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 ∼ 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 [59]. 12 Fig. 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 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 H = 1.4 × 10 4 cm −3 and T = 200K) 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, ∼ 17 times longer than the previous calculation [60]. 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 th = 10 15 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 [60]. 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 th . This trend has already been known previously [60], 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.
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 ∼ a few 10 AU, which is too wide to end up with BBHs that coalesce within a We assume ρ ad is 2 × 10 −5 (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 [60]. Results in the literature [61][62][63][64][65][66][67][68][69][70][71] are also shown.  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.
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 ∼ 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 ∼ 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 [72], 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.

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 [7,73]. Those works suggested a local merger rate density of ∼ 5 Gpc −3 yr −1 [73,74]. They also predicted BBHs formed in GCs have several distinct features from those formed under isolated environments. The GC-originated BBHs have more isotropic spinorbit misalignment angles [75]. Some of those BBHs may contain massive (∼ 100M ) BHs locating in the so-called "the pair-instability mass gap" [76], or may have finite eccentricities at a GW frequency of ∼ 10 Hz [77,78]. Such Monte-Carlo simulations, however, need careful calibrations of some parameters from comparison with the N -body results [79,80]. On the other hand, direct N -body simulations with a realistic number of particles (i.e., N ∼ 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 [81], although the situation may change in near future by adoption of tree-direct hybrid methods [82,83] and novel schemes for binary integration [84].
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 [85]. As an initial condition, we adopted a Plummer model [86] with a total mass of 2.5 × 10 3 M with 0.1 solar metallicity. We calculated 360 realizations per a given set of the parameters and obtained ∼ 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 [87]. 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 [88]. This value corresponds to 20-50 % of those formed in GCs citeRo-driguez16. This means BBHs formed in OCs contribute to GW events much more than previously expected.
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 [89] and a cosmic star formation rate density [90], 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 ∼ 70 yr −1 Gpc −3 [91], comparable to that estimated from LIGO/Virgo O2 and O3, 64.0 +73.5 −33.0 yr −1 Gpc −3 [92]. Our simulations also give a prediction on BH-MS binary distribution observable with Gaia [93]. In the next data release of Gaia (Gaia DR3), binaries will be included. Using the results of our simulations, we estimated that ∼ 10 BH-MS binaries would be detectable with Gaia [94], even if we assume relatively stringent observational constraints [95]. In contrast to BH-MS binary formed in galactic fields, those formed in OCs would have smaller mass of the MS companion (< 5M ), 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., [96]). Coupled with star cluster simulations, BPS can also predict the merger rate of BBHs formed in dense stellar clusters (e.g., [74] 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 ( [98,99], respectively). The BSE code is based on several BPS calculations (e.g., [96,97]), and the NBODY6 code and its derivatives are widely used for obtaining BBH merger rates in dense star clusters (e.g., [88,[100][101][102][103][104][105]). Here, we describe these fitting formulae briefly (see [106], 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 [107][108][109][110]. We follow the time evolution of stars with 8 ≤ M/M ≤ 160 for −8 ≤ log(Z/Z ) ≤ −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 17/89  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 ∼ 20M ) 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.
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 ) = −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 18/89 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.

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 [97], 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.

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 [111][112][113][114][115][116][117][118][119][120][121][122][123][124][125][126]. The tidal deformability λ is defined through Q ij = −λ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 λ = 2k 2 R 5 . Therefore, the behavior of λ 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. λ 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 [127,149].
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 [117,124]. Furthermore, the waveforms employing the effective-one-body (EOB) formalism are derived to incorporate higher-order PN effects [111][112][113][114][115][128][129][130][131][132]. 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 [116,119,125,126]. 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 [128,[133][134][135][136][137][138][139][140][141][142]. The recent progress in numerical techniques 19/89 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 [140,143,144], 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, η, in the range of 0.247-0.25, the binary with the total mass of ≈ 2.7 and 2.5 M , and 5 different EOSs which cover the binary tidal deformability in the range of ≈ 300-1900. Here, the binary tidal deformability,Λ, is defined bỹ where Λ i = λ i /m 5 i (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 [140,143,144].
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 (SEOB-NRv2T [129,132]) waveforms can be accurate even up to ≈ 3 ms before the onset of merger [140] (see Fig. 9 for an example). However, the phase difference between the tidal-EOB waveforms and the NR results is still larger than ≈ 1 rad after two NSs come into contact for the case that the NS radii are larger than ≈ 13 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 [145] for convenience in data analysis. We decompose a frequency-domain gravitational waveform, h (f ), into the frequency-domain amplitude, A (f ), and phase, Ψ (f ), (with an ambiguity in the origin of the phase) byh and we define the corrections of the NS tidal deformation to GW amplitude and phase by and  Fig. 9 Comparison between NR and SEOBNRv2T waveforms for the models with the total mass of ≈ 2.5 M 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. respectively. Here, A BBH (f ) and Ψ BBH (f ) are the GW amplitude and phase of a binary black hole with the same mass as the BNS, respectively (hereafter referred to as the pointparticle parts). In this work, we employ the SEOBNRv2 waveforms [146] 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. [115] and [139] for the tidal-part amplitude and phase correction: for the phase correction and for the amplitude correction. Here, a, p, b, and q are the free parameters of the models.  Fig. 10 (Left) Difference in the tidal-part phase between the hybrid waveforms and our analytical model for the various binary configurations (see [140] and [143] 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.
Since our NR waveforms only contain the waveforms of which frequency is larger than ≈ 400 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. [129,132] (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. [139,143].
To focus on the inspiral-phase waveform and to avoid the contamination from the postmerger waveforms, which would have large uncertainties, we restrict the GW frequency range in 10-1000 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 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,F , defined bȳ where (·|·) and || · || are defined by  Fig. 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Λ = 1000 (dot-dashed, blue), 400 (dashed, blue), and 100 (dotted, blue), the NRTidal model (solid, red [135]) and the NRTidalv2 model (solid, cyan [136]) are shown in the plot. Note that latter two models are independent ofΛ when normalized by the leading tidal phase.
Here S 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 [147] for it. We found that the mismatch between our phenomenological waveform model and the hybrid waveforms is always within 10 −5 [143,144]. Gravitational waveform models for BNSs based on NR waveforms are also derived in [135] and [136] 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 [135] and [136] normalized by the leading-order PN term (see also [148] and [136]). 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 [136]).

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 [152]. The point-particle [153,154], and spin parts [155][156][157] 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 [115]. Since the tidal effects 23/89 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. [135,136] (NRTidal) and Kawaguchi et al. [143] (KyotoTidal) independently construct new waveform models in which the tidal part is improved by the calibration using NR waveforms. In KyotoTidal, nonlinear effects ofΛ 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 [146,158], 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 [159][160][161][162][163][164][165][166][167][168][169].
Narikawa et al. separately analyze the GW data of a binary-neutron-star merger, GW170817, by each of the Advanced LIGO twin detectors [170]. 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Λ for GW170817 between NR calibrated waveform models (the KyotoTidal and NRTidalv2 models) as shown in Fig. 13 [148]. Here, NRTidalv2 model is an upgrade of the NRTidal model [136]. The order of peak values ofΛ for the different waveform models in Fig. 13 can be understood by the difference in the magnitude of the tidal phase at around Λ ≈ 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 [171]. The total mass of the system was about 3.4M , 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Λ > 1100 was disfavored, which is consistent with the several analyses done for GW170817. For the low-spin prior (component spin is enforced as χ 1,2 ∼ U [−0.05, 0.05]), the 90% upper limit of GW190425 becomes much smaller than that of GW170817,Λ ≤ 600. Since the total mass of the system was large, the possibility of BH-NS system is discussed in [171].

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; [172]) near the detection limit [173]. The emission with a short hard pulse and a soft tail 24  lasted for ∼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 [174].
Amazingly, the binary coalescence detected by the Laser Interferometer Gravitationalwave Observatory (LIGO) occurred ∼1.7s before the Fermi/GBM trigger [9] and the sky position determined by the LIGO GW observation is consistent with the Fermi/GBM one. 25  The gamma-ray spectrum for the initial hard pulse is well fitted by a power-law function with a photon index of α = -1.6±0.4 and an exponential cutoff at E peak = 185 ± 62keV. For the soft tail emission, a blackbody with k B T = 10.3 ± 1.5keV well represents the spectrum. The fluence over the total duration is ∼3 × 10 −7 erg cm −2 , which falls in the ∼50th percentile of the fluence distribution obtained by Fermi/GBM [173]. Assuming a distance to the host galaxy NGC 4993 as 43Mpc, the isotropic equivalent energy and peak luminosity in the 1keV -10MeV band are ∼3 × 10 46 erg and ∼2 × 10 47 erg 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; [175]) 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 [17]. Fermi/LAT resumed scientific operation ∼10 3 s after the merger time and put a flux upper limit of 4.5 × 10 −10 erg 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 [176][177][178][179]. In particular, the nearby short GRB 150101B detected by Fermi/GBM (z = 0.134 [180]), which could accompany the kilonova and the off-axis jet emission suggested by the optical and X-ray observations [177], has a hard spike and a soft tail. Thus, the two-component sigunature could be a common feature for short GRBs [176]. The isotropic equivalent energy and luminosity in the 1keV -10MeV band are ∼2 × 10 49 erg and ∼4 × 10 50 erg 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. 26/89 MAXI.
Monitor of All-sky X-ray Image (MAXI; [181]) 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; [182]) GSC consists of 12 position sensitive proportional counters. Six of them covers 3 degree × 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 [183].
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 [184]. Fig. 15 shows the luminosity and upper limits for early (< 10 days) X-ray afterglows [185,186] 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 ∼ 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.

Swift.
The Neil Gehrels Swift Observatory (Swift; [187]) consists of three scientific instruments. The Burst Alert Telescope (BAT; [188]) is a wide field of view telescope monitoring about 1/6 of the sky in the 15-150keV band. The X-ray Telescope (XRT; [189]) 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; [190]) 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 [185,191]. 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 (∼ 1 • 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 27 [184]. 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 [185,186]. The dotted, dashed, and dash-dotted lines are models for luminosity L as a function of the time t, L ∝ t α , and α 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 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 [185].

CALET.
The Calorimetric Electron Telescope (CALET; [192,193]) is the scientific instrument attached on the International Space Station (ISS). The CALET Gamma-ray Burst Monitor (CGBM; [194]) 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 ∼60 • from the boresight) and from 40keV to 20MeV (Soft Gamma-ray Monitor; SGM, FoV of ∼110 • 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 × 10 −7 erg cm −2 s −1 in the 10-1000keV band assuming no shielding by the structure of ISS [36].
In addition, the main CALET instrument, the high-energy Calorimeter (CAL), observes gamma-rays from ∼1GeV up to 10TeV with a field of view of nearly 2sr. The field of view  Fig. 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 × 10 20 cm −2 and the photon index of 1.59 [198].
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 × 10 −10 erg cm −2 s −1 for gamma-rays above 1GeV and 4.0 × 10 −10 erg cm −2 s −1 above 10GeV [195].

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 × 10 38 erg s −1 at a distance of 40Mpc [196].
According to the year-long monitoring observations by Chandra, the X-ray flux continues to increase up to ∼160 days after the merger time by a power-law index of 0.9 [197,198]. After the peak, the flux shows the decline by a power-law index of ∼ −2 (Fig. 16). There is a hint of an X-ray flaring feature around 155 days after the merger time [199]. 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.

Theoretical
Interpretation. The GW event GW170817 [9] and the associated short GRB 170817A [36,173,174] 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 29 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.

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 [11,12] 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 [198,200,201], excluding a failed-jet scenario that the jet is choked by the ejecta from the neutron star merger [202][203][204].
The jet should be strong enough for successfully breaking out the merger ejecta in this event [205,206]. The breakout time since the jet launch (t = t 0 ) is approximately estimated as where L 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. [207]). 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 ∼ 1.7s after the merger, the jet luminosity should be large enough L iso,0 3 × 10 49 erg s −1 and the delay time should be short enough t 0 − t m 1.3 s [207]. The observed isotropic energy of the gamma-ray emission is very small, several orders of magnitude smaller than typical [36,173,174]. Therefore the jet should be off-axis for a reasonable radiative efficiency [13,14]. An off-axis observer receives photons emitted outside 30/89 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 Xray, are not consistent with a top-hat jet [208], but indicate a structured jet [12,198,201,[209][210][211][212][213][214]. These are the first confirmation of previous suspicion that realistic GRB jets are structured [215,216].
However, there remains confusion about the jet structure as different authors give different structures, such as Gaussian [198,201,211] and power-law profiles [12,209] (see Fig. 1 of Ref. [13]). 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 [217]. 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.
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 [217].
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 31/89 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 [14], 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 [218,219] (because the off-center is different from the usually-observed jet core) and causes the compactness problem [220,221] (because the off-centre jet is much less energetic and much closer to the line of sight than the jet core).
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 [222]. Future observations of off-axis events will reveal the jet structure, where roughly ∼ 10% events are expected to be brighter at smaller viewing angles than GRB 170817A.

Optical/Infrared Counterparts and Heavy Element Nucleosynthesis
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 [223][224][225][226][227][228][229][230][231]. However, there has been no observational evidence of r-process by NS mergers. 32 Fig. 20 Optical, infrared, and radio telescopes in the J-GEM network.
When NS mergers occur, a small fraction of mass of NSs is ejected into interstellar space [228,230,[232][233][234][235]. 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 [236][237][238][239][240][241]. 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.   [265]. Right: three-color image of the counterpart of GW170817 [267]. [242]. 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, [1]). Within 4 days from the GW detection (2 days from the initial alert), about 24deg 2 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 [242].
More intensive observing campaign has been performed for the second GW source GW151226 (BBH, [243]). For this event, 238 nearby galaxies were observed by galaxy targeted surveys [244], which is great improvement compared with the observations of GW150914. Also, Kiso Schmidt telescope and MOA-II telescopes covered 778deg 2 and 145deg 2 area, respectively [244]. Furthermore, Subaru/HSC covered 63.5deg 2 area with unprecedented sensitivity (limiting magnitudes of 24.6 and 23.8 for the i and z bands, respectively, [245]). The total area covered by the wide-field surveys was 986.5deg 2 , which corresponds to ∼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.8deg 2 ) is the largest among 8-10m class telescopes. Table 3 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. 3.3.2. GW170817. The first GW detection from a NS merger has been achieved for GW170817 [9]. The detection triggered EM follow-up observations covering entire wavelength range [10]. Thanks to the observations with three GW detectors, the position of the GW source is localized to ∼ 30deg 2 . In optical/infrared wavelengths, a counterpart AT2017gfo was identified by several groups .
We have performed wide-field survey using Subaru/HSC, covering 23.6deg 2 [265]. 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 [267]. 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 ∼ 0.03M of ejecta including lanthanide elements, which have a high opacity [269]. Also, the blue optical component, revealed by observations with Subaru/HSC, B&C/Tripol5, and 35/89  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.

36/89
3.3.3. Advances of Theoretical Models. The observed properties of AT2017gfo are consistent with the expectation from numerical relativity simulations [270]. 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 [271]. We find that post-merger ejecta (> 10 −2 M ) dominate over the dynamical ejecta (< 10 −2 M ). 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 [272]. 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 [273],   [274] 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.

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 [3,4]. 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, 37 [277] 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 deg 2 for the signals such as GW150914, [275]). Figure 24 shows the expected event number per year as a function of the distance. Once all the sensitivity of GW detectors reach ∼ 200Mpc, we expect > 10 events every year. Kilonova associated with such events is expected to be faint, ∼ 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 [276], will be important for firm identification of kilonova. Such multi-wavelength observations will provide us with unique information of r-process nucleosynthesis.
The most generic GW emission feature seen in recent self-consistent 3D CCSN models is the one from the PNS oscillation [298][299][300][301]. The characteristic GW frequency increases almost monotonically with time due to an accumulating accretion to the PNS, which ranges approximately from ∼ 100 to 1000Hz. On the other hand, the typical frequency of the SASIinduced GW signals is concentrated in the lower frequency range of ∼ 100 to 260Hz and persists when the SASI dominates over neutrino-driven convection [299][300][301][302]. 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 [299].  [300]). The time is measured after core bounce (:t = 0). A source distance of 10kpc is assumed. The waveform was extracted via a standard quadrupole formula with GR corrections [285]. Since the waveform is weakly dependent on the source orientation, an unbiased direction is chosen (e.g., from the north pole ((θ, φ) = (0, 0)). This figure is taken from [304]. Shown Fig. 25 is an example waveform taken from a 3D CCSN simulation in [300]. 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 star [303], through core bounce, up to ∼ 350ms after bounce. As consistent with the outcomes from recent 3D models (e.g., [299,301]), the hydrodynamic evolution and the associated GW waveform is characterized by the prompt convection phase shortly after bounce (T pb 20ms with T 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 pb 140ms) when the vigorous SASI activity (as well as the growing GW amplitudes) was observed for the model (see [300] for more details). The dominance of the SASI over neutrino-driven convection persists over 140ms T pb 300ms, 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 39/89 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., [305][306][307]).
Recently, the S-method utilizing the Wigner-Ville distribution was applied in [304], 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 [300]), 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 [304] for more details), five modes (A B, C, 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 [297]. 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 ∼ 130Hz) is originated from SASI, the frequency of which is twice of the SASI frequency (f SASI ∼ 65Hz) because of the quadrupole nature of the GW emission (see also [302]). Mode D is the overtone of the mode B (f D ≈ 2f B ), which is most likely to come from the PNS core oscillation (see Fig. 8 in [304]). Mode C intersects mode A at T pb ∼ 180ms, after which this mode is denoted as C # . Although yet to be investigated elaborately, modes C and C # are likely to arise in the innermost region (∼ 10km) [308]. Note that the typical frequency of the SASI-induced GWs is in the range of ∼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 ∼ 4 -10 for the advanced LIGO [299]. 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.

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 [309] in the context of rapidly rotating core-collapse. More recently,   [310] presented analysis of the GW circular polarization using results from core-collapse of a non-rotating 15M star [300], the waveform of which is already shown in Fig. 25. Figure 27 visualizes the time evolution of the GW circular polarization (hereafter, CP) of models SFHx (left panel) and TM1 (right panel) from [310]. Note that SFHx [311] and TM1 denotes the name of the nuclear equation of state (EOS) employed in the 3D-GR simulation [300]. 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 40

Fig. 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 [304] 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 [297]. Note that the deviation of mode A from the white dashed line (especially in the later postbounce phase) was previously observed in the literature [287,299]. These figures are taken from [304]. The source distance is assumed at 10kpc. These figures are taken from [310].
(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 [310]). The non-axisymmetric flows associated with spiral SASI give rise to the circular polarization. Before T pb ∼ 140ms, the CP is small because the SASI activity is still weak. But, after 41/89 T pb ∼ 140ms when the (spiral) SASI activity begins to be vigorous (e.g., [300]), the righthanded CP firstly emerged, which is followed by the left-handed CP (the blue line) until T pb ∼ 320ms, after which neutrino-driven convection dominates over the SASI.
To explore the detectability of the CP, Monte Carlo simulations of the coherent network analysis [312] 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 (SNR CP ) relative to that of the GW signal itself (SNR 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 [310] for more detail).

Rapid
Rotation: GW and Neutrino signals from an exploding 27M star. Rapid rotation in the iron core (with the initial rotation frequency typically greater than 0.5rad/s) leads to significant rotational flattening of the collapsing and bouncing core, which produces a time-dependent quadrupole (or higher) GW emission (e.g., [313] for a review). For the bounce signals having a strong and characteristic signature, the iron core must rotate enough rapidly 1 . The GW frequency associated with the rapidly rotating collapse and bounce is in the range of ∼ 600 − 1000Hz [280]. 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 ∼ 20kpc for the most rapidly rotating model ( [312], see also [315,316]). Figure 28 shows the neutrino and GW signals [317] obtained in 3D core-collapse supernova simulation of a rapidly rotating 27M star that is exploding due to the growth of the socalled low-T/|W | instability [290]. 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 [288] 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 [320,321]).   Fig. 28 The left panel shows detection rates ofν e at 10kpc for a rapidly rotating 27M model as a function of time after bounce [317]. 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 ∼ 120Hz (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 ( [318]) and KAGRA ( [319]). The peak GW frequency (the vertical red band) is about twice of the neutrino modulation frequency.
scenarios is a binary stellar evolution in a low-metallicity environment (see [6] for a review). It has been proposed that two massive stars in the approximate range of 40 to 100M lead to the formation of a massive helium core after experiencing the Roche-lobe overflow and common envelope phase (e.g., [322][323][324] for collective references therein). The gravitational collapse of the massive core (∼ 30M ) 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., [325]).
In order to clarify the formation process of the BH, one requires GR neutrino radiationhydrodynamics 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., [25,27] and collective references therein). In the context of multi-D simulations with multi-energy neutrino transport, Ref. [326] reported 1D and twodimensional (2D) core-collapse simulations of a solar-metallicity 40M star using two-flavor IDSA scheme ( [327]) and a post-Newtonian gravity to include GR effects. More recently, a BH-forming 3D-GR simulation of a zero-metallicity 40M star with an approximate neutrino transport (FMT) scheme was reported in [328]. It is only recently that the first 3D-GR simulation (using a 70M star [329]) with detailed neutrino transport [330] following the dynamics up to BH formation was published [331]. In the simulation, the evolution equations of metric, hydrodynamics are solved based on the Baumgarte-Shibata-Shapiro-Nakamura formalism [332,333] and the evolution of neutrino radiation field based on the multi-energy 43/89 M1 scheme [331]. Casting the helium core of ∼ 31M , the zero-metallicity 70M star was chosen as one of the possible progenitor candidates of the BBHs.
The left panel of Fig. 29 shows temporal evolution of the maximum density ρ max (solid lines) and the minimum lapse α min (dotted lines) for the 70M star (denoted as Z70, red line) and the 40M star (as S40, black line), respectively. Note that S40 is shown for comparing with previous results. The compactness parameter of S40 is much smaller (ξ 2.5 ∼ 0.26) than that of Z70. The maximum density at bounce (ρ max ∼ 4 × 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 fin = 293ms after bounce, it shows a drastic drop to α 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 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 b at some representative snapshots near T fin for Z70 (red lines) and at T pb = 400ms for S40 (green line). As a BHformation diagnostics, the ratio of r s /R is shown, where r s and R denotes the Schwarzschild radius and the radial coordinate, respectively. One can see that the maximum r s /R is ∼ 0.3 at 2.359ms before T fin (the thin red line labelled by "-2.359 [ms]") and rapidly increases with time, approaching to unity (precisely, 0.932) at T fin (thickest red line), which was judged as the epoch of the BH formation in [331]. 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., [334]) 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 b(g),BH ∼ 2.60(2.51)M and R iso ∼ 4km, respectively. By contrast, S40 shows significantly less compact structure 44/89 Fig. 30 A snapshot of the entropy distribution (k B baryon −1 ) for Z70 at T pb ∼ 294ms just before the BH formation ( [331]). The sheet represents the lapse function (α) on the z = 0 plane. This figure is taken from [331].
(green line) at the final simulation time (T pb = 400ms). The BH formation should occur much later, possibly when the mass shell at R(M b = 2.6M ) ∼ 10 9 cm accretes to the stalled shock. Using the same EOS (LS220 [335]), this expectation is in line with [326], which reported the BH formation at T pb ∼ 700ms and with [328] at T pb ∼ 1 s.
The time evolution of the (angle-averaged) shock radii R s , the gain radius R g , the ratio of the advection timescale to the neutrino-heating timescale in the gain region τ adv /τ heat , and the mass in the gain region M gain are presented in Fig. 2 of [331], respectively. For model Z70, the shock revival was obtained after T pb 260ms. At this time, the maximum temperature in the core becomes as high as T ∼ 100 MeV at a slightly off-center region at R iso ∼ 10km (equivalently at M b ∼ 1.0M ). Subsequently the high temperature region propagates outward in the mass coordinate, although spatially inward, due to the continuous mass accretion. The maximum temperature reaches ∼ 170 MeV at R iso ∼ 1km (M b ∼ 1.4M ) just before T 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 pb 260ms (τ adv /τ heat ≥ 1). This also results in the increase in the gain mass (see the blue line) due to the shock expansion. Figure 30 visualizes a 3D hydrodynamics feature near at the BH formation. During the first ∼ 160ms after bounce, the neutrino heating is still weak and high entropy bubbles are yet to appear. Only after T pb 230ms, the formation of high entropy plumes (s 15k 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 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 45/89 Fig. 31 A gravitational waveform for the BH-forming core-collapse of a 70M star. Note that h + and D denote the GW amplitude of the + polarization and the distance to the source, respectively. 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 = 1, 2) is less than ∼ 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 10cm before the shock expansion occurs (T pb ∼ 260ms). The strong GW emission thereafter mainly originates from strong convection motion behind the shock. Just before the final simulation time, the GW amplitude reaches ∼ 100cm, 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) ( [433,434]) for the source distance less than ∼ Mpc. Regarding the neutrino signals, both electron type (anti-)neutrinos show decreasing and plateau trend for T pb 260ms, whereas heavylepton 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 ( [336]), which is due to rapid contraction of the PNS to the forming BH (see also, [337,338] ). The detection of the short-live (∼ 300ms after bounce) neutrino signals are basically limited to Galactic events (see [32] 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., [339,340]).

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 -1000Hz 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 ∝ M/R 2 ) of the fundamental (f /g) modes is predominantly dependent on the PNS mass (M ) and the radius (R) [341]. 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 46/89 asteroseismology of the PNS has just started [342][343][344], 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 ∼ 4-10 by the most recent 3D CCSN models [299]. 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 ∼ 20kpc for the most rapidly rotating model ( [312]; see also Refs. [315,316]). See a recent review by Ref. [345] 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 [346][347][348]), 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.

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 starthey 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 47/89 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.

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 (ν e p → e + n), then backgrounds would be greatly reduced. As a result, supernova neutrinos from explosions 35,000 times more distant than SN1987A could be seen by an improved SK, covering 50% of the entire Universe. Perhaps more important for less distant explosions, neutron tagging of inverse beta events would facilitate the de-convolution of a galactic supernova's various signals, allowing a much more complete interpretation of the physics of the burst. In coincidence with a GW signal, detailed information concerning the explosion's dynamics would become even more valuable.
The key is to add 0.2% by mass of a soluble gadolinium compound like gadolinium sulfate, Gd 2 (SO 4 ) 3 , to the water. Doing so will make >90% of the neutrons visible as a consequence of the gamma rays released by gadolinium's capture of thermalized neutrons. This technique and the various new scientific advances it would make possible was first proposed in the Physical Review Letters article, "GADZOOKS! Anti-neutrino spectroscopy with large water Cherenkov detectors" [356]. The publication of this paper introduced the concept of gadolinium-enhanced water Cherenkov detectors, a transformational technology with a strong impact on the physics community. It provided a clear and cost-effective path to extend and build upon the pioneering Kamioka work of years past.
In order to test this new technology, a dedicated gadolinium R&D facility was built in 2009 in the Kamioka mine near SK. Called EGADS (Evaluating Gadolinium's Action on Detector Systems) [357], it consists of a gadolinium-loaded 200-ton scale model of SK, complete with 240 50-cm photomultiplier tubes, its own DAQ and readout electronics, a novel selective water filtration system, and water transparency evaluation equipment. Based on its successful operation, in 2015 the SK Collaboration formally approved the plan to add gadolinium to SK, designating this new phase of the experiment "SK-Gd". The T2K Collaboration, a 48/89 long-baseline neutrino oscillation experiment based at J-PARC which uses SK as their far detector, formally approved the gadolinium-loading plan in 2016.
In order to prepare SK for the addition of gadolinium, the detector had to be opened and refurbished for the first time in twelve years. There were four main tasks: (1) Replace the photomultiplier tubes that had failed (a few hundred out of 13,000) since the previous in-tank refurbishment in 2006. (2) Fix a small water leak in the SK tank.
(3) Clean any rust and other dirt that had accumulated in the detector since its original completion in 1996. (4) Install additional water piping to increase the total water flow for increased water purification, and to enable better control of the flow direction in the tank.
This in-tank work, which required over 3000 person-days of effort, was successfully carried out between May 2018 and January 2019. By February 2019 the detector had been refilled with pure water, and was ready for the addition of Gd 2 (SO 4 ) 3 . After taking the T2K beam schedule into account, the first gadolinium is now expected to go into the tank in mid-2020. By searching for the small but constant, diffuse flux of neutrinos produced by all core collapse explosions since the onset of stellar formation in the early Universe, it is expected that by 2022 we will have collected the world's first additional supernova neutrinos since SN1987A.

Theoretical studies of supernova explosion and neutrino emission.
In the rest of this section, we report the progress of theoretical study of supernova neutrinos for the detection of SK with the emphasis on detailed description of neutrino transfer and microphysics. In order to provide the prediction of supernova neutrino detections at SK with gadolinium loading, it is mandatory to get rid of the uncertainties in numerical simulations to describe the explosion mechanism and neutrino emission. One of the remaining uncertainties of supernova simulations is the neutrino transfer, which has been routinely treated with approximate ways in multi-dimensional simulations. We perform the sophisticated numerical simulations of neutrino-radiation hydrodynamics based on the direct solver of Boltzmann equation. The first principle-type calculations enable us to determine the final outcome of explosion, provide the solid prediction of neutrino emissions and evaluate the uncertainties of microphysics such as equation of state.

Boltzmann-radiation-hydrodynamics simulation of the core-collapse supernovae.
We have performed several simulations using the Boltzmann-radiation-hydrodynamics code under 2D/3D. The details of the code are illustrated in [358][359][360]. Various nuclear equations of states (EOSs), progenitor models, and the rotational velocities are employed. Here, we show the important features of the obtained results.
First, we examined the effect of the nuclear EOSs [361]. We employ the Lattimer-Swesty (LS) EOS [362] and the Furusawa-Shen (FS) EOS [363]. The former (latter) EOS model employs the soft (hard) nuclear force model and the single nuclear approximation (nuclear statistical equilibrium treatment) for the nuclear composition. The progenitor model is the 11.2 M progenitor model taken from [364]. The simulations with the different EOS models show a significant difference: the LS model shows shock revival while the FS model not (the left panel of Fig. 32). A useful quantity to understand the difference is the timescale ratio, which is defined as τ adv /τ heat , where τ adv := M gain /Ṁ and τ heat := E gain /Q gain are the 49 advection timescale and the heating timescale, respectively. The advection timescale is the timescale with which a fluid element flows through the gain region. The heating timescale is the timescale with which a fluid element is heated and become gravitationally unbound. The gain region is the region where the neutrino heating exceeds cooling. The symbols M gain ,Ṁ , E gain , and Q gain are the mass in the gain region, the mass accretion rate, the total energy including the gravitational binding energy in the gain region, and the net neutrino heating rate in the gain region, respectively. If the timescale ratio exceeds unity, the neutrino heating proceeds sufficiently fast and the explosion succeeds. The right panel of Fig. 32 shows that the timescale ratio of the LS model exceeds unity, while that of the FS model not. Althougḣ M , E gain , and Q gain are similar between the LS and FS models, M gain for the LS model is larger than that of the FS model. This is because the turbulence is stronger for the LS model. Indeed, the prompt convection is stronger for the LS model, and it enhances the neutrino driven convection at the later stage. The stronger prompt convection in the LS model is originated from the stronger photodissociation of the accreting nuclei: the accretion flow outside the shock in the LS model contains more heavy nuclei, and hence the energy consumed by the photodissociation of these nuclei is larger for the LS model. Therefore the proper treatment of the nuclear composition of an EOS is important to assess the influence of EOSs on the CCSNe. Second, we investigated the effect of rotation [365]. The employed EOS and progenitor model are the FS EOS and the 11.2 M model from [364]. The rotational velocity profile is so-called the shelluler rotation profile: Ω(r) = 1 rad/s/(1 + (r/1000 km) 2 ). However, the employed rotational velocity is too slow to affect the postbounce dynamics even though the velocity is almost the highest according to the current stellar evolution theory [366]: the shock radii and other dynamical features are similar between the rotating and non-rotating models. What is more interesting here is the momentum space distributions of neutrinos. Figure 33 shows the neutrino angular distributions in the laboratory frame outside the shock on the equator. Especially, the distribution of the high-energy neutrinos is tilted to the φ-direction. This is because matter rotates and drags the neutrino distribution to the rotational direction. This detailed angular distribution is only accessible by the Boltzmann solver. Since we have such information, we can assess the accuracy of the M1-closure scheme, 50/89  one of the approximate neutrino transport method. The second angular moment of the distribution divided by the zero-th moment is called the Eddington tensor. The M1-closure scheme estimates the Eddington tensor from the energy flux and density of neutrinos. In Fig. 34, we compare the eigenvalues of the Eddington tensor calculated from the distribution function directly and the M1-closure scheme. The largest eigenvalue, or the Eddington factor, is ∼ 20% larger for the M1-closure scheme. This difference originates from so-called raycollision between outward and inward rays. If we can get some information about these rays from neighboring matter, we may improve the accuracy of the M1-closure scheme. 51/89 Third, we suggested a new mechanism for the proto-neutron star kick motion [367]. We performed the core-collapse simulation with the 15 M model of [364] and Togashi-Furusawa (TF) EOS [368]. The TF EOS is one of the most realistic EOSs. The most interesting feature of this model is the proto-neutron star kick. Thanks to the exact treatment of the protoneutron star in our code, the proto-neutron star moves with the velocity of O(10) km/s. So far, the driving force of the proto-neutron star motion is considered to be the gravity of the asymmetric ejecta [369]. However, the driving force observed in our simulation seems the recoil of the asymmetric neutrino emission. With the proto-neutron star motion, the asymmetric distribution of neutrinos are sustained, and hence the driving force by the neutrinos persists. This result suggests the different mechanism of the proto-neutron star motion.
Finally, we have performed the 3D simulation with the Boltzmann-radiationhydrodynamics code. Due to the limited computational resources, the dynamics and the neutrino distributions until ∼ 20 ms after the core bounce are investigated. With this simulation, the neutrino distribution function with no spatial symmetry is obtained for the first time in the world. This result provides us the important clues to understand the behavior of the neutrino in the supernovae.

Light curve of supernova neutrinos.
In order to extract the information as much as possible from the detection of supernova neutrinos from the next supernova in future, it is necessary to prepare the templates of neutrino signals by systematically covering the variation of progenitors and microphysics. One of such systematic sets is the supernova neutrino database provided by [370]. This data is now routinely used to evaluate the event rates at SK, replacing the rather classic data set by [371].
The supernova neutrino database is constructed by the combination of supernova simulations for the set of progenitors with different masses and metallicities. The numerical code of general relativistic neutrino-radiation hydrodynamics under the spherical symmetry [372,373], which is the first-principle calculation, is utilized to study the series of phenomena; gravitational collapse, core bounce, formation of central object, the shock propagation and the associated neutrino emission. The numerical code for the thermal evolution of proto-neutron star with the flux-limited neutrino diffusion in general relativity [374,375] is utilized to study the neutrino emission from the cooling phase of the proto-neutron star over 20 s. In order to connect the two stages, we take out the profile of a central object from the core-collapse simulation and continue the cooling simulation to simulate the shock revival and the cease of accretion in the explosion.
We demonstrated the method to extract the physics information from neutrino emissions of the Galactic supernova based on the supernova neutrino database and additional simulations [376]. We explored the time profile of neutrino detection events at the SK for the set of progenitors with different masses and shock revival time. We provided the basic feature of the rise and fall of neutrino burst to extract the information of the core bounce of massive stars in the early phase within 0.3 s. We explored also the long term behavior of neutrino burst over 20 s up to the phase of fading out. We performed additional long simulations of the proto-neutron star cooling in the case of light and massive neutron stars. We found that the SK and Hyper-Kamiokande are able to detect the long term evolution over 100 s for the massive neutron stars in Galactic supernovae. We proposed a method to determine the mass of central object by plotting the cumulative event number summed from the last 52/89 event detection in time-backward manner. These results are used as the basis to prepare for the future detection of the next supernova in order to determine the details of supernova mechanism and the neutron star. The dependence of the neutrino signals on the equation of state (EOS) is now in progress under the development of the EOS table for supernova simulations as described below.
4.2.5. Information of supernova mechanism from supernova neutrino detectors. The accumulation of predicted signals of supernova neutrinos will enable us to extract the information of supernovae from the detection of neutrino bursts in the future [349,350]. SK will detect ∼10,000 neutrino events from a supernovae in our Galaxy and provide valuable and unprecedentedly detailed information regarding the supernova mechanism. For nearby explosions, pre-supernova neutrinos from the silicon burning stage just before explosion will be studied by both KamLAND and the newly gadolinium-loaded SK.
Events from the dynamical phase will reveal hydrodynamical instabilities and rotation, which are augmented by the detection at IceCube [317,351]. Observation of events from the cooling phase will clarify the birth of compact object with mass and radius, which may constrain the equation of state of dense matter [376].
A combined detection of GWs and neutrinos at KAGRA and SK would serve to unveil the supernova mechanism by the correlation between them [31,352]. The planned next generation of big detectors such as DUNE, JUNO, and Hyper-Kamiokande [353,354] with their extended supernova neutrino detection ranges can be expected to provide additional information (See [355] for further references).
Regarding multimessenger astronomy, in the case of a galactic supernova an early alert is possible since neutrinos are generated earlier in the explosion and therefore arrive before the electromagnetic radiation. Such an alert could be generated and disseminated either through large, directional, high-confidence detectors like the gadolinium-loaded SK acting independently to announce a burst, or a network of detectors which detect lower-confidence signals acting in coincidence [355] to reduce the chance of individual false alerts. 4.2.6. Neutrino bursts from black hole formation. Studying the neutrino burst from the black hole formation is an interesting target of the SK among the variety of gravitational collapse of massive stars [377]. In the case of stars of more massive than those for ordinary supernovae, 40-50M for example, with the intense accretion of matter from outer layers, the retreat of shock wave is inevitable and there is no chance to have the explosion. The mass of central object increases monotonically and attains the maximum mass, which can be supported by the equation of state. The black hole is formed due to the re-collapse of proto-neutron star at the critical mass. The neutrino signal in this case has a characteristic signature with a short duration, typically about 1 s, and increasing average energies and luminosities. This is cased by the increasing mass due to the accretion and the associated increase of density and temperature. This type of short neutrino burst can be a hallmark of black hole formation and can be used to constrain the equation of state of hot and dense matter [378,379].
The detection of GWs and neutrino bursts from the black hole formation may provide the information on the central object and the equation of state in addition to the case of corecollapse supernovae (see sec. 4.1). We revealed the characteristic feature in the frequencies of 53/89 GWs from the dynamical evolution of proto-neutron star toward the black hole formation. We analyzed the time evolution of accreting proto-neutron stars to determine the fundamental and gravity modes of GWs [380] by utilizing the results of the numerical simulations of massive stars with a set of equation of state [377]. The density increase toward the black hole mass leads to the rise of frequencies of GWs as a function of average density of the proto-neutron star. The ratio of the frequencies of the two modes can be characterized by the compactness of the proto-neutron star, therefore, the information of mass and radius. Moreover, the termination of neutrino signal can provide the timing of the black hole formation and the information of maximum density through the combined analysis of GWs and neutrino signals.

4.2.7.
Diffuse supernova neutrino background. The numerical simulations of the supernova explosions and the black hole formation from various massive stars are essential for the detection of diffuse supernova neutrino background, which is the main target of the SK with gadolinium loading. The numerical simulations of the core-collapse supernovae in multidimensions (see also sec. 4.1) and the black hole formation as discussed above have been applied to provide the integrated energy spectra of neutrino emission from the progenitors in the wide range of stellar mass. The feature of neutrino burst depends on the compactness of massive stars, which in turn determines the accretion rate of matter and affects the neutrino emission. In the wide coverage of massive stars with different compactness, the case of black hole formation contribute to hard energy spectra due to the energetic bursts. Hence the contribution of black hole formation for large compactness increases the event rates of diffuse supernova neutrino background. The long term observation at SK and Hyper-Kamiokande can provide the constraint on the critical compactness and the ratio of black hole formation in various stellar collapse [381].

Equation of state for supernovae and neutron star mergers.
The equation of state is one of the important ingredients in core-collapse supernovae as well as neutron star mergers. It is also one of remaining uncertainties in nuclear physics to determine the outcome of explosions, the birth of compact object and the neutrino emission. In addition to the progress of numerical simulations of supernovae, developments of the data table of supernova EOS for simulations have been made over the decades. One of the popular sets of supernova EOS is the Shen EOS based on the relativistic mean field theory [382][383][384], which has been applied to provide extended versions of data table of supernova EOS [385]. The Shen EOS table is publicly available on the web and has set the standard format to provide thermodynamical quantities for the usage of numerical simulations. The Shen EOS has been used in the numerical simulation for the supernova neutrino database.
Recently, the Shen EOS has been revised to improve the density dependence of the symmetry energy [386]. This modification is motivated by the recent progress of observational data on neutron stars from the GW detection and the X-ray observations. The symmetry energy of the original Shen EOS, which was determined by the fitting to the nuclear masses and radii, has been claimed to be rather large as compared with other frameworks. Recent nuclear experiments also provide constraints on the behavior of symmetry energy and help to improve the description of symmetry energy. A new term of density-dependent iso-vector interaction in the relativistic mean field theory provides a smaller symmetry energy while 54/89 keeping good properties of symmetric nuclei and matter. It provides also smaller radii of neutron stars within the observational constraints. The updated Shen EOS is now applied to numerical simulations of core-collapse supernovae and proto-neutron star cooling. It has been recently shown that the influence of the updated Shen EOS at high densities mainly appears in the evolution of proto-neutron stars when the matter becomes neutron-rich and remains minor in the dynamics around the core bounce. It is interesting to see the influence of the full data table of the updated Shen EOS in core-collapse supernovae and neutron star mergers.
In addition to the revisions of the nuclear interaction, the improvement of the supernova EOS of hot and dense matter has been made to describe the mixture of nuclei under the nuclear statistical equilibrium [387,388]. There is recent progress of the construction of EOS tables by the microscopic nuclear many body theories such as the variational method and the Dirac Brueckner Hartree-Fock theory [389,390] based on the nucleon-nucleon interactions, being different from the effective nuclear many body theories. In addition, the sets of EOS tables with systematic coverage of EOS parameters are under the development and will be applied to the study of core-collapse supernovae and neutron star mergers. These development will help to study the signal of neutrinos and GWs and to probe the dense matter in these astrophysical events.

Supernovae simulations on GPUs.
The dynamics of the supernovae is described by hydrodynamic equations for dense matter and the Boltzmann equation for neutrino transport under the gravitational effect described by the theory of general relativity. In addition to these coupled equations, data of physics processes such as sets of equation of state and reaction rates for neutrinos are needed to be implemented in the numerical simulations. Since the neutrino distribution is a function in six dimensions (three dimensions in space and three dimensions in neutrino momentum space), the required computational resources quickly increases with grid resolution of these degrees of freedom. Therefore the numerical simulations of supernovae tend to be challenges in large-scale computational techniques and algorithms.
So far most of large-scale simulations have been performed on massively parallel clusters, such as the K computer. Recently another type of architecture, which makes use of arithmetic accelerators such as GPUs, has become popular in high performance computing. Although the accelerator device has an advantage in cost performance, the code implementation becomes much more involved to achieve desired performance. Furthermore, whether accelerators work efficiently depends strongly on the structure of numerical algorithms.
Presumably for these reasons, application of GPUs to the simulations of core-collapse supernovae has been restrictive. As for the simulation code including the neutrino transport, the VERTEX code was ported to GPUs by employing CUDA [391]. In the VERTEX code, its most time consuming part is calculation of the collision term of the Boltzmann equation that exhausts almost half the simulation time. Offloading one dominant reaction process to GPUs achieved 1.8 times acceleration of whole simulation time.
We develop an efficient scheme to exploit accelerator devices such as GPUs in the numerical simulations of the core-collapse supernovae [392,393]. As the first step, we apply the offloading of simulation bottlenecks to the computations of neutrino-radiation hydrodynamics under spherical symmetry [373]. By adopting the implicit scheme for the evolution equation, an 55/89 iterative linear equation solver for the coefficient matrix is the most time consuming part. To offload this part to the GPU devices, we employ OpenACC as a framework of implementation as well as make use of cuBLAS library that is available on NVIDIA's CUDA environment. We change the data layout to maximize the efficiency of data transfer between the device memory and cores, in accord with so-called coalesced access. With both the OpenACC and cuBLAS, significant acceleration is achieved so that the linear equation solver is no longer a primary bottleneck [392].
The secondary bottlenecks of the original code are the computation of collision term in the Boltzmann equation and the inversion of block diagonal part of the coefficient matrix. The latter is required in the weighted Jacobi-type preconditioner for the iterative linear equation solver [394]. Offloading to GPU devices successfully accelerates these tasks [393]. For the inversion of block matrices, we employed a blocked version of the Gauss-Jordan algorithms [395] that is suitable to the many-core architecture. With these improvements, systematic simulations with better resolution have become feasible. Further optimization as well as extension to 2D and 3D simulations are underway. We are also implementing a code for PEZY-SC processor, another kind of accelerator device, by applying the technique developed for GPUs.

Probes of new physics
5.1. Test of gravity theories using gravitational wave data 5.1.1. Varieties of gravitational waveform within general relativity. The GW signals from coalescing binaries can be decomposed into inspiral, merger and ringdown waveforms.
Here we focus on the inspiral phase, in which the orbital frequency and hence the GW frequency gradually increases as the binary separation decreases. In the case that the binary with aligned spins evolves along a sequence of quasi-circular orbit, the inspiral waveform in the frequency domain can be approximately represented as where the overall amplitude A is ∝ M 5/6 , and (πMf ) −5/3 1 + 20 9 in the post-Newtonian (PN) approximation [153]. Here, η = µ/M is the symmetric mass ratio, v = (πM f ) 1/3 is the orbital velocity, and The second term in the right hand side of the above equation becomes a small contribution because η √ 1 − 4η ≤ √ 3/18 = 0.096 · · · , but we should note that this term can be dominant for spins antialigned with each other.
The first term in the square brackets of Eq. (14) is the term obtained by using the Newtonian dynamics and the quadrupole energy loss formula due to GW emission. So, when we discuss the phase evolution, we count the PN order relative to this leading term. The second term is relatively O(v 2 ), which we refer to as 1 PN order. In the following, we summarize additional effects which should be carefully discriminated from modifications due to extension of general relativity (GR).

56/89
The effect of spins: The lowest PN effect of spins on the phase evolution appears at the 1.5 PN order. The effect depends only on the spin components projected to the direction parallel to the orbital angular momentum at this order, which is encapsulated in the coefficient β. When the spins are not aligned with the direction of the orbital angular momentum, spin precession occurs. As a result, the amplitude modulation occurs, which helps to solve the degeneracy between spins and other parameters [396]. The GW templates for precessing binaries are ready (see, e.g., Ref. [397] and references therein). Eccentric binaries: The orbital eccentricity decays like e ∝ a 19/12 at the late stage of binary evolution, where a is the orbital semi-major axis [398]. Therefore, binaries are circularized well before they merge unless they are borne with a close separation and a large eccentricity. One possible scenario to form binaries that maintain a sufficiently large eccentricity at the coalescence time is the binary formation by capture in star clusters. Hence, detection of eccentricity would be useful to distinguish the formation channels of binaries that become GW sources. In the GW waveform, the effect appears most significantly in the phase evolution, unless the eccentricity is not extremely large. The correction appears in the phase evolution at −19/6 PN order. The minus sign of the PN order reflects the fact that the orbital eccentricity decays rapidly. Dark matter cloud: In the presence of dark matter, the waveform of GWs can be affected [399]. In the scenario of the existence of the dark matter mini-spike around an IMBH [400,401], the waveform of the GWs from an EMRI/IMRI system is modified due to the effects of the gravitational pull, the dynamical friction, and the accretion of dark matter [402][403][404].
They have shown a possibility of the measurement of the power-law index of the dark matter mini-spike radial profile with space-borne GW detectors such as LISA. (See also Ref. [405].) On the other hand, when the mass ratio is around O(10 3 ) ∼ O(10 4 ), energy dissipation due to the dynamical friction is comparable to the total binding energy of dark matter. If one considers a modification of dark matter distribution due to the energy dissipation, the differences of the waveform from the vacuum case becomes much smaller, even the spike has a large power-law index [406].

5.1.2.
Model independent test of waveform consistency with general relativity. From the inspiral waveform, we can determine the binary parameters. At the same time, we can test if there is no deviation from the standard waveform predicted for non-spinning quasi-circular binaries in GR. We can test the presence of spin and eccentricity of each binary component BH. At the same time we can also test the deviation from the prediction of GR.
Here we summarize the test of GR performed by LIGO/Virgo collaborations [35]. The first test is the one to test if the residual data which is obtained by subtracting the best-fit template is consistent with the Gaussian noise. There is no inconsistency observed.
They also have made a consistency test between the estimations of parameters obtained from the inspiral part and from the merger-ring down part. Small perturbation modes of a BH have its own frequencies and damping rates determined by its mass and spin. Those damped oscillation modes, which are called quasinormal modes (QNMs), are excited at the epoch of formation of the remnant BH after coalescence. See, e.g., Ref. [407] as a review of 57/89 QNMs including the history. Among various QNMs, the dominant mode is the fundamental (n = 0) mode with the harmonic ( , m) = (2, 2), and the GW waveform is written as where f R is the oscillation frequency, and τ is the damping time, and t 0 and φ 0 are the starting time and its phase, respectively. A is the amplitude at the starting time. We can find the information of f R and τ which are related to the mass M rem and spin a rem of the remnant BH, in Refs. [408,409]. As for the initial phase, we can maximize the SNR with an analytical method presented in Ref. [410,411] which is based on Ref. [412]. Since it is difficult to determine the starting time (e.g., see Ref. [413] for the estimation), the fractional differences between the quantities for the inspiral and the post-inspiral phases have been used as an inspiral-merger-ringdown consistency test for BBHs in Ref. [35]. In the GR case, we will have within the measurement error range. The analysis result for GW events in the LIGO-Virgo Catalog GWTC-1 is shown in Fig. 2 and Table III of Ref. [35]. Here, we should also note that the SNR only for the ringdown phase is not so high in the current observations. Therefore, the merger-ringdown phase with a combination of phenomenological and analytical BH perturbation theory parameters [414], is treated in the parameterized test of GWs. This analysis result for GW events in the LIGO-Virgo Catalog GWTC-1 is also shown in Fig. 3 of Ref. [35]. On the other hand, when we estimate the mass and spin of the remnant BH (e.g., the right panel of Fig. 4 in Ref. [5]) without testing gravity, we use the inspiral-merger-ringdown waveform from numerical relativity [415][416][417].
Using only the ringdown phase, one of the best way to test gravity is "black hole spectroscopy" by observing multiple QNMs [418] (see also Refs. [419,420] for testing the no-hair theorem with black hole ringdowns, and another simple method is proposed by Ref. [421]). For example, we may consider Eq. (17) with substitution, if the ( = 3, m = 3) QNM is the next dominant mode. In order to do so, it is necessary to extract f R and τ of the weak ringdown signal from the noisy data accurately for each QNM. In Ref. [422]. As the first step, mock data of GWs which include some deviation from the GR prediction was prepared, and it was applied to the following five methods to extract the dominant QNM: (1) plain matched filtering with ringdown part method, (2) matched filtering with both merger and ringdown parts method, (3) Hilbert-Huang transformation method, (4) autoregressive modeling method, and (5) neural network method. It was found that the determination of f R is much easier than that of τ although the accuracy depends on the analysis methods. Interestingly, it turned out that the standard matched filtering does not give the optimal inference in this problem.

58/89
As a unified framework to test possible modifications of gravity, one can use parametrized post-Einstein (PPE) approach [423]. In this approach, we consider the modification of the functions that appear in the GW waveform in GR (13) without specifying the model as where α i and β i specify the amplitude of modification, while a i and b i specify the PN order of modification. The point is that the leading order correction tends to be given by a particular PN order. In many cases the leading order effects appear in the phase evolution. Of course, the test targeting at a particular modification can be more sensitive. LIGO/Virgo collaborations did not give the constraints on the PPE parameters directly. Instead, they examined the constraint on the modification to the model parameters contained in the waveform model IMRPHENOMPv2, following the method proposed in Ref. [424] and implemented by Ref. [425]. For the inspiral part, they test the modification to the existing PN coefficients in the phase evolution up to 3.5 PN order in GR as well as the −1 PN and 0.5 PN coefficients which are not present in GR templates. Also, the modification of the parameters that characterize the merger and ring-down phases is also tested.
For the modification of GW propagation, what they did is similar to PPE, but the constraints are derived for the amplitude of modification per propagation distance by stacking the data. In all cases mentioned above, significant deviation from GR prediction has not been obtained.

Model independent test of extra polarizations.
The GW interferometers basically detect the tidal deformation of the physical distances among nearby test masses, which can be described by the geodesic deviation equation for slowly moving objects Namely, what we can observe is the rank 2 tidal tensor E ij = R i0j0 . The absolute distance cannot be measured by the current GW interferometer, and they measure the difference of the distance changes in different directions. Hence, the trace part of E ij is irrelevant. In generic dynamics of geometry, therefore, there are five components which can be decomposed into three dimensional 1 scalar, 2 transverse vector, 2 transverse traceless tensor components. It is not easy to think of gravity models in which scalar or vector components of the tidal waves are emitted without modifying the GW waveform significantly. For example, in the case of scalar-tensor theories, the effect of scalar wave emission, which also contributes to the excitation of the scalar component of the tidal tensor, appears most significantly in the orbital evolution of the binary due to the radiation reaction that it causes. If we consider such models that can excite larger magnitude of tidal tensor in the scalar or vector mode without losing much energy, the coupling of matter field to those modes tends to be too large to hide the influence on the tests of gravity in solar system.
For the reason mentioned above, sensible models would be a conversion of the excitation energy from the tensor modes, i.e., standard GW modes, to the other exotic modes during the propagation. In that case, the waveform can be identical to the original GW waveform although the phase velocity of the scalar and vector waves can be different from the speed of light. In that case, it is difficult to put some constraint to the conversion of the energy during the propagation. So far, only comparison of likeliness between pure tensor modes and pure 59/89 scalar mode or between pure tensor modes and pure vector modes has been done, finding that pure tensor modes are largely preferred [35,426,427].
In a general framework, a GW signal can be a superposition of scalar and/or vector modes in addition to the ordinary tensor modes: Here F P a is the antenna pattern function of the a-th detector for the polarization P , and h P is the GW waveform for the polarization P , where P = +, ×, V , W , S and L [428][429][430]. V and W modes are vector modes and S and L modes are scalar modes, which satisfy F S a = −F L a [428]. The antenna pattern functions are dependent on a GW source location θ and φ, which are the latitude and longitude, respectively. The separability of the polarization modes of GWs from a point source has been investigated by solving an inverse problem for the first time in [431]. It has been shown that in a nontensorial polarization search, the necessary number of detectors is at least the same as the number of polarization modes to be searched, e.g., at least three detectors are necessary to search for a superposition of two tensor and one scalar modes. The mode separability in the case of compact binary coalescences is subtle because their GW waveforms are well modeled with multiple source parameters. It is shown in [432] that even with correlations and degeneracies among the parameters, a mixture of the polarization modes are separable with the same number of detectors as the number of the modes. Then the scalar or vector modes with similar amplitudes to the tensor modes can be detected. For a binary neutron stars (BNS) observed with the third-generation detectors such as Einstein telescope [433] or cosmic explorer [434], the mode separation with smaller number of detectors is possible due to the long duration of a signal ( 1 hour) and the time-evolving detector antenna patterns [435].
With the help of an electromagnetic counterpart, more polarization modes than mentioned above can be searched once the sky location is fixed. We will discuss this possibility in the next.
Direct polarization test with the help of an electromagnetic counterpart.
We consider to analyse the data from a GW source, S a (t,Ω) + n a (t), where n a is the detector noise. Equation (22) means that tests of all polarization modes with a network of interferometers require at least five detectors. However, Hagihara et al. pointed out that there exist particular sky positions that allow for a vector mode test, because scalar modes, even if they exist, can be perfectly eliminated from a certain combination of the strain outputs at the four detectors only for these sky directions [436] (See Fig. 35). They also found that a vector mode test can be done by using four detectors [437]. They put also a direct upper bound on a vector component from GW170817 event, since the scalar modes can be largely suppressed for this event [437,438].
Very recently, they proposed a general formulation for directly testing extra GW polarization and a certain condition that allows for a scalar test. First, we focus on the network composed of the four ground-based interferometers. For a given source, we know its sky position, as is the case of GW events with an electromagnetic counterpart such as GW170817. Then, one knows exactly how to shift the arrival time of the GW from detector to detector. 60/89 Fig. 35 Seventy sky positions in the earth-centered coordinates. At these points for a GW source, the spin-1 test can be done in principle, because the spin-0 and spin-2 modes can be perfectly eliminated. The numerical calculations were done in Reference [436].
For example, the projection operator is defined to eliminate W , + and × modes. It is denoted as Π aW +× . For this example, h W , h + and h × in the strain outputs {S a } are eliminated as If the coefficient of h V in Eq. (23) vanishes in a certain sky region, there remains only the spin-0 part in the null steam. Therefore, the spin-0 polarization test is possible, if a GW source is found in this sky region. The vanishing coefficient condition is Note that this is invariant for choosing a reference axis for the polarization angle. D 4 = 0 (or sufficiently small D 4 practically) is a condition for directly testing scalar modes separately from the other modes only by four detectors [438]. A possible straightforward procedure of the GW data analysis along this direction is as follows. At first, we determine the sky location of the GW/EM source from multimessenger astronomy, especially by optical and VLBI observations. Secondly, by using the sky location, the arrival time shift for each GW detector is taken into account. Next, the (time-shifted) strain output at each detector is substituted into S a in the left-hand side of Eq. (23). If the left-hand side of Eq. (23) is above the noise level (the last term of the right-hand side), then, the existence of extra GW polarizations could be suggested. If it is comparable to (or lower than) the noise level, a certain direct upper bound can be placed on extra GW polarizations (through the first and second terms the right-hand side). In this procedure, explicit templates of GW waveforms are not needed. In this sense, this test of extra GW polarizations is robust. pipelines would be important for real data analysis, because the nature of detector noises is quite complicated. With four detectors, consistency tests with GR can be done by examining if the wave in the fourth detector, say KAGRA, is consistent with the three-detector network. Equation (22) for strain outputs (a = 1, · · · , 5) from the five detectors including LIGO-India can be always solved for h S − h L , h V , h W , h + and h × . In principle, adding the fifth detector will thus allow for a direct search of extra GW polarizations for any sky region. This would be interesting in probing new gravitational physics beyond GR.

5.1.4.
Gravitational waveform corresponding to various modifications of gravity. The sources of GWs are composed of the objects of extremely strong gravity in the sense that the deviation from the Minkowski spacetime is significant. The spacetime curvature radius is comparable to the size of the system in the case of BH and/or NS binaries.
In the case of BBHs, the emission of electromagnetic counterpart is not expected much. In fact, no detection of electromagnetic counterpart has been reported so far. In the case of binaries that include a NS, we have a chance to detect the electromagnetic counterpart. In fact, various follow-up observations succeeded in finding the counterpart for GW170817, the first coalescing BNS event detected by GWs [9] (see Sec. 3). However, even in that case the emission of the electromagnetic counterpart in the region of truly strong gravity cannot be observed directly because of the large opacity of the high density matter.
In contrast, GWs are emitted from the bulk motion of the objects where the strong gravity is really at work. Hence, GWs are thought to be an indispensable probe of strong gravity. Here, we summarize various modifications on the GW waveform by considering representative extended theories of gravity.

Scalar tensor gravity.
As the most typical modification of gravity, one can consider scalartensor theories. A representative theory will be Brans-Dicke theory [439], defined by the where ω BD is the so-called Brans-Dicke parameter. When we send ω BD to ∞, Einstein gravity is recovered. m a (φ) is the effective mass of a star labeled by a, and it depends on the local value of the scalar field around the star. The effective gravitational coupling depends on ω BD as G eff = (4 + 2ω BD )/(φ(3 + 2ω BD )), and the dimensionless scalar charge is defined by s a := ∂ ln m a (φ)/∂ ln φ, which can be interpreted as the dependence of the gravitational binding energy of the star on G eff . The modification to the GW waveform appears at −1 PN order like ψ(f ) = · · · + (3/128)(πMf ) 5/3 αv −2/3 + 1 + · · · in Eq. (14), where the coefficient α is given by α = −5(s 1 − s 2 ) 2 /(64ω BD ) [440]. This is caused by the dipole radiation of the scalar field. If the scalar charge of two stars are proportional to their mass, i.e., if the dimensionless scalar charges are identical, the dipole radiation does not occur. Hence, the dipole radiation is not expected from equal-mass BNSs.
In this model, numerical approach is also possible, but the effect is larger for lower frequencies. Therefore, the template based on the PN approach mentioned above would be sufficient, unless positive detection of the presence of scalar charge is reported.
One scenario to give a scalar charge to compact stars is spontaneous scalarization in the context of generalized Brans-Dicke model, in which ω BD is promoted to a function of φ [441]. Rather intuitive understanding of this process can be obtained by introducing the canonically normalized field ϕ = 2ω BD (φ)/φdφ. In terms of this new field, the gravitational action would be rewritten as From this expression, one can find that the ϕ-dependent part of the energy of a static body would be given approximately by where, after the curvature is replaced with the energy density ρ using the Einstein equation as an approximation, the deviation from the Minkowski metric is neglected. In the last line, the size of the system is set to R and ϕ is assumed to vanish outside the compact object. If the function φ(ϕ) is convex upward at the origin and the density becomes high enough or the size of the compact object is large enough, non-trivial configuration of the scalar field develops inside the star, which leads to giving a scalar charge to the star. The point is that this happens only in the situation where ρR 2 ≈ M/R is sufficiently large. Within the category of scalar-tensor theory, more generalization is possible. The Horndeski theory [442] is the most general scalar-tensor theory having second-order equations of motion. This theory can, however, be generalized further in a healthy way to the so-called degenerate higher-order scalar-tensor (DHOST) theories [443][444][445], which inevitably have higher-order equations of motion, but still propagate the same number of degrees of freedom as the Horndeski theory does (i.e., one scalar and two tensor modes). The Horndeski and DHOST theories offer us a powerful unifying framework to study the physics of dark energy and modified gravity (see, e.g., [446] for a review). 63/89 Measuring the speed of GWs can constrain the Horndeski and DHOST theories. The Lagrangian for the DHOST theories satisfying c GW = 1 is given by [37][38][39][40] where X := −φ µ φ µ /2 and We thus have 4 free functions of φ and X. For A 3 = 0 and f = f (φ), this reduces to a subclass of the Horndeski theory. Solar System and astrophysical experiments/observations are potentially able to put further bounds on the functions in the above Lagrangian. For Solar System and astrophysical tests of the Horndeski and DHOST theories, it is important to understand the Vainshtein screening mechanism [458], which suppresses the scalar-mediated force and typically is implemented in scalar-tensor theories whose Lagrangian depends on the second derivatives of the scalar field. In the Horndeski subclass with A 3 = 0 and f = f (φ), the Vainshtein mechanism is known to work perfectly for a static spherically symmetric body, Φ = Ψ = −G N M/r, inside a certain radius (which is sufficiently larger than the size of stars and the Solar System), where Φ and Ψ are two metric potentials defined by δg 00 = −2Φ and δg ij = −2Ψδ ij , G N := 1/16πf is the effective Newton constant, and M is the mass contained within r. Therefore, it would be difficult to place Solar System/astrophysical constraints on this subclass.
In generic DHOST theories, the behavior of weak gravitational fields is more interesting. It was found in [41][42][43] where G N and Υ i are written in terms of f and A 3 , and the prime denotes differentiation with respect to r. Since M = constant outside the source body, Eq. (30) shows that gravity is modified only inside astrophysical bodies in generic DHOST theories. (In other words, Vainshtein screening operates nicely in the vacuum exterior region.) This partial breaking of Vainshtein screening was first discovered in [447]. A number of astrophysical tests of Vainshtein-breaking theories have been proposed so far. Going beyond the weak gravity regime, it was shown in [448] that the structure of relativistic stars is quite sensitive to the Vainshtein-breaking effect, which would potentially give tight constraints on DHOST theories in the strong gravity regime. It was pointed out in [449,450] that gravitons can decay into φ in DHOST theories. To avoid this, it is required that A 3 = 0. Otherwise, GWs would not be observed. Upon imposing this constraint we have A 4 = 3f 2 X /2f and A 5 = 0. This yields a special subclass of DHOST theories characterized by 3 free functions, in which the behavior of weak gravitational fields is completely different from (30), as shown in [451,452]. (In fact, Eq. (30) is derived assuming A 3 = 0.) The main result of [451,452] is summarized as follows: (i) a kind of fine-tuning is required even in the vacuum exterior region so that Solar System tests are evaded; (ii) in the interior region the metric potentials Φ and Ψ obey the standard inverse power law, but the two do not coincide; (iii) the interior and exterior values of the effective Newton constant 64/89 are different, All of these features are in sharp contrast with the case of generic DHOST theories. A stringent bound on this particular class of DHOST theories with A 3 = 0 comes from the fact that the measured value of the Newton constant (which is supposed to be G ext ) is different from the gravitational coupling of GWs by a factor of 1 − Xf X /f . This can be seen from the quadratic Lagrangian for GWs, .
Assuming that the scalar radiation does not take part in the energy loss, one may naively replace the Newton constant in the standard quadrupole formula with G GW . Then, from the orbital decay of the Hulse-Taylor pulsar we find that |Xf X /f | 10 −3 . It would be interesting to explore other astrophysical constraints on this special class of Vainshteinbreaking modified gravity.
As a low energy effective theory of the dynamics of the metric, the term in Lagrangian at the lowest order of mass dimension is the cosmological constant and the next leading order is the Einstein-Hilbert action. In this sense, possible natural extension of gravity is to consider to add terms quadratic in curvature to the Lagrangian.
There are three parity preserving terms, R 2 , R µν R µν and R µνρσ R µνρσ . R 2 term is classified as an alternative representation of scalar-tensor theories [453]. One combination, R GB = R 2 − 4R µν R µν + R µνρσ R µνρσ , is know as the Gauss-Bonnet term, which does not contribute to the equations of motion as it can be written as a total divergence. Hence, truly remaining new possibility is, e.g., adding R µν R µν . However, this theory has a ghost degree of freedom, i.e., it contains infinitely many negative norm/energy states, which leads to an absence of stable vacua. This is caused by the presence of higher derivative terms in R µν R µν . Hence, we need to treat such a model simply as a low energy effective theory. Then, the higher derivative terms can be replaced by using the lower order equations of motion, which is now the Einstein equations, up to an appropriate field redefinition [454]. Once we apply this procedure to replace R µν R µν with T µν T µν , the model is reduced to just Einstein gravity with matter fields described by a little complicated Lagrangian.
As a parity violating term quadratic in the curvature, we can think of the Chern-Simons term, RR := αβ σχ R σχ µν R µν αβ /2, where ε µνρσ is the Levi-Civitá tensor. However, this term is also written as a total divergence, and hence it does not contribute to the equations of motion.
One interesting extension is to introduce a scalar field coupled to these higher order curvature terms. The two cases with a scalar field coupling to the Gauss-Bonnet term and the Chern-Simons term have been extensively studied. The Lagrangian is given by where is the length scale that characterizes the coupling of the axion field to the Chern-Simons term. The deviation from GR in such theories can be easily hidden in the weak 65/89 gravity regime, An interesting phenomena in this setup is that BHs can have a scalar hairs. The constraint coming from the GW observations is discussed in Sec. 5.1.5.

Massive gravity.
If we simply add a mass to the graviton, the equations for the GW propagation in flat spacetime h µν = 0 would be expected to be modified to ( − m 2 )h µν = 0. However, this simple-minded extension is not theoretically very attractive. One consistent way that preserves the Lorentz symmetry on Minkowski background is the dRGT ghost-free massive gravity [455]. However, it turned out to be difficult to construct even the homogeneous isotropic universe model in this framework, and hence various generalizations were proposed. In any case, the dRGT ghost-free massive gravity model requires auxiliary metric or tetrad that do not transform as dynamical fields. As a result, the general covariance is not maintained.
In this class of models of massive gravity, linear perturbation predict a large excitation of a scalar graviton around a star, even if we send the graviton mass to a small value, which is known as the van Dam-Veltman-Zakharov discontinuity [456,457]. However, in this limit, the linear analysis cannot be trusted any further and the non-linear effects become large, strongly suppressing the excitation of the scalar graviton. This is the very same mechanism as the Vainshtein mechanism introduced above [458]. As a result, the gravity around massive objects is thought to behave very similarly to GR.

Bigravity.
A natural way to recover the general covariance in the theory with massive graviton is to promote the auxiliary metric field in the dRGT ghost-free massive gravity to a dynamical one. Such a model is know as dRGT bi-gravity theory, which is also shown to be ghost free by Hassan and Rosen [459]. However, the general covariance guarantees the presence of a massless graviton. In this model the field that becomes massive is the second graviton, which corresponds to the difference of two metrices, roughly speaking.
The generation and propagation of GWs from binaries were studied in Ref. [460], in which the possibility of graviton oscillation between two graviton modes: one is massless while the other is massive, was pointed out. In the ordinary situation, only the metric directly coupled to the matter field is curved, while the other remains to be almost flat. In a particular choice of the model parameters, however, the massive graviton becomes very heavy around a star, and hence the difference between the two metrices is suppressed. In this case, the graviton oscillation becomes possible, without violating the test of gravity, say, tests in the solar system.
The Einstein-aether theory, which for hypersurface-orthogonal configurations of the aether vector can be considered as the low energy limit of a candidate theory of quantum gravity called non-projectable Hořava-Lifshitz gravity and which also serves as a testing ground of GR. We thus review the theoretical and observational constraints on the Einstein-aether theory The basic variables of the gravity sector are [461] the fourdimensional metric g µν , the aether vector u µ and a Lagrange multiplier λ, where µ, ν = 0, · · · , 3 and we adopt the signatures (−, +, +, +) for the metric. The total action is [462] S = S ae + S m , where Here ψ denotes matter fields, R and g are the Ricci scalar and the determinant of g µν , and where D µ denotes the covariant derivative compatible with g µν , and It is convenient to introduce linear combinations of c i 's (i = 1, . . . , 4) as c ij ≡ c i + c j and c ijk = c i + c j + c k . The Newton constant G N is related to G ae as [463] . The Minkowski spacetime with u µ = δ 0 µ is a solution of the Einstein-aether theory. It is then straightforward to analyze linear perturbations around the Minkowski background and investigate properties of spin-0, -1 and -2 modes. The coefficients of the time kinetic term of each excitation must be positive, In addition to the ghost-free condition, we must also require the absence of gradient instabilities by demanding that the squared speeds of propagation be non-negative. Moreover, the almost simultaneous detection of GW170817 [9] and GRB 170817A [36] sets a stringent constraint on the speed of the spin-2 mode as −3 × 10 −15 < c T − 1 < 7 × 10 −16 , where c 13 ). This implies |c 13 | < 10 −15 .
Applying the theory to cosmology, one can find that the gravitational constant appearing in the effective Friedman equation is given by [463], One then needs to impose the nucleosynthesis bound as Among the 10 parameterized post-Newtonian (PPN) parameters [465], the only two parameters that deviate from GR in this context are In the weak-field regime, Solar System observations set the constraints [465] |α 1 | ≤ 10 −4 , |α 2 | ≤ 10 −7 .
where (α 1 ,α 2 ) denotes the strong-field generalization of (α 1 , α 2 ) [470] given by [471], Here, σ ae is the sensitivity. Unfortunately, the sensitivities σ ae of a neutron star, which depend on c i 's and the equation of state of nuclear matter [471], are not known for the new ranges of the parameters. We thus leave the analysis of these two constraints to a future work. Putting all theoretical and observational constraints together, we have (38), as well as in the (c 2 , c 14 )-plane. In Fig. 37, we show the constraints (47). All the constraints except (38) are divided into two classes, those in the (c 1 , c 14 )-plane and those in the (c 2 , c 14 )-plane.
For further discussions about the constraints such as the comparison with those on the hypersurface-orthogonal theory [472], see [473].

Dedicated test for particular models.
To optimally constrain particular models, it would be better to employ dedicated test. One can imagine various extensions of the simple GW waveform adopted in the PPE framework.
One simple example is to add a mass to the field that extracts energy from the system through the dipole radiation. Since the radiation is emitted only in the frequency range above the mass scale, the modifications in the orbital evolution and hence in the GW phase evolution show up only at high frequencies [474,475]. Therefore, even if a NS has a scalar hair, such scalar-tensor theories can evade the constraint from pulsar 68/89 observations due to the mass of the scalar field. A simple model of scalar-tensor theories with a massive scalar field has been constrained by several observations: (1) m 10 −27 eV for stability on cosmological scales [476], (2) m 10 −16 eV from the observations of a pulsarwhite dwarf binary [477], and (3) 10 −13 eV m 10 −11 eV [478,479], which relies on the measurements of high spins of stellar mass BHs [480].
As mentioned earlier, Einstein-dilaton-Gauss-Bonnet (EdGB) theory, is an interesting candidate of modification of gravity obtained by adding a quadratic term in curvature. We need to introduce a scalar field, which we call dilaton here, to make the Gauss-Bonnet term relevant. In this model a BH can possess a scalar charge and hence dipole radiation is expected. The constraint from observed GW data on the model in which the scalar field is massive was examined in Ref. [481]. In the massless limit, on the other hand, the current constraint on EdGB theory is √ α EdGB 1.9 km, which is obtained by using low-mass X-ray binaries [482].
Activation of dipole radiation due to massive fields and modifications to GW waveform has been analyzed for the GW events in the catalog GWTC-1 [483] to constrain the magnitude of the modification nd [481]. Since the ground based detectors, LIGO, Virgo, and KAGRA, are sensitive at frequencies around 10-1000 Hz, GW signals for those are useful to test dipole radiation of the massive field in the range 10 −14 eV m 10 −13 eV. For m 10 −14 eV the modification of GWs effectively reduces to the massless limit.
Taking into account the dominant correction to the energy flux due to the dipole radiation, which appears at −1 PN order, the GW waveform in the frequency domain,h(f ), can be expressed as [484,485] The corrections to the amplitude and the phase are, respectively, expressed as [481] with whereω = πf /m, m is the mass of the field, Θ is the Heaviside step function, A is a parameter denoting the amplitude of the dipole radiation relative to the quadrupole one. The modified waveform, Eq. (48) with Eqs. (49) and (50), can model the modification due to the vector field as well as the scalar one. In this parametrization, there are two free parameters, i.e., the mass of the field m and the relative amplitude of dipole radiation A. Note that A is a function of the system parameters, such as the masses and spins, and the parameters of the modified gravity theory that one considers. Using the above waveform as a template, 90% confidence level (CL) constraints on A for each LVC GW event has been derived [481]. Since the modified waveform is valid only for the inspiral phase, the calculations were terminated when the mass scale reaches an appropriately chosen threshold frequency for each event. We should recall that in general a strong constraint on A does not always imply a strong constraint on the modification to theory. This is because A is proportional to the squared difference of the scalar charges of 69/89 constituents of the binary as we will see in Eq. (52). Therefore, GW events with vanishing of this difference are not sensitive to the modification, even for a quite large coupling constant.
The constraints on A can be translated to a constraint on the coupling strength, using the relation [486][487][488] A = 5π 6 where s EdGB describes the spin-dependence of the BH's scalar charges. After combining five BBH events possessing relatively high SNR, i.e., GW150914, GW170814, GW170608, GW170104, GW151226.
√ α EdGB is constrained as √ α EdGB 2.47 km for all mass below 6 × 10 −14 eV for the first time with 90% CL including √ α EdGB 1.85 km in the massless limit. This constraint in the massless limit is much stronger than the results of in Ref. [489,490], while it is accidentally almost the same with that by low-mass X-ray binaries [482].
Since the number of events with a lighter chirp mass, such as BNSs, must increase in the near future, it is interesting to stack those events by assuming various theories in which charged NSs are allowed consistently. Furthermore, the NS-BH binary is also interesting [491], because, in addition to their light chirp mass, A can be large in a theory such that only either NS or BH has a charge, as in the case of EdGB theory. Moreover, future multiband observations will improve the current upper bounds of various modified theories [492].
Propagation of GWs in parity violating gravity.
The nearly simultaneous detection of GWs and the gamma-ray burst from the merger of neutron stars, GW170817/GRB170817 A [9,36], offered us an unprecedented opportunity to measure the speed of GWs, c GW , at the level of one part in 10 15 . This puts a tight constraint on modified gravity as an alternative to dark energy. We now discuss the impact of this constraint on parity-violating sector of gravity [493] Let us consider parity-violating gravity whose action is of the form where L PV is the parity-violating Lagrangian and L φ is the Lagrangian for a scalar field φ which may be coupled nonminimally to gravity. The most frequently studied example is Chern-Simons (CS) gravity, for which L PV is given by Recently, it was found that one can construct other parity-violating Lagrangians for the gravity sector such as [494] 70/89 with L 1 := ε µναβ R αβρσ R ρ µν λ φ σ φ λ , L 2 := ε µναβ R αβρσ R ρσ µλ φ ν φ λ , L 3 := ε µναβ R αβρσ R σ ν φ ρ φ µ , L 4 := φ λ φ λ P , and φ µ := ∇ µ φ. If 4a 1 + 2a 2 + a 3 + 8a 4 = 0 is satisfied, the dangerous Ostrogradsky ghosts are removed in the unitary gauge. 2 We consider GWs h ij (t, x) propagating on a cosmological background. The above examples yield the quadratic Lagrangian of the form where ijk is the antisymmetric symbol, Λ is some energy scale, and α and β are dimensionless functions of time determined essentially by the cosmological evolution of the scalar field. In the most familiar case of CS gravity, we identically have α = β. However, we may have α = β in general. The Lagrangian (57) can thus be used to describe the GW propagation in parity-violating gravity in a unifying way [493] (see also [496] for an effective-field-theory approach to arrive at the same Lagrangian).
In the presence of the parity-violating interactions (57), the + and × polarization modes are coupled in the equations of motion, but one can obtain two decoupled equations by employing the circular polarization basis. In general parity-violating gravity, the equations of motion for GWs in Fourier space can be written as where the prime denotes differentiation with respect to the conformal time, H := a /a, and Here, A labels left and right circular polarization modes and λ A = ±1. It can be seen that CS gravity corresponds to the special case with (c A GW ) 2 = 1. Therefore, from the propagation speed of gravity no constraint can be imposed on CS gravity. However, parity violation in the gravity sector in general leads to (c A GW ) 2 = 1 and thus is tightly constrained by GW170817. For k/a ∼ 100 Hz, one has Λ −1 |α − β| 10 −11 km .
This is the new limit on parity violation in the gravity sector derived in [493]. By parametrizing the modification of GW waveforms during propagation in this theory, the data in the O1/O2 catalog has been reanalyzed by a grid search with the matched filtering method [497]. The results imply that the above constraint improves by 7 digits since the matched filtering is sensitive to the phase modifications, while the constraints on the amplitude modification are similar to the previous result by Yagi and Yang [498].
For a face-on binary, the effects of parity violation of gravity cannot be constrained because only one of the circular polarization modes can be observed. Therefore, obtaining a more stringent constraint or find the signature of parity-violating gravity requires the detection of GW signals from a nearly edge-on binary. 2 The ghost degrees of freedom might not be problematic if the theory is treated as a low-energy truncation of a fundamental theory. Although these ghost modes do cause instabilities if the theory is regarded as a complete theory, it is argued in [495] that the would-be ghost mode that is not manifest in the unitary gauge may be an instantaneous mode which does not in fact propagate and hence is not dangerous.
Exotic compact objects are proposed as alternative models to BHs. They are assumed to be compact enough to possess the light ring, but without the event horizon. In such models, the feature of GWs after the ringdown phase of compact binary coalescences is considered to be different from the BH case. The infalling waves generated during merger-ringdown phase might be reflected at the surface of the object, which would be just absorbed if the horizon exists. The reflected waves will be partly reflected back again at the angular momentum barrier but partly transmitted. This process occurs repeatedly, which phenomenon is dubbed as GW "echoes" [499,500] (see also Ref. [501] for a simple explanation with the causal structure of the Green's function of the wave equation with two potentials having disjoint bounded support). Therefore, the detection of echoes implies the existence of an exotic compact objects. The simplest model of the compact objects is that the horizon is replaced by a reflecting membrane at ∼ Planck proper length from the horizon radius in Kerr spacetime. Abedi et al. have analyzed three O1 BBH events and reported that they found echo signals at 3σ significance level (0.011 in p-value) [502]. In their study, they assume the GW waveform in the frequency domain as where γ is the overall reflection rate, ∆t echo is the interval between neighboring echoes, N is the number of echoes, and φ(f ) is the phase shift. The functionh 0 (f ) is a mergerringdown waveform extracted by cutting off the inspiral part smoothly from the inspiralmerger-ringdown waveform. The parameters γ and ∆t echo are the most relevant parameters to characterize the echo waveforms on which the significance of the signal depends most. Since Abedi et al. used only 32 seconds for the background estimation, Westerweck et al. [503] have improved the background estimation using 4096 seconds and shown lower significance, 0.032 in p-value, with the same template used in the paper by Abedi et al. [502]. The Bayesian analysis is also performed with the same template and small Bayes factors are reported [504,505]. The spacetime considered in Abedi et al. is exactly Kerr spacetime outside the surface, while the assumption of reflection rate in their study is not physically reasonable. Nakano et al. have derived the reflection rate at the angular momentum potential barrier numerically, and calculated a fitting formula for the reflection rate R(f ) [506]. The reflection rate becomes frequency dependent, unlike in the assumption in Abedi et al., and behaves as a high pass filter, i.e., R(f ) → 0 [R(f ) → 1] for larger (smaller) f . As a result, the lower frequency part of the waveform is strongly suppressed in the template given by Nakano et al.
In this setup, Uchikata et al. have shown that the significance of echo signals becomes much lower, p-value is larger than 0.9 for eight BBH events in O1 and O2 [507], although neither reduction nor increase of the significance was observed by increasing the number of samples when the original waveform proposed by Abedi et al. is employed. So far, we have focused on a model dependent search. A morphology independent search was done by Tsang et al. [508,509], where they assumed that echo waveforms for any model can be effectively described as superpositions of sine-Gaussian waveforms. In their 72/89 study, no significant echo signal is found in O1 and O2 events, including the BNS merger event. The reflection rate at the surface can also vary due to the model of the object. In the analytical template proposed by Refs. [510,511], the reflection rate at the object's surface is also assumed as a parameter. The reflection rate at the potential barrier is given from BH perturbations. Further details of testing exotic compact objects includes the echo phenomenon are reviewed in Ref. [512].

Dark matter
Arguably, one of the most important questions in science is to elucidate the nature of dark matter. Identifying the nature of dark matter definitely brings huge impact on cosmology, particle physics, and astronomy. Many dark matter candidates have been proposed in the literature. Since different dark matter candidates yield different type of observational signals, various experiments targeted for particular dark matter have been conducted. In this subsection, we give a brief overview of how the GW experiments can shed light on some dark matter candidates.

5.2.1.
Axion. An axion is a strong candidate for the dark matter [513]. To probe the axion dark matter, it is important to study electromagnetic waves [514] and GWs propagating in the axion dark matter [515,516]. In particular, the gravitational Chern-Simons term allows us to probe the parity violation in the gravity sector. We consider the model specified by the Lagrangian (33) with V = m 2 φ 2 /2.
It is instructive to write down the quadratic action deduced from Eq. (57) using the circular polarization modes as To avoid a ghost, due to the signature flip of the factor 1 − kλ A 2φ /M p , we need the cutoff scale below which we can use this effective action k < k g = 6.7 × 10 9 Hz 10 8 km 2 0.3GeV/cm 3 ρ .
Here, we used the observed local dark matter density ρ. The coherent oscillation induces the parametric resonance of GWs with a frequency corresponding to the axion mass. Using the current upperbound ≤ 10 8 km [517], we can estimate the length scale R ×10 at which the GW grows ten times larger as follows Since the Jeans length r J of the axion dark matter, which is the coherence length of the axion oscillation is much larger than the growth length, there occurs a huge enhancement of the amplitude of GWs, which gives rise to the strong constraint on the Chern-Simons coupling constant or the abundance of the axion dark matter. We refer the reader to the original paper [515] for the details. For the detection strategy to probe this effect is under investigation.

Primordial black holes.
Primordial black holes (PBHs) are a dark matter candidate that is not an elementary particle. The idea of PBHs dates back to 1971 when Hawking proposed a possibility that BHs are produced directly by the gravitational collapse of an overdense region in the primordial universe much prior to the recombination era. PBHs heavier than about 10 15 g do not lose their mass significantly due to the Hawking radiation over the age of the Universe. Since such PBHs are non-relativistic objects, interact with baryonic matter only gravitationally, and do not emit light, they have been considered as a dark matter candidate. Since the original proposal of PBHs, various observational searches for PBHs have been continuously conducted for a wide PBH mass range. So far, there are no observational evidence of PBHs, and upper limits on the PBH abundance have been placed, which, for some mass range, excluded the possibility that PBHs comprise all dark matter under the assumption that all PBHs have the same mass. Yet, there is still an open window where PBHs can provide all dark matter (see Ref. [518] for the most recent summary).
LIGO/Virgo detection of binary BHs provoked explosive research activities on PBHs (see Ref. [8] for a review in this respect). An exciting possibility is that these BHs (or a part of them) are PBHs. Indeed, features of the LIGO BHs such as comparatively large masses and small spins do not contradict with PBHs. It is known that a dominant formation channel of the PBH binaries is by a tidal force by the surrounding PBHs in the radiation dominated era [519]. This channel can explain the merger rate of the BH binaries derived by the LIGO/Virgo observations if the PBHs constitute about 10 −3 fraction of dark matter. Alternatively, a conservative assumption that LIGO BHs were produced by the standard astrophysical processes yields an upper limit on the abundance of stellar mass PBHs. Irrespective of whether LIGO BHs are PBHs or not, that PBHs can comprise at most ∼ 10 −3 of dark matter is the most severe constraint for PBHs in the stellar mass range. This vividly demonstrates the powerfulness of the GW observations to test PBHs as dark matter. It is worth to mention that LIGO/Virgo Collaboration extended searches of the PBH mergers to sub-solar mass range [520].
Although stellar mass PBHs are unlikely to explain all the dark matter, a question still remains that such PBHs constitute a tiny fraction of dark matter. Answering to that question definitely gives a profound consequence to cosmology. To this end, it is essential to discriminate between PBHs and the astrophysical BHs from GW observations. In what follows, we give a short review of the several ideas along this direction proposed in the literature.

Mergers at high redshifts.
Contrary to astrophysical BH binaries which formed long after the Big Bang, PBH binaries exist since the Universe was still dominated by radiation. Due to a broad distribution of orbital radius and eccentricity of the PBH binaries, merger rate of the PBH binaries continuously extends to high redshift even beyond z ∼ 100. On the other hand, the merger rate of the astrophysical BHs sharply drops beyond z ∼ 10. Thus, searches for mergers beyond z ∼ 10 provide a clean test of the existence of PBHs. Reference [521] 74 /89 has shown that the proposed space-interferometer called B-DECIGO is able to achieve this goal.

Mass distribution.
Mass distribution is another observable that can be used to test PBHs. Let m 1 and m 2 be masses of the individual BHs in a binary and R(m 1 , m 2 , t) be the distribution of the merger rate of the BH binaries at cosmic time t. Notice that R includes all the mergers and the observational bias, which depends on the detectors, should be taken into account in translating R to directly measurable merger distribution. Since this bias is possible to handle, we here regard R as an observable quantity.
Here C is a constant independent of m 1 , m 2 , ψ(m) is a function which depends on the PBH mass function, and α ≈ 1. Quite interestingly, R almost linearly depends on the total mass m 1 + m 2 irrespective of the PBH mass function which differs for each inflation model and is unknown 3 . It is worth mentioning that other physical mechanisms leading to BH mergers predict different α values. For instance, the PBH binaries formed by two-body encounters in the low-redshift universe yield α ≈ 1.43 [524]. Astrophysical BH binaries formed in dense star clusters give α ∼ 4 [525]. BH binaries formed in mass-segregated environments such as galactic nuclei give α that vary with the total binary mass m t [526]. Thus measurement of α provides a unique test of the PBH scenario which is independent of the assumption for the PBH mass function. Alternative approach is to directly compare the observed mass distribution with the theoretical ones in the PBH hypothesis computed for several representative PBH mass functions. This program has been already applied to the observational data, and some PBH functions are excluded [527].

Spin of PBHs.
Another potentially powerful GW observable to test the PBH scenario is spin distribution. GW observations are primarily sensitive to a quantity defined by Eq. (2). PBHs form binaries dynamically and there is no correlation between the spin of individual PBHs and the orbital angular momentum. Thus, probability density of χ eff is an even function, which is a robust prediction of the PBH scenario. Concrete shape of the probability density of χ eff is determined by the spin distribution of the individual PBHs W (J). Formally, the spin distribution of PBHs at cosmic time t is given by where P (δ M , J) is the probability distribution of the density contrast δ M and the angular momentum J of an overdense region that collapses to PBH if δ M ≥ δ th . Q(J, J , t) represents the evolution of the PBH spin from its initial value J to J at time t. P (δ M , J) is determined once the statistical properties of the primordial density perturbations are specified.
All these quantities must be determined in order to derive W (J, t). In Refs. [528][529][530], analyses related to P (δ M , J) and Q(J, J , t) have been performed. For the Gaussian primordial density perturbation, the angular momentum of PBHs coming from P (δ M , J) is estimated to be at most a few percent [529,530]. In Ref. [531], δ th (J) for PBHs formed in the radiation dominated epoch was obtained as where a K is the Kerr parameter. This result shows that the threshold increases as the angular momentum is increased and the formation of highly spinning PBHs are suppressed compared to slowly spinning PBHs. On the other hand, PBHs formed in the matter-dominated phase, which could happen prior to the radiation dominated epoch, likely to have significant spin [532].
So far, we have considered GWs from PBHs in the LIGO frequency band. This band corresponds to the merger phase of the stellar mass PBH binaries. However, GWs in different frequency bands could be also generated. The PBH binaries in the inspiral phase, due to their enormous number and weak signals, constitute the low-frequency stochastic GW background. Furthermore, primordial density perturbations as a seed of the PBHs provide another stochastic GW background by the non-linear mode couplings. These low-frequency GWs are not probed by the LIGO-like detectors, but could be detected by future detectors such as LISA.
In what follows, we introduce possible constraints on the abundance of PBHs expected in future experiments by evaluating energy-density spectra of stochastic GW backgrounds in two independent ways. We show that the experiments are sensitive to constrain the fraction to cold dark matter (CDM) for 10 −5 f PBH 1 (10 −13 f PBH 1) in case of the GWs from coalescing events (curvature perturbations) [533].
By using the merger rate of PBH binaries [534], we can calculate the energy-density spectrum of stochastic gravitational wave background (SGWG), which is discussed in Ref. [535] to be where dEGW dνs (ν s ) is the energy spectrum of GW produced by a BBH coalescence [536,537], ν s = (1 + z)ν, and ν cut is the cutoff frequency.
On the other hand, according to Appendix D of Ref. [533], we obtain the dimensionless energy-density spectrum of the induced SGWB to be Ω IGW ν = k 2π = Ω r,0 g * (T (k)) g * (T eq ) g * ,s (T (k)) g * ,s (T eq ) where ν = k/2π denotes the frequency of GW, and the dimensionless wavenumberk = k/k 0 is introduced for simplicity. We used notations of physical quantities written in Ref. [533]. Assuming a null detection of the GW signals from coalescing events in Eq. (71) or curvature perturbations in Eq. (72), we can constrain the abundance of PBHs, which is plotted in Fig. 38. The solid lines give the expected upper limits on the abundance by the non-detection of the GW signals from the coalescing events, and the shaded regions mean the parameter space excluded by the non-detection of the GW signals from curvature perturbations.

Conclusion
Gravitational wave physics and astronomy have finally become science that can be compared with the actual observational data after a long period of preparation. In this paper, we introduced important findings such as the discovery of massive binary black holes and binary neutron stars, which bring us a new implication about physical properties of nuclear matter, heavy element synthesis in the Universe, and the origin of short gamma-ray bursts. It is expected that this trend will be further accelerated as the international network of gravitational wave detectors on the ground expands and space gravitational wave antennas extends. Further improvements in sensitivity are planned for ground-based gravitational wave interferometers, including the construction of the third generation detectors [433,434]. An increase in sensitivity by an order of magnitude increases the chance to detect an event having a very high signal-to-noise ratio, and at the same time increases the event rate by a factor of 1000 or more, which enables highly statistically accurate analyses of gravitational waves. It is certain that future gravitational wave observations will revolutionize the 77/89 understanding of the physical universe that we have accumulated so far, powered by the improvement of the accuracy of gravitational wave source position and multi-messenger follow-up observations for gravitational wave events. It is clear that we have not seen the great potential of this new probe of physics, gravitational waves, yet. As the observation technology advances, it becomes more important to analyze and interpret the accumulated data correctly. The real benefit from gravitational wave observations can be extracted only when we can construct an appropriate physical picture to explain the observed signal, by combining the data from other observational means to follow-up gravitational wave detections. To that end, we need to further promote multi-messenger astronomy, and advance the theoretical research for understanding gravitational wave sources. In addition, in order to detect gravitational wave signals that have not been discovered yet to open up a new paradigm, it is necessary to further develop gravitational wave data analysis methods, making full use of the accumulated theoretical knowledge.