* Your assessment is very important for improving the workof artificial intelligence, which forms the content of this project

# Download Merger of binary neutron stars with realistic equations of state in full

Survey

Document related concepts

Astrophysical X-ray source wikipedia , lookup

Standard solar model wikipedia , lookup

Photon polarization wikipedia , lookup

Hawking radiation wikipedia , lookup

Astronomical spectroscopy wikipedia , lookup

Kerr metric wikipedia , lookup

Main sequence wikipedia , lookup

Gravitational lens wikipedia , lookup

Accretion disk wikipedia , lookup

Nuclear drip line wikipedia , lookup

Gravitational wave wikipedia , lookup

Stellar evolution wikipedia , lookup

Transcript

PHYSICAL REVIEW D 71, 084021 (2005) Merger of binary neutron stars with realistic equations of state in full general relativity Masaru Shibata,1 Keisuke Taniguchi,2 and Kōji Uryū3 1 Graduate School of Arts and Sciences, University of Tokyo, Komaba, Meguro, Tokyo 153-8902, Japan 2 Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA 3 Astrophysical Sector, SISSA, via Beirut 2/4, Trieste 34013, Italy (Received 30 January 2005; published 18 April 2005) We present numerical results of three-dimensional simulations for the merger of binary neutron stars in full general relativity. Hybrid equations of state are adopted to mimic realistic nuclear equations of state. In this approach, we divide the equations of state into two parts as P Pcold Pth . Pcold is the cold part for which we assign a fitting formula for realistic equations of state of cold nuclear matter slightly modifying the formula developed by Haensel and Potekhin. We adopt the SLy and FPS equations of state for which the maximum allowed Arnowitt-Deser-Misner (ADM) mass of cold and spherical neutron stars is 2:04M and 1:80M , respectively. Pth denotes the thermal part which is written as Pth th 1" "cold , where , ", "cold , and th are the baryon rest-mass density, total specific internal energy, specific internal energy of the cold part, and the adiabatic constant, respectively. Simulations are performed for binary neutron stars of the total ADM mass in the range between 2:4M and 2:8M with the rest-mass ratio QM to be in the range 0:9 & QM 1. It is found that if the total ADM mass of the system is larger than a threshold Mthr , a black hole is promptly formed in the merger irrespective of the mass ratios. In the other case, the outcome is a hypermassive neutron star of a large ellipticity, which results from the large adiabatic index of the realistic equations of state adopted. The value of Mthr depends on the equation of state: Mthr 2:7M and 2:5M for the SLy and FPS equations of state, respectively. Gravitational waves are computed in terms of a gauge-invariant wave extraction technique. In the formation of the hypermassive neutron star, quasiperiodic gravitational waves of a large amplitude and of frequency between 3 and 4 kHz are emitted. The estimated emission time scale is & 100 ms, after which the hypermassive neutron star collapses to a black hole. Because of the long emission time, the effective amplitude may be large enough to be detected by advanced laser interferometric gravitational wave detectors if the distance to the source is smaller than 100 Mpc. Thermal properties of the outcome formed after the merger are also analyzed to approximately estimate the neutrino emission energy. DOI: 10.1103/PhysRevD.71.084021 PACS numbers: 04.25.Dm, 04.30.–w, 04.40.Dg I. INTRODUCTION Binary neutron stars [1,2] inspiral as a result of the radiation reaction of gravitational waves, and eventually merge. The most optimistic scenario based mainly on a recent discovery of binary system PSRJ0737-3039 [3] suggests that such mergers may occur approximately once per year within a distance of about 50 Mpc [4]. Even the most conservative scenario predicts an event rate approximately once per year within a distance of about 100 Mpc [4]. This indicates that the detection rate of gravitational waves by the advanced LIGO will be 40–600 yr1 . Thus, the merger of binary neutron stars is one of the most promising sources for kilometer-size laser interferometric detectors [5,6]. Hydrodynamic simulations employing full general relativity provide the best approach for studying the merger of binary neutron stars. Over the last few years, numerical methods for solving coupled equations of the Einstein and hydrodynamic equations have been developed [7–15] and now such simulations are feasible with an accuracy high enough for yielding scientific results [12]. With the current implementation, radiation reaction of gravitational waves in the merger of binary neutron stars can be taken into account within 1% error in an appropriate computational 1550-7998= 2005=71(8)=084021(26)$23.00 setting [12]. This fact illustrates that it is now a robust tool for the detailed theoretical study of astrophysical phenomena and gravitational waves emitted. So far, all the simulations for the merger of binary neutron stars in full general relativity have been performed adopting an ideal equation of state [8,9,12 –14]. (But, see [16] for a simulation in an approximately general relativistic gravity.) For making better models of the merger which can be used for quantitative comparison with observational data, it is necessary to adopt more realistic equations of state as the next step. Since the lifetime (from the birth to the merger) of observed binary neutron stars is longer than 100 million yrs [2], the thermal energy per nucleon in each neutron star will be much lower than the Fermi energy of neutrons [17,18] at the onset of the merger. This implies that for modeling the binary neutron stars just before the merger, it is appropriate to use cold nuclear equations of state. During the merger, shocks will be formed and the kinetic energy will be converted to the thermal energy. However, previous studies have indicated that the shocks are not very strong in the merger of binary neutron stars, in contrast to those in iron core collapse of massive stars. The reason is that the approaching velocity at the first contact of two neutron stars is much smaller than the orbital velocity and the 084021-1 2005 The American Physical Society MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ sound speed of nuclear matter (e.g., [12]). This implies that the pressure and the internal energy associated with the finite thermal energy (temperature) are not still as large as those of the cold part. From this reason, we adopt a hybrid equation of state in which the finite-temperature part generated by shocks is added as a correction using a simple prescription (see Sec. II B). On the other hand, realistic equations of state are assigned to the cold part [19–21]. The motivation which stimulates us to perform new simulations is that the stiffness and adiabatic index of the realistic equations of state are quite different from those in the -law equation of state with 2 (hereafter referred to as the 2 equation of state) which has been widely adopted so far (e.g., [12]).1 It can be expected that these differences will modify the properties of the merger quantitatively as described in the following. Since the realistic equations of state are softer than the 2 one, each neutron star becomes more compact (cf. Figure 2). This implies that the merger will set in at a more compact state which is reached after more energy and angular momentum are already dissipated by gravitational radiation. Namely, compactness of the system at the onset of the merger is larger. This will modify the dynamics of the merger, and accordingly, the threshold mass for prompt black hole formation (hereafter Mthr ) will be changed. The adiabatic index of the equations of state is also different from that for the 2 equation of state. This will modify the shape of the hypermassive neutron stars,2 which are formed after the merger in the case that the total mass is smaller than Mthr . Previous Newtonian and postNewtonian studies [23–25] have indicated that for smaller adiabatic index of the equations of state, the degree of the nonaxial symmetry of the formed neutron star becomes smaller. However, if its value is sufficiently large, the formed neutron star can be ellipsoidal. As a result of this change, the amplitude of gravitational waves emitted from the formed neutron star is significantly changed. Since the adiabatic index of the realistic equations of state is much larger than that of the 2 equation of state for supranuclear density [19–21], the significant modification in the shape of the hypermassive neutron stars and in the amplitude of gravitational waves emitted from them is expected. The paper is organized as follows. In Sec. II A, II B, and II C, basic equations, gauge conditions, methods for extracting gravitational waves, and quantities used in the analysis for numerical results are reviewed. Then, the 1 In this paper, we distinguish the stiffness and the magnitude of the adiabatic index clearly. We mention ‘‘the equation of state is softer (stiffer)’’ when the pressure at a given density is smaller (larger) than another. Thus, even if the adiabatic index is larger for the supranuclear density, the equation of state may be softer in the case that the pressure is smaller. 2 The hypermassive neutron star is defined as a differentially rotating neutron star for which the total baryon rest-mass is larger than the maximum allowed value of rigidly rotating neutron stars for a given equation of state: See [22] for definition. PHYSICAL REVIEW D 71, 084021 (2005) hybrid equations of state adopted in this paper are described in Sec. II D. In Sec. III, after briefly describing the computational setting and the method for computation of initial condition, the numerical results are presented. We pay particular attention to the merger process, the outcome, and gravitational waveforms. Section IV is devoted to a summary. Throughout this paper, we adopt the geometrical units in which G c 1 where G and c are the gravitational constant and the speed of light. Latin and Greek indices denote spatial components (x; y; z) and space-time p components (t; x; y; z), respectively: r x2 y2 z2 . ij ij denotes the Kronecker delta. II. FORMULATION A. Summary of formulation Our formulation and numerical scheme for fully general relativistic simulations in three spatial dimensions are the same as in [12], to which the reader may refer for details of basic equations and successful numerical results. The fundamental variables for the hydrodynamics are : rest-mass density, " : specific internal energy, P : pressure, u : four-velocity, and dxi ui vi t; (1) dt u where subscripts i; j; k; denote x; y and z, and the space-time components. The fundamental variables for geometry are : lapse function, k : shift vector, ij : metric in three-dimensional spatial hypersurface, e12 ~ij e4 ij : conformal three-metric, and Kij : detij , extrinsic curvature. For a numerical implementation of the hydrodynamic equations, we define a weighted density, a weighted fourvelocity, and a specific energy defined, respectively, by ut e6 ; (2) u^ i hui ; (3) P ; (4) ut where h 1 " P= denotes the specific enthalpy. General relativistic hydrodynamic equations are written into the conservative form for variables , u^ i , and ^ and solved using a high-resolution shock-capturing e, scheme [11]. In our approach, the transport terms such as @i are computed by an approximate Riemann solver with third-order (piecewise parabolic) spatial interpolation with a Roe-type averaging [26]. At each time step, ut is determined by solving an algebraic equation derived from the normalization u u 1, and then, the primitive variables such as , ", and vi are updated. An atmosphere of small density 109 g=cm3 is added uniformly outside neutron stars at t 0, since the vacuum is not allowed in the shock-capturing scheme. The integrated mass of the 084021-2 e^ hut MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) atmosphere is at most 1% of the total mass in the present simulation. Furthermore, we add a friction term for a matter of low density 109 g=cm3 to avoid infall of such atmosphere toward the central region. Hence, the effect of the atmosphere for the evolution of binary neutron stars is very small. The Einstein evolution equations are solved using a version of the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism following previous papers [7,9,12,27]: ~ij , , A~ij e4 Kij ij Kk k , and the trace We evolve of the extrinsic curvature Kk k together with three auxiliary ~ik using an unconstrained free evofunctions Fi jk @j lution code. The latest version of our formulation and numerical method is described in [12]. The point worthy to note is that the equation for is written to a conservative form similar to the continuity equation, and solving this improves the accuracy of the conservation of the ArnowittDeser-Misner (ADM) mass and angular momentum significantly. As the time slicing condition, an approximate maximal slice (AMS) condition Kk k 0 is adopted following previous papers [9]. As the spatial gauge condition, we adopt a hyperbolic gauge condition as in [12,28]. Successful numerical results for the merger of binary neutron stars in these gauge conditions are presented in [12]. In the presence of a black hole, the location is determined using an apparent horizon finder for which the method is described in [29]. Following previous works, we adopt binary neutron stars in quasiequilibrium circular orbits as the initial condition. In computing the quasiequilibrium state, we use the socalled conformally flat formalism for the Einstein equation [30]. A solution in this formalism satisfies the constraint equations in general relativity, and hence, it can be used for the initial condition. The irrotational velocity field is assumed since it is considered to be a good approximation for coalescing binary neutron stars in nature [31]. The coupled equations of the field and hydrostatic equations [32] are solved by a pseudospectral method developed by Bonazzola, Gourgoulhon, and Marck [33]. Detailed numerical calculations have been done by Taniguchi and part of the numerical results are presented in [34]. B. Extracting gravitational waves Gravitational waves are computed in terms of the gaugeinvariant Moncrief variables in a flat space-time [35] as we have been carried out in our series of paper (e.g., [12,36,37]). The detailed equations are described in [12,37] to which the reader may refer. In this method, we split the metric in the wave zone into the flat background and linear perturbation. Then, the linear part is decomposed using the tensor spherical harmonics and gaugeinvariant variables are constructed for each mode of eigenvalues l; m. The gauge-invariant variables of l 2 can be regarded as gravitational waves in the wave zone, and hence, we focus on such mode. In the merger of binary neutron stars of nearly equal mass, the even-parity mode of l; jmj 2; 2 is much larger than other modes. Thus, in the following, we pay attention only to this mode. Using the gauge-invariant variables, the luminosity and the angular momentum flux of gravitational waves can be defined by r2 X dE 2 [email protected] RE j2 [email protected] RO (5) lm j ; 32% l;m t lm dt r2 X dJ O [email protected] RElm RElm j [email protected] RO lm Rlm j; 32% l;m dt (6) where RElm and RO lm are the gauge-invariant variables of even and odd parities. The total radiated energy and angular momentum are obtained by the time integration of dE=dt and dJ=dt. To search for the characteristic frequencies of gravitational waves, the Fourier spectra are computed by Z R lm f e2%ift Rlm tdt; (7) where f denotes a frequency of gravitational waves. Using the Fourier spectrum, the energy power spectrum is defined as X dE % r2 jR ffj2 f > 0; (8) df 4 l2;m0 lm where for m 0, we define q R lm f jR lm fj2 jR lm fj2 m > 0; (9) and use jR lm fj jR lm fj for deriving Eq. (8). We also use a quadrupole formula which is described in [37–39]. As shown in [38], a kind of quadrupole formula can provide approximate gravitational waveforms from oscillating compact stars. In this paper, the applicability is tested for the merger of binary neutron stars. In quadrupole formulas, gravitational waves are computed from 2 1 2 d I kl hij Pi k Pj l Pij Pkl ; (10) 2 r dt2 where I ij and Pij ij ni nj (ni xi =r) denote a tracefree quadrupole moment and a projection tensor. In fully general relativistic and dynamical space-times, there is no unique definition for the quadrupole moment Iij . Following [37–39], we choose the formula as Z Iij xi xj d3 x: (11) Then, using the continuity equation, the first time derivative is computed as Z (12) I_ ij vi xj xi vj d3 x: 084021-3 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ PHYSICAL REVIEW D 71, 084021 (2005) To compute I"ij , we carry out the finite differencing of the numerical result for I_ ij . In this paper, we focus only on l 2 mass quadrupole modes. Then, the gravitational waveforms are written as 2s 14 5 h fR 1 cos2 ,cos2’ r 64% 22 3 s 15 R sin2 ,5; (13) R22 1 cos2 ,sin2’g 64% 20 and provide the amplitude of a given mode measured by an observer located in the most optimistic direction. s 2 5 R22 cos, sin2’ R22 cos, cos2’; h r 64% (14) in the gauge-invariant wave extraction technique, and 1 I"xx I"yy h 1 cos2 , cos2’ I"xy 1 cos2 , r 2 I"xx I"yy 2 " sin , ; sin2’ Izz (15) 2 " Ixx I"yy 2 cos, sin2’ I"xy cos, cos2’ ; h r 2 (16) in the quadrupole formula. In Eqs. (13) and (14), we use the variables defined by M R22 RE22 RE22 p r; 2 (17) (18) R20 RE20 r: For the derivation of h and h , we assume that the wave part of the spatial metric in the wave zone is written as dl2 dr2 r2 1 h d,2 sin2 ,1 h d’2 2 sin,h d,d’; (19) RE21 0 and Ixz Iyz 0 since we assume the and set reflection symmetry with respect to the equatorial plane. In the following, we present s 5 R ; R (20) 16% 22 s 5 R ; R 16% 22 (21) C. Definitions of quantities and methods for calibration In numerical simulations, we refer to the total baryon rest-mass, the ADM mass, and the angular momentum of the system, which are given by Z M d3 x; (24) 1 I @ dSi 2% r!1 i Z e5 ~ ~ij 2 k 4 5 2 k ~ A A Kk Rk e H e d3 x; 3 16% ij (25) 1 I ’i A~i j e6 dSj 8% r!1 Z 1 ~j 1 6 i ~ij Ai @j ’i A~ij ’k @k e Ji ’ 8% 2 2 ’j @j Kk k d3 x; (26) 3 where dSj r2 @j rdcos,d’, ’j [email protected] j [email protected] j , ^ Ji u^ i , and R~k k denotes the e , H ut e, ~ij . To derive the expressions Ricci scalar with respect to for M and J in the form of volume integral, the Gauss law is used. Here, M is a conserved quantity. We also use the notations M1 and M2 which denote the baryon rest-mass of the primary and secondary neutron stars, respectively. In terms of them, the baryon rest-mass ratio is defined by QM M2 =M1 1. In numerical simulation, M and J are computed using the volume integral shown in Eqs. (25) and (26). Since the computational domain is finite, they are not constant and decrease after gravitational waves propagate to the outside of the computational domain during time evolution. Therefore, in the following, they are referred to as the ADM mass and the angular momentum computed in the finite domain (or simply as M and J, which decrease with time). The decrease rates of M and J should be equal to the emission rates of the energy and the angular momentum by gravitational radiation according to the conservation law. Denoting the radiated energy and angular momentum from the beginning of the simulation to the time t as &Et and &Jt, the conservation relations are written as J in the gauge-invariant wave extraction method, and as the corresponding variables, Mt &Et M0 ; (27) A I"xx I"yy ; (22) Jt &Jt J0 ; (28) A 2I"xy ; (23) where M0 and J0 are the initial values of M and J. We check if these conservation laws hold during the simulation. in the quadrupole formula. These have the unit of length 084021-4 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) Significant violation of the conservation laws indicates that the radiation reaction of gravitational waves is not taken into account accurately. During the merger of binary neutron stars, the angular momentum is dissipated by several 10%, and thus, the dissipation effect plays an important role in the evolution of the system. Therefore, it is required to confirm that the radiation reaction is computed accurately. The violation of the Hamiltonian constraint is locally measured by the equation as 5 2 ij k k 5 2 ~ ~ ~ ~ f Aij A Kk & Rk 2%H 8 3 8 ~ j R~ k j& j2%H 5 j 8 k 1 5 2 jA~ij A~ij j Kk k 2 : (29) 3 8 Following [26], we define and monitor a global quantity as 1 Z f d3 x: (30) H M Hereafter, this quantity will be referred to as the averaged violation of the Hamiltonian constraint. D. Equations of state Since the lifetime of binary neutron stars from the birth to the merger is longer than 100 million yrs for the observed systems [2], the temperature of each neutron star will be very low ( & 105 K) [17,18] at the onset of merger; i.e., the thermal energy per nucleon is much smaller than the Fermi energy of neutrons. This implies that for modeling the binary neutron stars just before the merger, it is appropriate to use cold nuclear equations of state. On the other hand, during the merger, shocks will be formed and the kinetic energy will be converted to the thermal energy to increase the temperature. However, previous studies have indicated that the shocks in the merger are not strong enough to increase the thermal energy to the level as large as the Fermi energy of neutrons, since the approaching velocity at the first contact of two neutron stars is much smaller than the orbital velocity and the sound speed of nuclear matter. This implies that the pressure and the internal energy associated with the finite temperature are not still as large as those of the cold part. From this reason, we adopt a hybrid equation of state. In this equation of state, we write the pressure and the specific internal energy in the form P Pcold Pth ; (31) " "cold "th ; (32) TABLE I. The values of pi we choose in units of c G M 1. i pi (SLy) pi (FPS) i pi (SLy) pi (FPS) 1 2 3 4 5 6 7 8 0.1037 0.1956 39264 1.9503 254.83 1.3823 1:234 1:2 105 0.15806 0.220 5956.4 1.633 170.68 1.1056 0:703 2 104 9 10 11 12 13 14 15 9 105 4 0.75 0.057 0.138 0.84 0.338 9 105 5 0.75 0.0627 0.1387 0.56 0.308 and " are computed from hydrodynamic variables and ^ Thus, "th is determined by " "cold . e. For the cold parts, we assign realistic equations of state for zero-temperature nuclear matter. In this paper, we adopt the SLy [20] and FPS equations of state [19]. These are tabulated as functions of the baryon rest-mass density for a wide density range from 10 g=cm3 to 1016 g=cm3 . To simplify numerical implementation for simulation, we make fitting formulas from the tables of equations of state, slightly modifying the original approach proposed in [21]. In our approach, we first make a fitting formula for "cold as "cold 1 p1 p2 p3 p4 1 p5 p6 p7 1 fp8 p10 p12 p13 fp8 p10 fp9 p11 p14 p15 fp9 p11 ; (33) where 1 : (34) ex 1 The coefficients pi i 1–15 denote constants, and are listed in Table I. In making the formula, we focus only on the density for 1010 g=cm3 in this work, since the matter of lower density does not play an important role in the merger. Then, the pressure is computed from the thermodynamic relation in the zero-temperature limit d" Pcold 2 cold : (35) d With this approach, the accuracy of the fitting for the pressure is not as good as that in [21]. However, the first law of the thermodynamics is completely satisfied in contrast to that in [21]. In Fig. 1, we compare Pcold and "cold calculated by the fitting formulas (solid curves) with the numerical data tabulated (dotted curves)3. It is found that two results agree approximately. The relative error between two is within 10% for > 1010 g=cm3 and & 2% for supranuclear density with * 2 1014 g=cm3 . fx 3 where Pcold and "cold are the cold (zero-temperature) parts, and are written as functions of . Pth and "th are the thermal (finite-temperature) parts. During the simulation, The tables for the SLy and FPS equations of state, which were involved in the LORENE library in Meudon group (http:// www.lorene.obspm.fr), were implemented by Haensel and Zdunik. 084021-5 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ PHYSICAL REVIEW D 71, 084021 (2005) (a) (b) FIG. 1. Pressure and specific internal energy as functions of baryon rest-mass density (a) for the SLy and (b) for the FPS equations of state. The solid and dotted curves denote the results by fitting formulas and numerical data tabulated, respectively. In Fig. 2, we show the relations among the ADM mass M, the total baryon rest-mass M , the central density c , and the circumferential radius R for cold and spherical neutron stars in the SLy and FPS equations of state. For comparison, we present the results for the 2 polytropic equation of state P Kp 2 which was adopted in [12]. In the polytropic equations of state, there exists a degree of freedom for the choice of the polytropic constant Kp . Here, for getting approximately the same value of the maximum ADM mass for cold and spherical neutron stars as that of the realistic equations of state, we set Kp 1:6 105 cgs unit: (36) In this case, the maximum ADM mass is about 1:72M . We note that for the 2 equation of state, the ADM mass M, the circumferential radius R, and the density can be rescaled by changing the value of Kp using the following rule: (a) M / Kp1=2 ; R / Kp1=2 ; and / Kp1 : (37) (b) FIG. 2. (a) ADM mass (solid curves) and total baryon rest-mass (dotted curves) as functions of central baryon rest-mass density c and (b) relation between the circumferential radius and the ADM mass for cold and spherical neutron stars in equilibrium. ’FPS’ and ’SLy’ denote the sequences for the FPS and SLy equations of state, respectively. The relations for the 2 polytropic equation of state with Kp 1:6 105 in the cgs unit are also drawn by the dashed curves. 084021-6 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) Hence, the mass and the radius are arbitrarily rescaled although the compactness M=R is invariant in the rescaling. Figure 2 shows that in the realistic equations of state, the central density and the circumferential radius are in a narrow range for the ADM mass between 0:8M and 1:5M . Also, it is found that neutron stars in the realistic equations of state are more compact than those in the 2 polytropic equation of state for a given mass. Namely, the realistic equations of state are softer than the 2 one. On the other hand, the adiabatic index d lnP=d ln for the realistic equations of state is much larger than 2 for the supranuclear density [19–21]. These properties result in quantitatively different results in the merger of two neutron stars from those found in the previous work [12]. The thermal part of the pressure Pth is related to the specific thermal energy "th " "cold as Pth th 1"th ; (38) where th is an adiabatic constant. As a default, we set th 2 taking into account the fact that the equations of state for high-density nuclear matter are fairly stiff. (We note that for the ideal nonrelativistic Fermi gas, th 5=3 [40]. For the nuclear matter, it is reasonable to consider that it is much larger than this value.) To investigate the dependence of the numerical results on the value, we also choose th 1:3 and 1.65. The thermal part of the pressure plays an important role when shocks are formed during the evolution. For the smaller value of th 1, local conversion rate of the kinetic energy to the thermal energy at the shocks should be smaller. III. NUMERICAL RESULTS A. Initial condition and computational setting Several quantities that characterize irrotational binary neutron stars in quasiequilibrium circular orbits used as initial conditions for the present simulations are summarized in Table II. We choose binaries of an orbital separation which is slightly larger than that for an innermost orbit. Here, the innermost orbit is defined as a close orbit for which Lagrange points appear at the inner edge of neutron stars [33,41]. If the orbital separation becomes smaller than that of the innermost orbit, mass transfer sets in and dumbbell-like structure will be formed. Until the innermost orbit is reached, the circular orbit is stable, and hence, the innermost stable circular orbit does not exist outside the innermost orbit for the present cases. However, we should note that the innermost stable circular orbit seems to be very close to the innermost orbit since the decrease rates of the energy and the angular momentum as functions of the orbital separation are very small near the innermost orbit. The ADM mass of each neutron star, when it is in isolation (i.e., when the orbital separation is infinity), is chosen in the range between 1:2M and 1:45M . Models SLy1212, SLy1313, SLy135135, SLy1414, FPS1212, FPS125125, FPS1313, and FPS1414 are equal-mass binaries, and SLy125135 and SLy135145 are unequal-mass ones. For the unequal-mass case, the mass ratio QM is chosen to be * 0:9 since all the observed binary neutron stars for which each mass is determined accurately indeed have such mass ratio [2]. Mass of each neutron star in model SLy125135 is approximately the same as that of PSRJ0737-3039 [3], while the mass in model SLy135145 is similar to that of PSRB1913 16 [1]. The total baryon rest-mass for models SLy1313 and SLy125135 and for models SLy1414 and SLy135145 are approximately identical, respectively. For all these binaries, the orbital period of the initial condition is about 2 ms. This implies that the frequency of emitted gravitational waves is about 1 kHz. The simulations were performed using a fixed uniform grid and assuming reflection symmetry with respect to the equatorial plane (here, the equatorial plane is chosen to be TABLE II. A list of several quantities for quasiequilibrium initial data. The ADM mass of each star when they are in isolation, the maximum density for each star, the baryon rest-mass ratio QM M2 =M1 , the total baryon rest-mass, the total ADM mass M0 , nondimensional spin parameter q0 J0 =M02 , orbital period P0 , the orbital compactness [C0 M0 *2=3 ], the ratio of the total baryon rest-mass to the maximum allowed mass for a spherical and cold star (Q M =Msph max ), and the frequency of gravitational waves f0 2=P0 . The density, mass, period, and wave frequency are shown in units of 1014 g=cm3 , M , ms, and kHz, respectively. Model Each ADM mass max QM M M0 q0 P0 C0 Q f0 SLy1212 SLy1313 SLy135135 SLy1414 SLy125135 SLy135145 1.20, 1.30, 1.35, 1.40, 1.25, 1.35, 1.20 1.30 1.35 1.40 1.35 1.45 8.03, 8.57, 8.86, 9.16, 8.29, 8.85, 8.03 8.57 8.86 9.16 8.86 9.48 1.00 1.00 1.00 1.00 0.9179 0.9226 2.605 2.847 2.969 3.093 2.847 3.094 2.373 2.568 2.666 2.763 2.568 2.763 0.946 0.922 0.913 0.902 0.921 0.901 2.218 2.110 2.083 2.012 2.110 2.013 0.103 0.112 0.116 0.122 0.112 0.122 1.075 1.175 1.225 1.277 1.175 1.277 0.902 0.948 0.960 0.994 0.948 0.994 FPS1212 FPS125125 FPS1313 FPS1414 1.20, 1.25, 1.30, 1.40, 1.20 1.25 1.30 1.40 9.93, 10.34, 10.79, 11.76, 9.93 10.34 10.79 11.76 1.00 1.00 1.00 1.00 2.624 2.746 2.869 3.120 2.371 2.469 2.566 2.760 0.925 0.914 0.903 0.882 1.980 1.935 1.882 1.750 0.111 0.116 0.121 0.134 1.251 1.309 1.368 1.487 1.010 1.034 1.063 1.143 084021-7 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ the orbital plane). The detailed simulations were performed with the SLy equation of state. In this equation of state, the used grid size is (633, 633, 317) or (377, 377, 189) for x; y; z. In the FPS equation of state, simulations were performed with the (377, 377, 189) grid size to save the computational time. The grid covers the region L x L, L y L, and 0 z L where L is a constant. The grid spacing is determined from the condition that the major diameter of each star is covered with about 50 grid points initially. We have shown that with this grid spacing, a convergent numerical result is obtained [12]. The circumferential radius of spherical neutron stars with the SLy and FPS equations of state is about 11.6 km and 10.7 km for M 1:4M , respectively (see Fig. 2). Thus, the grid spacing is 0:4 km. Accuracy in the computation of gravitational waveforms and the radiation reaction depends on the location of the outer boundaries if the wavelength, 5, is larger than L [12]. For L & 0:45, the amplitude and the radiation reaction of gravitational waves are significantly overestimated [12,42]. Because of the restriction of the computational power, it is difficult to take a huge grid size in which L is much larger than 5. As a consequence of the present restricted computational resources, L has to be chosen as 0:450 where 50 denotes 5 of the l m 2 mode at t 0. Hence, the error associated with the small value of L is inevitable, and thus, the amplitude and radiation reaction of gravitational waves are overestimated in the early phase of the simulation. However, the typical wavelength of gravitational waves becomes shorter and shorter in the late inspiral phase, and hence, the accuracy of the wave PHYSICAL REVIEW D 71, 084021 (2005) extraction is improved with the evolution of the system. This point will be confirmed in Sec. III. The wavelength of quasiperiodic gravitational waves emitted from the formed hypermassive neutron star (denoted by 5merger ) is much shorter than 50 and as large as L (see Table III), so that the waveforms in the merger stage are computed accurately (within 10% error) in the case of neutron star formation irrespective of the grid size. We performed simulations for models SLy1313, SLy125135, SLy1414, and SLy135145 with the two grid sizes, and confirmed that this is indeed the case. We demonstrate this fact in Sec. III by comparing the results for models SLy1313a and SLy1313b. From the numerical results for four models, we have also confirmed that the outcome in the merger does not depend on the grid size. Thus, when we are interested in the outcome or in gravitational waves emitted by the hypermassive neutron stars, simulations may be performed in a small grid size such as (377, 377, 189). With the (633, 633, 317) grid size, about 240 GBytes computational memory is required. For the case of the hypermassive neutron star formation, the simulations are performed for about 30 000 time steps (until t 10 ms) and then stopped to save the computational time. The computational time for one model in such a simulation is about 180 CPU hours using 32 processors on FACOM VPP5000 in the data processing center of National Astronomical Observatory of Japan (NAOJ). For the case of the black hole formation, the simulations crash soon after the formation of apparent horizon because of the socalled grid stretching around the black hole formation TABLE III. A list of computational setting. th , L, &, 50 , and fmerger denote the adiabatic index for the thermal part, the location of outer boundaries along each axis, the grid spacing, the wave length of gravitational waves at t 0, and the frequency of gravitational waves from the formed hypermassive neutron stars, respectively. 5merger denotes the wave length of gravitational waves from the formed hypermassive neutron stars 5merger c=fmerger . The length and the frequency are shown in units of km and kHz. In the last column, the outcome is shown. NS implies that a hypermassive neutron star is produced and remains stable at t 10 ms. NS ! BH implies that a hypermassive neutron star is formed first, but as a result of gravitational radiation reaction, it collapses to a black hole in t & 10 ms. BH implies that a black hole is promptly formed. Grid number Model th SLy1212b SLy1313a SLy1313b SLy1313c SLy1313d SLy135135b SLy1414a SLy125135a SLy135145a 2 2 2 1.3 1.65 2 2 2 2 (377, (633, (377, (377, (377, (377, (633, (633, (633, 377, 633, 377, 377, 377, 377, 633, 633, 633, FPS1212b FPS125125b FPS1313b FPS1414b 2 2 2 2 (377, (377, (377, (377, 377, 377, 377, 377, L & 50 fmerger 5merger Product 189) 317) 189) 189) 189) 189) 317) 317) 317) 77.8 130.8 77.8 77.8 77.8 77.8 130.8 130.8 130.8 0.414 0.414 0.414 0.414 0.414 0.414 0.414 0.414 0.414 333 316 316 316 316 316 302 316 302 3.1 3.2 3.2 3.7 3.4 3.6 3.2 97 94 94 81 88 83 94 NS NS NS NS ! BH NS NS ! BH BH NS BH 189) 189) 189) 189) 69.5 69.5 69.5 69.5 0.370 0.370 0.370 0.370 297 297 282 262 3.5 86 NS ! BH NS ! BH BH BH 084021-8 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) region. In this case, the computational time is about 60 CPU hours for about 10 000 time steps. B. Characteristics of the merger 1. General feature In Figs. 3–5, we display the snapshots of the density contour curves and the velocity vectors in the equatorial plane at selected time steps for models SLy1414a, SLy1313a, and SLy125135a, respectively. Figure 6 displays the density contour curves and the velocity vectors in the y 0 plane at a late time for SLy1414a and SLy1313a. Figs. 3 and 6(a) indicate typical evolution of the density contour curves in the case of prompt black hole formation. On the other hand, Figs. 4, 5, and 6(b) show those in the formation of hypermassive neutron stars. Figure 7 displays the evolution of the maximum value of (hereafter max ) and the central value of (hereafter c ) for models SLy1414a, SLy135145a, SLy135135b, SLy1313a, SLy125135a, and SLy1212b [Fig. 7(a)] and for models FPS1414b, FPS1313b, FPS125125b, and FPS1212b [Fig. 7(b)]. This shows that in the prompt black hole formation, c monotonically decreases toward zero. On the other hand, c and max settle down to certain values in the hypermassive neutron star formation. For model SLy135135b, a hypermassive neutron star is formed first, but after several quasiradial oscillations of high amplitude, it collapses to a black hole due to dissipation of the angular momentum by gravitational radiation. The large oscillation amplitude results from the fact that the selfgravity is large enough to deeply shrink surmounting the centrifugal force. These indicate that the total ADM mass of this model (M 2:7M ) is only slightly smaller than the threshold value for the prompt black hole formation. The quasiradial oscillation of the large amplitude induces a characteristic feature in gravitational waveforms and the Fourier spectrum (cf. Sec. III C). In the case of black hole formation (models SLy1414, SLy135145, SLy135135, FPS1414, FPS1313, FPS125125, and FPS1212), the computation crashed soon after the formation of apparent horizons since the region around the apparent horizon of the formed black hole was stretched significantly and the grid resolution became too poor to resolve such region. On the other hand, we stopped the simulations for other cases to save the computational time, after the evolution of the formed massive neutron stars was followed for a sufficiently long time. At the termination of these simulations, the averaged violation of the Hamiltonian constraint H remains of order 0.01 (cf. Figure 18). We expect that the simulations could be continued for a much longer time than 10 ms if we could have sufficient computational time. In every model, the binary orbit is stable at t 0 and the orbital separation gradually decreases due to the radiation reaction of gravitational waves for which the emission time scale is longer than the orbital period. If the orbital separation becomes sufficiently small, each star is elongated by tidal effects. As a result, the attraction force due to the tidal interaction between two stars becomes strong enough to make the orbit unstable to merger. The merger starts after about one orbit at t 2 ms irrespective of models. Since the orbital separation at t 0 is very close to that for a FIG. 3 (color online). Snapshots of the density contour curves for in the equatorial plane for model SLy1414a. The solid contour curves are drawn for 2 1014 i g=cm3 i 2 10 and for 2 1014 100:5i g=cm3 i 1 7. The dotted curves denote 2 1014 g=cm3 . The number in the upper left-hand side denotes the elapsed time from the beginning of the simulation in units of ms. The initial orbital period in this case is 2.012 ms. Vectors indicate the local velocity field vx ; vy , and the scale is shown in the upper right-hand corner. The thick circle in the last panel of radius r 2 km denotes the location of the apparent horizon. 084021-9 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ FIG. 4. FIG. 5. PHYSICAL REVIEW D 71, 084021 (2005) The same as Fig. 3 but for model SLy1313a. The initial orbital period is 2.110 ms in this case. The same as Fig. 3 but for model SLy125135a. The initial orbital period is 2.110 ms in this case. 084021-10 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) (a) (b) FIG. 6 (color online). Snapshots of the density contour curves for and the local velocity field vx ; vz in the y 0 plane (a) at t 2:991 ms for model SLy1414a and (b) at t 8:621 ms for model SLy1313a. The method for drawing the contour curves and the velocity vectors is the same as that in Fig. 3. marginally stable orbit, a small decrease of the angular momentum and energy is sufficient to induce the merger in the present simulations. If the total mass of the system is high enough, a black hole is directly formed within about 1 ms after the merger sets in. On the other hand, for models with mass smaller than a threshold mass Mthr , a hypermassive neutron star is formed at least temporarily. The hypermassive neutron star is stable against gravitational collapse for a while after its formation, but it will collapse to a black hole eventually due to radiation reaction of gravitational waves or due to outward angular momentum transfer (see discussion later). In the formation of the hypermassive neutron stars, a double-core structure is first formed, and then, it relaxes to a highly nonaxisymmetric ellipsoid [cf. Figs. 4, 5, and 6(b)]. The contour plots drawn for a high-density region (a) (b) FIG. 7 (color online). Evolution of the maximum values of and the central value of (a) for models SLy1414a (long-dashed curves), SLy135145a (dotted long-dashed curves which approximately coincide with long-dashed curve for c ), SLy135135b (dotted curves), SLy1313a (solid curves), SLy125135a (dotted-dashed curves), and SLy1212b (dashed curves), and (b) for models FPS1414b (long-dashed curves), FPS1313b (solid curves), FPS125125 (dotted-dashed curves), and FPS1212b (dashed curves). The dotted horizontal lines denote the central values of the lapse and density of the marginally stable and spherical star in equilibrium for given cold equations of state. The reason that the maximum density decreases in the final stage for the black hole formation case is as follows: We choose as a fundamental variable to be evolved and compute from =ut e6 . In the final stage, is very large ( > 1) and, hence, a small error in results in a large error in . Note that the maximum value of increases monotonically by many orders of magnitude. 084021-11 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ 14 3 with 4 10 g=cm show that the axial ratio of the bar measured in the equatorial plane is 0:5; the axial lengths of the semi major and minor axes are 20 and 10 km, respectively. Figure 6(b) also shows that the axial length along the z axis is about 10 km. Namely, a highly elliptical rotating ellipsoid is formed. This outcome is significantly different from the previous ones found with the 2 equation of state [12], in which nearly axisymmetric spheroidal neutron stars are formed. The reason is that the adiabatic index of the realistic equations of state adopted in this paper is much larger than 2 that is adopted in the previous one. According to a Newtonian study [43], a uniformly rotating ellipsoid (Jacobi-like ellipsoid) exists only for * 2:25. This fact suggests that rapidly rotating stars with a large adiabatic index are only subject to the ellipsoidal deformation. Note that similar results have been already reported in Newtonian and post-Newtonian simulations [23–25]. The rotating hypermassive neutron stars also oscillate in a quasiradial manner (cf. Figure 7). Such oscillation is induced by the approaching velocity at the collision of two stars. By the radial motion, shocks are formed at the outer region of the hypermassive neutron stars to convert the kinetic energy to the thermal energy. The shocks are also generated when the spiral arms hit the oscillating hypermassive neutron stars. These shocks heat up the outer region of the hypermassive neutron stars for many times, and as a result, the thermal energy of the envelope increases fairly uniformly. The further detail of these thermal properties is discussed in Sec. III B 3. Since the degree of the nonaxial symmetry is sufficiently large, the hypermassive neutron star found in this paper is a stronger emitter of gravitational waves than that found in [12]. The significant radiation decreases the angular momentum of the hypermassive neutron stars. The nonaxisymmetric structure also induces the angular momentum transfer from the inner region to the outer one due to the hydrodynamic interaction. As a result of these effects, the rotational angular velocity * v’ decreases and its profile is modified. In Fig. 8, we show the evolution of * of the hypermassive neutron stars along x and y axes at t 3:897, 6.069, and 8.621 ms for models SLy1313a and SLy125135a. At its formation, the hypermassive neutron stars are strongly differentially and rapidly rotating. The strong differential rotation yields the strong centrifugal force, which plays an important role for sustaining the large selfgravity of the hypermassive neutron stars [22,44]. Since the angular momentum is dissipated by the gravitational radiation and redistributed by the hydrodynamic interaction, * decreases significantly in the central region, and hence, the steepness of the differential rotation near the center decreases with time. This effect eventually induces the collapse to a black hole. It should be also noted that * along two axes is significantly different near the origin. The reason at t 3:897 ms PHYSICAL REVIEW D 71, 084021 (2005) is that the formed hypermassive neutron stars have a double-core structure (cf. Figures 4 and 5) and the angular velocity of the cores is larger than the low density region surrounding them. The reason for t * 6 ms is that the hypermassive neutron star is not a spheroid but an ellipsoid of high ellipticity and the angular velocity depends strongly on the coordinate ’. Figures 4 and 5 show that even after the emission of gravitational waves for 10 ms, the hypermassive neutron star is still highly nonaxisymmetric. This indicates that gravitational waves will be emitted for much longer time scale than 10 ms. Thus, the rotational kinetic energy and the angular momentum will be subsequently dissipated by a large factor, eventually inducing the collapse to a black hole. We here estimate the lifetime of the hypermassive neutron stars using Fig. 7 which shows that the value of c decreases gradually with time. It is reasonable to expect that the collapse to a black hole sets in when the value of c becomes smaller than a critical value. Since the angular momentum should be sufficiently dissipated before the collapse sets in, the threshold value of c for the onset of the collapse will be approximately equal to that of marginally stable spherical stars (i.e., the dotted horizontal line in Fig. 7). One should keep in mind that the threshold value depends on the slicing condition, and thus, this method can work only when the same slicing is used for computation of the spherical star and for simulation. In this paper, the (approximate) maximal slicing is adopted both in the simulation and in computation of spherical equilibria so that this method can be used. The results for models SLy135135b and FPS1212b indeed illustrate that the prediction by this method is appropriate. For models SLy1313a, SLy125135a and SLy1212b, the decrease rate of the value of c estimated from the data for 5 ms & t & 10 ms is 0:005 ms1 . Extrapolating this result suggests that the hypermassive neutron stars will collapse to a black hole at t 30 ms for models SLy1313a and SLy125135a and at t 50 ms for model SLy1212b. These time scales are much shorter than the dissipation time scale by viscosities or the redistribution time scale of the angular momentum by the magnetic field [22]. Therefore, the gravitational radiation or the outward angular momentum transfer by the hydrodynamic interaction plays the most important role in the evolution of the hypermassive neutron stars. In the prompt formation of a black hole, most of the fluid elements are swallowed into the black hole in 1 ms after the merger sets in. Thus, the final outcome is a system of a rotating black hole and a surrounding disk of very small mass [cf. Figure 6(a)]. In Fig. 9, we plot the evolution of the total baryon rest-mass located outside a radius r, M r, for models SLy1414a and SLy135145a. M r is defined by 084021-12 M r Z rr0 rmax d3 x0 ; (39) MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) (a) (b) FIG. 8. Evolution of the angular velocity * along x (solid curves) and y (dashed curves) axes (a) for models SLy1313a and (b) SLy125135a. The time is shown in the upper right corner of each panel in units of ms. * along x and y axes is computed by vy =x and vx =y, and hence, it is not determined at the origin. Since the formed hypermassive neutron star is not spheroidal, * along two axes is significantly different. where rmax is introduced to exclude the contribution from the small-density atmosphere. In the present work we choose as rmax 7:5M0 31 km within which the integrated mass of the atmosphere is negligible ( < 0:01% of the total mass). The results are plotted for r=M0 1:5, 3, and 4.5. Note that the apparent horizon is located for r 0:5M0 at t 3:0 ms for models SLy1414a and SLy135145a, and inside the horizon about 99% and 98% of the initial mass are enclosed for these cases, respectively. Figure 9 indicates that the fluid elements still continue to fall into the black hole at the end of the simulation. This suggests that the final disk mass will be smaller than 1% of the total baryon rest-mass. In [12], we found that even the small mass difference with QM 0:9 increases the fraction of disk mass around the black hole significantly. However, in the present equations of state, QM 0:9 is not small enough to significantly increase the disk mass. This results from the difference in the equations of state. The detailed reason is discussed in Sec. III B 4. The area of the apparent horizons A is determined in the black hole formation cases. We find that A 0:85; 16%M02 2 fined by q JBH =MBH , where JBH and MBH are the angular momentum and the mass of the black hole, are computed from A 1 1 1 q2 1=2 : 2 16%MBH 2 (41) 2 0:85, q 0:7. Equation (41) implies that for A=16%MBH (40) for models SLy1414a and SLy135145a. Since most of the fluid elements are swallowed into the black hole and also the energy carried out by gravitational radiation is at most 0:01M0 (cf. Figure 18), the mass of the formed black hole is approximately 0:99M0 . Assuming that the area of the apparent horizon is equal to that of the event horizon, the nondimensional spin parameter of the black hole de- FIG. 9 (color online). Evolution of M r=M for r 1:5, 3, and 4:5M0 for models SLy1414a (curves except for the dotted curve) and SLy135145a (dotted curves). M0 and M denote the ADM mass and total baryon rest-mass of binary neutron stars at t 0. Definition for M r is shown in Eq. (39). 084021-13 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ On the other hand, we can estimate the value of q in the following manner. As shown in Sec. III C, the angular momentum is dissipated by 10%–15% by gravitational radiation, while the ADM mass decreases by 1%. As listed in Table II, the initial value of q is 0:9. Therefore, the value of q in the final stage should be 0:75– 0.8. The values of q derived by two independent methods agree with each other within 10% error. This indicates that the location and the area of the black holes are determined within 10% error. For q 0:7– 0.8 and MBH 2:8M , the frequency of the quasinormal mode for the black hole oscillation is about 6.5–72:8M =MBH kHz [45]. This value is rather high and far out of the best sensitive frequency range of the laser interferometric gravitational wave detectors [5]. Thus, in the following, we do not touch on gravitational waveforms in the prompt black hole formation. 2. Threshold mass for black hole formation The threshold value of the total ADM mass above which a black hole is promptly formed is Mthr 2:7M for the SLy equation of state and Mthr 2:5M for the FPS one with th 2. For the SLy case, we find that the value does not depend on the mass ratio for 0:9 & QM 1. The maximum allowed mass for the stable and spherical neutron stars is 2:04M and 1:80M for the SLy and FPS equations of state, respectively. This implies that if the total mass is by 30%–40% larger than the maximum allowed mass for stable and spherical stars, a black hole is promptly formed. In a previous study with the 2 equation of state [12], we found that threshold mass is by about 70% larger than the maximum allowed mass for stable and spherical neutron stars. Thus, comparing the threshold value of the total ADM mass, we can say that a black hole is more subject to be formed promptly with the realistic equations of state. In [46,47], the maximum mass of differentially rotating stars in axisymmetric equilibrium (hereafter Mmax:dif ) is studied for various equations of state. The authors compare Mmax:dif with the maximum mass of spherical stars (hereafter Mmax:sph ) for given equations of state. They find that the ratio Mmax:dif =Mmax:sph for FPS and APR equations of state (APR is similar to SLy equation of state) is much smaller than that for 2 equation of state. Their study is carried out for axisymmetric rotating stars in equilibrium and with a particular rotational law, and hence, their results cannot be simply compared with our results obtained for dynamical and nonaxisymmetric spacetime. However, their results suggest that the merged object may be more susceptible to collapse to a black hole with the realistic equations of state. This tendency agrees with our conclusion. The compactness in each neutron star of no rotation in isolation is defined by C GMsph =Rsph c2 where Msph and Rsph denote the ADM mass and the circumferential radius PHYSICAL REVIEW D 71, 084021 (2005) of the spherical stars. For the SLy equation of state, C 0:151, 0.165, 0.172, and 0.178 for Msph 1:2, 1.3, 1.35, and 1:4M , respectively. For FPS one, C 0:162, 0.169, 0.177, and 0.192 for Msph 1:2, 1.25, 1.3, and 1:4M , respectively. This indicates that a black hole is promptly formed for C * 0:17 after merger of two (nearly) identical neutron stars. In the 2 equation of state, the threshold value of C is 0:15– 0.16 [12]. Thus, comparing the threshold value of the compactness, we can say that a black hole is less subject to be formed with the realistic equations of state. The reason that the threshold mass for the prompt black hole formation is smaller with the realistic equations of state may be mentioned in the following manner: In the realistic equations of state, the compactness of each neutron star is larger than that with the 2 equation of state for a given mass. Accordingly, for a given total mass, the binary system at the onset of the merger is more compact. This implies that the angular momentum is dissipated more before the merger sets in with the realistic equations of state. In the case of the hypermassive neutron star formation, the centrifugal force plays the most important role for sustaining the large selfgravity. Thus, the large dissipation of the angular momentum before the merger helps the prompt black hole formation. Therefore, a black hole is more subject to be formed in the realistic equations of state. 3. Thermal properties In Fig. 10(a)–10(c), we show profiles of " and "th as well as that of along x and y axes at t 2:991 ms for model SLy1414a, at t 8:621 ms for model SLy1313a, and at t 8:621 ms for model SLy125135a, respectively. The density contour curves at the corresponding time steps are displayed in the last panel of Figs. 3–5. Figure 10(d) shows the evolution of the total internal energy and thermal energy defined by Z U "d3 x; (42) Uth Z "th d3 x; (43) for models SLy1313a, SLy125135a, and SLy1212b. Note that in the absence of shock heating, " should be equal to "cold . Thus, "th denotes the specific thermal energy generated by the shock heating. First, we focus on the thermal property for models SLy1313a and SLy125135a which are the representative models of hypermassive neutron star formation. In these cases, the heating is not very important in the central region. This is reasonable because the shocks generated at the collision of two stars are not very strong, and thus, the central part of the hypermassive neutron stars is formed without experiencing the strong shock heating. On the 084021-14 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) other hand, the shock heating plays an important role in the outer region of the hypermassive neutron stars and in the surrounding accretion disk since the spiral arms hit the hypermassive neutron stars for many times. The typical value of "th is 0:01c2 –0:02c2 . Here, we recover c for making the unit clear. In the following, we assume that the components of the hypermassive neutron stars and surrounding disks are neutrons. Then, the value of "th 0:01c2 implies that the thermal energy per nucleon is "th 9:37 MeV=nucleon: (44) 0:01c2 Since the typical value of "th is 0:01c2 –0:02c2 , the typical thermal energy is 10–20 MeV. This value agrees approximately with that computed in [48,49]. Figure 10(d) shows that the total internal energy and thermal energy are relaxed to be U 0:1M c2 5 1053 erg; (45) Uth 0:025M c2 1 1053 erg; (46) for both models SLy1313a and SLy125135a. Thus, the thermal energy increases up to 25% of the total internal energy. We note that these values are approximately identical between models SLy1313a and SLy125135a. This implies that the mass ratio of QM 0:9 does not significantly modify the thermal properties of the hypermassive neutron stars in the realistic equations of state. The region of "th 0:01c2 will be cooled via the emission of neutrinos [48,49]. According to [48,49], the emis- (b) (a) (c) (d) FIG. 10. Profiles of " (solid curves) and "th " "cold (dashed curves) as well as that of along x and y axes (a) at t 2:991 ms for model SLy1414a, (b) at t 8:621 ms for model SLy1313a, and (c) at t 8:621 ms for model SLy125135a. Note that the region of r & 2 km for panel (a) is inside the apparent horizon [see Figs. 3 and 6(a)]. (d) Evolution of U and Uth in unit of the rest-mass energy M c2 for models SLy1313a (solid curves), SLy125135a (dotted curves), and SLy1212b (dashed curves). 084021-15 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ sion rate in the hypermassive neutron star with the averaged value of "th 10–30 MeV is 1052 –1053 erg=s. Thus, if all the amounts of the thermal energy are assumed to be dissipated by the neutrino cooling, the time scale for the emission of the neutrinos will be 1–10 s. This is much longer than the lifetime of the hypermassive neutron stars & 100 ms. Therefore, the cooling does not play an important role in their evolution. Since the lifetime of the hypermassive neutron stars & 100 ms is nearly equal to the time duration of the short gamma-ray bursts [50], it is interesting to ask if they could generate the typical energy of the bursts. In a model for central engines of the gamma-ray bursts, a fireball of the electron-positron pair and photon is produced by the pair annihilation of the neutrino and antineutrino [49,50]. In [49], Janka and Ruffert estimate the efficiency of the annihilation as several 103 for the neutrino luminosity 1052 erg=s, the mean energy of neutrino 10 MeV, and the radius of the hypermassive neutron star 10 km (see Eq. (1) of [49]). This suggests that the energy generation rate of the electron-positron pair is 1050 erg=s. Since the lifetime of the hypermassive neutron stars is & 100 ms, the energy available for the fireball will be at most 1049 erg. This value is not large enough to explain typical cosmological gamma-ray bursts. Furthermore, as Janka and Ruffert found [49], the pair annihilation of the neutrino and antineutrino is most efficient in a region near the hypermassive neutron star, for which the baryon density is large enough [cf. Figure 6(b)] to convert the energy of the fireball to the kinetic energy of the baryon. Therefore, it is not very likely that the hypermassive neutron stars are the central engines of the typical short gamma-ray bursts. Now, we focus on model SLy1414a in which a black hole is promptly formed after the merger. Comparing Fig. 10(a) with the last panel of Fig. 3, the region of high thermal energy is located along the spiral arms of the accretion disk surrounding the central black hole. (Note that the region of r & 2 km is inside the apparent horizon, and hence, we do not consider such region.) The part of the matter in the spiral arms with small orbital radius r & 5 km is likely to be inside the radius of an innermost stable circular orbit around the black hole, and hence, be swallowed into the black hole. Otherwise, the matter in the spiral arms will form an accretion disk surrounding the black hole. Thus, eventually a hot accretion disk will be formed. However, the region of high thermal energy for r * 10 km is of low density with & 1012 g=cm3 , and the total mass of the disk will be 103 M & M & 0:01M (see Fig. 9). The total thermal energy of the accretion disk is estimated as Mdisk " th Uth Mdisk " th 1:8 1050 erg: (47) 0:01M 0:01c2 where Mdisk and " th denote the mass of the accretion disk and the averaged value of the specific thermal energy. Hence, even if all the amounts of the thermal energy are PHYSICAL REVIEW D 71, 084021 (2005) dissipated by the emission of neutrinos, the total energy of the radiated neutrinos will be at most several 1050 erg. According to [49], the efficiency of the annihilation of the neutrino and antineutrino is several 105 for the neutrino luminosity 1050 erg=s, the mean energy of neutrino 10 MeV, and the disk radius 10 km. This indicates that the energy of the fireball is at most 1046 erg. Although the density of the baryon at the region that the pair annihilation is likely to happen is small enough to avoid the baryon loading problem, this energy is too small to explain cosmological gamma-ray bursts. 4. Effects of mass difference Comparing the evolution of the contour curves, the maximum density, and the central value of the lapse function for models SLy1313a and SLy125135a [see Figs. 4, 5, and 7(a)], it is found that the mass difference plays a minor role in the formation of a hypermassive neutron star as far as QM is in the range between 0.9 and 1. Figs. 7(a) and 9 also illustrate that the evolution of the system to a black hole is very similar for models SLy1414a and SLy135145a. In [12] in which simulations were performed using the 2 equation of state, we found that the mass difference with QM 0:9 significantly induces an asymmetry in the merger which contributes to formation of large spiral arms and the outward angular momentum transfer, which are not very outstanding in the present results. The reason seems to be as follows. In the previous equation of state, the mass difference with QM 0:9 results in a relatively large ( 15%) difference of the compactness between two stars. On the other hand, the difference in the compactness between two stars with the present equations of state is 10% for QM 0:9. This is due to the fact that the stellar radius depends weakly on the mass in the range 0:8M & M & 1:5M (see Fig. 2). As a result of the smaller difference in the compactness, the tidal effect from the more massive star to the companion becomes smaller, and therefore, the asymmetry is suppressed. To yield a system of a black hole and a massive disk, smaller mass ratio with QM < 0:9 will be necessary in the realistic equations of state. Another possible reason is that neutron stars in the realistic equations of state are more compact than those in the 2 equation of state. Because of this fact, at the merger, the system is more compact, and hence, even in the formation of the asymmetric spiral arms, they cannot spread outward extensively but wind around the formed neutron star quickly. Consequently, the mass of the disk around the central object is suppressed to be small and also the asymmetric density configuration does not become very outstanding. 5. Dependence of dynamics on the grid size and th For model SLy1313, we performed simulations changing the value of th with the grid size (377, 377, 189). 084021-16 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) In Fig. 11, evolution of c and max is shown for models SLy1313b–SLy1313d. Note that the grid size and grid spacing are identical for these models. The results for model SLy1313a are shown together for comparison with those for model SLy1313b for which the parameters are identical but for the grid size. By comparing the results for models SLy1313a and SLy1313b, the magnitude of the error associated with the small size of L is investigated. Figure 11 shows that two results are approximately identical besides a systematic phase shift of the oscillation. This shift is caused by the inappropriate computation of the radiation reaction in the late inspiral stage for t & 2 ms: For model SLy1313b, L is smaller, and hence, the radiation reaction in the inspiral stage is significantly overestimated to spuriously accelerate the merger resulting in the phase shift. However, besides the phase shift, the results are approximately identical. In particular, the results agree well in the merging phase. This indicates that even with the smaller grid size (377, 377, 189), the formation and evolution of the hypermassive neutron star can be followed within a small error. Comparison of the results among models SLy1313b– SLy1313d tells that for the smaller value of th , the maximum density (central lapse) of a hypermassive neutron star formed during the merger is larger (smaller). This is due to the fact that the strength of the shock formed at the collision of two stars, which provides the thermal energy in the outer region of the formed hypermassive neutron stars to expand, is proportional to the value of th 1. Figure 11 also indicates that the evolution of the system does not FIG. 11 (color online). The same as Fig. 7 but for models SLy1313a (solid curves), SLy1313b (dashed curves), SLy1313c (dotted-dashed curves), and SLy1313d (long-dashed curves). Note that the simulation for SLy1313b was stopped at t 7:5 ms to save computational time. depend strongly on the value of th for th * 1:65. However, for th 1:3, the formed hypermassive neutron star is very compact at its birth, and hence, collapses to a black hole in a short time scale (at t 8 ms) after the angular momentum is dissipated by gravitational radiation. This time scale for black hole formation is much shorter than that for models SLy1313b and SLy1313d for which it would be 30 ms. This indicates that for small values of th 1 & 0:3, the collapse to a black hole is significantly enhanced. In reality, in the presence of cooling processes such as neutrino cooling, the adiabatic index th decreases effectively. Such cooling mechanism may accelerate the formation of a black hole. However, the emission time scale of the neutrino is 1–10 s as mentioned in Sec. III B 3. Thus, the effect does not seem to be very strong. C. Gravitational waveforms 1. Waveforms in the formation of hypermassive neutron stars In Figs. 12 –15, we present gravitational waveforms in the formation of hypermassive neutron stars for several models. Figure 12(a) shows R and R for model SLy1313a. For comparison, gravitational waves computed in terms of a quadrupole formula (A and A ) defined in Sec. II B are shown together by the dashed curves. The amplitude of gravitational waves, h, observed at a distance of r along the optimistic direction (, 0) is written as h 1022 R; 100 Mpc : r 0:31 km (48) Thus, the maximum amplitude observed along the most optimistic direction is 2 1022 at a distance of 100 Mpc. In the real data analysis of gravitational waves, a matched filtering technique [5] is employed. In this method, the signal of the identical frequency can be accumulated using appropriate templates, and as a result, the effective amplitude increases by a factor of N 1=2 where N denotes the number of the cycle of gravitational waves for a given frequency. We determine such effective amplitude in Sec. III C 3 (cf. Equation (50)). The waveforms shown in Fig. 12(a) are typical ones in the formation of a hypermassive neutron star. In the early phase (tret & 2 ms), gravitational waves associated with the inspiral motion are emitted, while for tret * 2 ms, those by the rotating and oscillating hypermassive neutron star are emitted. In the following, we focus only on the waveforms for tret * 2 ms. Gravitational waves from the hypermassive neutron stars are characterized by quasiperiodic waves for which the amplitude and the frequency decrease slowly. The amplitude decreases with the ellipticity, which is decreased by the effects that the angular momentum decreases due to 084021-17 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ PHYSICAL REVIEW D 71, 084021 (2005) (a) (b) FIG. 12 (color online). (a) Gravitational waves for model SLy1313a. R and R (solid curves) and A and A (dashed curves) as functions of the retarded time are shown. (b) R and R as functions of the retarded time for model SLy125135a (solid curves). For comparison, those for SLy1313a are shown by the dashed curves. the radiation reaction and is transferred from the inner region to the outer one by the hydrodynamic interaction associated with the nonaxisymmetric structure. However, the time scale for the decrease appears to be much longer than 10 ms as illustrated in Figs. 12 –15. The oscillation frequency varies even more slowly. The reason seems to be FIG. 13 (color online). R and R for models SLy1313a (dashed curves) and SLy1313b (solid curves) as functions of the retarded time. It is found that gravitational waveforms in the merger stage agree well. Note that the simulation for SLy1313b was stopped at t 7:5 ms to save computational time. that the following two effects approximately cancel each other: (i) with the decrease of the angular momentum of the hypermassive neutron stars due to the radiation reaction as well as the angular momentum transfer by the hydrodynamic interaction with outer envelope, the characteristic frequency of the figure rotation decreases, while (ii) with the decrease of the angular momentum, the centrifugal force is weakened to reduce the characteristic radius for a spin-up. (We note that the radiation reaction alone may FIG. 14 (color online). R and R for models SLy135135b (solid curves) and SLy1313a (dashed curves) as functions of the retarded time. 084021-18 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) (b) (a) FIG. 15 (color online). R and R (a) for models SLy1313b (dotted curves), SLy1313c (solid curves), and SLy1313d (dashed curves), and (b) for models SLy1212b (solid curves) and FPS1212b (dashed curves) as functions of the retarded time. increase the frequency of the figure rotation [51]. In the hypermassive neutron stars formed after the merger, the angular momentum transfer due to the hydrodynamic interaction is likely to play an important role for the decrease of the frequency.) In gravitational waveforms computed in terms of the quadrupole formula (the dashed curves in Fig. 12), the amplitude is systematically underestimated by a factor of 30%– 40%. This value is nearly equal to the magnitude of the compactness of the hypermassive neutron star, GM=Rc2 , where M and R denote the characteristic mass and radius. Since the quadrupole formula is derived ignoring the terms of order GM=Rc2 , this magnitude for the error is quite reasonable. In simulations with Newtonian, post-Newtonian, and approximately relativistic frameworks, gravitational waves are computed in the quadrupole formula (e.g., [16,24,25]). The results here indicate that the amplitudes for quasiperiodic gravitational waves from hypermassive neutron stars presented in those simulations are significantly underestimated.4 Although the error in the amplitude is large, the wave phase is computed accurately except for a slight systematic phase shift. From the point of view of the data analysis, the wave phase is the most important information on gravitational waves. This suggests that a quadrupole formula may be a useful tool for computing approximate gravitational waves. We note that 4 Besides systematic underestimation of the wave amplitude, rather quick damping of quasiperiodic gravitational waves is seen in the results of these references. This quick damping also seems to be due to a systematic error, which is likely to result from a relatively dissipative numerical method (smoothed particle hydrodynamics method) used in these references. these features have been already found for oscillating neutron stars [38] and for nonaxisymmetric protoneutron stars formed after stellar core collapse [37]. Here, we reconfirm that the same feature holds for the merger of binary neutron stars. In Fig. 12(b), we display gravitational waveforms for model SLy125135a. For comparison, those for SLy1313a are shown together (dashed curves). It is found that two waveforms coincide each other very well. As mentioned in Sec. III B 4, the mass difference with QM 0:9 does not induce any outstanding change in the merger dynamics. This fact is also reflected in the gravitational waveforms. In Fig. 13, we compare gravitational waveforms for models SLy1313b (solid curves) and SLy1313a (dashed curves). For SLy1313b, the simulation was performed with a smaller grid size and gravitational waves were extracted in a near zone with L=50 0:25 and L=5merger 0:83 (cf. for model SLy1313a, L=50 0:41 and L=5merger 1:39). This implies that the waveforms for model SLy1313b are less accurately computed than those for SLy1313a. Indeed, the wave amplitude for tret & 2 ms is badly overestimated. However, the waveforms from the formed hypermassive neutron stars for two models agree very well except for a systematic phase shift, which is caused by the overestimation for the radiation reaction in the early phase (tret & 2 ms). Thus, for computation of gravitational waves from the hypermassive neutron stars, we may choose the small grid size. Making use of this fact, we compare gravitational waveforms among several models computed with the small grid size in the following. In Fig. 14, we compare gravitational waves from the hypermassive neutron stars for models SLy1313a (dashed 084021-19 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ curves) and SLy135135b (solid curves). As shown in Figs. 12 and 13, quasiperiodic waves for which the frequency is approximately constant are emitted for model SLy1313a. On the other hand, the frequency is not constant but modulates with time for model SLy135135b (e.g., see the waveforms at tret 3:6, 4.6, 5.6, and 6.6 ms for which the wavelength is relatively short). The reason is that the formed hypermassive neutron star quasiradially oscillates with a large amplitude and the frequency of gravitational waves varies with the change of the characteristic radius. Because of this, the Fourier spectra for models SLy1313 and SLy135135 are significantly different although the difference of the total mass is very small [cf. Figure 17(b)]. In Fig. 15(a), we compare gravitational waveforms for models SLy1313b (dotted curves), SLy1313c (solid curves), and SLy1313d (dashed curves). For these models, the cold part of the equation of state is identical but the value of th is different. As mentioned in Sec. III B 5, with the smaller values of th , the shock heating is less efficient, and as a result, the formed hypermassive neutron star becomes more compact. Since the characteristic radius decreases, the amplitude of gravitational waves decreases and the frequency increases. This shows that the strength of the shock heating affects the amplitude and the characteristic frequency of gravitational waves. In Fig. 15(b), we compare gravitational waveforms for models SLy1212b (solid curves) and FPS1212b (dashed curves). For these models, the equations of state are different, but the total ADM mass is approximately identical. Since the FPS equation of state is slightly softer than the SLy one, the compactness of each neutron star is larger by a factor of 5%–10% (cf. Figure 1) and so is for the formed hypermassive neutron star. As a result, the frequency of gravitational waves for the FPS equation of state is slightly ( 15%) higher [cf. Figure 17(d)]. On the other hand, the amplitude of gravitational waves is not very different. This is due to the fact that with increasing the compactness, the radius of the hypermassive neutron star decreases while the angular velocity increases. These two effects approximately cancel each other, and as a result, dependence of the amplitude is not remarkable between two models. 2. Emission rate of the energy and the angular momentum In Fig. 16(a), the emission rates of the energy and the angular momentum by gravitational radiation are shown for models SLy1313a (solid curves) and SLy125135a (dashed curves). In the inspiral phase for tret & 2 ms, they increase with time since the amplitude and the frequency of the chirp signal increase. After the peak is reached, the emission rates quickly decrease by about 1 order of magnitude since the merged object becomes a fairly axisymmetric transient object. However, because of its large angular momentum, the formed hypermassive neutron star soon changes to a highly ellipsoidal object PHYSICAL REVIEW D 71, 084021 (2005) which emits gravitational waves significantly. The luminosity from the ellipsoidal neutron star is as high as the first peak at tret 2:2 ms. This is in contrast with the results obtained with the 2 equation of state in which the magnitude of the second peak is 30%–50% as large as that of the first peak [9]. This reflects the fact that the degree of the ellipticity of the formed hypermassive neutron star is much higher than that found in [9] because of the large adiabatic index for the realistic equations of state. The emission rates of the energy and the angular momentum via gravitational waves gradually decrease with time, since the degree of the nonaxial symmetry decreases. However, the decrease rates are not very large and the emission rates at tret 10 ms remain to be as high as that in the late inspiral phase as dE=dt 7 1054 erg=s and dJ=dt 7 1050 g cm2 =s2 . The angular momentum at t 10 ms is J 0:7J0 4 1049 g cm2 =s. Assuming that the emission rate of the angular momentum does not change and remains 7 1050 g cm2 =s, the emission time scale is evaluated as J=dJ=dt 50 ms. For more accurate estimation, we should compute J Jmin =dJ=dt where Jmin denotes the minimum allowed angular momentum for sustaining the hypermassive neutron star. Since Jmin is not clear, we set Jmin 0. Thus, the estimated value presented here is an approximate upper limit for the emission time scale (see discussion below), and hence, the hypermassive neutron star will collapse to a black hole within 50 ms. This estimate agrees with the value 30 ms obtained in terms of the change rate of c (cf. Sec. III B 1). Therefore, we conclude that the lifetime of the hypermassive neutron stars and hence the time duration of the emission of quasiperiodic gravitational waves are as short as 30–50 ms for models SLy1313a and SLy125135a. Figure 16(b) displays the emission rates of the energy and the angular momentum for models SLy1212b (solid curves) and FPS1212b (dashed curves). For these models, the value of L is not large enough to accurately compute gravitational waves in the inspiral phase for tret & 2 ms. Thus, we only present the results for the merger phase. The emission rates for SLy1212b are slightly smaller than those for model SLy1313a. This results from the fact that the total mass of the system is smaller. On the other hand, the emission rates for FPS1212b are slightly larger than those for SLy1212b. The reason is that the FPS equation of state is softer than the SLy one, and as a result, the formed hypermassive neutron star is more compact and the rotational angular velocity is larger. The hypermassive neutron star formed for FPS1212b collapses to a black hole at t 10 ms. This is induced by the emission of the angular momentum by gravitational waves. However, the collapse time is shorter than the emission time scale evaluated by J=dJ=dt. This is reasonable because the hypermassive neutron star formed for model FPS1212 is close to the marginally stable configuration, and hence, a small amount 084021-20 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) (b) (a) FIG. 16 (color online). dE=dt and dJ=dt of gravitational radiation (a) for models SLy1313a (solid curves) and SLy125135a (dashed curves), and (b) for models SLy1212b (solid curves) and FPS1212b (dashed curves). of the dissipation leads to the collapse. This illustrates that the time scale J=dJ=dt should be regarded as the approximate upper limit for the collapse time scale. 3. Fourier power spectrum To determine the characteristic frequency of gravitational waves, we carried out the Fourier analysis. In Fig. 17, the power spectrum dE=df is presented (a) for models SLy1313a and SLy125135a, (b) for SLy1313a and SLy135135b, (c) for SLy1313a, SLy1313c, SLy1313d, and (d) for SLy1212b and FPS1212b, respectively. Since the simulations were started with the initial condition of the orbital period 2 ms (i.e., frequency of gravitational waves is 1 kHz), the spectrum of inspiraling binary neutron stars for f < 1 kHz cannot be taken into account. Thus, only the spectrum for f * 1 kHz should be paid attention. In the panel (a), we plot the following Fourier power spectrum of two point particles in circular orbits in the second post-Newtonian approximation [52]: dE M 3 x 1 x df 3f 2 6M 27 19 2 3 (49) x2 : 8 8M 24M2 Here, and M denote the reduced mass and the total mass of the binary, and x M%f2=3 . We note that the third post-Newtonian terms does not significantly modify the spectrum since their magnitude is 0:01 of the leadingorder term. Thus, the dotted curve should be regarded as the plausible Fourier power spectrum for f & 1 kHz. Figure 17 shows that a sharp characteristic peak is present at f 3– 4 kHz irrespective of models in which hypermassive neutron stars with a long lifetime ( * 10 ms) are formed (see also Table III for the list of the characteristic frequency). This is associated with quasiperiodic gravitational waves emitted by the formed hypermassive neutron stars. The amplitude of the peak is much higher than that in the 2 equation of state [12]. The reason is that with the realistic equations of state, the ellipticity of the formed hypermassive neutron stars is much larger, and as a result, quasiperiodic gravitational waves of a higher amplitude are emitted. Also, the elliptic structure of the hypermassive neutron stars is preserved for a long time duration. These effects amplify the peak amplitude in the Fourier power spectrum. The energy power spectra for models SLy1313a and SLy125135a are very similar reflecting the fact that the waveforms for these two models are very similar [Fig. 17(a)]. This indicates that the spectral shape depends very weakly on the mass ratio QM as far as it is in the range between 0.9 and 1. On the other hand, three peaks are present at f 2:6, 3.6, and 4.5 kHz in the energy power spectrum for model SLy135135b [Fig. 17(b)]. Thus, the spectral shape is quite different from that for model SLy1313a although the total mass is only slightly different between two models. The reason is that the amplitude of the quasiradial oscillation of the hypermassive neutron star is very large and the characteristic radius varies for a wide range for model SLy135135b, inducing the modulation of the wave frequency. Indeed, the difference of the frequencies for the peaks is approximately equal to that of the quasiradial oscillation 1 kHz. As a result, the intensity of the power spectrum is dispersively distributed to multi peaks in this case, and the amplitude for the major peak at f 3:6 kHz is suppressed. The similar feature is also 084021-21 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ PHYSICAL REVIEW D 71, 084021 (2005) (a) (b) (d) (c) FIG. 17 (color online). Fourier power spectrum of gravitational waves dE=df (a) for models SLy1313a (solid curve) and SLy125135a (dashed curve), (b) for models SLy1313a (solid curve) and SLy135135b (dashed curve), (c) for models SLy1313a (dashed curve), SLy1313c (solid curve), and SLy1313d (long-dashed curve), and (d) for models SLy1212b (solid curve) and FPS1212b (dashed curve). Since the simulations are started when the frequency of gravitational waves is 1 kHz, the spectrum for f < 1 kHz is not correct. The dotted curve in the panel (a) denotes the analytical result of dE=df in the second post-Newtonian and point-particle approximation. The real spectrum for f & 1 kHz is approximated by the dotted curves. found for models SLy1313c and FPS1212b for which the hypermassive neutron stars collapse to a black hole within 10 ms. Figs. 17(c) and 17(d) illustrate that the amplitude and the frequency for the peak around f 3– 4 kHz depend on the total mass, the value of th , and the equations of state as in the case of gravitational waveforms. Figure 17(c) indicates that for the larger total mass (but with M < Mthr ), the peak frequency becomes higher. Also, with the increase of the value of th , the peak frequency is decreased since the formed hypermassive neutron star becomes less compact. As Fig. 17(c) shows, the peak frequency is larger for the FPS equation of state than for the SLy one for the same value of the total mass. This is also due to the fact that the hypermassive neutron star in the FPS equation of state is more compact. The effective amplitude of gravitational waves observed from the most optimistic direction (which is parallel to the p axis of the angular momentum) is proportional to dE=df in the manner q hf jR j2 jR j2 f 1=2 dE=df 100 Mpc 21 ; (50) 1:8 10 r 1051 erg=Hz where r denotes the distance from the source, and R ; are 084021-22 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) the Fourier spectrum of R; . Equation (50) implies that the effective amplitude of the peak is about 4 –5 times larger than that at 1 kHz. Furthermore, the amplitude of the peak in reality should be larger than that presented here, since we stopped simulations at t 10 ms to save the computational time, and hence, the integration time 10 ms is much shorter than the realistic value. Extrapolating the decrease rate of the angular momentum, the hypermassive neutron star will dissipate sufficient angular momentum by gravitational radiation until a black hole is formed. As indicated in Secs. III B 1 and III C 2 , the lifetime would be 30–50 ms for models SLy1313a and SLy125135a and 50–100 ms for model SLy1212b. Thus, we may expect that the emission will continue for such time scale and the effective amplitude of the peak of f 3– 4 kHz would be in reality amplified by a factor of 31=2 –101=2 2–3 to be 3–5 1021 at a distance of 100 Mpc. Although the sensitivity of laser interferometric gravitational wave detectors for f > 1 kHz is limited by the shot noise of the laser, this value is larger than the planned noise level of the advanced laser interferometer 1021:5 f=1 kHz3=2 [5]. It will be interesting to search for such quasiperiodic signal of high frequency if the chirp signal of gravitational waves from inspiraling binary neutron stars of distance r & 100 Mpc are detected in the near future. Detection of the quasiperiodic gravitational waves will demonstrate that a hypermassive neutron star of a lifetime much longer than 10 ms is formed after the merger. Since the total mass of the binary should be determined by the data analysis for the chirp signal emitted in the inspiral phase [53], the detection of the quasiperiodic gravitational waves will provide the lower bound of the binary mass for the prompt formation of a black hole Mthr . As found in this paper, the value of Mthr depends sensitively on the equations of state. Furthermore, the values of Mthr (Mthr 2:7M and 2:5M for the SLy and FPS equations of state, respectively) are very close to the total mass of the binary neutron stars observed so far [2]. Therefore, the merge of mass Mthr is likely to happen frequently, and thus, the detection of gravitational waves from hypermassive neutron stars will lead to constraining the equations of state for neutron stars. For example, if quasiperiodic gravitational waves are detected from a hypermassive neutron star formed after the merger of a binary neutron star of mass M 2:6M , the FPS equation of state should be rejected. As this example shows, the merit of this method is that only one detection will significantly constrain the equations of state. The further detail about this method is described in [54]. 4. Calibration of radiation reaction Figure 18 shows evolution of the ADM mass and the angular momentum computed in a finite domain by Eqs. (25) and (26) as well as the violation of the Hamiltonian constraint H defined in Eq. (30) for models SLy1313a and SLy1414a. The solid curves in the upper two panels denote M and J while the dashed curves are M0 &E and J0 &J which are computed from the emitted energy and angular momentum of gravitational waves. The ADM mass and angular momentum computed by two methods should be identical because of the pres- (b) (a) FIG. 18. Evolution of ADM mass and angular momentum in units of their initial values M0 and J0 , and violation of the Hamiltonian constraint (a) for model SLy1313a and (b) for model SLy1414a. In the upper two panels, the solid curves denote M=M0 and J=J0 computed from Eqs. (25) and (26), while the long-dashed curves denote 1 &Et=M0 and 1 &Jt=J0 , respectively, (see Eqs. (27) and (28)). 084021-23 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ ence of the conservation laws. The figure indicates that the conservation holds within 2% error for the ADM mass and angular momentum (except for the case that a black hole is present). This implies that radiation reaction of gravitational waves is taken into account within 2% in our numerical simulation. The error in the angular momentum conservation is generated mainly in the late inspiral phase with t & 2 ms in which L is smaller than the wavelength of gravitational waves and the radiation reaction cannot be evaluated accurately. To improve the accuracy for the conservation in this phase, it is required to take a sufficiently large value of L that is larger than the wavelength. On the other hand, the magnitude of the error does not change much after the formation of the hypermassive neutron stars for t * 2 ms as found in Fig. 18(a). This implies that the radiation reaction of gravitational waves to the angular momentum for the formed hypermassive neutron stars is computed within 1% error. The bottom panels show that the violation of the Hamiltonian constraint is of order 0.01 in the absence of black holes. Also noteworthy is that the violation does not grow but remain small in the absence of black holes. This strongly indicates that simulations will be continued for an arbitrarily long duration for space-times of no black hole. On the other hand, the computation crashes soon after the formation of a black hole for model SLy1414a. This is mainly due to the fact that the resolution around the black hole is too poor. If one is interested in the longterm evolution of the formed black hole, it is obviously necessary to improve the resolution around the black hole to overcome this problem. In the current simulation, the radius of the apparent horizon is covered only by 5 grid points. The axisymmetric simulations for black hole formation (e.g., [26,28]) have experimentally shown that more than 10 grid points for the radius of the apparent horizon will be necessary to follow evolution of the formed black hole for 30M 0:4 ms. To perform such a betterresolved and longterm simulation with L * 0:550 , the grid size more than (1500, 1500, 750) is required, implying that a powerful supercomputer, in which the memory and the computational speed are by a factor of * 10 as large as those of the present computational resources, is necessary. If the computational resources are not improved in the near future, adopting the adaptive mesh refinement technique will be inevitable for following the evolution of the black hole [55]. IV. SUMMARY We performed fully general relativistic simulations for the merger of binary neutron stars adopting realistic equations of state. Since the stiffness is significantly different from that in the 2 equation of state adopted in the previous works (e.g., [12]), several new features have PHYSICAL REVIEW D 71, 084021 (2005) emerged. The following is the summary of the results obtained in this paper: (1) If the ADM mass of the system is larger than 2:7M ( 2:5M ), a black hole is promptly formed in the SLy (FPS) equation of state. Otherwise, a hypermassive neutron star of ellipsoidal shape is formed. This indicates that the threshold mass depends on the equations of state and the values are very close to those for observed binary neutron stars. (2) In the black hole formation, most of mass elements are swallowed into the horizon, and hence, the disk mass around the black hole is much smaller than 1% of the total baryon rest-mass as far as the mass ratio QM is larger than 0.9. Although the disk is hot with the thermal energy 10–20 MeV, the total thermal energy which is available for the neutrino emission is expected to be at most 1050 erg. Since the pair annihilation of the neutrino and antineutrino to the electron-positron pair would be <104 [48], it seems to be very difficult to generate cosmological gamma-ray bursts in this system. (3) The nondimensional angular momentum parameter (J=M2 ) of the formed Kerr black hole is in the range between 0.7 and 0.8. Then, for the system of mass 2:8M , the frequency of gravitational waves associated with quasinormal mode ringing of l m 2 modes would be 6:5–7 kHz, which is too high for gravitational waves to be detected by laser interferometric detectors. (4) The hypermassive neutron stars formed after the merger have a large ellipticity with the axial ratio 0:5. They rotate with the period of 0:5–1 ms, and thus, become strong emitters of quasiperiodic gravitational waves of a rather high frequency f 3–4 kHz. Although the frequency is far out of the best sensitive frequency range of the laser interferometric gravitational wave detectors, the effective amplitude of gravitational waves is very high as several 1021 at a distance of r 100 Mpc. Thus, if the merger happens for r < 100 Mpc, such gravitational waves may be detectable by advanced laser interferometers. The detection of these quasiperiodic gravitational waves will be used for constraining the equations of state for nuclear matter. (5) Because of the larger emission rate of gravitational waves, the angular momentum of the hypermassive neutron star is dissipated in a fairly short time scale & 100 ms for the mass M 2:4–2:7M . Also, due to a high degree of nonaxial symmetry, the angular momentum is transferred outward by the hydrodynamic interaction. As a result of these effects, the hypermassive neutron stars collapse to a black hole within 100 ms. This time scale is much shorter than the viscous dissipation time scale and the transport 084021-24 MERGER OF BINARY NEUTRON STARS WITH . . . PHYSICAL REVIEW D 71, 084021 (2005) time scale of the angular momentum by magnetic fields [22]. Therefore, the gravitational radiation or the outward angular momentum transfer by the hydrodynamic interaction plays the most important role. (6) The thermal energy of the outer region of the hypermassive neutron stars is high as 10–20 MeV, and the total emission rate of the neutrino energy is estimated as 1053 erg=s. The thermal energy is generated by the shocks due to the multiple collisions between the spiral arms and the oscillating hypermassive neutron star. Thus, the hypermassive neutron star will be a strong emitter of neutrinos. However, the emission time scale is 1–10 s which is much shorter than the lifetime <100 ms. This implies that the neutrino cooling plays a minor role in the evolution of the hypermassive neutron star. [1] R. A. Hulse and J. H. Taylor, Astrophys. J. 201, L55 (1975); J. H. Taylor and J. M. Weisberg, Astrophys. J. 345, 434 (1989); J. M. Weisberg and J. H. Taylor, Astrophys. J. 576, 942 (2002). [2] I. H. Stairs, Science 304, 547 (2004). [3] M. Burgey et al., Nature (London) 426, 531 (2003). [4] V. Kalogera et al., Astrophys. J. 601, L179 (2004); 614, L137 (2004). [5] C. Cutler and K. S. Thorne, in Proceedings of the 16th International Conference on General Relativity and Gravitation, edited by N. T. Bishop and S. D. Maharaj (World Scientific, Singapore, 2002), p. 72, and references therein. [6] TAMA Collaboration, M. Ando et al., Phys. Rev. Lett. 86, 3950 (2001). [7] M. Shibata, Phys. Rev. D 60, 104052 (1999). [8] M. Shibata and K. Uryū, Phys. Rev. D 61, 064001 (2000). [9] M. Shibata and K. Uryū, Prog. Theor. Phys. 107, 265 (2002). [10] J. A. Font et al., Phys. Rev. D 65, 084024 (2002). [11] J. A. Font, Living Rev. Relativity 6, 4 (2003). [12] M. Shibata, K. Taniguchi, and K. Uryū, Phys. Rev. D 68, 084020 (2003). [13] M. Miller, P. Gressman, and W.-M. Suen, Phys. Rev. D 69, 064026 (2004). [14] M. D. Duez, P. Marronetti, T. W. Baumgarte, and S. L. Shapiro, Phys. Rev. D 67, 024004 (2003). [15] L. Baiotti, I. Hawke, P. J. Montero, F. Löffer, L. Rezzolla, N. Stergioulas, J. A. Font, and E. Seidel, Phys. Rev. D 71, 024035 (2005). [16] R. Oechslin, K. Uryū, G. Poghosyan, and F. K. Thielemann, Mon. Not. R. Astron. Soc. 349, 1469 (2004). [17] S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars (Wiley Interscience, New York, 1983). (7) The mass difference with the mass ratio QM 0:9 does not modify the dynamics of the merger and the outcome after the merger significantly from that with QM 1. This disagrees with the previous result which was obtained in the simulations performed with the 2 equation of state [12]. The reason is that with the realistic equations of state, the radius of neutron stars is small as 11–12 km depending weakly on the mass in contrast to that in the 2 equations of state. ACKNOWLEDGMENTS Numerical computations were performed on the FACOM VPP5000 machines at the data processing center of NAOJ. This work was in part supported by Monbukagakusho Grant (Nos. 15037204, 15740142, and 16029202). [18] S. Tsuruta, Phys. Rep. 292, 1 (1998). [19] V. R. Pandharipande and D. G. Ravenhall, in Hot Nuclear Matter, edited by M. Soyeur, H. Flocard, B. Tamain, and M. Porneuf NATO Advanced Study Institutes Ser. B, Vol. 205, (Dordrecht, Reidel, 1989), p. 103. [20] F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001). [21] P. Haensel and A. Y. Potekhin, Astron. Astrophys. 428, 191 (2004). [22] T. W. Baumgarte, S. L. Shapiro, and M. Shibata, Astrophys. J. Lett. 528, L29 (2000). [23] F. Rasio and S. L. Shapiro, Astrophys. J. 432, 242 (1994). [24] X. Zhuge, J. M. Centrella, and S. L. W. McMillan, Phys. Rev. D 54, 7261 (1996). [25] J. A. Faber and F. A. Rasio, Phys. Rev. D 62, 064012 (2000); 63, 044012 (2001); 65, 084042 (2002). [26] M. Shibata, Phys. Rev. D 67, 024033 (2003). [27] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 (1995). [28] M. Shibata, Astrophys. J. 595, 992 (2003). [29] M. Shibata, Phys. Rev. D 55, 2002 (1997); M. Shibata and K. Uryū, Phys. Rev. D 62, 087501 (2000). [30] J. R. Wilson and G. J. Mathews, Phys. Rev. Lett. 75, 4161 (1995). [31] C. S. Kochanek, Astrophys. J. 398, 234 (1992); L. Bildsten and C. Cutler, Astrophys. J. 400, 175 (1992). [32] M. Shibata, Phys. Rev. D 58, 024012 (1998); S. A. Teukolsky, Astrophys. J. 504, 442 (1998). [33] E. Gourgoulhon et al., Phys. Rev. D 63, 064029 (2001). [34] K. Taniguchi and E. Gourgoulhon, Phys. Rev. D 66, 104019 (2002); 68, 124025 (2003). [35] V. Moncrief, Ann. Phys. (N.Y.) 88, 323 (1974). The Moncrief formalism was originally derived for the Schwarzschild space-time. We here apply his formalism in a flat space-time. 084021-25 MASARU SHIBATA, KEISUKE TANIGUCHI, AND KŌJI URYŪ [36] M. Shibata, Prog. Theor. Phys. 101, 1199 (1999). [37] M. Shibata and Y. Sekiguchi, Phys. Rev. D 71, 024014 (2005). [38] M. Shibata and Y. Sekiguchi, Phys. Rev. D 68, 104020 (2003). [39] M. Shibata and Y. Sekiguchi, Phys. Rev. D 69, 084024 (2004). [40] S. Chandrasekhar, Stellar Structure (Dover, New York, 1967), Chap. 10. [41] K. Uryū, M. Shibata, and Y. Eriguchi, Phys. Rev. D 62, 104015 (2000). [42] M. Shibata and K. Uryū, Phys. Rev. D 64, 104017 (2001). [43] R. A. James, Astrophys. J. 140, 552 (1964); J.-L. Tassoul, Theory of Rotating Stars (Princeton University Press, Princeton, 1978). [44] M. Shibata, T. W. Baumgarte, and S. L. Shapiro, Astrophys. J. 542, 453 (2000). [45] E. W. Leaver, Proc. R. Soc. London A 402, 285 (1985). T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. Phys. Suppl. 90, 1 (1987). [46] N. D. Lyford, T. W. Baumgrate, and S. L. Shapiro, Astrophys. J. 583, 410 (2003). PHYSICAL REVIEW D 71, 084021 (2005) [47] I. A. Morrison, T. W. Baumgrate, and S. L. Shapiro, Astrophys. J. 610, 941 (2004). [48] M. Ruffert, H.-Th. Janka, and G. Schäfer, Astron. Astrophys. 311, 532 (1996); M. Ruffert et al., ibid. 319, 122 (1997); M. Ruffert and H.-Th. Janka, ibid. 380, 544 (2001). [49] H.-Th. Janka and M. Ruffert, Astron. Astrophys. 307, L33 (1996). [50] T. Piran, Phys. Rep. 314, 575 (1999); 333, 529 (2000); Rev. Mod. Phys. 76, 1143 (2005). [51] D. Lai and S. L. Shapiro, Astrophys. J. 442, 259 (1995). [52] L. Blanchet, G. Faye, B. R. Iyer, and B. Joguet, Phys. Rev. D 65, 061501 (2002); L. Blanchet, T. Damour, G. Esposito-Farese, and B. R. Iyer, Phys. Rev. Lett. 93, 091101 (2004 ). [53] C. Cutler and E. E. Flanagan, Phys. Rev. D 49, 2658 (1994). [54] M. Shibata, in preparation. [55] M. Berger and J. Oliger, J. Comput. Phys. 53, 484 (1984); E. Enans, S. Iyer, E. Schnetter, W.-M. Suen, J. Tao, R. Wolfmeyer, and H.-M. Zhang, gr-qc/0501066; B. Zink, N. Stergioulas, I. Hawke, C. D. Ott, E. Schnetter, and E. Müller, gr-qc/0501080. 084021-26