Black hole-neutron star binary merger: Dependence on black hole spin orientation and equations of state
Kyohei Kawaguchi, Koutarou Kyutoku, Hiroyuki Nakano, Hirotada Okawa, Masaru Shibata, Keisuke Taniguchi
BBlack hole-neutron star binary merger: Dependence on black hole spin orientationand equations of state
Kyohei Kawaguchi, Koutarou Kyutoku, Hiroyuki Nakano,
3, 4
Hirotada Okawa,
1, 5
Masaru Shibata, and Keisuke Taniguchi
6, 7 Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Interdisciplinary Theoretical Science (iTHES) Research Group, RIKEN, Wako, Saitama 351-0198, Japan Department of Physics, Kyoto University, Kyoto 606-8502, Japan Center for Computational Relativity and Gravitation,and School of Mathematical Sciences, Rochester Institute of Technology,85 Lomb Memorial Drive, Rochester, New York 14623 Advanced Research Institute for Science and Engineering,Waseda University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan Department of Physics, University of the Ryukyus,1 Senbaru, Nishihara, Okinawa 903-0213, Japan Graduate School of Arts and Sciences, University of Tokyo, Komaba, Meguro, Tokyo 153-8902, Japan (Dated: October 8, 2018)We systematically performed numerical-relativity simulations for black hole (BH) - neutron star(NS) binary mergers with a variety of the BH spin orientation and nuclear-theory-based equations ofstate (EOS) of the NS. The initial misalignment angles of the BH spin measured from the directionof the orbital angular momentum are chosen in the range of i tilt , ≈ ◦ –90 ◦ . We employed fourmodels of nuclear-theory-based zero-temperature EOS for the NS with which the compactness ofthe NS is in the range of C = M NS /R NS = 0 . . M NS and R NS are the mass andthe radius of the NS, respectively. The mass ratio of the BH to the NS, Q = M BH /M NS , and thedimensionless spin parameter of the BH, χ , are chosen to be Q = 5 and χ = 0 .
75, together with M NS = 1 . M (cid:12) so that the BH spin misalignment has a significant effect on tidal disruption of theNS. We obtain the following results: (i) The inclination angle of i tilt , < ◦ and i tilt , < ◦ arerequired for the formation of a remnant disk with its mass larger than 0 . M (cid:12) for the case C = 0 . C = 0 . . M (cid:12) for C (cid:38) . . M (cid:12) is obtained for i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . g / cm is approximately aligned with theremnant BH spin for i tilt , ≈ ◦ . On the other hand, the disk axis is misaligned initially with ∼ ◦ for i tilt , ≈ ◦ , and the alignment with the remnant BH spin is achieved at ∼ ∼
100 ms and dependsonly weakly on the misalignment angle and the EOS. (iii) The ejecta velocity is typically ∼ . . c and depends only weakly on the misalignment angle and the EOS of the NS, while the morphologyof the ejecta depends on its mass. (iv) The gravitational-wave spectra contains the information ofthe NS compactness in the cutoff frequency for i tilt , (cid:46) ◦ . PACS numbers:
I. INTRODUCTION
The ground-based gravitational-wave detectors, suchas Advanced LIGO [1], Advanced VIRGO [2], and KA-GRA [3], will detect signals in the next decade. Blackhole (BH) - neutron star (NS) binary mergers are one ofthe most promising gravitational-wave sources for thesedetectors [4]. Since NSs and BHs are compact objectswith strong gravity, detection of gravitational waves fromcompact binary mergers, including BH-NS mergers, canbe a touchstone of the theory of gravity [5]. The gravi-tational waveforms from BH-NS mergers depend on theBH mass, the BH spin, the NS mass, and the NS ra-dius, and thus they will carry information of these bi-nary parameters [6, 7]. In particular, the information ofthe NS radius could be used for constraining the equa-tion of state (EOS) and makes a big contribution notonly to astrophysics but also to nuclear and fundamental physics [8–11]. Since the expected signal-to-noise ratioin the gravitational-wave detection is not very high, tak-ing a cross correlation between the observational dataand theoretical template is an important method to ex-tract the information from the detected signal at leastfor the next generation gravitational-wave detectors [12],and hence waveforms for a variety of binary parametershave to be derived.BH-NS mergers are also proposed as a potential pro-genitor of short-hard gamma-ray-bursts (sGRBs) [13].If the NS is tidally disrupted during the merger, a hotand massive disk with mass (cid:38) . M (cid:12) could be formedaround the remnant spinning BH. This BH-disk systemcould launch a relativistic jet by releasing its gravita-tional energy through the neutrino emission or an elec-tromagnetic energy flow in a short time scale (cid:46) a r X i v : . [ a s t r o - ph . H E ] J un for explaining its duration and the estimated event rates(see Refs. [14, 15] for reviews and references therein).For extracting physical information of the central engineof sGRB from electromagnetic observations, we have toclarify the formation process of the BH-disk system andthe relation between this system and sGRB.A fraction of the NS material would be ejected dur-ing tidal disruption [16, 17]. Since the NS consists ofhighly neutron-rich matter, r-process nucleosynthesis isexpected to take place in the ejecta [18], and the emis-sion powered by decay of the radioactive nuclei would oc-cur (kilonova/macronova) [19–21]. This electromagneticcounterpart of the binary merger is useful for determin-ing the position of the source when gravitational wavesare detected. Also, their light curves will reflect the bi-nary parameters, and could be useful for extracting thephysical information of the binary. Furthermore, the r-process nuclei produced in the ejecta are considered tocontribute to the chemical evolution history of the uni-verse [18, 22, 23].The property of disks and ejecta, particularly theirmass, is essential for predicting the electromagnetic coun-terparts, such as sGRB and kilonova. The mass of thematter that remains outside the BH after the merger de-pends strongly on whether and at which orbital separa-tion the NS is disrupted by the tidal force of the com-panion BH [24]. If the tidal disruption occurs at a suf-ficiently distant orbit from the innermost stable circularorbit (ISCO) of the BH, an appreciable amount of masswould remain outside the remnant BH. On the otherhand, if the tidal disruption does not occur or occursnear the ISCO of the BH, entire NS material would beswallowed by the BH. The orbital separation at whichthe tidal disruption occurs depends on binary parame-ters, particularly on the mass ratio of the binary, the NScompactness, and the BH spin. Thus, we have to clar-ify the dependence of tidal disruption processes on theseparameters in order to predict the mass of the remnantdisk and the ejecta.The orbital angular momentum and BH spin can bemisaligned. A population synthesis study suggests thatabout a half of BH-NS binaries can have an initial con-figuration in which the misalignment angle between theorbital angular momentum and the BH spin is larger than45 ◦ [25]. A post-Newtonian (PN) study shows that theorbit precesses due to spin-orbit coupling (or draggingof inertial frames) in the presence of misalignment [26].This could affect dynamics of BH-NS binaries not onlyin the inspiral phase but also in the merger phase.Recently, a variety of numerical-relativity (NR) sim-ulations have been performed for BH-NS binaries [6, 7,16, 17, 27–42], and quantitative dependence of the mergerprocess on binary parameters have been revealed. Someworks were done taking into account magnetic fields [34,37, 38], nuclear-theory-based EOS [16, 17, 35, 40, 42],and neutrino cooling [40, 42]. However, most of themwere done for the case that the BH spin is aligned withthe orbital angular momentum and there are only two studies for the misaligned case [36, 41]. Moreover, massejection has not yet been studied for the misaligned-spincase based on NR simulations (there is a study with anapproximate treatment of general relativity [43]).The primary purpose of this paper is to clarify thequantitative dependence of the disk formation and themass ejection from BH-NS mergers on BH spin mis-alignment and nuclear-theory-based EOS by perform-ing NR simulations systematically. We employed fournuclear-theory-based EOSs described by piecewise poly-trope with four pieces [44], and varied the misalignmentangle between the orbital angular momentum and theBH spin, while fixing the NS mass, the magnitude of theBH spin, and the mass ratio of the binary as 1 . M (cid:12) ,0 .
75, and 5, respectively.This paper is organized as follows. In Sec. II, we de-scribe methods for computing initial conditions, piece-wise polytropic EOS, and models of BH-NS binaries weemployed in this paper. In Sec. III, the formulation andthe methods of numerical simulations are summarized. InSec. IV, we present the numerical results for misaligned-spin BH-NS mergers. Finally, summaries and discussionsof this work are presented in Sec. V. Throughout this pa-per, we adopt the geometrical units in which G = c = 1,where G and c are the gravitational constant and thespeed of light, respectively. Exceptionally, c is some-times inserted for clarity when we discuss the velocityof ejecta. Our convention of notation for physically im-portant quantities is summarized in Table I. The dimen-sionless spin parameter of the BH, total mass of the sys-tem at infinite separation, mass ratio, and compactnessof the NS are defined as χ = S BH /M , Q = M BH /M NS , m = M BH + M NS , and C = M NS /R NS , respectively. J = ( J x , J y , J z ) is the total angular momentum of thesystem (see Sec. III B for its definition). The angle be-tween the BH spin and the orbital angular momentum, i tilt , is defined by i tilt := cos − (cid:18) S BH · L S BH L (cid:19) , (1)where L is the orbital angular momentum of the binary,which is defined using the total angular momentum, J ,and the BH spin, S BH , as L := J − S BH with S BH = | S BH | and L = | L | . We use i tilt , to describe i tilt at the initialcondition. Latin and Greek indices denote spatial andspacetime components, respectively. II. INITIAL CONDITIONA. Formulation and methods
We prepare quasiequilibrium states of BH-NS binariesas initial conditions of numerical simulations in a similarmanner to our previous works for aligned-spin binaries[6, 7, 47]. Numerical computations are performed by us-ing multidomain spectral method library
LORENE [46]. As
TABLE I. Our convention of notation for physically important quantities.Symbol M irr The irreducible mass of the BH M BH The gravitational mass of the BH in isolation S BH The magnitude of the BH spin angular momentum M NS The gravitational mass of the NS in isolation R NS The circumferential radius of the NS in isolation m The total mass of the system at the infinite separation M The Arnowitt-Deser-Misner mass of the system J The total angular momentum of the system Q The mass ratio M BH /M NS C The compactness parameter of the NS M NS /R NS χ The dimensionless spin parameter of the BH S BH /M i tilt The misalignment angle of the BH spin i tilt , i tilt at the initial condition gravitational radiation reaction reduces the orbital eccen-tricity [48], typical BH-NS binaries settle to a quasicir-cular orbit. For the case that the orbital separation islarge, the time scale of radiation reaction, τ GW , is muchlonger than the orbital period, P orb , because their ratiois given by [24] τ GW P orb ≈ . Q ) Q (cid:18) r m (cid:19) / . (2)For the case that the BH has a spin angular momentuminclined to the orbital angular momentum, the orbitalplane precesses due to the spin-orbit coupling effect ofgeneral relativity, and a closed orbit is not obtained evenin the absence of radiation reaction [26]. This impliesthat even in one orbital cycle, the gravitational interac-tion between two objects varies depending on the anglebetween the BH spin and the line connecting two centersof mass, and hence, the definition of the quasicircular or-bit is not trivial. At a large orbital separation, however,the precession time scale, P prec , is also longer than P orb as [26] P prec P orb ≈ . Q ) Q + 3) r m , (3)where we take the small-spin limit. Thus, we can neglectthe orbital precession as well as the gravitational radi-ation reaction for computing initial data at a large or-bital separation, and the binary can be regarded approx-imately as an equilibrium configuration in the comovingframe.We then compute quasiequilibrium states assuming thepresence of an instantaneous helical Killing vector fieldwith the orbital angular velocity Ω, ξ µ := ( ∂ t ) µ + Ω ( ∂ ϕ ) µ . (4)For nonspinning or aligned-spin BHs, this reduces to agenuine helical Killing vector as far as the radiation re-action is neglected. Accordingly, the orbital plane can be taken to be a plane perpendicular to the rotational axis.This does not hold for misaligned-spin BHs, and we donot restrict the orbital plane to be a plane perpendicularto the rotational axis. Instead, we compute initial databy requiring that neither the BH nor NS has the velocitycomponent along the rotational axis [36]. This impliesthat the binary is located at the extrema of coordinateseparation along the rotational axis, where the helicalsymmetry should hold instantaneously.The line element in the 3 + 1 form is written as ds = g µν dx µ dx ν = − α dt + γ ij (cid:0) dx i + β i dt (cid:1) (cid:0) dx j + β j dt (cid:1) , (5)where α , β i , and γ ij are the lapse function, shift vector,and three-diemensional spatial metric, respectively. Ini-tial data of the gravitational field consist of γ ij and theextrinsic curvature defined by K ij := − L n γ ij , (6)where n µ is the future-oriented timelike unit normal vec-tor field to the initial hypersurface.We employ the extended conformal thin-sandwich for-malism [49, 50] with the conformal transverse-tracelessdecomposition of Einstein’s equation [51] to compute α , β i , γ ij , and K ij . We assume the conformal flatness ofthe spatial metric, γ ij = ψ ˜ γ ij = ψ f ij , the stationar-ity of the conformal metric, ∂ t ˜ γ ij = 0, and the maximalslicing condition, K = 0 = ∂ t K . Here, f ij and K arethe flat spatial metric and the trace part of the extrinsiccurvature, K := γ ij K ij , respectively. To handle a coordi-nate singularity associated with the BH, we decomposethe conformal factor ψ and a weighted lapse functionΦ := αψ into singular and regular parts using constants M P and M Φ as [27, 47] ψ = 1 + M P r BH + φ, (7)Φ = 1 − M Φ r BH + η, (8)where r BH := | x − x P | is the coordinate distance fromthe puncture located at x P . We also decompose ˆ A ij := ψ − K ij into regular and singular parts asˆ A ij = ˜ D i W j + ˜ D j W i − f ij ˜ D k W k + K P ij , (9)where ˜ D i denotes a covariant derivative associated with f ij and the index of W i is raised/lowered by f ij . Thesingular part is given by [52] K P ij := 32 r (cid:2) ˆ x i P BH j + ˆ x j P BH i − ( f ij − ˆ x i ˆ x j )ˆ x k P BH k (cid:3) + 3 r (cid:2) (cid:15) ikl S l P ˆ x k ˆ x j + (cid:15) kjl S l P ˆ x k ˆ x k (cid:3) , (10)where ˆ x i := ( x i − x i P ) /r BH . The index of ˆ x i is alsoraised/lowered by f ij . The parameters P BH i and S i P areconstants associated with the linear momentum and spinangular momentum of the puncture, respectively (see be-low).The equations to determine φ , β i , η , and W i are de-rived by combining the Hamiltonian constatint, momen-tum constraint, ∂ t K = 0, and ∂ t ˜ γ ij = 0 as∆ φ = − πψ ρ H − ψ − ˆ A ij ˆ A ij , (11)∆ β i + 13 ˜ D i ˜ D j β j = 16 π Φ ψ j i + 2 ˆ A ij ˜ D j (cid:0) Φ ψ − (cid:1) , (12)∆ η = 2 π Φ ψ ( ρ H + 2 S ) + 78 Φ ψ − ˆ A ij ˆ A ij , (13)∆ W i + 13 ˜ D i ˜ D j W j = 8 πψ j i . (14)where ∆ := f ij ˜ D i ˜ D j . The matter source terms are de-fined by ρ H := T µν n µ n ν , (15) j i := − T µν n µ γ νi , (16) S ij := T µν γ µi γ νj , (17)with S := γ ij S ij . The asymptotic flatness gives the outerboundary conditions as φ | ∞ = β i | ∞ = η | ∞ = W i | ∞ = 0 . (18)In contrast to initial data of BH-NS binaries computedin the excision framework [53–57], we do not have to giveinner boundary conditions at the horizon. We also do nothave to give nonzero boost contributions to the shift vec-tor [36], because the Arnowitt-Deser-Misner linear mo-mentum of the system can be set to zero by choosing P BH i appropriately. Hence, BH-NS binaries with mis-aligned spins do not exhibit the center-of-mass motion inthe puncture framework.Free parameters associated with the puncture are de-termined as follows. The so-called puncture mass, M P ,is adjusted for obtaining a desired BH mass, M BH . The other mass parameter, M Φ , is determined by the con-dition that the Arnowitt-Deser-Misner mass and theKomar mass agree with each other for stationary andasymptotically flat spacetime [58, 59], that is, (cid:73) r →∞ ∂ i Φ dS i = − (cid:73) r →∞ ∂ i ψdS i = 2 πM . (19)The linear momentum of the puncture, P i BH , is deter-mined by the condition that the Arnowitt-Deser-Misnerlinear momentum of the system vanishes as P BH i = − (cid:90) j i ψ d x. (20)The spin parameter of the puncture, S i P , is given inCartesian coordinates: S i P = S P (sin i (cid:48) tilt , , cos i (cid:48) tilt ) . (21)The magnitude S P is adjusted for obtaining a desiredvalue of BH spin, S BH , measured on the horizon. Theinclination angle i (cid:48) tilt is chosen to be 30 ◦ , 60 ◦ , and 90 ◦ inthis study. Because the direction of orbital angular mo-mentum does not agree with the axis of helical symmetry, i (cid:48) tilt does not alway agree with i tilt , which is defined asthe angle between the orbital and BH spin angular mo-menta. The typical difference is 3 ◦ , and we do not adjustvalues of i (cid:48) tilt to control i tilt in this study.The NS matter is assumed to be composed of an idealfluid. The energy-momentum tensor is given by T µν := ρhu µ u ν + P g µν , (22)where ρ is the rest-mass density, P is the pressure, h = 1 + ε + P/ρ is the specific enthalpy with ε the spe-cific internal energy, and u µ is the four-velocity of thefluid. The velocity field of the fluid is expected to be ir-rotational, because the viscosity of the NS matter is low[60, 61] and the rotational periods of observed NS in com-pact binaries are not very short (see, e.g., [62]). The zerorelativistic vorticity condition, or irrotationaliry condi-tion, is written as ω µν := ( δ αµ + u α u µ )( δ βν + u β u ν ) ( ∇ α u β − ∇ β u α )= h − [ ∇ µ ( hu ν ) − ∇ ν ( hu µ )]= 0 , (23)where the energy-momentum conservation and adiabac-ity condition are used for deriving the second-line expres-sion [63, 64]. This implies the presence of a velocity po-tential Ψ such that hu µ = ∇ µ Ψ, and the elliptic equationto determine Ψ is derived from the continuity equation ∇ µ ( ρu µ ) = 0 together with the helical symmetry. The ir-rotational conditions and helically symmetric conditions Although we distinguish S P from S BH , the difference betweentwo is at most O (10 − ). of the specific momentum, L ξ ( hu µ ) = 0, are combined togive hξ µ u µ = − C (= const . ) , (24)which is used to determine h . The integration constant C is determined by the condition that the baryon restmass of the NS takes a desired value. As we explain inthe next section, the specific enthalpy determines all theother thermodynamical quantities in the computation ofinitial data.The relative location of each component of the binaryis determined as follows. We fix the binary separation inthe direction perpendicular to the rotational axis, whichis chosen to be the z -axis. The centers of BH and NSare put on the xz -plane, but now they are not limited tothe x -axis. The orbital angular velocity of the binaryis determined by the force-balance condition that the NSdoes not move perpendicular to the rotational axis, andthis amounts to requiring dh/dx = 0 in our coordinates.The location of the rotational axis with respect to the bi-nary components is determined by the condition that themagnitude of the orbital angular momentum agrees withthe value derived by the third-and-a-half post-Newtonianformulas for a given value of the orbital angular velocity(see Appendix D of [65]). Finally, the binary separationalong the rotational axis is determined by the conditionthat the NS has no velocity component along the rota-tional axis, and this amounts to requiring dh/dz = 0 inour coordinates [36].
B. Piecewise polytropic equations of state
Since the cooling time scale of NS is much shorter thanthe lifetime of typical compact binaries, we can employa zero-temperature EOS for the NS just before binarymergers [66]. Employing a zero-temperature EOS, thethermodynamical quantities, such as P , ε , and h , can bedescribed as functions of ρ as P = P ( ρ ) , ε = ε ( ρ ) , h = h ( ρ ) . (25)From the first-law of thermodynamics, these quantitiessatisfy relations, dε = Pρ dρ, (26) dh = 1 ρ dP, (27)which determine ε and h from given P ( ρ ), respectively. We can also fix the locations of both components by forcing themto be on the x -axis. This method respects the instantaneoushelical symmetry as well as the method adopted in this study.Taking the fact that no criteria are available to determine whichcondition gives superior initial data, we simply follow the methodadopted in [36]. In this work, we employ a piecewise polytropicEOS [44] to describe a zero-temperature EOS of theNS. This is a phenomenologically parametrized EOS,which reproduces a zero-temperature nuclear-theory-based EOS at high density only with a small numberof polytropic constants and indices as P ( ρ ) = κ i ρ Γ i for ρ i − ≤ ρ < ρ i (1 ≤ i ≤ n ) , (28)where n is the number of the pieces used to parametrizean EOS. ρ i is the rest-mass density at a boundary of twoneighboring i -th and ( i + 1)-th pieces, κ i is the i -th poly-tropic constant, and Γ i is the i -th adiabatic index. Note,here, ρ = 0 and ρ n → ∞ . Requiring the continuityof the pressure, κ i ρ Γ i i = κ i +1 ρ Γ i +1 i (1 ≤ i ≤ n − κ , ρ i and Γ i . ε is determined by integrating Eq. (26) with the integra-tion constant, ε (0) = 0. It was shown that the piecewisepolytropic EOS with four pieces reproduces the nuclear-theory-based EOSs within ∼
5% errors in pressure forthe nuclear density range [44].Table II lists the EOSs which we employ in our study.We employ the models of the NS EOS which can realizethe NS with M NS ≈ M (cid:12) which satisfies the recent ob-servational constraint [67, 68]. For these models, the NSradius is in the range ∼ −
15 km for M NS = 1 . M (cid:12) ,which is largely consistent with the recent theoretical andobservational suggestion [69, 70]. Following [6, 7, 44] ,we always fix the parameters of EOS in the subnuclear-density region asΓ = 1 . , (29) κ /c = 3 . × − (cid:0) g / cm (cid:1) − Γ , (30)and we set ρ = 10 . g / cm and ρ = 10 g / cm . Here,we insert c for clarity. Instead of giving ρ , we give P inTable II for each EOS, which is the pressure at ρ = ρ . C. Models
As we already mentioned, we choose i tilt , ≈ ◦ , ◦ ,and 90 ◦ . We employ four different piecewise polytropicEOSs, APR4, ALF2, H4, and MS1, for each value of i tilt , . On the other hand, we set the NS mass, M NS , themass ratio, Q , and dimensionless spin paramter, χ , to befixed values ( M NS , Q, χ ) = (1 . , , . normalized by the total mass is set to be m Ω = 0 . i tilt , .Specifically, “i30”, “i60”, and “i90” denote the modelswith i tilt , ≈ ◦ , ◦ , and 90 ◦ , respectively. For all themodels, we rotate the initial data before we start thesimulation so that the initial direction of the total angularmomentum agrees with the direction of the z -axis. TABLE II. The key quantities for piecewise polytropic EOSs [44] which we employ in this paper. P is the pressure at ρ = ρ shown in the unit of dyne / cm , Γ i is the adiabatic index for each piecewise polytrope, and M max is the maximum mass of thespherical NS for a given EOS. R . , ρ . , M ∗ , . , and C . are the radius, the central rest-mass density, the baryon rest mass,and the compactness parameter for the NS with M NS = 1 . M (cid:12) , respectively.Model log P Γ Γ Γ M max [ M (cid:12) ] R . [km] ρ . [g / cm ] M ∗ , . [ M (cid:12) ] C . APR4 34.269 2.830 3.445 3.348 2.20 11.1 8.9 × × × × i tilt , ), the ADM mass ( M ), and the total angu-lar momentum ( J ) , respectively. Note that M NS = 1 . M (cid:12) and m = 8 . M (cid:12) Model EOS i tilt , [ ◦ ] M [ M (cid:12) ] J [ GM (cid:12) /c ]APR4i30 APR4 33 8.04 63APR4i60 APR4 63 8.05 57APR4i90 APR4 94 8.05 47ALF2i30 ALF2 33 8.04 63ALF2i60 ALF2 63 8.05 57ALF2i90 ALF2 94 8.05 47H4i30 H4 33 8.04 63H4i60 H4 63 8.05 57H4i90 H4 94 8.05 47MS1i30 MS1 32 8.04 63MS1i60 MS1 63 8.05 57MS1i90 MS1 93 8.05 48 III. METHODS OF SIMULATIONS
Numerical simulations are performed using an adap-tive mesh refinement (AMR) code
SACRA [72]. Here, weemploy a Baumgarte-Shapiro-Shibata-Nakamura (BSSN)formulation partially incorporating Z4c prescription [73].The gauge conditions, the numerical scheme, and thediagnostics are essentially the same as those describedin [6, 7].
A. Formulation and Numerical Methods
We solve Einstein’s evolution equation partially incor-porating the Z4c formulation (see [74] for our prescrip-tion) with a moving puncture gauge. The Z4c formula-tion is a modified version of the BSSN-puncture formu-lation, introducing a new variable. In the BSSN formu-lation, the long-term simulation causes a slow accumula-tion of the numerical error that leads to gradual violationof the constraints, which should be zero if the evolutionequations are solved exactly. This accumulation comesfrom the fact that in the BSSN formulation, the evolu- tion equations of the constraints have a non-propagatingmode. In the Z4c formulation, by introducing a new vari-able, Θ, the evolution equation of the constraints changesentirely to a wave equation. Therefore, the violation ofthe constraints propagates away, and the local accumu-lation of the numerical error can be suppressed.In the Z4c formulation, we evolve the conformal factor, W := γ − / , the conformal three-metric, ˜ γ ij = γ − / γ ij ,a variable slightly modified from the trace of the extrinsiccurvature, ˆ K := K − A ij = γ − / ( K ij − Kγ ij / i , and the new variable, Θ. Here, γ := det γ ij . The evolution equations are written as, (cid:0) ∂ t − β i ∂ i (cid:1) W = 13 W (cid:104) α (cid:16) ˆ K + 2Θ (cid:17) − ∂ i β i (cid:105) , (31) (cid:0) ∂ t − β k ∂ k (cid:1) ˜ γ ij = − α ˜ A ij + ˜ γ ik ∂ j β k + ˜ γ jk ∂ i β k −
23 ˜ γ ij ∂ k β k , (32) (cid:0) ∂ t − β i ∂ i (cid:1) ˆ K = − D i D i α + α (cid:20) ˜ A ij ˜ A ij + 13 (cid:16) ˆ K + 2Θ (cid:17) + 4 π ( ρ H + S )]+ ακ (1 − κ ) Θ , (33) (cid:0) ∂ t − β j ∂ j (cid:1) Θ = (cid:20) α (cid:26) R − ˜ A ij ˜ A ij + 23 (cid:16) ˆ K + 2Θ (cid:17) (cid:27) − α { πρ H + κ (2 + κ ) Θ } (cid:21) e − ( r/r ) , (34) (cid:0) ∂ t − β k ∂ k (cid:1) ˜ A ij = − W (cid:18) D i D j α − γ ij D k D k α (cid:19) + W α (cid:18) R ij − γ ij R (cid:19) + α (cid:104)(cid:16) ˆ K + 2Θ (cid:17) ˜ A ij − A ik ˜ A kj (cid:105) − πW α (cid:18) S ij − γ ij S (cid:19) + ˜ A ik ∂ j β k + ˜ A jk ∂ i β k −
23 ˜ A ij ∂ k β k , (35) (cid:0) ∂ t − β j ∂ j (cid:1) ˜Γ i = − A ij ∂ j α + 2 α (cid:20) ˜Γ ijk ˜ A jk − W ˜ A ij ∂ j W −
13 ˜ γ ij ∂ j (cid:16) K + Θ (cid:17) − π ˜ γ ij j j (cid:21) + ˜ γ jk ∂ j ∂ k β i + 13 ˜ γ ij ∂ j ∂ k β k − ˜Γ k d ∂ k β i + 23 ˜Γ id ∂ k β k − ακ (cid:16) ˜Γ i − ˜Γ i d (cid:17) , (36)where D i denotes a covariant derivative associated with γ ij , ˜Γ i d = − ∂ j ˜ γ ij , and κ and κ are coefficients asso-ciated with the constraint damping. An overall factor, e − ( r/r ) , is multiplied in the right-hand side of Eq. (34)so that Θ plays a role only in the inner region of thesimulation box. In our simulation, we set κ = κ = 0,and r = L/
2, where L is the size of the computationaldomain on one side (see Table. IV). The spatial deriva-tives in the evolution equations are evaluated by fourth-order centered finite differencing except for the advectionterms, which are evaluated by fourth-order upwind finitedifferencing. A fourth-order Runge-Kutta scheme is em-ployed for the time evolution.Following [77], we employ a moving-puncture gauge inthe form (cid:0) ∂ t − β i ∂ i (cid:1) α = − αK, (37) (cid:0) ∂ t − β j ∂ j (cid:1) β i = 34 B i , (38) (cid:0) ∂ t − β j ∂ j (cid:1) B i = (cid:0) ∂ t − β j ∂ j (cid:1) ˜Γ i − η B B i , (39)where B i is an auxiliary variable, and η B is a coefficientintroduced to suppress a strong oscillation of the shiftvector. In this work, we set η B = 0 . /M (cid:12) .The EOS is basically the same as those described in [6,7]: We decompose the pressure and the specific internalenergy into a zero-temperature part and a thermal partas P = P cold + P th , ε = ε cold + ε th . (40)Here, P cold and ε cold are functions determined by thepiecewise polytropic EOS. Then, the thermal part of thespecific internal energy is calculated by ε th = ε − ε cold ,where ε is given from the hydrodynamics. Finally, wedetermine the thermal part of the pressure using a simpleΓ-law, ideal-gas EOS as P th = (Γ th − ρε th , (41)where Γ th is an adiabatic index for the thermal part,which we choose Γ th = 1 .
8. As is discussed in AppendixA. of [17], the difference of the disk mass and the ejectamass among the different values of Γ th is expected tobe small compared to the numerical errors due to finitegridding, thus, we do not study the dependence on Γ th in this paper.Since the vacuum is not allowed in any conservativehydrodynamics scheme (see [72] for details), we put anartificial atmosphere of small density outside the NS. The atmosphere density is set to be ρ atm = 10 − ρ max ∼ g / cm for the inner region, r < R crit , and ρ atm =10 − ρ max ( R crit /r ) for the outer region, r ≥ R crit , with R crit ≈ L/
16. The total rest mass of the atmosphereis always less than 10 − M (cid:12) , and hence we can safelyneglect the effect of the artificial atmosphere as far asappreciable tidal disruption occurs. B. Diagnostics
We estimate the mass and the spin angular momentumof BH with misaligned spin assuming that the deviationfrom Kerr spacetime is negligible in the vicinity of theBH, at least for the case that the separation of the binaryis large or the system is approximately regarded as asteady state. In a steady state with stationary slicing, theevent horizon (EH) agrees with apparent horizon (AH).Thus, the irreducible mass of the BH is determined bythe area of the AH, A AH , as M irr =: (cid:114) A AH π . (42)In Kerr spacetime, a relation M = 2 M (cid:16) − (cid:112) − χ (cid:17) χ (43)holds between the gravitational and irreducible massesof the BH. If the irreducible mass and the BH spin areknown, we can calculate the BH mass using this relation.To determine the spin angular momentum, we use the re-lation between the spin and the intrinsic scalar curvatureof the horizon, R (2)EH , following [78]. In Kerr spacetime, R (2)EH is written as R (2)EH ( θ ) = 2 (cid:0) ˆ r + χ (cid:1) (cid:0) ˆ r − χ cos θ (cid:1) M (cid:0) ˆ r + χ cos θ (cid:1) . (44)Here, ˆ r + = 1 + (cid:112) − χ is a normalized radius of theEH, and θ is the latitude in Boyer-Lindquist coordinates.The minimum and maximum values of R (2)EH at θ = 0and π/
2, i.e., at the pole and the equator of the BH are,respectively, R (2)min = − (cid:112) − χ M , (45) R (2)max = − (cid:16) − χ + 2 (cid:112) − χ (cid:17) M χ . (46)Solving these equations with respect to χ , we have χ = 1 − (cid:18)
12 + M R (2)min (cid:19) , (47) χ = − (cid:113) M R (2)max M R (2)max . (48) TABLE IV. Setups of the grid structure for the simulationwith our AMR algorithm. ∆ x is the grid spacing at the finest-resolution domain. R diam is the semimajor diameter of the NSin the direction perpendicular to the axis of helical symmetry. L is a half of the edge length of the largest domain. λ = π/ Ω is the gravitational wavelength of the initial configuration.Model ∆ x/M R diam / ∆ x L/λ L (km)APR4i30 0.0134 101 2.357 2444APR4i60 0.0133 102 2.343 2429APR4i90 0.0136 100 2.398 2486ALF2i30 0.0153 102 2.687 2786ALF2i60 0.0153 102 2.687 2786ALF2i90 0.0156 101 2.729 2829H4i30 0.0172 102 3.032 3144H4i60 0.0172 102 3.032 3144H4i90 0.0173 101 3.046 3158MS1i30 0.0186 102 3.273 3394MS1i60 0.0188 101 3.308 3430MS1i90 0.0188 101 3.308 3430 Using these relations, we can estimate the BH spin ap-proximately from the intrinsic scalar curvature of the AH.The direction of the BH spin can be determined from thelocation of R (2)min . We note that the direction of the BHspin is not gauge invariant in this definition. However,we expect that it gives a reasonable measure of the spindirection if the tidal forces are negligible [78, 79]. For allthe models, the values of χ min and χ max agree with eachother up to significant ( ∼
3) digits . Thus, we use thevalue of χ max for calculating M BH and S BH in this paper.We define the total angular momentum of the system, J = ( J x , J y , J z ), from a rotational invariance of the grav-itational Hamiltonian at spatial infinity as [45], J i := 18 π (cid:15) ijk (cid:73) r →∞ x j (cid:0) K lk − Kγ lk (cid:1) dS l , (49)where (cid:15) ijk is the spatial Levi-Civita tensor.Finally, we define the time at which the binary merges.For this purpose, we define the rest mass inside the AHas M ≤ AH := (cid:90) r ≤ r AH ρ ∗ d x, (50)where, r AH = r AH ( θ, ϕ ) is the radius of the AH as a func-tion of the angular coordinates, and ρ ∗ = ραu t √ γ is theconserved rest-mass density. Then, we define the mergertime, t merge , as the time at which M ≤ AH ≥ − M (cid:12) isachieved. C. Setups for AMR grids In SACRA , the Einstein and hydrodynamical equationsare solved in an AMR algorithm described in [72]. Here,we briefly describe the settings for AMR grids, and thedetails are found in [72]. In this work, we prepare nine refinement levels with different grid resolutions anddomain sizes. Each domain is composed of the uni-form vertex-centered cubic grid with the grid number(2 N + 1 , N + 1 , N + 1) for ( x, y, z ). We always chose N = 60 for the best resolved runs in this work. Wealso performed simulations with N = 40 and 48 to checkthe convergence of the result. As described in [72], theAMR domains are classified into two categories: One isthe coarser domains which cover wider regions with theirorigin fixed approximately at the center of the mass of thesystem. The other is the finer domains which cover theregions around the BH or the NS and comove with it. Weset four coarser domains and five pairs of finer domainsfor all the simulations which we performed in this paper.The grid spacing for each domain is h l = L/ (2 l N ), where2 L is the edge length of the largest cubic domain and l is the depth of the domain.Table IV summarizes the parameters of the grid struc-ture for the simulations. In all the simulations, the semi-major diameter of the NS in the direction perpendicularto the axis of helical symmetry is covered with ≈
100 gridpoints. For N = 60, the total memory required for thesimulation with 14 domains is about 35 GBytes. We per-form all the simulations using personal computers of 64–128 GBytes memory and 6–24 processors with OpenMPlibrary. The typical computational time required to per-form one simulation is 9 weeks for the 24 processors case. IV. NUMERICAL RESULTS
We present numerical results of our simulations in thissection.
A. Orbital Evolution
Figure 1 plots bird’s-eye views for the evolution ofthe coordinate separation x i sep := x i NS − x i BH for mod-els H4i30, H4i60, and H4i90. Here x i NS denotes the lo-cation of the maximum rest-mass density, and x i BH isthe location of the puncture. This figure shows that thenumber of the orbits before the merger decreases withthe increase of i tilt , as ≈ . , .
5, and 5 .
5, respectivelyfor models H4i30, H4i60, and H4i90. This dependencestems primarily from the general relativistic spin-orbitinteraction, which is well known for inducing an “orbitalhung up” effect [26, 75–77]. The additional energy of thespin-orbit interaction is written, in the leading order, as E SO = 1 r M NS M BH L · S BH , (51)and hence the spin-orbit interaction essentially weak-ens the attractive force of gravity if L · S BH >
0. Forthis situation, the angular velocity of the binary, Ω,decreases, and so does the luminosity of gravitationalwaves, which is proportional to Ω / . The reduction of -90-60-30 0 30 60 90-90-60-30 0 30 60 90-90-60-30 0 30 60 90z(km) H4i30 x(km)y(km)z(km) -90-60-30 0 30 60 90-90-60-30 0 30 60 90-90-60-30 0 30 60 90z(km) H4i60 x(km)y(km)z(km) -90-60-30 0 30 60 90-90-60-30 0 30 60 90-90-60-30 0 30 60 90z(km) H4i90 x(km)y(km)z(km) FIG. 1. Evolution of the orbital separation x i sep := x i NS − x i BH of binaries with ( Q, M NS , χ ) = (5 , . M (cid:12) , . i t il t [ deg r ee ] t-t merge [ms]H4i30H4i60H4i90 FIG. 2. Evolution of i tilt for the models with H4. the gravitational-wave luminosity makes the approachingvelocity smaller, and thus the time to merger becomeslonger. Since L · S BH ∝ cos i tilt , this effect can be signifi-cant when the BH spin is aligned with the orbital angularmomentum, and thus, the binary with a small value of i tilt , merges later.The figure also illustrates that the orbits of the bina-ries are precessing. This is also primarily due to the spin-orbit interaction. For models H4i30, H4i60, and H4i90,the elevation angles of the orbits measured from the xy -plane are always ≈ ◦ , ◦ , and 45 ◦ , respectively. Wenote that gravitational waves are radiated primarily tothe direction of the orbital angular momentum, and thedirection of the total angular momentum J changes dur-ing the inspiral phase due to the gravitational radiationreaction. However, the angle between J and the z -axis isalways smaller than 5 ◦ .It is known that, at least approximately, i tilt is a con-stant of motion [26] for the case that the spin of the NSis absent. This feature is also seen in our simulation:see Fig. 2, in which we plot the time evolution of i tilt . φ p r e c [ deg r ee ] t[ms]H4i30H4i60H4i90H4i60(1.5PN) FIG. 3. Evolution of the precession angle, φ prec , for the mod-els with H4. Irrespective of i tilt , and EOS, we indeed find that i tilt approximately keeps their initial values and their fluctu-ation is smaller than ≈ ◦ irrespective of models. There-fore i tilt would be regarded approximately as the valuewhich is determined when the BH-NS binary was born,even just before the merger. Also this property ensuresthat setting of the simulation models is well-defined; themodels with different values of i tilt , describes entirelydifferent physical systems.Next, we analyze the evolution of the precession angle.We define the precession angle of the orbit as φ prec := cos − (cid:20) ( L × S BH ) · ¯ y | L × S BH | (cid:21) . (52)Here, ¯ y is the coordinate basis of the y -axis. The evolu-tion of φ prec is shown in Fig. 3 up to the time of merger.We also plot the value obtained by integrating the leadingpost-Newtonian formula of the precession angular veloc-0ity [26], ω PNprec = | J | r (cid:0) Q − + 4 (cid:1) , (53)using instantaneous values of J and r for the simulationof the model H4i60.For models H4i30, H4i60 and H4i90, the final valuesof φ prec are ≈ ◦ , ◦ , and 145 ◦ , respectively, whilethe number of the orbits are ∼ . .
5. We find that theprecession angular velocity ω prec computed by the timederivative of φ prec is always smaller than orbital angu-lar velocity by an order of magnitude. The evolution of φ prec agrees quantitatively with the one calculated withthe leading-order post-Newtonian formula, despite thedifference in the gauge condition. Figure 3 shows thatthe final value of φ prec is smaller for the case that i tilt , is larger. It is simply because a longer inspiral phase isrealized for binaries of a smaller value of i tilt , . B. Tidal disruption
Figure 4 shows the rest-mass density together withthe location of apparent horizons for models H4i30 andH4i60. The images are generated using a volume ren-dering method truncating the density below 10 g / cm ,and the color on the AH surfaces describes the value of2-dimensional Ricci scalar on it (cf. Sec. III B).For the model H4i30, the NS is tidally disrupted form-ing a one-armed tidal tail around the BH. An efficientangular momentum transport process works, and a frac-tion of the NS material becomes gravitationally unboundduring this phase. At a few milliseconds after the onsetof the tidal disruption, a large fraction of the NS materialis swallowed by the BH, and ≈
20% of the NS materialremains outside the BH. The inner part of the tidal tail issubsequently wound around the BH, and a disk with itsradius ≈
150 km is formed. Also, some material, whichwas not able to get enough kinetic energy to escape fromthe system, falls back to the disk continuously.Initially, the tidal tail around the BH is tilted withthe angle less than ≈ ◦ , which reflects the elevationangle of the orbit. Also, the tidal tail is slightly warp-ing due to the spin-orbit interaction. However, the diskformed finally is nearly aligned with the BH spin, and themorphology of the disk for this model resembles the diskformed with aligned-spin BH-NS mergers. The remnantBH spin axis agrees approximately with the direction ofthe initial total angular momentum, i.e., the z -axis.For the model H4i60, general features of the tidal dis-ruption of the NS is similar to the model H4i30, whilethe precession of the tidal tail is more appreciable forthis model. More than ≈
90% of the NS material fallsinto the BH in a few milliseconds after the tidal disrup-tion for this model. This is a result of the fact thatthe tidal disruption occurs in the vicinity of the ISCO ofthe BH. The elevation angle of the tidal tail measuredfrom the xy -plane is different for each part, and it is pointed out in [41] that this may prevent the collision ofthe tidal tail. However, the tidal tail still collides withitself to form a disk or torus for this model. This mightbe due to the difference of the initial parameters of thebinary, that a binary with larger mass ratio and dimen-sionless spin parameter are employed in [41], for whichlarger elevation angle is expected to be achieved. More-over, we performed the simulation for a longer time afterthe merger than in [41], and this also makes more chancesfor the tidal tail to collide with itself. The disk appearsto be misaligned with the BH spin by ≈ ◦ at least at ≈
10 ms after the tidal disruption, while the BH spin isaligned approximately with the z -axis.For the model H4i90, the NS is tidally disrupted veryweakly by the BH and a tiny tidal tail is formed. Sincethe tidal disruption occurs at a close orbit to the ISCOor perhaps inside the ISCO of the BH, most of the NSmaterial is swallowed by the BH. Thus, only a tiny ac-cretion disk is formed and the total amount of the ejectais not appreciable for this model.Figure 5 plots the time evolution of the total rest massoutside the BH defined by M > AH := (cid:90) r>r AH ρ ∗ d x, (54)for a variety of EOSs and i tilt , . From the comparisonamong the models with the same value of i tilt , , we findthat the values of M > AH increases in the order of APR4,ALF2, H4, and MS1, i.e., in the order of the compactness.This result shows that the tidal disruption occurs at amore distant orbit for the case that the compactness ofthe NS is small. This dependence of the tidal disruptionon the compactness is the same as the dependence whichwas found in the study on BH-NS mergers with alignedspins [6, 7].For a fixed EOS, the values of M > AH after the mergerdecrease as the values of i tilt , increase, and thus theincrease of i tilt , prevents the tidal disruption of theNS. This result agrees qualitatively with the result ofprevious study on the misaligned-spin BH-NS binarymerger [36, 41], and is explained primarily by the reduc-tion of the spin-orbit coupling. According to the study onthe aligned-spin BH-NS merger [7], a larger mass remainsoutside the BH after the merger for the case that the BHspin is larger and parallel with the orbital angular mo-mentum. This is because the spin-orbit interaction worksas a repulsive force and the ISCO radius of the BH be-comes small for this situation. Since the spin-orbit inter-action energy is proportional to L · S BH = LS BH cos i tilt , at the leading PN order, the spin-orbit interaction isweakened as i tilt , increases. Thus, the increase of i tilt , enlarges the ISCO radius of the BH effectively and re-duces the remnant mass after the merger.We summarize the value of M > AH at ≈
10 ms afterthe onset of merger in Table V. We compare these resultswith the fitting formula for the aligned-spin case obtainedin [80]. Since the spin-orbit interaction is proportionalprimarily to L · S BH , we compare our numerical results1 FIG. 4. Snapshots of the volume rendered density map as well as the location of AHs at ≈ ≈
300 km areshown. M > A H [ M s un ] t-t merge [ms]APR4i60ALF2i60H4i60MS1i60 M > A H [ M s un ] t-t merge [ms]H4i0H4i30H4i60H4i90 FIG. 5. Evolution of the rest mass outside the apparent horizon M > AH for the models with i tilt , ≈ ◦ (left figure) and withH4 (right figure). The result for the aligned-spin case with Q = 5 are picked up from [17]. with the values derived by the fitting formula using theeffective spin parameter defined by χ eff = χ cos i tilt , . (55)The deviations of the value calculated by fitting formula∆ fit = | M fit − M > AH | /M fit are within 50% for M > AH (cid:38) . M (cid:12) and within 30% for M > AH (cid:38) . M (cid:12) . As far asthe value of M > AH is larger than 0 . M (cid:12) , M fit gives areasonable estimate for M > AH . In [41, 81, 82], the effective spin parameter is defined by a differ-ent form.
C. Disk formation and mass ejection
1. Mass of ejecta and disk
Next we evaluate the mass of disk and ejecta, whichare the key quantities for the electromagnetic emissionfrom the remnant of BH-NS mergers. We calculate theejecta mass by M eje := (cid:90) r>r AH ,u t < − ρ ∗ d x. (56)Here, we assume that the contribution of the internal en-ergy of the ejecta is negligible for estimating the unboundmaterial. Then, we define the rest mass of remnant disksby M disk := M > AH − M eje . (57)2 TABLE V. The list of M > AH , M disk , M eje , v ave , and v eje . The subscripts 10 ms and 20 ms denote the values evaluated at ≈
10 msand ≈
20 ms after the onset of merger, respectively. The results for the aligned-spin case are obtained from [17]. – implies thatwe were not able to take the data for them. For the model APR4i0, the simulation was stopped before t − t merge ≈
20 ms. Forthe model APR4i90, the mass of the disk and ejecta are so small that accurate values cannot be derived for them. For the i90models, the data for P eje ,i was not output.Model M > AH , [ M (cid:12) ] M disk , [ M (cid:12) ] M disk , [ M (cid:12) ] M eje , [ M (cid:12) ] v ave , [ c ] v eje , [ c ]APR4i0 0.068 0.059 – 8 × − × − .
30 0.057APR4i60 4 × − × − × − × − .
27 0.078APR4i90 – – – < − .
24 –ALF2i0 0.24 0.20 0.16 0.046 0.21 0.15ALF2i30 0.16 0.13 0.10 0.033 0 .
27 0.17ALF2i60 0.026 0.016 0.013 0.010 0 .
28 0.048ALF2i90 2 × − × − < − < − .
26 –H4i0 0.32 0.27 0.21 0.050 0.22 0.18H4i30 0.25 0.21 0.19 0.042 0 .
27 0.21H4i60 0.084 0.072 0.061 0.012 0 .
25 0.14H4i90 3 × − × − × − × − .
28 -MS1i0 0.36 0.28 0.22 0.079 0.24 0.19MS1i30 0.30 0.23 0.19 0.070 0 .
28 0.23MS1i60 0.18 0.14 0.11 0.041 0 .
27 0.21MS1i90 0.022 0.012 0.011 0.010 0 .
27 – M disk [M sun ] 0.2 0.15 0.1 0.05 0.01 0.14 0.15 0.16 0.17 0.18Compactness 0 10 20 30 40 50 60 70 80 90 i t il t [ deg r ee ] M e j e [ M s un ] . . . . . . . . . . C o m pa c t ne ss i t il t [ deg r ee ] . . . . . FIG. 6. The contour for M disk (left panel) and M eje (right panel) evaluated at ≈
10 ms after the onset of merger in the planeof the NS compactness C and initial value of i tilt . We note that M disk is described as M bd in [17].We list the values of M disk and M eje at ≈
10 ms and ≈
20 ms after the onset of merger in Table V. This showsthat M disk and M eje monotonically decrease with the in-crease of i tilt , . This reflects the fact that the effectiveISCO radius of the BH increases with the increase of i tilt , . We find that M disk and M eje with i tilt , (cid:38) ◦ areappreciably smaller than those for i tilt , = 0 ◦ . In partic- ular, little amount of disk and ejecta are produced for allthe EOS with i tilt , ≈ ◦ . For the moderate misalign-ment angle, i tilt , ≈ ◦ , the values of M disk and M eje are sensitive to the EOS: For the MS1 EOS, the diskwith M disk > . M (cid:12) and ejecta with M eje > − M (cid:12) are produced. On the other hand, M disk < − M (cid:12) and M eje < − M (cid:12) for the APR4 EOS.To clarify the dependence of M disk and M eje on i tilt , M d i sk [ M s un ] χ eff APR4ALF2H4MS1APR4a5ALF2a5H4a5MS1a5H4a375 M e j e [ M s un ] χ eff APR4ALF2H4MS1APR4a5ALF2a5H4a5MS1a5H4a375
FIG. 7. M disk (left panel) and M eje (right panel) evaluated at ≈
10 ms after the onset of merger as a function of an effectivespin parameter χ eff = χ cos i tilt , . The plotted lines are the linear interpolation for the results of misaligned BH-NS mergersobtained in this paper, and the plotted points are the results for aligned-spin BH-NS mergers with Q = 5 obtained in [17].“a5” and “a375” denote the aligned models with χ = 0 . χ = 0 . and EOS, we plot contours of M disk and M eje at ≈
10 msafter the onset of merger as functions of i tilt , and thecompactness parameter C of the NS in Fig. 6. The depen-dence of M disk and M eje is clear: Both of them decreasemonotonically with the increases of C and i tilt , . For amoderate value of compactness C = 0 . i tilt , shouldbe smaller than 50 ◦ for M disk to be larger than 0 . M (cid:12) .On the other hand, M disk (cid:38) . M (cid:12) even if i tilt , ≈ ◦ fora stiff EOS that realizes C = 0 . C (cid:38) . M disk is smaller than 0 . M (cid:12) for anyvalue of i tilt , . M eje > . M (cid:12) is possible for i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . M disk and M eje obtained by numeri-cal simulations for aligned-spin BH-NS mergers [17] withthose for the misaligned-spin cases. Each line describesthe results of M disk and M eje for the misaligned-spin BH-NS mergers interpolated linearly for χ eff . Each point inFig. 7 shows the results of the aligned-spin BH-NS merg-ers with the same mass ratio ( Q = 5) and the same EOSas we employed in this paper, but with smaller BH spin χ = 0 .
5. We also plot a new result for the model with Q = 5, H4 EOS, and χ = 0 . M disk and M eje , the results of the aligned-spin case agree approx-imately with the interpolated line in the error margindue to the finite grid resolution (see Appendix A), whileslightly larger mass is realized for the misaligned-spincase. The slope between H4a375 and H4a5 also agreeswith the slope of the interpolated line for H4. Excep-tionally for models with ALF2, results obtained by thealigned-spin BH-NS mergers deviate from the interpo-lated plot lines by ∼ M disk and ≈
30% for M eje , which might be a little bit larger than the errormargin due to the finite grid resolution. We note thatthe deviation depends on the interpolation method anda systematic error associated with employing χ eff .
2. Disk morphology and accretion
Figure 8 shows snapshots of the density profiles at ≈
20 ms after the onset of merger for selected models(see figure caption). We find that the size of the re-gion with ρ > g / cm is always ∼
150 km, whilethat with ρ > g / cm depends on the mass ofthe disk: For the larger disk mass, the size of the re-gion with ρ > g / cm becomes larger. The maxi-mum density of the disk also becomes larger for largerdisk mass: For instance, the maximum density exceeds10 g / cm for the model H4i30, while there is no regionwith ρ > g / cm for the model H4i60. The ratio ofthe disk height to disk radius depends only weakly on thebinary parameters, and it is always ∼ . . i tilt , ≈ ◦ , the rotational axesof the disks are approximately aligned with the z -axis.However, the density distribution on the xy cross sec-tion is not axisymmetric and an approximately station-ary spiral-shape shock wave is seen in the disks. On theother hand, for i tilt , ≈ ◦ , the disks are misaligned withthe z -axis. Also, in a similar manner to i30 models, thedisks have a non-axisymmetric structure for i tilt , ≈ ◦ .To quantify the misaligned structure of the disk, wedefine the total angular momentum of the disk J i disk as J i disk = (cid:90) r>r AH ,u t > − ρ ∗ (cid:15) ijk x j ˆ u k d x, (58)where ˆ u i = hu i , and we plot the evolution of the tilt angleof J i disk measured from the direction of the BH spin afterthe merger in Fig 9. This shows that the tilt angle of J i disk at ≈ ≈ ◦ and ≈ ◦ formodels with i tilt , ≈ ◦ and i tilt , ≈ ◦ , respectively.These tilt angles of J i disk reflect the elevation angle ofthe orbits just before the merger and gradually decrease4 (a):(b):(c):(d): FIG. 8. The density profiles of the accretion disk at ≈
20 ms after the onset of merger for selected models. The left, middle,and right columns show the plots for the xy , xz , and yz -planes, respectively. The top, second top, third top, and bottom rawsshow the plots for models ALF2i30, H4i30, H4i60, and MS1i60, respectively. d i sk t il t[ deg r ee ] t-t merge [ms] ALF2i30H4i30MS1i30APR4i60ALF2i60H4i60MS1i60 FIG. 9. Evolution of the tilt angle between J i disk and the BHspin for the models with i tilt , ≈ ◦ and 60 ◦ .FIG. 10. The density profiles on the yz cross section for themodel MS1i30 at ≈ as the system evolves. Figure 9 shows that these tiltangles appear to be slightly larger than those expectedfrom Fig. 8: We should note that the tilt angle shown inFig. 9 does not really describe representative values forthe tilt angle of the dense part of the disk but indicatesthe tilt angle of the remnant matter including the tidaltail with larger orbital radii (see below).Figure 10 shows a larger-scale density profile ( ≈ yz -plane for the model MS1i60 at ≈
40 ms after the onset of merger. While there is arelatively dense torus in the central region ( ≈
200 km),the tidal tail is widely spreading with an elevation angle ≈ ◦ . This tidal tail is not unbound although it has alarge orbital angular momentum, and thus it contributesto the tilt of J i disk which we defined above.Figure 11 shows the density profiles of the disk on the yz cross section for the model MS1i60 at ≈ ,
40, and70 ms after the onset of merger. While the dense part with ρ > g / cm is tilted by ≈ ◦ from the xy -plane at ≈
10 ms after the onset of merger, the tilt angleof the disk decreases gradually as the system evolves,and its axis is approximately aligned with the z -axis at ≈
70 ms after the onset of merger. The same featurecan also be seen on the plot of the xz -plane, and thuswe conclude that the disk has a tendency to align withthe BH spin during the evolution. The time scale forthe disk to align is ≈
50 ms and is comparable to oreven bit shorter than the time scale of the disk preces-sion t prec ∼
100 ms ( r disk /
150 km) . In the presence offluid viscosity, the so-called Bardeen-Petterson effect [87]is known as a mechanism for the disk to be aligned withthe BH spin. Since any effect of viscosity is not takenaccount, the Bardeen-Petterson effect cannot play a rolein our simulation. However, we suspect that Bardeen-Petterson-like effect induced by a purely hydrodynami-cal mechanism, such as angular momentum redistribu-tion due to a shock wave excited in a non-axisymmetricmanner of the disk, should work in the disk for the align-ment. To summarize, a disk with tilt angle of ≈ ◦ –30 ◦ would be formed for models with i tilt , ≈ ◦ . However,the dense part of the disk is subsequently aligned withthe BH spin in ≈
50 ms while the tidal tail with largeorbital radii keeps its elevation angle.Figure 12 plots the time evolution of M disk for the mod-els with the MS1 EOS and i tilt , ≈ ◦ . This figure showsthat M disk gradually decreases after a steep initial de-crease. This shows that the infall of the material intothe BH continues for a long time scale. This is inducedprimarily by a hydrodynamical process associated withthe non-axisymmetric torque in the disk, as we can inferfrom Fig. 8. Also, the fallback of the matter could givean impact to the disk material. To check that numericalviscosity is not a main source of the angular momentumredistribution, we plot the mass accretion rate of the diskmaterial into the BH for three different grid resolutionsin Fig. 13. Here the accretion rate is defined as the timederivative of M ≤ AH . We find that the value of ˙ M ≤ AH de-pends only very weakly on the grid resolution. Since thenumerical viscosity should depend on the grid resolution,this result shows that the contribution of the numericalviscosity to the mass accretion is negligible, and we cansafely consider that the accretion is induced by a physicalprocess. The non-axisymmetric and dynamical feature ofthe disk would be responsible for this.Lower panels of Fig. 12 plot the time scale of the accre-tion, t acc := M disk / ˙ M ≤ AH , for the models with the MS1EOS and for i tilt , ≈ ◦ . We find that it depends onlyweakly on i tilt , and EOS and t acc ≈ t acc is asshort as ≈
200 ms at 50 ms after the onset of merger.In other words, the mass accretion rate is as large as ∼ . M (cid:12) / s for M disk = 0 . M (cid:12) even at ≈
50 ms afterthe onset of merger. This shows the importance of thenon-axisymmetric structure of the disk that is preservedfor a long time scale and governs the angular momentum6
FIG. 11. The density profiles on the yz cross section for the model MS1i60 at ≈ ,
40, and 70 ms after the onset of merger. M d i sk [ M s un ] MS1i0MS1i30MS1i60MS1i90 0 50 100 150 200 250 300 0 10 20 30 40 50 t a cc [ m s ] t-t merge [ms] M d i sk [ M s un ] APR4i30ALF2i30H4i30MS1i30 0 50 100 150 200 250 300 0 10 20 30 40 50 t a cc [ m s ] t-t merge [ms] FIG. 12. Evolution of M disk and the accretion time scale t acc = M disk / ˙ M ≤ AH for the models with MS1 (left figure) and with i tilt , ≈ ◦ (right figure). transport.
3. Ejecta morphology and velocity
In Fig. 14, we plot the image of the ejecta with volumerendering for models H4i30, H4i60, and ALF2i60. Wefind that for the case that the ejecta is massive (H4i30),it has a crescent-like shape with its opening angle ≈ ◦ .On the other hand, for the case that the ejecta has smallmass (ALF2i60), the opening angle becomes larger than360 ◦ winding around the center of mass. We also findthat the ejecta is warped due to the orbital precessionfor this case. In particular, the rear-end collision for theejecta with its opening angle larger than 360 ◦ is less pro-nounced in the misaligned-spin case than in the aligned-spin case [17].After the matter becomes gravitationally unbound, theejecta expand in an approximately self-similar manner.To analyze the feature of the ejecta, we define an average velocity of the ejecta as, v ave = (cid:115) T eje M eje , (59)where T eje is kinetic energy of the ejecta, defined follow-ing [83]. We also compute the linear momentum of theejecta defined by P eje ,i = (cid:90) ρ ∗ ˆ u i d x, (60)and calculate its magnitude by P eje = (cid:115)(cid:88) i ( P eje ,i ) . (61)Using P eje , we define the bulk velocity of the ejecta v eje by v eje = P eje M eje . (62)7 d M < A H / d t[ M s un / s ] t-t merge [ms] MS1i30N40MS1i30N48MS1i30N60 FIG. 13. Evolution of the accretion rate, ˙ M ≤ AH , to the AHfor MS1i30 but with different grid resolutions. Here, note that v eje may reflect the morphology of theejecta: The linear momentum of the ejecta vanishes ifits morphology is isotropic, while the value of v eje be-comes close to v ave if the mass is ejected coherently to aparticular direction.We summarize the values of v ave and v eje measured at10 ms after the onset of merger in Table V. This showsthat irrespective of the models, v ave ∼ . c and theirdependence on the misalignment angle and the EOS isweak. Although the magnitude of v ave is quite universal, v eje varies from model to model. v eje becomes large forthe case that the ejecta mass is large ( ≈ . M (cid:12) ). Thisis consistent with the result in Fig. 14 that the massejection proceeds in an anisotropic manner for the casethat M eje is large. On the other hand, v eje becomes smallfor the ejecta with small mass ( (cid:46) . M (cid:12) ). This is alsoconsistent with a quasi-axisymmetric ejection shown inFig. 14. D. The properties of remnant BH
We list the values of M irr , M BH , χ , S BH , and the tiltangle of the BH spin, cos − ( S z BH /S BH ), at 10 ms afterthe onset of merger in Table VI. The mass of the remnantBH increases monotonically as the value of i tilt , increasesand the EOS becomes stiff. The reason for this is thatthe tidal disruption is suppressed for larger values of i tilt , and stiffer EOSs, and more matter of the NS is swallowedby the BH.The dimensionless spin parameter decreases with theincrease of i tilt , , and it depends weakly on the EOS.In particular, for the models with i tilt , ≈ ◦ , the fi-nal value of the dimensionless spin parameter becomessmaller than the intial value of 0.75. This is because thedimensionless spin parameter is defined by S BH /M ,and the increase of M is larger than S BH for thesemodels. In fact, we can see the increase of S BH . Thevalue of cos − ( S z BH /S BH ) shows that the BH spin be- comes approximately parallel with the z -axis after themerger. This is because the direction of the total angularmomentum approximately preserves its initial direction,and the BH swallows nearly entire angular momentumof the system. The tilt angle of S BH in the final statebecomes larger for larger values of i tilt , . However, it isalways smaller than 5 ◦ .The BH spin increases because the BH swallows theNS matter of positive angular momentum. From theincrements of the BH mass ∆ M BH := M BH , f − M BH , i and the BH spin ∆ S BH := | S BH , f − S BH , pm | , we calcu-late the mean value of the specific angular momentum∆ l := ∆ S BH / ∆ M BH that the BH gained due to thefalling material. Here, M BH , i is the initial BH mass and S BH , pm is the value of the BH spin evaluated just beforethe merger, t ≈ t merge − l is expected to reflect the specific angularmomentum at the ISCO. Thus, we compare ∆ l with thespecific angular momentum of the ISCO of the BH withan effective spin parameter which we introduced in theprevious section [see Eq. (55)].In Fig. 15, we plot ∆ l/M BH as a function of χ eff . Wealso plot the specific orbital angular momentum at theISCO of the aligned-spin BH as a function of χ eff . Thevalues of ∆ l/M BH are approximately the same as the spe-cific angular momentum at an effective ISCO. This resultimplies that the separation at which the orbital motion ofthe binary becomes unstable and the NS falls into the BHis given effectively by the ISCO in the equatorial motionaround the BH with χ eff . Taking a closer look, ∆ l/M BH tends to be smaller than the value for the effective ISCO.This is likely to stem from the gravitational-wave emis-sion that dissipates the orbital angular momentum whilethe matter falls into the BH, or the redistribution of thespecific angular momentum due to the tidal torque. E. The gravitational waveform
The misalignment between the orbital angular momen-tum and the BH spin causes the precession of the orbitand induces the modulation in gravitational waves. Also,the misalignment angle of the BH spin and the EOS ofthe NS affect the tidal-disruption process, and as a re-sult, gravitational waveforms are modified by them. In
SACRA , we extract the out-going component of the com-plex Weyl scalar Ψ at finite radii and project it ontothe spin-weighted spherical harmonic functions. Here, wetook the axis of the spherical harmonics to be z -axis: theinitial direction of the total angular momentum. Thento obtain gravitational waveforms, we integrate Ψ twicein time as h ( t ) = h + ( t ) − ih × ( t ) = (cid:90) t dt (cid:48) (cid:90) t (cid:48) dt (cid:48)(cid:48) Ψ ( t (cid:48)(cid:48) ) . (63)In the following, we plot the normalized amplitude Dh/m or the amplitude observed at a hypothetical dis-8 FIG. 14. The volume-rendered density map of the ejecta at the time that the ejecta expands to ≈ ≈
10 ms after the onset of merger. The irreducible mass of the BH( M irr , f ), the mass of the BH ( M BH , f ), the dimensionless spin parameter ( χ f ), the BH spin ( S BH , f ), the tilt angle of the BH spin(cos − ( S z BH /S BH )), and the dominant QNM frequency derived by the fitting formula Eq. (67) using the result for M BH and χ ( f QNM ), respectively. The results for the aligned-spin case are taken from [17]. For the model APR4i90, we failed to find thelocation of the AH because of inappropriate setting of the AMR domains.Model M irr , f [ M (cid:12) ] M BH , f [ M (cid:12) ] χ f S BH , f [ GM (cid:12) /c ] cos − ( S z BH /S BH )[ ◦ ] f QNM [kHz]APR4i0 6.83 7.82 0.85 52 0 ◦ < ◦ ◦ ◦ < ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ ◦ tance D = 100 Mpc as a function of approximate re-tarded time defined by t ret = t − D − M ln DM . (64)The Fourier spectrum of the gravitational waveformcould reflect more quantitative information. In this pa-per, we define the Fourier power spectrum of gravita-tional waves as the root mean square of two independentpolarizations as˜ h ( f ) = (cid:115) | ˜ h + ( f ) | + | ˜ h × ( f ) | , (65)˜ h A ( f ) = (cid:90) h A ( t ) e πift dt, (66)( A = + , × ) . We will plot a dimensionless Fourier spectrum ˜ h eff ( f ) := f ˜ h ( f ) observed at a hypothetical distance D = 100 Mpc as a function of the frequency f , or a normalized spec-trum D ˜ h eff ( f ) /m as a function of dimensionless fre-quency f m .Figure 16 shows plus-mode gravitational waveformsof ( l, m ) = (2 , , ,
0) for models APR4i0and APR4i60. While the amplitude of ( l, m ) = (2 , ,
0) is smaller than that of ( l, m ) = (2 ,
2) for thealigned-spin case, they could have a significant contribu-tion to gravitational waveforms for the misaligned-spincase. This is because the direction of the orbital an-gular momentum does not always agree with the axisof the spin-weighted spherical harmonic function for themisaligned-spin case.As we show in Sec. IV A, the angular velocity of theorbital precession is always smaller than the orbital an-gular velocity by an order of magnitude. Thus, gravita-tional waves for the misaligned-spin case have the featuresimilar to gravitational waves from the aligned-spin caseobserved from an inclined direction with respective to L no r m a li z ed s pe c i f i c angu l a r m o m en t u m χ eff APR4ALF2H4MS1l
ISCO /M BH FIG. 15. The average specific angular momentum of materialbrought into BH normalized by M BH , ∆ l/M BH , as a functionof an effective dimensionless spin parameter χ eff = χ cos i tilt , .The specific angular momentum at the ISCO of the BH is alsoplotted. for each instant. Indeed, it has already shown for pre-cessing binary BH cases (see, e.g., [84]) that the wave-forms take a far simpler form in the quadrupole align-ment (QA) frame: the frame in which z -axis agrees withthe instantaneous direction of L . If we project gravita-tional waves onto the spherical-harmonic function in theQA frame, and describe these expansion coefficients by( l (cid:48) , m (cid:48) ), ( l (cid:48) , m (cid:48) ) = (2 , ±
2) modes are the dominant modes.Under the rotational transformation, these componentsmixes not only into ( l, m ) = (2 , ±
2) modes but also intodifferent m modes with l = 2. This is consistent withthe fact that the dominant frequency of ( l, m ) = (2 , ,
0) modes agrees with the frequency of (2 ,
2) moderather than the half (this fact is also pointed out in [36]).We note that different l modes do not mix under therotation of the axis of spherical harmonics.Other m (cid:48) modes, such as ( l (cid:48) , m (cid:48) ) = (2 , ± , ± , ±
3) modes, also contribute to the gravitationalwaveform. Because the phase of the m (cid:48) mode is m (cid:48) Φ,where Φ is the phase of the orbit, the mixing amongdifferent m (cid:48) modes causes modulation in the amplitudeof the waveforms. For example, when the ( l (cid:48) , m (cid:48) ) = (2 , , ±
2) modes, the amplitudeexhibits modulation with the periods of 2 π and 6 π interms of the orbital phase, Φ. Indeed, Fig. 17 showsthat the amplitude observed along the z -axis modulatesprimarily with the period of ≈ π in terms of Φ.Obviously, the mixing of several m (cid:48) modes in gravi-tational waves can occur for the aligned-spin case if wechoose the axis of spherical harmonics which disagreeswith the orbital angular momentum. One thing to benoted is that because the orbital angular momentum pre-cesses for the misaligned-spin case, we cannot avoid thesituation that the orbital angular momentum disagreeswith the axis of spherical harmonics. This implies thatthe mixing among several m (cid:48) modes is unavoidable forthe misaligned-spin case. Figure 18 shows plus-mode gravitational waveforms(the left panels) and gravitational-wave spectra (the rightpanels) observed from different inclinations, θ = 0 ◦ , 45 ◦ ,and 90 ◦ with respect to the z -axis. As we show inSec. IV A, J is always aligned approximately with the z -axis, and hence θ is regarded approximately as the an-gle between J and the direction of the observer. We takeinto account the contributions from l = 2 − θ = 0 ◦ has asimilar feature to those for the aligned-spin case; the am-plitude and the frequency approximately monotonicallyincrease as the system evolves. Only slight modulationin the amplitude and the frequency is found in the wave-forms. There are much more appreciable modulationin the waveforms for θ = 45 ◦ and 90 ◦ . The amplitudeof the waveforms becomes smaller for larger inclinationangles. This dependence stems from the fact that for( l (cid:48) , m (cid:48) ) = (2 , ±
2) modes the amplitude of the plus-modewaveform is proportional to (1 + cos θ ) /
2. The phasesof the waveforms agree among different values of θ .For MS1i90, the modulation in the amplitude and thefrequency are larger than that for MS1i30. There is nomonotonic dependence of the amplitude of the waveformson θ . While the amplitude of the waveform for θ = 0 ◦ increases as the time evolves, that of θ = 90 ◦ approxi-mately decreases until t ret ≈
20 ms. The phase evolutionof the waveforms is also different among different valuesof θ . This is in particular significant for the last ≈
10 msof the inspiral phase. These features of the waveformsare due to the precession of the orbital plane. Althoughthe orbital plane precesses, the angle between the z -axisand L is always approximately constant until the onsetof merger for simulations in this study. Therefore, theangle between the direction of the observer and L is ap-proximately unchanged for θ = 0 ◦ , and the effect of theorbital precession to the waveform is small. On the otherhand, for the observer at θ = 90 ◦ , the orbital precessionchanges the angle between L and the direction of the ob-server, and affect the waveforms. In particular, the am-plitude of the waveform is strongly suppressed when theorbital plane is edge-on to the direction of the observer.In the top-right and the bottom-right panels in Fig. 18,we find that the amplitude of the spectra for MS1i30decreases approximately monotonically for larger valuesof θ , while there is no simple dependence of the amplitudeon θ for MS1i90. While the spectrum for the aligned-spinmodel has a flat shape (see Fig. 19 or [6, 7] for the spectrafor the aligned-spin case), there are some bumps in thespectra for the misaligned-spin models. In particular,bumps at ≈ θ = 45 ◦ and 90 ◦ change the feature of the cutoff. These bumps are theconsequences of the mixing among the different modesof gravitational waves. We find that l = 4 modes arenegligible for the bumps in the spectra. Since m (cid:48) = 2and 0 do not contribute to the modulation, ( l (cid:48) , m (cid:48) ) =(2 , ± , ± , ±
3) are primarily responsible forthe bumps.0 -5e-22-4e-22-3e-22-2e-22-1e-22 0 1e-22 2e-22 3e-22 4e-22 5e-22 5 10 15 20 25 30 -0.15-0.1-0.05 0 0.05 0.1 0.15 h [ M p c ] D h / m t ret [ms]APR4i0(l,m)=(2,2)(l,m)=(2,1)(l,m)=(2,0) -5e-22-4e-22-3e-22-2e-22-1e-22 0 1e-22 2e-22 3e-22 4e-22 5e-22 5 10 15 20 25 30 35 -0.15-0.1-0.05 0 0.05 0.1 0.15 h [ M p c ] D h / m t ret [ms]APR4i60(l,m)=(2,2)(l,m)=(2,1)(l,m)=(2,0) FIG. 16. The plus-mode gravitational waveforms for ( l, m ) = (2 , , (2 , ,
0) for APR4i0 (the left panel) and APR4i60(the right panel). ( h + + h x ) / D / m Φ /2 π APR4i30APR4i60APR4i90
FIG. 17. The square of the gravitational-wave amplitude ob-served along the z -axis as a function of the orbital phase de-vided by 2 π , Φ / π , for models with the MS1 EOS. The spectra of gravitational waves reflects the fate ofthe merger. In particular, as already clarified by [6, 7, 32],the location of the cutoff in the spectra has a correlationwith the compactness of the NS. We expect that this alsoholds for the misaligned-spin case. However, the bumpsof the spectra for the misaligned-spin models may make itdifficult to determine the cutoff frequency, particular forthe case that gravitational waves are observed from largeinclination angles. Thus, in the following, we check thebehavior of the bump in the spectra and the dependenceof the cutoff in the spectra on the models qualitatively asa first step to understand the gravitational-wave spectrafrom the misaligned-spin BHNS. Here, we only considerthe waveforms for an observer along the z -axis, the direc-tion of the initial total angular momentum, which wouldbe the simplest case.We show the spectra of the waveforms for different val- ues of i tilt , with APR4 in the top-left panel in Fig. 19.The bumps in the spectra become deeper as the valueof i tilt , increases. This is because the amplitude of( l (cid:48) , m (cid:48) ) = (2 , ± , ± , ±
3) modes depend onthe angle between L and the direction of the observer,and θ (cid:48) becomes large for large values of i tilt , . Forthe models with APR4, the quasi-normal mode (QNM)is excited, and the cutoff frequency of the spectra re-flects the frequency of the QNM. The frequency of the( l, m ) = (2 ,
2) least-damped mode of the QNM is givenapproximately by the fitting fomula [85] as f QNM = 12 πM BH ,f (cid:2) . − . − χ f ) . (cid:3) . (67)The values of the QNM frequency are summarized inTable VI. Figure 19 shows that the spectra have a cutoffaround the QNM frequency, although the bumps of thespectra for the misaligned-spin models make it difficultto determine the cutoff frequency accurately.Next, we compare the spectra of the waveforms withdifferent EOSs. The top right, bottom left, and bottomright panels of Fig. 19 show the spectra for i tilt , ≈ ◦ ,60 ◦ , and 90 ◦ , respectively with four different EOSs. Weagain find the bumps in the spectra, and they have ap-proximately the same shape in f ≈ f > i tilt , ≈ ◦ than with i tilt , ≈ ◦ , we still find that thecutoff frequency should be, at least in principle, differ-ent among the models with i tilt , ≈ ◦ . For example,1 -5e-22-4e-22-3e-22-2e-22-1e-22 0 1e-22 2e-22 3e-22 4e-22 5e-22 5 10 15 20 25 30 35 40 -0.15-0.1-0.05 0 0.05 0.1 0.15 h [ M p c ] D h / m t ret [ms]MS1i30 θ =0[deg.] θ =45[deg.] θ =90[deg.] h e ff [ M p c ] D h e ff / m f[Hz]MS1i30fm θ =0[deg.] θ =45[deg.] θ =90[deg.] -5e-22-4e-22-3e-22-2e-22-1e-22 0 1e-22 2e-22 3e-22 4e-22 5e-22 5 10 15 20 25 30 -0.15-0.1-0.05 0 0.05 0.1 0.15 h [ M p c ] D h / m t ret [ms]MS1i90 θ =0[deg.] θ =45[deg.] θ =90[deg.] h e ff [ M p c ] D h e ff / m f[Hz]MS1i90fm θ =0[deg.] θ =45[deg.] θ =90[deg.] FIG. 18. Plus-mode gravitational waveforms (left panels) and gravitational-wave spectra (right panels) observed from differentinclination angles with respect to the z -axis. The top and bottom figures are the results for models MS1i30 and MS1i90,respectively. Gravitational waveforms are plotted as functions of retarded time, t ret . The left axes in the plots of the waveformsdenote the amplitude observed at a hypothetical distance D = 100 Mpc, and the right does the normalized amplitude Dh/m .The upper axis in the plots of the spectra denotes the dimensionless frequency, fm , and the right axis denotes the normalizedamplitude D ˜ h eff ( f ) /m . The bottom axis in the plots of spectra denotes the frequency of the gravitational waveform, in Hz,and the left axis denotes the amplitude observed at a hypothetical distance D = 100 Mpc. if we say the cutoff frequency as the frequency at which h eff ≈ × − is achieved, the difference of the cutofffrequency between models with APR4 and H4 is ≈ i tilt , ≈ ◦ and i tilt , ≈ ◦ cases, and the dif-ference among the EOSs is always larger than 10% for i tilt , ≈ ◦ and 6% for i tilt , ≈ ◦ . By contrast to thecases with i tilt , ≈
30 and 60 ◦ , the difference in the spec-tra is small among different EOSs for i tilt , ≈ ◦ , andhence it might be difficult to distinguish the EOS fromthe spectra. The reason for this is that tidal disruptionoccurs so weakly that the difference in the EOS is notappreciable. The cutoff frequency at f ≈ F. Emitted energy, linear momentum, angularmomentum by gravitational waves
A binary loses its orbital energy and orbital angularmomentum by the gravitational radiation. The amountof energy and angular momentum emitted by gravita-tional waves depends on the binary parameter, and thatis one of the interests in studying the binary merger.Gravitational waves could also carry linear momentum ofthe system. Non-zero linear momentum radiation causesrecoil of the system. We evaluate these energy, linearmomentum, and angular momentum emitted by gravita-tional waves, using the formula of [86].In Table VII, we list the total energy, ∆ E (and its ra-tio to the initial ADM mass), the total linear momentumnormalized by the initial ADM mass, P GW /M , and the2 h e ff [ M p c ] D h e ff / m f[Hz]fm APR4i0APR4i30APR4i60APR4i90 h e ff [ M p c ] D h e ff / m f[Hz]fm APR4i30ALF2i30H4i30MS1i30 4e-23 5e-23 6e-23 8e-23 1e-22 2e-22 3e-22 4e-22 5e-22 300 500 700 1000 2000 3000 0.01 0.02 0.03 0.04 0.06 0.08 0.1 0.2 0.01 0.02 0.03 0.04 0.06 0.08 0.1 0.2 h e ff [ M p c ] D h e ff / m f[Hz]fm APR4i60ALF2i60H4i60MS1i60 h e ff [ M p c ] D h e ff / m f[Hz]fm APR4i90ALF2i90H4i90MS1i90
FIG. 19. Gravitational wave spectra observed along the z -axis of the simulation. The top-left plot shows the spectra for themodels with APR4 and with different initial values of i tilt , . The top-right, bottom-left, and bottom-right plots compare thespectra among four different EOSs for i tilt , ≈ ◦ , 60 ◦ , and 90 ◦ , respectively. The vertical lines at the upper axis show theQNM frequency of the remnant BH for each models (see Table VI).TABLE VII. The list of total radiated energy, ∆ E , its ratio to initial ADM mass ∆ E/M , total radiated angular momentum∆ J , total radiated linear momentum normalized by the initial ADM mass, P GW /M , and the linear momentum of the ejectanormalized by the initial ADM mass, P eje /M , respectively. For the i90 models, the data for P eje ,i was not output.Model ∆ E [ M (cid:12) c ] (∆ E/M [%]) ∆ J [ GM (cid:12) /c ] P GW /M [km / s] P eje /M [km / s]APR4i30 0.15(1.9) 11 8 . × . . × . × − APR4i90 0.098(1.2) 6.2 5 . × − ALF2i30 0.12(1.4) 9.4 7 . × . × ALF2i60 0.11(1.4) 7.8 2 . × . × ALF2i90 0.092(1.1) 6.1 3 . × − H4i30 0.093(1.2) 8.5 7 . × . × H4i60 0.093(1.2) 7.2 6 . × . × H4i90 0.085(1.1) 5.8 3 . × − MS1i30 0.074(0.93) 7.6 6 . × . × MS1i60 0.075(0.93) 6.6 9 . × . × MS1i90 0.073(0.90) 5.5 1 . × − J , emitted by gravitational waves.We also list the recoil velocity caused by the mass ejec-tion, P eje /M . We take into account the contributionsfrom l = 2 − l = 2 modes contribute to ∆ E and ∆ J by more than 83%, while l = 3 by ≈ l = 4 by ≈ i tilt , .Table VII shows that ∆ E and ∆ J decrease monoton-ically with the decrease of the compactness of the NS.The same dependence of ∆ E and ∆ J on the compact-ness of the NS is found for the results obtained by thenon-spinning and aligned-spin BH-NS mergers [6, 7], andit is due to the fact that a longer inspiral phase (i.e.,longer gravitational-wave emission phase) is realized fora softer EOS. Table VII also shows that, for a fixed EOS,∆ E and ∆ J monotonically decrease with the increase of i tilt , . The dependence of ∆ E on i tilt , is weaker thanthat on the EOS as far as tidal disruption is appreciable,while ∆ J depends appreciably on i tilt , .The recoil velocity induced by the gravitational-waveemission, P GW /M , decreases as the compactness of theNS becomes small. This is because the smaller compact-ness of the NS results in earlier tidal disruption duringthe insprial phase, resulting in an earlier shutdown of thegravitational-wave emission. The recoil velocity causedby the mass ejection is larger for the models with smallercompactness of the NS because the ejected mass becomeslarger for these models. This opposite dependence of therecoil velocity on the compactness for P GW and P eje re-verses the dominant component for the recoil [17]. Whilethe recoil due to the gravitational-wave emission is dom-inant for models with a large compactness, the recoil in-duced by the mass ejection becomes dominant for mod-els with a small compactness. These two components arecomparable for the case that M eje ∼ . M (cid:12) . V. SUMMARY AND DISCUSSION
We performed numerical-relativity simulations for themerger of BH-NS binaries with various BH spin misalign-ment angles, employing four models of nuclear-theory-based EOSs described by a piecewise polytrope. We in-vestigated the dependence of the orbital evolution in thelate inspiral phase, tidal-disruption process of the NS,properties and structures of the remnant disk and ejecta,properties of the remnant BH, gravitational waveformsand their spectra on the BH spin misalignment angle andthe EOS of the NS.We showed that a large BH spin misalignment anglesuppresses the NS tidal disruption event by the reductionof the spin-orbit interaction. The remnant mass of thematerial outside the BH decreases as the misalignmentangle of the BH increases. This dependence agrees withthe previous results for misaligned-spin BH-NS merg-ers [36, 41]. Also we reconfirm the findings in [6, 7] thatthe remnant mass increases as the compactness of the NS decreases. In our study, this dependence on the com-pactness of the NS is shown irrespective of the BH mis-alignment angle. The deviation of the result from theprediction of the fitting formula [80] is within 50% for M > AH (cid:38) . M (cid:12) and within 30% for M > AH (cid:38) . M (cid:12) ,even though we employ a simple definition of the effec-tive spin parameter. This reconfirms the argument of [41]that the mass of the material outside the remnant BHcan be modeled with a good accuracy by considering theresult of the aligned-spin cases.Effects of the orbital precession are reflected in thetidal tail. The elevation angle of the tidal tail measuredfrom the xy -plane is different for each part. Although itis pointed out in [41] that the elevation angle may preventits elements to collide, the material of the tidal tail stillcollides with each other, and a weakly inclined torus iseventually formed at least for the model with i tilt , (cid:46) ◦ .Monotonic dependence of the remnant disk mass andthe ejecta mass on the orbital misalignment angle wasshown. Both M disk and M eje decrease as the misalign-ment angle increases. However, we still found that ifthe compactness of the NS is moderate C = 0 .
160 oreven small C = 0 . i tilt , (cid:46) ◦ and i tilt , (cid:46) ◦ can produce disks larger than 0 . M (cid:12) .Such a system could be a candidate for the progenitorsof sGRB. M eje > . M (cid:12) is achieved for i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . i tilt , < ◦ with C = 0 . i tilt , ≈ ◦ , the structure of thedisk is similar to the disk formed for the aligned-spinBH-NS mergers. In particular, the rotational axis of thedense part ( (cid:38) g / cm ) of the disk is aligned approx-imately with the remnant BH spin. On the other hand,for the models with i tilt , ≈ ◦ , we found that the axisof the disk is misaligned with the direction of the rem-nant BH spin initially with ≈ ◦ , although the mis-alignment angle of the dense part of the disk approacheszero in ≈ ≈ ◦ and ≈ ◦ for models with i tilt , ≈ ◦ and ≈ ◦ ,respectively. This reflects the orbital elevation duringthe inspiral phase. It is pointed out in [82] that themisalignment of the disk may affect the light curve ofthe sGRB. However, since the high-density part of thedisk with ρ > g / cm would play main roles, theeffect of the BH spin misalignment may not be observ-able in the sGRBs, because the dense part of the diskbecomes aligned with the BH spin in a relatively shorttime scale.We suspect that Bardeen-Petterson-like effectinduced by a purely hydrodynamical mechanism, such asangular momentum redistribution due to a shock waveexcited in a non-axisymmetric manner of the disk, shouldwork in the disk for the alignment.We found that the accretion time scale of the matterin the disk to the BH is typically ≈
100 ms, and depends4weakly on the binary parameters. The main mechanismof the accretion in the present context is the redistribu-tion of the angular momentum due to the torque exertedby non-axisymmetric structure of the disk, which is seenin Fig. 8. In reality, the viscosity induced by the magne-torotational instability turbulence could play an impor-tant role for this phase [88]. Since we did not take thoseeffects into account, the accretion rate for the late phasemight not be very quantitative. However, the present re-sult shows that the purely hydrodynamical effect is im-portant for the accretion of the matter, and this effectshould be considered whenever we study the evolution ofthe accretion disk formed by BH-NS.We found that the velocity of the ejecta is typically0 . . c , and has only weak dependence on the misalign-ment angle of the BH spin and the EOS of the NS. Wealso found that the morphology of the ejecta changes de-pending on the ejecta mass: Crescent-like-shaped ejectawith its opening angle ≈ ◦ is formed for relativelymassive ejecta ( (cid:38) . M (cid:12) ), while spiral-shaped ejectawith its opening angle larger than 360 ◦ is formed for rel-atively less massive ejecta ( (cid:46) . M (cid:12) ). In particular,the spiral shape reflects the orbital precession. This de-pendence of the ejecta morphology was also found forthe aligned-spin case [17], and might be explained by theperiastron advance in general relativity.We found that the dimensionless spin parameter of theremnant BH depends only weakly on the EOS of the NS,but it depends strongly on the misalignment angle, i tilt , .The final direction of the BH spin becomes aligned ap-proximately with the initial direction of the total angularmomentum. We also found an approximate relation be-tween the misalignment angle and the increase of the BHspin, and that the specific angular momentum that theBH gained during the merger approximately agrees witha specific angular momentum at the ISCO of the BH with χ eff = χ cos i tilt , .We showed that the mixing among the components ofspherical harmonics occurs and causes the modulation ingravitational waveforms for the misaligned-spin case. Inparticular, we found that the period in the modulationis primarily ≈ π in terms of the orbital phase. We alsostudied the dependence of waveforms on the direction ofthe observer, and found that the modulation due to theorbital precession becomes significant for the case thatthe observer is located along the direction perpendicularto the total angular momentum. The bump-shape modu-lation in the power spectrum of gravitational waveformsis found, and the depth of the bump becomes large as i tilt , becomes large.In the presence of the bumps in the spectra, the lo-cation of the cutoff frequency becomes obscured. Nev-ertheless, for the case that gravitational waves are ob-served along the axis of total angular momentum, thedifferences of the location of the cutoff in the spectraamong the EOSs are seen for i tilt , (cid:46) ◦ , while they arehardly found for i tilt , ≈ ◦ . This result shows that, inprinciple, gravitational waves from BH-NS binaries with Q = 5 and χ = 0 .
75 contain the information of the EOSof the NS even if the misalignment angle of the BH spinis large up to ≈ ◦ . To discuss whether we can extractthe information of the NS EOS from the waveform bythe observation, we need to define a quantitative indica-tor which reflects the information of the EOS, such as acutoff frequency in the spectra, in an appropriate mannereven in the presence of the orbital precession, and discussthe detectability considering the noise in the signal. Weleave these tasks for our future study.The dependence of the energy, linear momentum, andangular momentum radiated by gravitational waves onthe misalignment angle and EOS was shown. We foundthat the recoil induced by the mass ejection dominatesthe total recoil velocity for the case that the ejecta massis larger than ≈ . M (cid:12) , while the recoil induced by thegravitational radiation is dominant for the case that theejecta mass is smaller.Finally, we list several issues to be explored in thefuture. In this paper we studied the models only with Q = 5 and χ = 0 .
75 to focus on the dependence on theBH spin misalignment and the EOS of the NS. As it isknown that the mass ratio and the BH spin magnitudeinfluences on the merger process, we also need to clarifythe dependence on these parameters with the spin mis-alignment systematically. In particular, the larger BHspin enhances the tidal disruption of the NS, and thuscharacteristic features of misaligned-spin BH-NS mergerscould be revealed more clearly. Also, we plan to performmore detailed analysis of gravitational waveforms for themisaligned-spin cases, because the waveforms may con-tain rich information on the misaligned BH-NS system.
ACKNOWLEDGMENTS
Kyohei Kawaguchi is grateful to Hiroki Nagakura,Kenta Hotokezaka, Kunihito Ioka, Kenta Kiuchi, Sho Fu-jibayashi, Takashi Yoshida, and Yuichiro Sekiguchi forvaluable discussions. Koutarou Kyutoku is grateful toFrancois Foucart for valuable discussions. This workwas supported by JSPS Grant-in Aid for Scientific Re-search (24244028) and RIKEN iTHES project. KyoheiKawaguchi is supported by JSPS Research Fellowship forYoung Scientists (DC1). Hiroyuki Nakano acknowledgessupport by JSPS Grant-in-Aid for Scientific Research(24103006). Keisuke Taniguchi acknowledges support byJSPS Grant-in-Aid for Scientific Research (26400267).
Appendix A: Convergence with respect to the gridresolution
Table VIII compares M > AH and M eje among differentgrid resolutions for selected models. If we assume thefirst-order convergence between N = 48 and N = 60,the errors with N = 60 results are always smaller than ≈
40% and ≈
32% for M > AH and M eje , respectively.5 Model
N M > AH [ M (cid:12) ] M eje [ M (cid:12) ]ALF2i30 40 0.172 3 . × −
48 0.166 3 . × −
60 0.156 3 . × − H4i30 40 0.264 4 . × −
48 0.257 4 . × −
60 0.248 4 . × − H4i60 40 0 .
103 1 . × −
48 9 . × − . × −
60 8 . × − . × − MS1i90 40 2 . × − . × −
48 2 . × − . × −
60 2 . × − . × − TABLE VIII. M > AH and M eje for runs with different gridresolutions for selected models. Errors become large for a smaller mass. In particular,the error for M > AH is ≈
16% for model H4i30, while theerror is ≈
28% for model MS1i90. , 044049 (2010).[7] K. Kyutoku, H. Okawa, M. Shibata, and K. Taniguchi,Phys. Rev. D , 064018 (2011).[8] L. Lindblom, Astrophys. J. , 569 (1992).[9] M. Vallisneri, Phys. Rev. Lett. , 3519 (2000).[10] J. S. Read, C. Markakis, M. Shibata et al. , Phys. Rev.D , 124033 (2009).[11] V. Ferrari, L. Gualtieri, and F. Pannarale, Phys. Rev. D , 064026 (2010).[12] B. S. Sathyaprakash, B. F. Schutz, Living Rev. Relativity , 2 (2011).[13] B. Paczynski, Acta Astronomica , 257 (1991).[14] E. Nakar, Phys. Rep. , 166 (2007).[15] E. Berger, Annu. Rev. Astron. Astrophys. , 43 (2014).[16] K. Kyutoku , K. Ioka, M. Shibata, Phys. Rev. D ,041503 (2013).[17] K. Kyutoku, K. Ioka, H. Okawa, M. Shibata, and K.Taniguchi, arXiv:1502.05402 (2015).[18] J. M. Lattimer, D. N. Schramm , Astrophys. J. Lett., , L145 (1974).[19] L. Li, B.Paczynski, Astrophys. J. , L59 (1998).[20] S. R. Kulkarni, arXiv:astro-ph/0510256 (2005).[21] B. D. Metzger, G. Martinez-Pinedo, S. Darbha et al.,Mon. Not. R. Astron. Soc. , 2650 (2010).[22] B. S. Meyer, Astrophys. J. , 254 (1989).[23] C. Freiburghaus, S. Rosswog, F.-K. Thielemann, Astro-phys. J. Lett. , L121 (1999).[24] M. Shibata and K. Taniguchi, Living Rev. Relativity ,6 (2011).[25] K. Belczynski, R. E. Taam, E. Rantsiou, M. van denSluys, Astrophys. J. , 474 (2008).[26] L. E. Kidder, Phys. Rev. D, , 821 (1995).[27] M. Shibata and K. Ury¯u, Phys. Rev. D , 121503(R)(2006). [28] M. Shibata and K. Ury¯u, Class. Quant. Grav. , S125(2007).[29] M. Shibata and K. Taniguchi, Phys. Rev. D , 084015(2008).[30] Z. B. Etienne, J. A. Faber, Y. T. Liu et al., Phys. Rev.D , 084002 (2008).[31] M. D. Duez, F. Foucart, L. E. Kidder et al., Phys. Rev.D , 104015 (2008).[32] M. Shibata, K. Kyutoku, T. Yamamoto, and K.Taniguchi, Phys. Rev. D , 044030 (2009).[33] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baum-garte, Phys. Rev. D , 044024 (2009).[34] S. Chawla, M. Anderson, M. Besselman et al., Phys. Rev.Lett. , 111101 (2010).[35] M. D. Duez, F. Foucart, L. E. Kidder et al., Class. Quant.Grav. , 114106 (2010).[36] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukol-sky, Phys. Rev. D , 024005 (2011).[37] Z. B. Etienne, Y. T. Liu, V. Paschalidis, and S. L.Shapiro, Phys. Rev. D , 064029 (2012).[38] Z. B. Etienne, V. Paschalidis, and S. L. Shapiro, Phys.Rev. D , 084026 (2012).[39] G. Lovelace, M. D. Duez, F. Foucart et al., Class. Quant.Grav. , 135004 (2013).[40] M. B. Deaton, M. D. Duez, F. Foucart et al., Astrophys.J. ,47(2013).[41] F. Foucart, M. B. Deaton, M. D. Duez et al., Phys. Rev.D , 084006 (2013).[42] F. Foucart et al., Phys. Rev. D , 024026 (2014).[43] E. Rantsiou, S. Kobayashi, P. Laguna, and F. A. Rasio,Astrophys. J. ,1326(2008).[44] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman,Phys. Rev. D , 124032 (2009).[45] T. Regge and C. Teitelboim, Phys. Lett. B , 124018 (2009).[48] P. C. Peters, Phys. Rev. , B1224 (1964).[49] J. W. York, Phys. Rev. Lett. , 1350 (1999).[50] H. P. Pfeiffer and J. W. York, Phys. Rev. D , 044022(2003). [51] G. B. Cook, Living Rev. Relativity , 5 (2000).[52] J. M. Bowen and J. W. York, Phys. Rev. D , 2047(1980).[53] K. Taniguchi, T. W. Baumgarte, J. A. Faber, and S. L.Shapiro, Phys. Rev. D 75, 084005 (2007).[54] K. Taniguchi, T. W. Baumgarte, J. A. Faber, and S. L.Shapiro, Phys. Rev. D 77, 044003 (2008).[55] P. Grandclement, Phys. Rev. D , 124002 (2006).[56] P. Grandclement, Phys. Rev. D , 129903(E) (2007).[57] F. Foucart, L. E. Kidder, H. P. Pfeiffer, S. A. Teukolsky,Phys. Rev. D , 124051 (2008).[58] R. Beig, Phys. Lett. A , 153 (1978).[59] A. Ashtekar and A. Magnon-Ashtekar, J. Math. Phys. , 793 (1979).[60] C. S. Kochanek, Astrophys. J. , 234 (1992).[61] L. Bildsten and C. Cutler, Astrophys. J. , 175 (1992).[62] D. R. Lorimer, Living Rev. Relativity , 8 (2008).[63] S. A. Teukolsky. Astrophys. J. , 442 (1998).[64] M. Shibata. Phys. Rev. D. , 024012 (1998).[65] A. Boh´e, S. Marsat, G. Faye, and L. Blancet, ClassicalQuantum Gravity , 075017 (2013).[66] J. M. Lattimer, M. Prakash, Science, 1090720 (2004).[67] P. B. Demorest, T. Pennucci, S. M. Ransom et al. , Na-ture , 1081 (2010).[68] J. Antoniadis , P. C. C. Freire, N. Wex et al. Science ,448 (2013).[69] K. Hebeler, J. M. Lattimer, C. J. Pethick and A.Schwenk, Astrophys. J. , 11 (2013).[70] A. W. Steiner and J. M. Lattimer, Astrophys. J. ,123 (2014).[71] L. Bildsten and C. Cutler, Astrophys. J. , 175 (1992).[72] T. Yamamoto, M. Shibata, K. Taniguchi, Phys. Rev. D, , 064054 (2008). [73] D. Hilditch, S. Bernuzzi, M. Thierfeder et al., Phys. Rev.D, , 084057 (2013).[74] K. Kyutoku, M. Shibata, and K. Taniguchi, Phys. Rev.D , 064006 (2014).[75] T. Damour, Phys. Rev. D , 124013 (2001).[76] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys.Rev. D , 041501(R) (2006).[77] B. Brugmann, J. A. Gonzalez, M. Hannam et al., Phys.Rev. D , 024027 (2008).[78] G. Lovelace, R. Owen, H. P. Pfeiffer, T. Chu, Phys. Rev.D , 084017 (2008).[79] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnanand D. Merritt, Phys.Rev.D , 064030 (2007).[80] F. Foucart, Phys. Rev. D , 124007 (2012).[81] P. C. Fragile, O. M. Blaes, P. Anninos, and J. D.Salmonson, Astrophys. J. , 417 (2007).[82] N. Stone, A. Loeb, E. Berger, Phys. Rev. D , 084053(2013).[83] K. Hotokezaka, K. Kiuchi, K. Kyutoku et al., Phys. Rev.D , 024001 (2013).[84] P. Schmidt, M. Hannam, and S. Husa, arXiv:1207.3088(2012).[85] E. Berti, V. Cardoso, A. O. Starinets, Class. Quant.Grav. , 163001 (2009).[86] M. Ruiz, M. Alcubierre, D. Nez, and R. Takahashi, Gen.Relativ. Gravit. , 1705 (2008).[87] J. M. Bardeen, J. A. Petterson, Astrophys. J. , L65(1975).[88] S. A. Balbus, J. F. Hawley, Astrophys. J.376