Table Of Content

The mass ejection from the merger of binary neutron stars Kenta Hotokezaka1, Kenta Kiuchi2, Koutarou Kyutoku3, Hirotada Okawa4, Yu-ichiro Sekiguchi2, Masaru Shibata2, and Keisuke Taniguchi5 1Department of Physics, Kyoto University, Kyoto 606-8502, Japan 2Yukawa Institute of Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan 3Theory Center, Institute of Particles and Nuclear Studies, KEK, Tsukuba, Ibaraki 305-0801, Japan 4CENTRA, Departamento de F´ısica, Instituto Superior Técnico, Universidade Técnica de Lisboa - UTL,Av. Rovisco Pais 1, 1049 Lisboa, Portugal 5Graduate School of Arts and Sciences, The University of Tokyo, Tokyo 153-8902, Japan Numerical-relativitysimulationsforthemergerofbinaryneutronstarsareperformedforavariety of equations of state (EOSs) and for a plausible range of the neutron-star mass, focusing primarily on the properties of the material ejected from the system. We find that a fraction of the material 3 is ejected as a mildly relativistic and mildly anisotropic outflow with the typical and maximum 1 velocities ∼0.15 – 0.25c and ∼0.5 – 0.8c (where c is the speed of light), respectively, and that the 0 totalejectedrestmassisinawiderange10−4 –10−2M ,whichdependsstronglyontheEOS,the 2 (cid:12) total mass, and the mass ratio. The total kinetic energy ejected is also in a wide range between n 1049 and 1051 ergs. The numerical results suggest that for a binary of canonical total mass 2.7M , (cid:12) a the outflow could generate an electromagnetic signal observable by the planned telescopes through J the production of heavy-element unstable nuclei via the r-process [1–3] or through the formation 1 of blast waves during the interaction with the interstellar matter [4], if the EOS and mass of the 3 binary are favorable ones. ] PACSnumbers: 04.25.Dm,04.30.-w,04.40.Dg E H . I. INTRODUCTION study for this subject, have been extensively performed h since the first success in 2000 [21] (see, e.g., [22, 23] for p a review of this field). However, most of the simulations - Coalescence of binary neutron stars is one of the o have focused on the studies of gravitational waveforms most promising sources for next-generation kilo-meter- r and the resulting product formed in the central region. t sizegravitational-wavedetectorssuchasadvancedLIGO, s Few attention has been paid to the study for the mate- a advanced VIRGO, and KAGRA (LCGT) [5]. These de- rial ejected (but see [13] for a study in an approximate [ tectors will detect gravitational waves in the next 5 – framework of general relativity, and see [14–17] for an 10 yrs. Statistical studies have predicted that the de- 2 early effort in the context of Newtonian gravity). v tection rate of gravitational waves emitted by binary 5 neutron stars for these detectors will be ∼ 1 – 100 per The material ejected from binary neutron star merg- 0 year [6, 7]. The typical signal-to-noise ratio for most ers may generate electromagnetic signals observable in 9 of these events will be ∼ 10 or less. Thus, it will be the current and future-planned telescopes. One possi- 0 quite helpful if electromagnetic or other signals observ- ble signal could be generated by the radioactive decay of . 2 able are associated with the gravitational-wave bursts unstable r-process nuclei, which are produced from the 1 and the gravitational-wave detection is accompanied by neutron-rich material in the ejecta [1–3, 12, 13, 17, 18]. 2 the detections of other signals. Short-hard gamma-ray A fraction of the unstable nuclei produced subsequently 1 bursts (SGRB) have been inferred to be associated with : decay in a short timescale and could heat up the ejecta, v the binary neutron star merger [8]. However, the jet which emits a UV and visible light that may be observ- i of SGRB would be highly collimated [9], and hence, it X ablebycurrentandfuture-plannedopticaltelescopes. In will not be always possible to detect SGRB associated thiscase, thetypicaldurationofapeakluminosityisex- r a with the binary neutron star mergers. Moreover, it is pected to be several hours to a day. Another possible not guaranteed that the telescopes for the observation of signal could be generated during the free expansion and SGRB will be in operation with the gravitational-wave the subsequent Sedov phase of the ejecta which sweeps detectors. Exploring other possible signals that could up the interstellar medium and forms blast waves [4]. In be detected is a very important subject in the fields of thisprocessturningon,theshockedmaterialattheblast gravitational-wave physics/astronomy [3, 4, 10–13, 18]. waves could generate magnetic fields and accelerate par- Thispaperpresentsourlatestresultsofnumericalsim- ticles that emit synchrotron radiation in the radio-wave ulations performed in the framework of numerical rela- band, for a hypothetical amplification of the electromag- tivity, focusing in particular on the exploration of the netic field and a hypothetical electron injection. It is material ejected from binary neutron star mergers. In also pointed out that the binary neutron star merger the past decade, numerical simulations for the merger of could drive ultra-relativistic outflows in every direction binary neutron stars in full general relativity, which is and emit synchrotron radiation in x-ray-to-radio bands probably the unique approach of the rigorous theoretical within a second-to-day timescale [19]. All these stud- 2 ies illustrate that exploring the process of the material A. Equations of state ejection from binary neutron star mergers in detail is an important subject. The exact EOS for the high-density nuclear matter is For the detailed numerical study of the ejected mate- stillunknown[28]. Thisimpliesthatanumericalsimula- rial, we have to be careful when following the motion of tionemployingasingleparticularEOS,whichmight not the materials in a low-density outer region. Most of the be correct, would not yield a scientific result. A study, numerical-relativity simulations of binary neutron star systematically employing a wide possible range of EOSs, mergers so far have been performed with a computa- is required for binary neutron star mergers. Neverthe- tional domain that was not wide enough for this pur- less,thelatestdiscoveryofahigh-massneutronstarPSR pose [22, 23]. We have to enlarge the computational do- J1614-2230 with mass 1.97 ± 0.04M(cid:12) [29] significantly main sufficiently widely to confirm that the outflowed constrains the model EOS to be chosen, because it sug- material is indeed ejected from the system (i.e., we have geststhatthemaximummassforsphericalneutronstars to confirm that the material is indeed unbound by the foragivenEOShastobelargerthan∼2M(cid:12). Thisindi- system by following the motion of the ejected material catesthattheEOSshouldberatherstiff, althoughthere for a long time). Another subtle issue in the hydrody- are still many candidate EOSs. namics simulations is that we have to put an artificial TomodelavarietyofthecandidateEOSs,specifically, atmosphere when employing a conservative shock cap- weemployapiecewisepolytropicEOSproposedbyRead turing scheme that is a standard one in this field [24]. In et al. [30]. This EOS is described assuming that neu- our previous simulations [25–27], we put an atmosphere tronstarsarecold(inazero-temperaturestate),i.e.,the with fairly large density (∼ 107 g/cm3) that did not af- rest-mass density, ρ, determines all other thermodynam- fect the motion of neutron stars but did for the motion ical quantities. To systematically model nuclear-theory- of the ejected material of low density which might es- based EOSs at high density with a small number of pa- cape to a far region. For the study of the mass ejection, rameters,thepressureiswritteninaparameterizedform we have to reduce the density of the atmosphere as low as as possible (which should be much lower than the density of the ejected material), and in addition, we have to P(ρ)=κiρΓi for ρi ≤ρ<ρi+1 (0≤i≤n), (1) carefully confirm that such an artificial atmosphere does not affect the properties of the ejected material. In the wherenisthenumberofthepiecesusedtoparameterize simulationreportedinthispaper, wesucceedinthesim- a high-density EOS, ρi is the rest-mass density at the ulation reducing the atmosphere density to a low level boundary of two neighboring (i−1)-th and i-th pieces, (<∼ 105 g/cm3) enough to obtain a scientifically quanti- κi is the polytropic constant for the i-th piece, and Γi tative result. is the adiabatic index for the i-th piece. Here, ρ0 = 0, ρ denotes a nuclear density ∼ 1014 g/cm3 determined 1 The paper is organized as follows: In Sec. II, we below, and ρ →∞. Other parameters (ρ ,κ ,Γ ) are n+1 i i i summarize the equations of state (EOSs) employed and determined by fitting with a nuclear-theory-based EOS. models of binary neutron stars. In Sec. III, we briefly Requiring the continuity of the pressure at each ρ , 2n i summarizeourformulationandnumericsforsolvingEin- free parameters, say (κ ,Γ ), determine the EOS com- i i stein’s equation and hydrodynamics equations as well as pletely. The specific internal energy, ε, and hence the the tools for diagnostics. In Sec. IV, numerical results specific enthalpy, h, are determined by the first law of are presented, focusing on the properties of the material thermodynamics and the continuity of each variable at ejected from the system. Section V is devoted to a sum- boundary densities, ρ . i mary and discussion. Throughout this paper, we employ Readetal.[30]showedthatapiecewisepolytropicEOS the geometrical units c = 1 = G where c and G are the with three pieces above the nuclear density (i.e., n = 3) speed of light and gravitational constant, respectively, approximatelyreproducesmostpropertiesofthenuclear- although we recover c when we need to clarify the units. theory-based EOS at high density, and they derived the fitted parameters for a large number of nuclear-theory- based EOSs. In this paper, thus, we employ this piecewise polytropic EOS, determining the free parameters basically following [31–33] (in which a piecewise poly- II. EQUATIONS OF STATE AND CHOSEN trope with n = 1 was used). First, the EOS below the MODELS nuclear density ρ is fixed by the following parameters 1 Γ = 1.35692395, (2) In this section, we summarize the model EOSs em- 0 ployed in this paper, and initial condition of binary neu- κ0/c2 = 3.99873692×10−8 (g/cm3)1−Γ0. (3) tron stars chosen for numerical simulations. As shown in Sec. IV, the properties of the material ejected from The EOS for the nuclear matter was determined in [30] binary neutron star mergers depend strongly on these as follows: ρ was fixed to be ρ = 1014.7g/cm3, and P 2 2 2 inputs. at ρ = ρ was chosen as a free parameter. The reason 2 3 TABLE I: Parameters and key quantities for four piecewise polytropic EOSs employed in this paper. P is shown in units of 2 dyn/cm2. M isthemaximummassalongthesequencesofsphericalneutronstars(cf. Fig.2). (R ,ρ )and(R ,ρ ) max 1.35 1.35 1.5 1.5 are the circumferential radius and the central density of 1.35M and 1.5M neutron stars, respectively. We note that the (cid:12) (cid:12) values of the mass, radius, and density listed are slightly different from those obtained inthe original tabulated EOSs (see the text for the reason). MS1 is referred to as this name in [30], but in other references (e.g., [28]), it is referred to as MS0. We follow [30] in this paper. EOS (log(P ),Γ ,Γ ,Γ ) M (M ) R (km) ρ (g/cm3) R (km) ρ (g/cm3) 2 1 2 3 max (cid:12) 1.35 1.35 1.5 1.5 APR4 (34.269,2.830,3.445,3.348) 2.20 11.1 8.9×1014 11.1 9.6×1014 ALF2 (34.616,4.070,2.411,1.890) 1.99 12.4 6.4×1014 12.4 7.2×1014 H4 (34.669,2.909,2.246,2.144) 2.03 13.6 5.5×1014 13.5 6.3×1014 MS1 (34.858,3.224,3.033,1.325) 2.77 14.4 4.2×1014 14.5 4.5×1014 is that P is closely related to the radius and deforma- 2 1e+37 bility of neutron stars [34]. Namely, P2 primarily deter- APR4 mines the stiffness of an EOS. Second, ρ was fixed to ALF2 3 be ρ = 1015.0g/cm3. With these choices, the set of free 1e+36 H4 3 MS1 parameters becomes (P2,Γ1,Γ2,Γ3). These four param- 2) eters are determined by a fitting procedure (see [30] for m 1e+35 c the fitting procedure). / n With the given values of Γ1 and P2, κ1 and ρ1 are dy 1e+34 subsequently determined by P ( κ = P ρ−Γ1, (4) 1e+33 1 2 2 ρ = (κ /κ )1/(Γ1−Γ0). (5) 1 0 1 1e+32 1e+14 1e+15 By the same method, κ and κ are determined from 2 3 ρ (g/cm3) κ ρΓ2 =κ ρΓ1, κ ρΓ3 =κ ρΓ2. (6) 2 2 1 2 3 3 2 3 FIG. 1: Pressure as a function of the rest-mass density for Table I lists the EOSs and their parameters which we APR4,ALF2,H4,andMS1EOSs(thesolid,dashed,dotted, employ in this study. We choose four types of the rep- and dash-dotted curves, respectively). resentative EOSs. APR4 was derived by a variational method with modern nuclear potentials [35] for the hypothetical components composed of neutrons, protons, EOS that is soft at ρ = ρ has to be in general stiff for 2 electrons, and muons; MS1 was derived by a mean-field ρ >∼ ρ3. By contrast, H4 and MS1 have pressure higher theoryforthehypotheticalcomponentscomposedofneu- thanAPR4forρ<∼ρ3,althoughtheEOSsbecomesofter trons,protons,electrons,andmuons,aswell[36];H4was for a high-density region ρ>∼ρ3. In particular, MS1 has derived by a relativistic mean-field theory including ef- extremelyhighpressure(i.e.,ahighervalueofP )among 2 fects of hyperons [37]; ALF2 is a hybrid EOS which de- many other EOSs for ρ <∼ ρ3, and thus, it is the stiffest scribes a nuclear matter for a low density and a quark EOS as far as the canonical neutron stars are concerned. matter for a high density with the transition density is ALF2 has small pressure for ρ ≤ ρ as in the case of 2 3ρnuc where ρnuc ≈2.8×1014 g/cm3 [38]. We note that APR4, but for ρ2 <∼ ρ ≤ ρ3, the pressure is higher than the piecewise polytropic EOSs are slightly different from that for APR4. For ρ ≥ ρ the pressure of ALF2 is as 2 theoriginaltabulatedones,becauseoftheirsimplefitting high as that for H4. All the properties mentioned above formula. Thisresultsinasmallerrorinthemassandra- are reflected in the radius, R , and central density, 1.35 diusofneutronstars. However,asshownin[30],theerror ρ , of spherical neutron stars with the canonical mass 1.35 is small (at most several percent), and the semiquantita- M = 1.35M where M is the gravitational (Arnowitt- (cid:12) tive properties of the original EOSs are well captured by Deser-Misner; ADM) mass of the spherical neutron stars these simple EOSs. in isolation: see Table I. The pressure at ρ = ρ (P ) is 2 2 Figure 1 plots the pressure as a function of the rest- correlated well with this radius and central density (see mass density for four EOSs. APR4 has relatively small below). pressureforρ1 ≤ρ<∼ρ3whileithashighpressureforρ>∼ Here, a word of caution is necessary for our APR4. ρ . Thus, for ρ < ρ , which neutron stars of canonical The pressure in this piecewise polytropic EOS is ex- 3 3 mass 1.3 – 1.4M have, this EOS is soft, and hence, the tremely (unphysically) high in the high-density region (cid:12) value of P2 is relatively small. We note that for a small with ρ >∼ 1016 g/cm3. This results pathologically in the value of P , Γ and/orΓ have tobe large(∼3) because situation that the sound velocity exceeds the speed of 2 2 3 the maximum mass of spherical neutron stars, M for light for the high-density state. In reality, such a high max a given EOS has to be larger than ∼ 2M . Thus, an density is achieved only in the formation of a black hole (cid:12) 4 (i.e., inside the horizon), and such a pathology may not APR4-120150. affect the evolution of the system for the outside of the horizon. However, this pathology could still break a numerical simulation after the formation of a black hole. III. FORMULATION AND NUMERICAL To avoid this happens, we artificially set the maximum METHODS density as 1016 g/cm3 when employing this EOS. Figure 2 plots the gravitational mass as a function of Numerical simulations were performed using an the central density and as a function of the circumfer- adaptive-mesh refinement (AMR) code SACRA [42] (see ential radius for spherical neutron stars for four EOSs. also [43] for the reliability of SACRA). The formulation, All the EOSs chosen are stiff enough that the maxi- the gauge conditions, and the numerical scheme are ba- mum mass is larger than 1.97M . Because the pres- (cid:12) sically the same as those described in [42], except for sure in a density region ρ <∼ 1015 g/cm3 is relatively the improvement in the treatment of the hydrodynamics small (i.e., P is small) for APR4 and ALF2, the ra- 2 code for a far region. Thus, we here only briefly review dius for these EOSs is relatively small as ∼ 11 km and themanddescribethepresentsetupofthecomputational 12.5 km, respectively, for the canonical mass of neutron domain for the AMR algorithm and grid resolution. stars 1.3 – 1.4M [40]. By contrast, for H4 and MS1 (cid:12) for which P is relatively large, the radius becomes a 2 relatively large value 13.5 – 14.5 km for the canonical A. Formulation and numerical methods mass. The radius has also the correlation with the cen- traldensityρ . ForAPR4andALF2withM =1.35M , c (cid:12) ρ ≈8.9×1014 g/cm3 andρ ≈6.4×1014 g/cm3. ForH4 SACRA solves Einstein’s evolution equations in the c c Baumgarte-Shapiro-Shibata-Nakamuraformalismwitha andMS1withM =1.35M ,thecentraldensityisrather (cid:12) low as ρ ≈5.5×1014 g/cm3 and ρ ≈4.1×1014 g/cm3, moving-puncture gauge [44]. It evolves a conformal fac- c c tor W := γ−1/6, the conformal three-metric γ˜ := respectively. AsweshowinSec.IV,thepropertiesofthe ij material ejected from the merger of binary neutron stars γ−1/3γij, the trace of the extrinsic curvature K, a depend strongly on the radius of the neutron stars or ρ . conformally-weightedtrace-freepartoftheextrinsiccur- c vatureA˜ :=γ−1/3(K −Kγ /3),andanauxiliaryvari- ij ij ij able Γ˜i := −∂ γ˜ij. Introducing an additional auxiliary j B. Initial conditions variable Bi and a parameter ηs, which we typically set to be ≈0.8/m in units of c=G=M =1, we employ a (cid:12) moving-puncture gauge in the form [45] We employ binary neutron stars in quasiequilibria for the initial condition of numerical simulations as in our (∂ −βj∂ )α = −2αK, (7) series of papers [25, 26]. The quasiequilibrium state is t j computed in the framework described in [39] to which (∂t−βj∂j)βi = (3/4)Bi, (8) the reader may refer. The computation of quasiequilib- (∂ −βj∂ )Bi = (∂ −βj∂ )Γ˜i−η Bi. (9) t j t j s rium states is performed using the spectral-method li- brary LORENE [41]. We evaluate the spatial derivative by a fourth-order cen- Numerical simulations were performed, systematically tralfinitedifferenceexceptfortheadvectionterms,which choosing wide ranges of the total mass and mass ratio of are evaluated by a fourth-order lopsided upwind finite binary neutron stars. Because the mass of each neutron differencing scheme, and employ a fourth-order Runge- star in the observed binary systems is in a narrow range Kutta method for the time integration. between ∼ 1.2 – 1.45M(cid:12) [40], we basically choose the To solve hydrodynamics equations, we evolve ρ∗ := neutron-star mass 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, and ραutW−3, uˆ :=hu , and e :=hαut−P/(ραut). Here, i i ∗ 1.5M(cid:12). Also, the massratioof the observedsystem q := uµ denotes the four velocity of the fluid. The advection m1/m2(≤ 1) where m1 and m2 are lighter and heavier terms are handled with a high-resolution central scheme masses, respectively, is in a narrow range ∼ 0.85 – 1. by Kurganov and Tadmor [46] with a third-order piece- Thus, we choose q as 0.8 ≤ q ≤ 1. Specifically, the wise parabolic interpolation for the cell reconstruction. simulations were performed for the initial data listed in We note that the total rest mass of the system is calcu- Table II. lated by The initial data were prepared so that the binary has about 3 – 4 quasicircular orbits before the onset of the (cid:90) M = ρ d3x. (10) merger. For four EOSs chosen, this requirement is ap- ∗ ∗ proximately satisfied with the choice of the initial angular velocity mΩ = 0.026 for APR4 and ALF2 and FortheEOSemployedinthenumericalsimulation,we 0 mΩ = 0.025 for H4 and MS1. Here, m = m +m . decompose the pressure and specific internal energy into 0 1 2 For the following, the model is referred to as the name cold and thermal parts as “EOS”-“m ”“m ”; e.g., the model employing APR4, 1 2 m = 1.2M , and m = 1.5M is referred to as model P =P +P , ε=ε +ε . (11) 1 (cid:12) 2 (cid:12) cold th cold th 5 3 3 APR4 2.5 2.5 ALF2 ass) 2 ass) 2 MHS41 m m ar 1.5 ar 1.5 ol ol M (s 1 AAPLRF42 M (s 1 0.5 H4 0.5 MS1 0 0 0 1e+15 2e+15 10 15 20 ρ (g/cm3) c R (km) FIG. 2: Left: The gravitational mass as a function of the central density ρ for spherical neutron stars in APR4, ALF2, H4, c andMS1EOSs(thesolid,dashed,dotted,anddash-dottedcurves). Right: Thesameastheleftpanelbutforthegravitational mass as a function of the circumferential radius. We calculate the cold parts of both variables using the theatmosphereisalways∼10−6M orless. Intestsim- (cid:12) piecewise polytropic EOS (see section IIA) from the ulations,wealsoadoptedn=2andf =10−10 –10−12, at primitivevariableρ,andthenthethermalpartofthespe- and found that the numerical results on the ejected ma- cificinternalenergyisdefinedfromεasε =ε−ε (ρ). terial such as its mass and its total energy depend only th cold Becauseε vanishesintheabsenceofshockheating, ε weakly on the values of n and f (e.g., the ejected mass th th at is regarded as the finite temperature part determined by increases by ∼ 10% if we change n from 3 to 2 (denser the shock heating in the present context. In this paper, one)forsomemodelsofAPR4andH4). Hence,wecould we adopt a Γ-law ideal gas EOS for the thermal part as safely conclude that the tenuous atmosphere chosen in this work does not significantly affect the properties of P =(Γ −1)ρε . (12) th th th the ejected material. Following the conclusion of a detailed study in [47], Γ We extracted l = |m| = 2 modes of gravitational th is chosen in the range 1.6 – 2.0 with the canonical value waves, h and h , by calculating the outgoing part + × 1.8. For several models, we performed simulations vary- of the complex Weyl scalar Ψ at finite coordinate radii 4 ingthevalueofΓ ,andexploredtheeffectsoftheshock r ≈200M –400M andbyintegratingΨ twiceintime th (cid:12) (cid:12) 4 heating; as shown in Sec. IV, numerical results depend as in [33], to which the reader may refer (see also [48]). fairly strongly on the value of Γ (although the depen- Wealsoanalyzedtheevolutionofgravitational-wavefre- th dence on Γ is not as strong as the dependence on the quency, which is determined by extracting the phase th EOS, P ). of Ψ , arg(Ψ ), and by taking the time derivative as cold 4 4 Becausethevacuumisnotallowedinanyconservative 2πf := d(arg(Ψ ))/dt. To find the characteristic fre- 4 hydrodynamicsscheme(e.g.,toderivethevelocitybydi- quencyofgravitationalwaves, wealsodefinetheaverage vidingthemomentumdensitybythedensity),weputan value of f by artificialatmosphereofsmalldensityoutsidetheneutron (cid:90) stars. The atmosphere has to be as tenuous as possible f|h|dt because a dense atmosphere may significantly affect the motionofthematerialejectedfrombinaryneutronstars. fave := (cid:90) , (14) |h|dt Specifically, we set the density of the atmosphere in the following simple rule (cid:26) where we used |h| = (h2 +h2)1/2 as the weight factor. f ρ (r ≤r ), + × ρ = at max uni (13) Then, we define the physical dispersion of f by at f ρ (r/r )−n (r ≥r ), at max uni uni (cid:90) where ρ denotes the maximum rest-mass density of max (f −f )2|h|dt tThaeblneeuI)t.roWnestatyrspiactaltlhyesientitfial s=tate10<∼−1100,15ng=/cm33, (asnede σf2 := (cid:90) ave . (15) at |h|dt r = 16L where 2L denotes the side length of uni min min the finest computational domain in the AMR algorithm (see Sec. IIIC and Table III). For MS1, a computational In the following, f and σ are calculated for gravi- ave f region is wider and we set f =10−11 to reduce the at- tational waves emitted by the remnant massive neutron at mosphere mass. In these settings, the total rest mass of stars. 6 TABLEII:Listoftheparametersoftheinitialconditionforbinarieschoseninnumericalsimulations: Totalmass,massratio, masses of two components, initial value of angular velocity, and initial frequency of gravitational waves (f =Ω /π). 0 0 Model m(M ) q m (M ) m (M ) mΩ f (Hz) (cid:12) 1 (cid:12) 2 (cid:12) 0 0 APR4-130160 2.90 0.813 1.30 1.60 0.026 579 APR4-140150 2.90 0.933 1.40 1.50 0.026 579 APR4-145145 2.90 1.000 1.45 1.45 0.026 579 APR4-130150 2.80 0.867 1.30 1.50 0.026 600 APR4-140140 2.80 1.000 1.30 1.50 0.026 600 APR4-120150 2.70 0.800 1.20 1.50 0.026 622 APR4-125145 2.70 0.862 1.25 1.45 0.026 622 APR4-130140 2.70 0.929 1.30 1.40 0.026 622 APR4-135135 2.70 1.000 1.35 1.35 0.026 622 APR4-120140 2.60 0.857 1.20 1.40 0.026 646 APR4-125135 2.60 0.926 1.25 1.35 0.026 646 APR4-130130 2.60 1.000 1.30 1.30 0.026 646 ALF2-140140 2.80 1.000 1.40 1.40 0.026 600 ALF2-120150 2.70 0.800 1.20 1.50 0.026 622 ALF2-125145 2.70 0.862 1.25 1.25 0.026 622 ALF2-130140 2.70 0.929 1.30 1.40 0.026 622 ALF2-135135 2.70 1.000 1.35 1.35 0.026 622 ALF2-130130 2.60 1.000 1.30 1.30 0.026 646 H4-130150 2.80 0.867 1.30 1.50 0.025 577 H4-140140 2.80 1.000 1.40 1.40 0.025 577 H4-120150 2.70 0.800 1.20 1.50 0.025 598 H4-125145 2.70 0.862 1.25 1.25 0.025 598 H4-130140 2.70 0.929 1.30 1.40 0.025 598 H4-135135 2.70 1.000 1.35 1.35 0.025 598 H4-120140 2.60 1.000 1.30 1.30 0.025 621 H4-125135 2.60 1.000 1.30 1.30 0.025 621 H4-130130 2.60 1.000 1.30 1.30 0.025 621 MS1-140140 2.80 1.000 1.40 1.40 0.025 577 MS1-120150 2.70 0.800 1.20 1.50 0.025 598 MS1-125145 2.70 0.862 1.25 1.25 0.025 598 MS1-130140 2.70 0.929 1.30 1.40 0.025 598 MS1-135135 2.70 1.000 1.35 1.35 0.025 598 MS1-130130 2.60 1.000 1.30 1.30 0.025 621 B. Analysis of the ejected material energy of the fluid element of |u |>1 by t (cid:90) M = ρ d3x, (16) ∗esc ∗ |ut|>1 (cid:90) √ E = T nµnν γd3x tot,esc µν In this section, we describe the method for analyzing |ut|>1 (cid:90) the material ejected from the merger of binary neutron = ρ e d3x, (17) ∗ ∗ stars. Here, the ejected material is composed of a fluid |ut|>1 element which is unbound by the gravitational potential (cid:90) U = ρ εd3x, (18) of binary neutron stars and an object formed after the esc ∗ merger. Thus, first of all, we have to determine which |ut|>1 fluid elements are unbound. To assess this point for all where Tµν is the stress-energy tensor, the fluid elements, we calculate u tµ = u at each grid µ t T =ρhu u +Pg , (19) point. Here, tµ is a timelike vector (1,0,0,0) which is a µν µ ν µν Killing vector at spatial infinity. If |u |>1, we consider andnµ istheunittimelikehypersurfacenormal. Wenote t that the fluid element there is unbound. that the total energy is not uniquely defined by E tot,esc for dynamical spacetimes, and thus, the total energy defined here should be considered as an approximate mea- sure for it. We here choose this expression for simplicity. We then define kinetic energy approximately by Thenwecalculatethetotalrestmass,totalenergy(ex- cludinggravitationalpotentialenergy),andtotalinternal T :=E −M −U . (20) ∗esc tot,esc ∗esc esc 7 TABLEIII:ThegridstructureforthesimulationinourAMRalgorithm. ∆xisthegridspacinginthefinest-resolutiondomain with L being the location of the outer boundaries along each axis and L = N∆x. R /∆x denotes the numbers of grid min diam assigned inside the semi-major diameter of the lighter and heavier neutron stars in the finest level. λ is the gravitational 0 wavelength for the initial configuration. The last column shows the values of Γ employed. th Model ∆x(km) R /∆x L (km) L (km) λ (km) Γ diam min 0 th APR4-130160 0.172 102, 96 2636 10.3 518 1.8 APR4-140150 0.167 102, 101 2572 10.0 518 1.8 APR4-145145 0.166 102, 102 2550 10.0 518 1.8 APR4-130150 0.172 102, 98 2636 10.3 500 1.8 APR4-140140 0.167 102, 102 2572 10.0 500 1.8 APR4-120150 0.172 103, 98 2644 10.3 482 1.6, 1.8, 2.0 APR4-125145 0.174 102, 100 2665 10.4 482 1.8 APR4-130140 0.170 103, 101 2609 10.2 482 1.8 APR4-135135 0.169 102, 102 2601 10.2 482 1.6, 1.8, 2.0 APR4-120140 0.174 102, 99 2679 10.5 464 1.8 APR4-125135 0.174 102, 100 2665 10.4 464 1.8 APR4-130130 0.171 102, 102 2629 10.3 464 1.8 ALF2-140140 0.195 102, 102 3001 11.7 500 1.8 ALF2-120150 0.200 102, 98 3065 12.0 482 1.8 ALF2-125145 0.199 102, 100 3054 11.9 482 1.8 ALF2-130140 0.198 102, 101 3044 11.9 482 1.8 ALF2-135135 0.195 103, 103 3001 11.7 482 1.8 ALF2-130130 0.199 102, 102 3054 11.9 464 1.8 H4-130150 0.222 102, 98 3429 13.4 480 1.8 H4-140140 0.219 102, 102 3358 13.1 480 1.8 H4-120150 0.228 102, 96 3501 13.7 463 1.6, 1.8, 2.0 H4-125145 0.226 102, 98 3465 13.5 463 1.8 H4-130140 0.223 102, 100 3430 13.4 463 1.8 H4-135135 0 221 102, 102 3393 13.3 463 1.6, 1.8, 2.0 H4-120140 0.230 101, 98 3537 13.8 446 1.8 H4-125135 0.227 102, 100 3494 13.6 446 1.8 H4-130130 0.223 103, 103 3430 13.4 446 1.8 MS1-140140 0.237 103, 103 3644 14.2 480 1.8 MS1-120150 0.249 101, 97 3823 14.9 463 1.8 MS1-125145 0.244 102, 99 3751 14.7 463 1.8 MS1-130140 0.244 101, 100 3751 14.7 463 1.8 MS1-135135 0.242 102, 102 3715 14.5 463 1.8 MS1-130130 0.244 102, 102 3751 14.7 446 1.8 We found irrespective of models that T is much (by anisotropic shell of uniform density as ∗esc about 1 – 2 orders of magnitude) larger than U . esc  To approximately analyze the configuration of the ρ π/2−θ ≤θ ≤π/2+θ  esc 0 0 ejected material, we also calculate the moments of in- ρ= and R ≤r ≤R , (23) − + ertia defined by 0 otherwise, (cid:90) Iii,esc = ρ∗(xi)2d3x, (no sum for i), (21) where ρesc, R±, and θ0 are time-varying parameters. In |ut|>1 this case, and then, define 4π M = ρ (R3 −R3)sinθ , (24) (cid:114) (cid:114) (cid:114) ∗esc 3 esc + − 0 I I I X¯ = xx,esc, Y¯ = yy,esc, Z¯ = zz,esc, 1R5 −R5 M M M R¯2 = + −(3−sin2θ ), (25) ∗esc ∗esc ∗esc 5R3 −R3 0 (22) + − and R¯ = √X¯2+Y¯2. From dR¯/dt and dZ¯/dt, we can Z¯2 = 51RR+35 −−RR−35 sin2θ0. (26) + − determine the typical (average) velocity of the ejected material, whichisdenotedbyV¯R andV¯Z inthefollow- Thus for an axial ratio, esc esc ing. We consider a model that the configuration of the Z¯ η = , (27) ejected material is approximated by an axisymmetric R R¯ 8 sinθ is calculated as diameter of each neutron star is covered approximately 0 by 100 grid points for N =60. 3η2 sin2θ = R . (28) 0 1+η2 R IV. NUMERICAL RESULTS Hence,fromtheaxialratiocalculatedforanumericalre- sult of the ejected material, we can approximately define the extent in the θ direction; e.g., for η = 0.4 and 0.5, Table IV summarizes the remnant formed, the rest R θ ≈40◦ and 51◦, respectively. massandkineticenergyoftheejectedmaterialmeasured 0 at10msaftertheonsetofthemergert=t ,andthe merge characteristic (average) frequency of gravitational waves C. Setup of AMR grids emitted by the hypermassive neutron star (HMNS) for N = 60 [61]. Here, t is chosen to be the time at merge An AMR algorithm implemented in SACRA can pre- which the amount of the rest mass of the ejected mate- pare a fine-resolution domain in the vicinity of compact rial steeply increases. In the following two subsections, objects as well as a sufficiently wide domain that covers we summarize the results for the formation of HMNSs a local wave zone. In the present study, we prepare ad- and black holes separately. ditional domains wider than those used in our previous studies [26, 32, 33], to follow the motion of the material ejected during the merger of binary neutron stars for a A. Properties of the merger and mass ejection: sufficiently long time (longer than 10 ms). HMNS case The chosen AMR grids consist of a number of computational domains, each of which has the uniform, vertex- Binaryneutronstarsinquasicircularorbitsevolvedue centeredCartesiangridswith(2N+1,2N+1,N+1)grid to the gravitational-wave emission. Their orbital separa- pointsfor(x,y,z)withtheequatorialplanesymmetryat tion decreases gradually, and eventually, the merger sets z =0. Since we chose that the grid spacing for three di- in. Previous studies (e.g., [26]) clarified that soon after rections is identical, the shape of each AMR domain is a the onset of the merger, either a long-lived HMNS or a half cube. We chose N = 60 for the best resolved runs black hole is formed. For most of the simulations in this in this work, and all the results shown in the following paper performed with stiff EOSs and with the canoni- were obtained with this resolution. We also performed cal total mass 2.6 – 2.8M , we found that a long-lived (cid:12) simulations with N = 40 and 50 (or 48) for several cho- HMNS is formed with its lifetime much longer than its sen models to check the convergence of the results (see dynamical timescale ∼ 0.1 ms and its rotation period Appendix A). ∼1 ms; the lifetime is longer than 10 ms for most of the We classify the domains of the AMR algorithm into models employed in this paper. In this section, we pay two categories: one is a coarser domain, which covers a attention to the case that such a HMNS is formed. wide region including both neutron stars with its origin Figures 3 – 5 display snapshots of the density profiles fixed at the approximate center of mass throughout the in the merger for models APR4-135135, APR4-120150, simulation. Theotherisafinerdomain,twosetsofwhich and H4-120150, respectively. Figure 6 also displays the comove with two neutron stars and cover the region in central density as a function of time for the models with their vicinity. We denote the side length of the largest m = m = 1.35M (left), and m = 1.2M and 1 2 (cid:12) 1 (cid:12) domain, number of the coarser domains, and number of m =1.5M (right). These figures show that a compact 2 (cid:12) the finer domains by 2L, l , and 2l , respectively. In and nonaxisymmetric object (proto HMNS) is formed in c f this work, l = 5 and l = 4 (in total, 13). The grid thecentralregionsoonaftertheonsetofthemerger. The c f spacing for each domain is h = L/(2lN), where l = 0 – shapeandcompactnessoftheHMNSdependstronglyon l l (= l +l −1) is the depth of each domain. In the the EOS and mass ratio; e.g., the presence of the asym- max c f following, we denote L/2lmax by Lmin and hlmax by ∆x. metricspiralarmsfoundinthetoppanelsofFigs.4and5 TableIIIsummarizestheparametersofthegridstruc- isthefeatureonlyfortheasymmetricbinaries;theampli- ture for the simulations. For all the simulations, L is set tude of the quasiradial oscillation is larger for the equal- tobeL/c>∼10ms. Thisimpliesthatthematerialcannot mass binaries; a high-amplitude quasiradial oscillation escape from the computational domain in ∼10 ms after is a unique property found only for models with APR4 the onset of the merger, even if it could move with the (seeFig.6). However,itisuniversalthattheHMNSsare speedoflight. Inreality,thespeedofmostoftheejected rapidlyrotatingandnonaxisymmetric,irrespectiveofthe material is smaller than ∼ 0.5c, and thus, the material EOS, total mass (m≤2.8M ), and mass ratio, as found (cid:12) staysinthesecondcoarsestlevelformorethan10ms. L in previous studies [21, 25, 26]. This rapid rotation to- isalsomuchlargerthanthegravitationalwavelengthsat gether with the nonaxisymmetric configuration not only the initial instant λ := π/Ω . This implies that a spu- results in the emission of strong gravitational waves but 0 0 rious effect caused by outer boundaries when extracting also is the key for an efficient mechanism of angular mo- gravitationalwavesisexcludedinthepresentworkmore mentum transport from the HMNS to the surrounding efficiently than in the previous works. The semi-major material because the HMNS exerts the torque. 9 FIG. 3: Snapshots of the density profile for the merger of binary neutron stars for an equal-mass model APR4-135135. The firstrowshowsthedensityprofilesintheequatorialplaneandinthecentralregion,andsecond–fourthonesshowthedensity profile for a wide region in the x-y, x-z, and y-z planes. t ≈11.3 ms for this model. merge 10 FIG. 4: The same as Fig. 3, but for unequal-mass model APR4-120150. t ≈10.3 ms for this model. merge