{'fraction_non_alphanumeric': 0.026650787349477307, 'fraction_numerical': 0.015588196374222575, 'mean_word_length': 4.321222363047458, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Formation of black hole and emission of gravitational waves.

Numerical simulations were performed for the formation process of rotating black holes. It is suggested that Kerr black holes are formed for wide ranges of initial parameters. The nature of gravitational waves from a test particle falling into a Kerr black hole as well as the development of 3D numerical relativity for the coalescing binary neutron stars are discussed.


The existence of black holes.
In 1967, the first pulsar(pulsating radio source) was found. At present the number of the observed pulsars is about 1700 and its period of the pulse ranges from 1.56 ms to about 10 s. 1) The pulsar is identified as the neutron star from the following argument. Let us assume that the period of the pulsar is due to the rotation of the unknown object of mass M and radius R. We do not know M and R at the moment but we know that the gravitational force should be larger than the centrifugal force for the equilibrium. This requires GM R 2 > RΩ 2 , [1] where G and Ω are the Newton's gravitational constant and the angular frequency of the pulsar, respectively. Equation [1] is rewritten as [2] whereρ is the mean density of the unknown object. For the shortest period of 1.56 ms, Ω = 4 × 10 3 s −1 . Thenρ is larger than 5.8 × 10 13 gcm −3 . In such extremely high density matter, an electron is easily captured by a proton to be a neutron so that the matter in the pulsar consists of mainly neutrons. This is the reason why the pulsar is identified as the neutron star. The radius of the neutron star should be small and is considered to be about 10 km.
Theoretically the mass of the neutron star can be determined as a function of the central density if the equation of state of the neutron star matter is known (i.e. the relation between the pressure and the density). In general, the mass first increases with the increase of the central density and then begins to decrease with a maximum value for each equation of state. Although the equation of state in such high density matter is not known well, various models of the equation of state give the maximum mass ranging from 0.7 to about two times the mass of the sun. Observationally masses of several neutron stars are determined and their values are about 1.4 times the mass of the sun. 1) The existence of the maximum mass of the neutron star is definite? Since we do not know the true equation of state of the neutron star matter, one may think that the maximum mass might be infinite in principle. However there exists a theoretical upper limit of the maximum mass of the neutron star irrespective of the equation of state of the high density matter. Assume that: 1) The Einstein equation which describes the equilibrium structure of the neutron star is correct. 2) The sound velocity of the matter is smaller than the light velocity.
3) The equation of state for low density matter is known. Then it is proved that the maximum mass of the neutron star is at most 3.2 times the mass of the sun irrespective of the equation of state. 2), 3) Therefore if a certain compact object has the mass larger than 3.2 times the mass of the sun, one can say definitely that the object is not a neutron star. The mass of the object is too large so that the gravitational collapse starts and would not stop. In short, the object would be a black hole. At present there exist about eighteen such black hole candidate X-ray stars 4) since their mass is larger than 3.2 times the mass of the sun.
More massive black hole candidates exist. For example by examining the orbits of the stars near the center of our galaxy, Ghez et al. (2003) 5) suggested the existence of the black hole of mass about 4 × 10 6 times the mass of the sun in the center of our galaxy. The other method to identify a black hole candidate is to measure the velocities of the stars near the central part of the extragalaxies. From the size of the spatial distribution of the stars there and the measured velocities, the existence of very massive black hole is suggested. The number of such black hole candidates at present is about 40 and the mass of the black holes ranges from 10 6 to 10 9 times the mass of the sun. 6) 2. The final state of the gravitational collapse. Let us consider a neutron star of mass near the maximum mass. If the extra mass falls into this neutron star from outside, the mass will exceed the maximum possible mass. Since the mass is too large for the equilibrium, the collapse of the neutron star begins and the size of the star decreases. The gravity increases so much that it becomes difficult for any matter to escape from the star. Finally the trapped surface will be formed. Here from the trapped surface even the outgoing light rays can not expand due to the strong gravity. Penrose and Hawking proved that under such a situation the space-time should have a singularity eventually under the plausible conditions on the pressure and the energy density of the matter. 7) Here singularity means that the space time is incomplete and the Einstein equations are violated there. Therefore the formation of the singularity is a crisis of the Einstein's theory of gravity since such a collapse could happen in reality. Penrose then proposed the "cosmic censorship hypothesis"; the singularity in the collapse should be covered by the event horizon in nature so that we can not observe the singularity of the space time in principle. 8) This means that inside the event horizon there exists a singularity where the Einstein equations are violated but the effect of the violation can not reach outside the event horizon where we are living. Therefore as far as the space time outside the event horizon is concerned, the Einstein theory is not violated and is a self-consistent effective theory.
What is the final ultimate fate of the gravitational collapse such as in the collapse of the neutron star of mass exceeding the maximum one? In such a collapse at first the matter density and the space time structure would rapidly change and the gravitational waves would be emitted. However eventually we expect that the space time would approach the stationary state. Under the cosmic censorship hypothesis (; the assumption that all singularities in space time are hidden behind the non-singular event horizon) Hawking proved that such a stationary final state is an axially symmetric solution to vacuum Einstein equations 7) and then Israel (1967) 9) and Carter (1971) 10) proved that they form discrete continuous families each depending on at most two parameters. Robinson (1975) proved that the Kerr solution is the unique one of the Israel-Carter theorem. 11) Here the Kerr solution is characterized by two parameters; the mass M and the angular momentum per unit mass a. Therefore two parameters of the Israel-Carter theorem are the mass and the angular momentum of the Kerr solution. Here we should notice that the Kerr solution has a black hole structure (i.e. the singularity is covered by the event horizon) only when a is smaller than GM/c where G and c are the Newton's gravitational constant and the light velocity, respectively.
However, if we do not adopt the cosmic censorship hypothesis, many other stationary solutions to vacuum Einstein equations are found. Especially Tomimatsu and Sato (1973) found a series of solutions characterized by the mass, the angular momentum, and the deformation parameter. 12) Their solutions have naked singularities, that is, they do not have the black hole structure. This means that the singularities of the space time exist outside the event horizon so that we can observe the singularities in principle. In the Kerr solution also, the singularity is naked if a is larger than GM/c. If the cosmic censorship is true, all the black holes in galactic X-ray sources and the center of the galaxies should be Kerr black holes with a < GM/c. They cannot be naked singularities like in Tomimatsu-Sato solutions. 12) Here the natural question arises; in the real gravitational collapse, the Kerr black hole is always formed irrespective of the initial conditions? If, for example, the final state were one of the Tomimatsu-Sato solution, the Einstein's theory of gravity should be modified.
Let us consider a star of mass M and angular momentum J. We can define a non-dimensional angular momentum q by Then if this star becomes a black hole as a result of gravitational collapse, q should be smaller than unity under the cosmic censorship hypothesis. Since the definition of q contains only the angular momentum and the mass of the system other than the physical constants, we can argue the value of q for any systems even if it will not evolve to a black hole. For example, the sun will not be a black hole while q is 0.18. For rapidly rotating massive stars which are believed to evolve to black holes, the value of q is usually greater than unity. Then what will happen when q > 1? Under the cosmic censorship, if all the matter of the star collapses to be a black hole, the value of q is too large. In this case, is the naked singularity formed finally? Is a Kerr black hole in reality formed irrespective of the initial conditions even when q < 1? To answer these questions, we needed the development of the method to obtain the numerical solutions of the Einstein equations for nonspherical dynamical space time.
3. Numerical relativity. Before 1976, it was difficult to compute the dynamical evolution of the non-spherical space time obeying the Einstein equations numerically due to the low power of the computer as well as the lack of the good formalism and the good numerical method. The basic variables in the Einstein equations consist of the density of the fluid and the velocity of the fluid as well as ten components of the metric tensor. Four components of the metric tensor are coordinate degrees of freedom. In general relativity, one can adopt any time and spatial coordinates so that these four components of the metric tensor express this freedom explicitly. However in numerical simulations, the choice of the coordinates is crucial for non-crushing computations. Secondly, we give initial data of the density and the velocity of the fluid as well as the gravitational waves. Then we have four equations to determine the initial data of the metric tensor. From Einstein equations we have four time dependent equations to determine the density and the velocity of the fluid as well as twelve time dependent equations to determine six components of the metric tensor and their * ) For details of the formalism and the numerical method, see refs. 15), 16). time derivatives. Usually the coordinate conditions are expressed by the elliptic type partial differential equations so that we have another four equations to determine four components of the metric tensor.
In 1976, Smarr first performed the numerical simulations of the collision of two black holes. 13) This is the first numerically generated non-spherical space time obeying the Einstein equations. He solved the vacuum space time outside the event horizon of two black holes. What our group tried to solve was the collapse of the axially symmetric rotating stars leading to the formation of black holes so that the matter density is not zero in our case. The first subject we attacked was the collapse of a spherically symmetric star. 14) In this first study we tried to find a numerical method which can be easily extended to axially symmetric cases. As a result we found the importance of numerical treatment of the regularity conditions at the origin. In the spherical polar coordinate (r, θ, φ), the origin is r = 0 with θ and φ being arbitrary. This means that the point in the real space does not have one to one correspondence to the coordinate point. The regularity condition guarantees that these apparently different coordinate points correspond to a single point in the real space. Since we are solving the general space time, if we lose the regularity condition numerically, the space time we are solving becomes a different one from what we like to simulate. Therefore to guarantee the regularity condition numerically is extremely important. This led us to the use of the regularized variables to ensure the regularity automatically. We developed these regularized variables to axially symmetric cases where we require the regularity condition on the symmetry axis. After many trials we completed a numerical code for the first time. * ) Using the formalism and the numerical method we developed, we performed many numerical simulations of the collapse of axially symmetric rotating stars. 15)- 19) We summarize here the main results. We use the cylindrical coordinates (R, φ, Z). The system is assumed to be axially symmetric and plane symmetric about Z = 0 plane. Units of mass, length and time are taken as where M i is the initial mass of the collapsing star. For example for M i = ten times the mass of the sun, L = 15 km and T = 5 × 10 −5 s.
We have simulated three types of collapsing rotating stars.
Type A (Formation of Super Massive Black Hole). The initial radius of the star r 0 is taken to be 14.5. At t = 0, the star is falling radially with the free fall velocity. The initial angular velocity (Ω) is given by where Ω 0 is a constant. The equation of state is where p, ρ and are the pressure, the density and the internal energy per unit mass, respectively. Type B (Formation of Super Massive Black Hole). The initial conditions are the same as Case A except that the initial angular velocity is more strongly centralized as Type C (M i =ten times the mass of the sun). The initial density distribution is taken as The equation of state is given as Typical results of the numerical simulations are shown in Fig. 1a and b. A50 in the caption of Fig. 1a means the simulation of Type A with the value of nondimensional angular momentum q = 0.5. In Fig. 1a the contour lines show the density and the arrows show the velocity vector. QMAX and VEL-MAX mean the maximum density and the maximum velocity, respectively. We are showing only Z ≥ 0. The density and the velocity for Z < 0 can be obtained from plane symmetry. The time in Fig. 1a is 0.287 in our unit explained before and is near the initial stage. The time in Fig. 1b is 18.6 and is near the final stage of A50. The dashed line shows the apparent horizon, that is, even the outgoing light rays can not expand from this surface. The general theorem tells us that the event horizon exists outside the apparent horizon. Even if one may emit the light ray We show only the density for Z > 0. The density for Z < 0 is obtained from the symmetry with respect to Z = 0 plane. Each line corresponds to ρ = QMAX · 10 − n 2 for n = 1, 2, · · ·, 11. QMAX is shown in the figure outward from the point inside the dashed line, the light will never reach us living on the earth. We see almost all the matter is inside the horizon so that the rotating black hole is formed in this model A50. Fig. 2a shows the density contour at Time=11.5 of model A105. We see some low density matter near the equator is expanding but the main part of the star is collapsing. In Fig. 2b, the dashed line shows the apparent horizon so that the black hole is formed. Some of the low density disk like matter is expanding but nothing peculiar is happening outside the apparent horizon. This suggests that the naked singularity is not formed. Fig. 3a-d show the evolution of Model A146. The result is completely different from two models before. The apparent horizon is not identified but the formation of the jet is seen in Fig. 3d. From these numerical results, we can say that there is a critical value of q c ∼ 1 for the formation of the apparent horizon in this Type A.
For Type B, there also exists q c ∼ 0.92 while q c ∼ 0.86 for Type C. These three types of the collapse are different from each other in their equation of state, the rotation law and the initial density distribution. However, the critical value of q c for the formation of the black hole is almost the same. In all the models in which the black hole is formed, nothing peculiar seems to happen outside the event horizon. If we remember that the singularity of the Kerr black hole is hidden by the event horizon for q < 1, our numerical simulations support the cosmic censorship hypothesis. 20) After our simulations, Piran (1985, 1986) performed five series of numerical simulations similar to ours and found that q c ranges from 0.8 to 1.2. 21), 22) Shibata (2000) also obtained the similar results 23) and Sekiguchi and Shibata (2004) generalized the criterion by defining q as the effective q parameter in the stellar central region 24) and confirmed that q c is ∼ 1. Therefore the results first obtained by our group are consistent with later simulations.
4. The gravitational radiation from the black hole. In the dynamical space time such as the collapse of the rotating star, the gravitational waves are emitted. Smarr (1976) first numerically computed the amount of the gravitational waves for the two black hole collision 13) while Stark and Piran (1985) did for their formation of rotating black holes. 21) In these two types of simulations the black holes were formed so that the variation of the space time was very large. We expected that the amount of the emitted gravitational waves was also large but their results showed that they were very small and at most ∼ 10 −3 M where M is the mass of the system. The key factor for this unexpected result is the symmetry of the system. Let us consider N particles rotating in the same circle with a regular equi-angle interval around the black hole. 25) It is found that the luminosity of the gravitational wave decreases exponentially as a function of N . For N → ∞, the system becomes an axially symmetric rotating ring and the luminosity is zero. It is known that a rotating ring does not emit the gravitational wave if its radius does not change. Physically each element of the ring does emit the gravitational waves but their phase cancels out due to the phase cancellation effect. Therefore the main reason for the weak emitter of the gravitational waves in Smarr's simulation as well as Stark and Piran's simulation seems to be the axial symmetry of the system. For non axially sym- Observationally, on the other hand, Hulse and Taylor found in 1974 the first binary pulsar PSR1913+16. 26) This binary pulsar is expected to coalesce in 300 million years from now on due to the emission of the gravitational waves. In the final phase of the coalescence, the strong gravitational waves will be emitted since the system is non axially symmetric. Theoretically the binary pulsar is formed as a result of the evolution of the binary massive stars. If the mass of these stars is large enough for the formation of black holes, it is possible that the binary black hole system is formed instead of the binary neutron star. 27) Therefore the coalescence of binary black holes is also possible theoretically although observationally its existence is not confirmed.
Since the full numerical simulation for such a non axially symmetric coalescence process was not possible at that time, we tried to mimic the process by a test particle of mass µ falling into a Kerr black hole of mass M and the specific angular momentum a. If µ M , we can consider the effect of the energy and momentum tensor of the test particle as the perturbation to the Kerr black hole space time. Considering = µ/M 1 as a smallness parameter, we can linearize the Einstein equations known as nonlinear differential equations. In this approximation the strong gravity and the fast motion of the source is taken into account but the non-linearity of the system is sacrificed.
There are two methods for the perturbation of the black hole space time. First we can perturb the metric tensor itself. Then using the tensor harmonics and the Fourier expansion with respect to time with the angular frequency ω, we have the radial equation with the source term due to the test particle orbiting the black hole. This radial equation is a second order ordinary differential equation with respect to the radial coordinate r so that it is easily solved under the appropriate boundary condition. From the solution, we can determine the amplitude of the gravitational waves for the angular frequency ω.
The above method is, however, applicable only for a = 0 case. In this case, the space time is spherically symmetric and is called as a Schwarzshild black hole. For a = 0, the space time is not spherically symmetric so that the tensor harmonics expansion of the metric tensor is not appropriate. However, for the perturbation of the Weyl tensor, * * ) Teukolsky (1973) showed that we can obtain the radial equation similar to the spherically symmetric case which is the second order ordinary differential equation with respect to the radial coordinate of the Kerr metric. 28) This equation is called the Teukolsly equation. The Teukolsky equation is similar to the Schroedinger equation and we can express its form symbolically as where we call V (r) and S(r) as the potential term and the source term, respectively. If we treat only the perturbation in the vacuum, the Teukolsky equation has no problem. However if we consider the space time perturbed by a test particle falling into a black hole, the Teukolsky equation has two difficulties: 1) The potential term has a long range nature so that the asymptotic behavior is not good for numerical computations. 2) The source term diverges like r 3.5 for r → ∞ so that it is difficult to obtain the accurate result numerically.
In the case of the metric perturbation, the radial equation does not have such difficulties so that if we can perform the metric perturbation of the Kerr metric such difficulties will not exist. Inspired by these considerations, Sasaki and Nakamura (1981) first considered the perturbation of the Weyl tensor for the spherically symmetric space time. 29) In the metric perturbation of the spherically symmetric space time, two kinds of perturbation equations are known depending on the parity of the tensor harmonics. The odd parity case is called the Regge-Wheeler equation (1957) 30) while the even parity case is called the Zerilli equation (1970). 31) Both the Regge-Wheeler and the Zerilli equations have short range potentials and converging source terms. 32) While the perturbation of the Weyl tensor for the spherically symmetric space time is known as * * ) The Weyl tensor is the traceless part of the Riemann tensor. While the Riemann tensor consists of the second derivative of the metric tensor and expresses the curvature of the space time geometrically and the tidal force physically. 33) which has a long range potential and a diverging source term. Sasaki and Nakamura noticed the transformation found by Chandrasekhar and Detweiler (1975) 34 Now we can compute the gravitational waves by a particle falling into a Kerr black hole. The important point here is that this computation may mimic the simulation of the coalescing binary black holes and give us the physical insight of the results. The gravitational waves emitted by a particle falling into a black hole is characterized by the energy E and the angular momentum L z per unit mass of the test particle as well as a of the black hole. We usually use the units of c = G = M = 1, so that E, L z and a are in units of c 2 , GM/c and GM/c, respectively. We consider E = 1 case so that the particle velocity is zero at r = ∞. First we show the results for a = 0. 37), 38) In Fig. 4 we show the wave pattern of the gravitational wave emitted by a test particle falling into a Schwarzshild black hole (a = 0) as a function of the retarded time t − r * where r * = r + 2M ln(r/2M − 1). h + is the plus polarization mode of the gravitational wave. Since the wave amplitude is in proportion to µ and decreases as 1/r we show rh + /µ in the figure. Fig. 4a-d cor- respond to L z = 0, 2, 3 and 3.9, respectively. * * * ) For L z = 3.9, the change of the angle ∆φ ∼ 2π, so that this may mimic the coalescence of the binary black hole. As shown in Table I, the energy radiated in this case is E GW = 0.5µ/M µc 2 . If we extrapolate this to µ → M , the efficiency E GW /M becomes 3.1%. Since this case is non axially symmetric, the efficiency is large as we expected. Fig. 5 shows the wave forms emitted by a ring of L z = 2.0 (a) and L z = 3.9 (b). The wave form for L z = 2.0 (a) resembles that of the wave form of the rotating collapse calculated by Stark and Piran (1985). 21) Moreover the very low efficiency shown as a = 0, L z = 2(ring) in Table I is * * * )

Bardeen-Press-Teukolsky equation (BPT equation)
For this case, Detweiler and Szedenits (1979) 39) also computed the energy using BPT equation with the integration by parts. also similar. Physically in the rotating collapse of the star to a Kerr black hole, the star becomes the oblate spheroidal shape due to the rotation. The central main part of the spheroid forms the black hole like space time first. This process does not contribute much to the gravitational waves. Then the remaining ring like part falls down to the newly formed Kerr black hole. This ring like part mainly contributes to the emission of the gravitational waves. This is the physical reason why the perturbation calculation and the full numerical simulation give almost the same wave form and the efficiency of the energy radiated as the gravitational waves.
In the case of the Kerr black hole (a = 0), even if the angular momentum is zero, the particle orbit falling in the equatorial plane is different from that falling along the symmetry z-axis. 40) This is due to the effect of dragging of the inertial frame of the Kerr black hole. In reality the wave forms from them are different. Fig. 6 shows the gravitational wave forms from a particle with zero angular momentum falling into the Kerr black hole a = 0.99 (the solid line) and the Schwarzshild black hole a = 0 (the dashed line). For a = 0.99, after the initial burst of the radiation we see the damped oscillation. This damped oscillation is the very characteristic of the black hole called the quasi normal mode of the black hole whose frequency is a complex number. The real part and the imaginary part of the quasi normal mode depend on the mass M and the angular momentum a of the black hole. From the comparison of the solid and the dashed lines, we see the real part increases and the imaginary part decreases with the increase of a.
If we observe the quasi normal mode, we can determine the mass and the angular momentum of the black hole. Let us consider the particle with angular momentum falling into the Kerr black hole in the equa- torial plane. 41)- 43) In this case, the sign of the angular momentum has the meaning. We call the particle is co-rotating (counter-rotating) if the angular momentum of the particle is parallel (anti-parallel) to that of the black hole. Fig. 7 shows the gravitational wave forms for a co-rotating and a counterrotating case. The solid line shows a = 0.85 and L z = 2.6 case while the dashed line does a = 0.85 and L z = −2.25 case. In Table I we show the energy radiated for a = 0.99, L z = 2 for which the change of the angle ∆φ ∼ 2π and mimics the coalescence of the binary spinning black hole. If we extrapolate this to µ → M , the efficiency E GW /M becomes 9.1% which is the strong gravitational wave source. In Table I, Dim means the dimension of the system and 2D means that the system is axially symmetric. As we have already discussed, we see that the efficiency for 2D case is very small. 5. 3D numerical relativity. Around 1988, in Japan a supercomputer with ∼GBytes (Giga bytes) memory and a few GFLOPS (Giga Floating point Operations Per Second) speed appeared. This means that non axially symmetric simulations were possible. For example, 3D post-Newtonian calculation with 150 × 150 × 150 Cartesian grids requires 1 GBytes memory. The estimated computation time was 100 hours up to 50,000 time steps. Therefore we started such simulations for final merging phase of Hulse-Taylor binary pulsar. Due to the emission of the gravitational waves, the separation of the binary neutron star decreases and finally they will coalesce. When the separation is ∼500 km, the frequency of the gravitational wave is ∼20 Hz. After three minutes, two neutron stars merge to be a single black hole. The wave for this so called "the Last Three Minutes" is one of the main gravitational wave sources for detectors such as TAMA300, 44) LIGO, 45) VIRGO, 46) GEO600 47) and proposed LCGT. 48) Theoretically we first need the accurate theoretical template of the wave form for "the Last Three Minutes". When the separation is large, we can regard each neutron star as a point particle and the analytic post Newtonian expansion calculations are adequate. The template has parameters such as the mass of each neutron star. Comparing the theoretical template with the observation in future, we can determine the mass of each neutron star, for example. This will open a new window to the universe and is called "gravitational wave astronomy".
However in the final merging phase, a point particle approximation is not good but numerical simulations are needed. For such a numerical simulation we first consider the radiation reaction due to the emission of the gravitational waves. In the case of the electro-magnetic wave, the radiation formula starts from the second time derivative of the dipole of the charge and the radiation reaction term is in proportion to the third time derivative of the dipole moment. However for the gravitational wave case, the mass dipole moment is constant since the mass dipole is in proportion to the center of mass. The radiation formula starts from the third time derivative of the mass quadrupole moment and the radiation reaction term is expressed by the fifth time derivative of the mass quadrupole moment. Therefore to compute the radiation reaction term, we must perform the fifth time derivative numerically but the result is highly noisy due to the truncation errors so that it was hopeless to include the radiation reaction terms.
A key idea to overcome this difficulty was as follows. 49  is obtained without the numerical time derivative.
To obtain the fifth time derivative, we perform the numerical time derivative of the fourth time derivative of the mass quadrupole moment which is not noisy.
Using this new formalism to compute the fifth time derivative of the mass quadrupole moment accurately, we performed many numerical simulations of coalescing binary neutron stars including the effects of the radiation reaction for the first time. Here we show a few of them. 51)-55) We show in Fig. 8a-i, 51) the density contour and the velocity vector (arrows) on the equatorial plane of the coalescing binary neutron star. The mass of each star is 1.5 times the mass of the sun. The time is shown on the upper right corner in units of milliseconds. The thick circle shows the radius of the black hole. Fig. 8a is the initial stage. Due to the emission of the gravitational waves, the merging of two neutron stars proceeds and finally the black hole is formed. (Fig. 8i) We show in Fig. 9 the gravitational wave form from the coalescence of the binary neutron star of Fig. 8 when it occurs at 30 million light years from the earth. 51) The amplitude of the wave is around 10 −21 and is the target sensitivity of such as LIGO. 45) The energy emitted in the gravitational wave is 3% of the mass energy. This is just the value expected in the perturbation calculation in section 4 (Table I).
We show in Fig. 10 the density contour and the velocity vector (arrows) on the equatorial plane of the coalescing unequal mass binary neutron star. 52) The masses are 1.7 and 1.28 times the mass of the sun, respectively. We see that the smaller mass neutron star is tidally disrupted and merged. The energy emitted in the gravitational wave is about 2.9% of the mass energy.
In 1992, Oohara and Nakamura first took into account the effect of the post-Newtonian terms to the hydrodynamics equations. 53) They compared the results with those without the post-Newtonian terms. In Fig. 11 we show the density and the velocity on the x-y plane. The left and right figures are for the Newtonian(N) and the post-Newtonian(PN) calculations, respectively. Fig. 12 shows wave forms of h + and h × observed on the z-axis at 30 million light years. The solid and dashed lines are for PN and N, respectively. We see the difference between the two cases. Therefore we should take into account the full details of the general relativity for the accurate wave forms. However in 1990's the power of the computer was not enough to simulate the full set of the Einstein equations in 3D. Looking back to 1987 or so again, however, the power of the computer was enough to simulate the propagation of the 3D gravitational waves without matter. I tried first to simulate the propagation of the gravitational waves whose solution is known analytically. 16) When we consider the low amplitude gravitational wave, we can linearize the Einstein equations. Then we can obtain the general solution analytically. Using this analytic solution, I started to construct a full 3D general relativistic code. I first omitted the matter terms and treated only the metric tensor. At t = 0, I put the low amplitude gravitational wave whose solution is known analytically. All the non-linear terms of the Einstein equation were included in my code. In each time step, I compared the numerical solution with the analytic one and found that the accuracy was not good. Since the non-linear term was included in my code, I first considered that the non-linear terms were the origin of the bad accuracy. Then I omitted the non-linear terms in the code but still the accuracy was bad. This difficulty opened a new formalism of the Einstein equations in 3D numerical relativity.
In the vacuum Einstein equations, the basic variables are ten components of the metric tensor. Four of them are coordinate degrees of freedom. Then true degrees of freedom seem to be six(= 10−4) spatial components of the metric tensor. However, the Einstein equations have ten components. Six components of the Einstein equations are essentially the time evolution equations of six spatial components of the metric tensor. The other four components of the Einstein equations are not the evolution equations but called the constraint equations of six spatial components of the metric tensor. This means that six spatial components of the metric tensor are not free although they obey six evolution equations. True degree of freedom is therefore two(= 6 − 4) which correspond to two degrees of the polarization of the gravitational waves. In 3D code, we first(t = 0) put the initial data of six spatial components of the metric tensor compatible with four constraint equations. Next we choose four components of the metric tensor which express the coordinate degrees of freedom. We then get the numerical value of six spatial components of the metric tensor at t = t + ∆t from the evolution equations. Analytically, the constraint equations are guaranteed in this formalism, but numerically they are not. The key idea is to use auxiliary variables in the evolution equations using the constraint equations. Then the accuracy of the numerical code becomes much better since the constraint equations are used (for details see refs. 16), 57), 56)). Fig. 13 shows one component of the metric tensor in the equatorial plane for the 3D time evolution of the gravitational waves. 16) In Fig. 14, the dotted line shows the accuracy of the numerical solution compared with the analytic one for low amplitude gravitational waves. 16) We see that the accuracy is a few %.
In 1993, I first performed the numerical simulation of the coalescence of the binary neutron stars. 59), 56) In Fig. 15, I show the first 3D full general relativistic simulations of coalescing binary neu-tron stars. The figure shows contours of the gravitational wave in x-y, y-z, z-x plane. t is the time in units of 5 × 10 −6 s. For t = 29.7 we see the spiral gravitational wave pattern in x-y plane while in y-z and z-x plane the pattern is expanding mainly along the z-axis. This wave pattern is typical for the quadrupole mode. However the lack of the power of super computers at that time did not allow me for full 3D accurate simulations for more detailed study.
The new formalism of the Einstein equations for numerical relativity proposed by me in 1987 was developed by Shibata and Nakamura (1995) 57) as well as Baumgarte and Shapiro (1999). 58) It is now called BSSN(Baumgarte-Shapiro-Shibata-Nakamura) formalism. The BSSN formalism is now widely adopted in 3D numerical codes. For example, very recently Baker et al. (2006) 60) and Campanelli et al. (2006) 61) succeeded in computing the gravitational waves from coalescing binary non spin-ning black holes in full general relativity adopting the BSSN formalism. Both of papers found that ∼3% of the rest energy is converted to the gravitational waves. It is surprising that this number is almost the same as the extrapolation from the perturbation calculations listed in Table I (line 2). As for the merger of binary neutron star, Shibata, Taniguchi and Uryu (2005) 62) as well as Shibata and Taniguchi (2006) 63) adopted the BSSN formalism and performed the detailed full general relativistic simulations with the detailed equation of state and found that the amplitude of the gravitational wave is ∼10 −21 at the distance of 30 million light year. This amplitude is very similar to the result of Nakamura and Oohara 51) with the post Newtonian approximation.