3.1 Numerical method

3.1.1 Initial condition


Table 2: Quasi-equilibrium study groups by which the initial data of each simulation group is supplied.
Simulation group Quasi-equilibrium study


AEI Grandclément [82Jump To The Next Citation Point, 83Jump To The Next Citation Point]
CCCW Foucart et al. [75Jump To The Next Citation Point]
KT Shibata & UryΕ« [202Jump To The Next Citation Point, 203Jump To The Next Citation Point] and Kyutoku et al. [106]
LBPLI Grandclément [82, 83]
UIUC Taniguchi et al. [209Jump To The Next Citation Point, 210Jump To The Next Citation Point]

As shown in Table 1, the studies of quasi-equilibrium and its sequences for BH-NS binaries have been performed by several groups. In the numerical simulations, quasi-equilibrium is employed as the initial condition (besides a minor modification of the quasi-equilibrium state, reviewed below). The five simulation groups, AEI, CCCW, KT, LBPLI, and UIUC, employ the quasi-equilibrium states derived by different groups as summarized in Table 2. Depending on the implementations on which each group relies, the type of the initial condition is different. Table 3 classifies the type of the initial data by the following six parameters:

(1)
The EOS of the NS

The notation of “Γ = 2” or “ Γ = 2.75” means the single polytropic EOS with an adiabatic index of Γ = 2 or Γ = 2.75.

The notations of “PwPoly” and “SH” mean the piecewise polytropic EOS and the Shen’s EOS, respectively.

(2)
The spatial background metric &tidle;γij and the extrinsic curvature K

The notation of “CF” means the conformally-flat condition &tidle;γ = f ij ij, and that of “M” means the maximal slicing K = 0.

The notation of “KS” indicates the Kerr–Schild metric. In this case, the extrinsic curvature K is also set to that derived from the Kerr–Schild metric.

(3)
The state of the fluid flow in the NS

The notation of “Co” means a corotating NS, and that of “Ir” an irrotational NS.

(4)
The method to treat the BH, i.e., the moving puncture or the excision

The notations of “Pu” and “Ex” mean that the moving-puncture and the excision approaches are used in the numerical simulation, respectively. “ExPu” means that the initial condition is prepared in the excision approach, but the simulation is done in the moving-puncture formulation (see below).

(5)
Presence or absence of the BH spin a and misalignment angle πœƒ between the spin axis and orbital angular momentum vector.
(6)
The mass ratio Q


Table 3: Summary of the initial data used in simulations performed so far. Note that LBPLI’s data includes magnetic fields.
Group & Ref. (1) (2) (3) (4) (5) (6)
KT [202] Γ = 2 CF,M Co Pu a = 0 Q ≈ 2.5
KT [203] Γ = 2 CF,M Co Pu a = 0 Q ≈ 2.5,3.1
KT [197] Γ = 2 CF,M Ir Pu a = 0 Q ≈ 2.5,3.1
KT [228Jump To The Next Citation Point] Γ = 2 CF,M Ir Pu a = 0 Q ≈ 3.1
KT [194Jump To The Next Citation Point] Γ = 2 CF,M Ir Pu a = 0 Q = 1.5,2,3,4,5
KT [107Jump To The Next Citation Point] PwPoly CF,M Ir Pu a = 0 Q = 2,3
KT [109Jump To The Next Citation Point] PwPoly CF,M Ir Pu a = 0,− 0.5,0.25,0.5,0.75 Q = 2,3,4,5
UIUC [62Jump To The Next Citation Point] Γ = 2 CF,M Ir ExPu a = 0 Q = 2,3
UIUC [63Jump To The Next Citation Point] Γ = 2 CF,M Ir ExPu a = 0 − 0.5,0.75 Q = 1,3,5
CCCW [58Jump To The Next Citation Point] Γ = 2 KS Ir Ex a = 0 Q = 1
CCCW [57Jump To The Next Citation Point] Γ = 2,2.75, SH KS Ir Ex a = 0.5 Q = 3
CCCW [74Jump To The Next Citation Point] Γ = 2 KS Ir Ex a = 0,0.5,0.9 Q = 3
          πœƒ = 0,20,40,60,80∘
LBPLI [41Jump To The Next Citation Point] Γ = 2 CF,M Ir Ex a = 0.5 Q = 5
AEI [154Jump To The Next Citation Point] Γ = 2 CF,M Ir ExPu a = 0 Q = 5

When a quasi-equilibrium configuration constructed by the excision approach is used as initial data for a simulation code based on the moving-puncture approach, one more step is needed before starting the simulations because the data inside the BH excision surface does not exist, and thus, one needs to fill the BH interior by extrapolating all field values radially from the excision surface toward the BH center. In the extrapolation procedure, one has to employ a method in which any constraint violation introduced inside the BH is not allowed to affect the exterior spacetime. Two groups proposed such a procedure [34, 61Jump To The Next Citation Point], and the UIUC group used the method they developed in their simulations [61].

Strictly speaking, the quasi-equilibrium derived in the framework of Section 2 is not realistic because an approaching radial velocity driven by gravitational radiation reaction is not taken into account. If we adopt such quasi-equilibrium data as the initial condition, the trajectories obtained in the numerical simulation usually result in a slightly eccentric orbit. One prescription to suppress this error is to choose an initial condition in which the binary separation is large enough that the effect of the radial velocity is not serious. Another method to reduce the eccentricity was proposed by the CCCW group [157, 75] as described in the following.

(1) At the beginning of their procedure, they prepare a quasi-equilibrium data, which has zero approaching velocity vr = aΛ™0r = 0 and orbital angular velocity Ω0. (2) The second step is to evolve the quasi-equilibrium initial data for at least one and a half orbits. Then one records the time derivative of the measured coordinate separation between the center of the compact objects, Λ™ d(t) and fits Λ™ d (t) to a function of the form

dΛ™= A0 + A1t + B sin(ωt + Ο•), (105 )
where the parameters A0, A1, B, ω, and Ο• are all determined by the fitting. The A0 + A1t part denotes the smooth inspiral, while the B sin(ωt + Ο•) part denotes the unwanted oscillations due to the eccentricity of the orbit. For a nearly circular Newtonian orbit, the eccentricity e of the orbit can be written as
B e = ----, (106 ) ωd0
where d0 = d(t = 0). This implies that reducing the orbital eccentricity is equivalent to reducing B. (3) The final step is to add the corrections to the approaching velocity Λ™a 0 and the orbital angular velocity Ω using the parameters, which appear in Equation (105View Equation). Such corrections should be chosen so that the eccentricity-induced initial radial velocity and radial acceleration can be removed. Namely, the initial radial velocity is changed as
δΛ™a = − B-sinΟ•-, (107 ) 0 d0
and the initial orbital angular velocity as
B ωcos Ο• δΩ0 = − -2d-Ω---. (108 ) 0 0
Using the corrected approaching radial velocity and the orbital angular velocity for the initial step (1), they iterated the procedure and obtained a better initial data. They showed that going twice through the iterative method described above, the orbital eccentricity is reduced by about an order of magnitude.

3.1.2 Evolving metric

There are currently two formalisms for evolving the spacetime metric. One is the Baumgarte–Shapiro–Shibata–Nakamura (BSSN) formalism [140Jump To The Next Citation Point, 195, 18] together with the moving-puncture gauge condition [39Jump To The Next Citation Point, 15Jump To The Next Citation Point, 35Jump To The Next Citation Point, 136], and the other is the generalized harmonic (GH) formalism [129Jump To The Next Citation Point] (see also [77, 86, 162Jump To The Next Citation Point, 163Jump To The Next Citation Point]). The AEI, KT, and UIUC groups employ the BSSN formalism, while the CCCW and LBPLI [6Jump To The Next Citation Point] groups use the GH formalism. These two formalisms are different from the standard 3+1 formalism [230Jump To The Next Citation Point], which was found to be inappropriate in numerical relativity because stable and longterm numerical simulations are not feasible due to the presence of an unphysical growth mode excited even by a tiny numerical truncation error.

The BSSN formalism, the original version of which was first proposed by Nakamura in 1987 [140], is in a sense a modified version of the 3+1 formalism [230]. The essence in this formalism is to adopt not only the three metric (γij) and the extrinsic curvature (Kij) but also a conformal factor of the three metric (ψ or Ο• = log (ψ) or W = ψ −2 or χ = ψ −4), the trace part of the extrinsic curvature (K), and the first spatial derivative of the three metric (F = γ&tidle; i ij,j or &tidle;Γ i = − &tidle;γij ,j where &tidle;γij is the conformal three metric) as new independent variables, and then, to appropriately rewrite the basic equations using these new variables. With this prescription, the evolution equations for the conformal three metric, &tidle;γij, and the conformal traceless extrinsic curvature, A&tidle;ij ≡ γ −1βˆ•3(Kij − K γijβˆ•3), are written in the form of simple wave equations. This prescription leads to avoiding the spurious growth of unphysical modes by numerical truncation errors and enables one to perform a stable and longterm evolution for a variety of systems. Typical basic equations are described in Appendix A.

The basic equations of the GH formalism are similar to those employed in the PN approximation [25Jump To The Next Citation Point]. Using the harmonic or generalized harmonic gauge condition,

∇ μ∇ μx ν = H ν, (109 )
where ν H is the arbitrary function, Einstein’s equation is written in a hyperbolic form of the spacetime metric gμν. The main modification, analogous to the introduction of Fi or &tidle;Γ i in the BSSN formalism, is treating Hν as independent functions. These may be either the functions themselves or determined by evolution equations. This is essential for the stable and longterm evolution of the system. The other modification is to introduce the first spatial derivative of the metric components as new independent variables Diμν = ∂igμν to make the GH formulation first-order [129Jump To The Next Citation Point].

In both formalisms, the number of constraint equations is increased as a result of defining the new variables, in addition to the Hamiltonian and momentum constraints. For the BSSN formalism, det(&tidle;γij) = 1, &tidle;γijA&tidle;ij = 0, and Fi − &tidle;γij,j = 0 or &tidle;Γ i + &tidle;γij,j = 0 are new constraints. For the GH formalism, D − ∂ g = 0 iμν iμν is a new constraint.

3.1.3 Evolving a black hole

For handling a BH, the original BSSN formalism was not satisfactory. Re-defining a spatial conformal factor [39Jump To The Next Citation Point], adopting more than fourth-order finite-differencing schemes, and employing an appropriate moving-puncture gauge condition [39, 15, 221, 35Jump To The Next Citation Point] are required for evolving BHs accurately and stably. In such a BSSN-puncture formulation, we choose a BH spacetime, which is free of the true singularity [32]. Although a coordinate singularity may appear in the center of the BH, this can be effectively (and, in a sense, fortunately) excised in the moving-puncture gauge condition [88]. Consequently, one can evolve the whole computational region without artificially excising inside the BH horizon (i.e., without special numerical treatments for the inside of the BH horizon).

On the other hand, in the code of the CCCW and LBPLI groups using GH formalism, the gauge, similar to the moving-puncture gauge, has not yet been developed. These groups employ the excision technique, i.e., a region inside the apparent horizon is cut out of the computation domain and replaced with an inner boundary. Pretorius has successfully simulated orbiting BH-BH binaries using the GH formalism with excision for the first time [162, 163]. Subsequently, the CCCW group has successfully simulated a longterm inspiral and merger of BH-BH binaries using a high-accuracy spectral method [181, 206], verifying that this technique is also robust for handling BH spacetimes.

3.1.4 Hydrodynamics

For the evolution of orbiting NS and a disk surrounding BHs formed after the merger, hydrodynamics (or magnetohydrodynamics) equations have to be solved. For a hydrodynamic simulation specifically, one has to solve the continuity equation for the rest-mass density ρ, the relativistic Euler and energy equations for the four velocity ui and internal energy πœ€ (or h), i.e.,

μ ∇ μ(ρu ) = 0, (110 ) γ ∇ Tμν = 0, (111 ) iν μ μν nν∇ μT = 0, (112 )
where Tμν is the stress-energy tensor. For the perfect fluid, it is written as
μν μ ν μν T = ρhu u + P g . (113 )
For a magnetohydrodynamics simulation, we have to solve Maxwell’s equation as well. In ideal magnetohydrodynamics, only the induction equation for magnetic-field components should be solved [73Jump To The Next Citation Point].

If one needs to follow the evolution of the lepton number densities, the following additional continuity equations for the leptons (the electron number density and/or the total lepton number density) have to be solved

∇ μ(ρYau μ) = ρSa, (114 )
where Ya denotes a lepton number fraction per baryon and Sa is the corresponding production/annihilation rates of the lepton in the baryon rest frame, which is associated with neutrino emission and absorption (e.g., [175, 176, 172, 173, 170, 146Jump To The Next Citation Point, 183Jump To The Next Citation Point, 184Jump To The Next Citation Point, 185Jump To The Next Citation Point]).

The hydrodynamic equation is schematically written in the form

i ∂U-- ∂F-- ∂t + ∂xi + SU = 0, (115 )
where U represents the set of the evolved variables and i F are the associated transport fluxes. The third term of Equation (115View Equation) denotes the source term, which is usually the coupling term between the hydrodynamic variables and geometric quantities or the first derivative of the metric. The third term is evaluated in a straightforward manner, which usually does not cause numerical instabilities in a straightforward evaluation (as long as an appropriate time step is chosen). A careful treatment is required in handling the transport term [73]. The numerical scheme for the transport term employed in all the groups is similar. The CCCW, LBPLI, and UIUC groups employ the scheme of Harten, Lax, and van Leer (HLL) [89], the AEI group employs several schemes prepared in their Whisky library [13, 14Jump To The Next Citation Point, 11, 12], and the KT group employs the scheme of Kurganov and Tadmor [105]. In handling the transport term, one further needs to determine the method of the interpolation by which the values of hydrodynamic quantities at the cell surfaces are specified. All these groups currently employ the third-order piecewise parabolic interpolation.

3.1.5 Equations of state

Employing realistic EOS for modeling NS is a key ingredient in accurately determining the final fate after the onset of the merger and in deriving realistic gravitational waveforms in the merger phase. However, the EOS of the nuclear matter beyond the normal nuclear density is still uncertain due to the lack of strong constraints obtained from experiments and observations (e.g., [115, 116Jump To The Next Citation Point, 117] and see [51Jump To The Next Citation Point] for a recent devlopment). Hence, we do not know exactly the detailed physical state of NS such as density and composition profiles as functions of radius for a given mass and the relation between the mass and the radius.

Gravitational waveforms emitted in the tidal-disruption phase as well as properties of the remnant BH-disk system such as the mass and the typical density of the disk depend strongly on the EOS, as shown in the following. The primary reason for this is that these depend strongly on the relation between the mass and the radius of the NS, the sensitivity of a NS to the tidal force by its companion BH depends primarily on its radius, e.g., a NS of larger radii (with a stiffer EOS) will be disrupted at a larger orbital separation (or a lower orbital frequency; see equations in Section 1); if the tidal disruption of a NS occurs at a larger distance, more material will be widely spread around its companion BH, and consequently, a high-mass remnant disk is likely to be formed; also, the gravitational-wave frequency at tidal disruption, which will be one of the characteristic frequencies, is lower for a NS of larger radius.

In the study of binaries composed of NS, rather than require one to prepare a “realistic” EOS, we should consider exploring the possibility of determining the NS EOS by gravitational-wave observation, as discussed in [127, 220, 168Jump To The Next Citation Point, 70]. For this purpose, we need to prepare theoretical templates of gravitational waves employing a wide variety of possible EOS for the NS matter, and systematically perform a wide variety of simulations employing a large number of plausible EOS. Then, the next task is to determine what kinds of EOS should be employed in their theoretical work.

Thermal energy per baryon inside NS, except for newly-born ones, is believed to be much lower than the Fermi energy of the constituent particles, because a young NS is quickly cooled down by neutrino emission (e.g., [116, 99]). This implies that we can safely neglect thermal effects on the NS in the inspiral and early merger phases, and can employ a cold EOS, for which the pressure, P, the specific internal energy, πœ€, and other thermodynamic quantities are written as functions of the rest-mass density ρ. In particular, employing the cold EOS is appropriate for the case in which the merger does not result in the formation of a disk surrounding the remnant BH, i.e., for the case in which the NS is simply swallowed by the companion BH, because a heating process, such as shock heating, never plays a role. By contrast, to follow a longterm evolution of tidally disrupted material, which forms a disk around the remnant BH, the finite-temperature effects of the EOS and the neutrino cooling have to be taken into account. This is because the temperature of the disk material is likely to be increased significantly by shock or a longterm viscous heating, while the density of the disk is lower than that of the NS, and as a result, the Fermi energy is lowered, resulting in a situation where the roles of the pressures by degenerate electrons and thermal gas become less important and increasingly important, respectively.

Because the simulation for BH-NS binaries in general relativity is still in an early phase (by 2010), much work has been done with a rather simple (not very realistic) EOS. The often-used EOS in the early phase of this study within the numerical-relativity community is the Γ-law EOS,

P = (Γ − 1)ρπœ€, (116 )
with Γ = 2. In this EOS, the initial condition for the NS is determined by the polytropic EOS, P = κρ Γ (see Section 2). However, this EOS is known to disagree with typical nuclear-theory-based EOS for a small value of Γ ∼ 2 (e.g., [107Jump To The Next Citation Point] and Figure 10View Image).

A better choice is to employ a piecewise polytropic EOS. This is a phenomenologically parameterized EOS, which approximately reproduces cold nuclear-theory-based EOS at a high density (above the nuclear density) only with a small number of polytropic constants and indices [167Jump To The Next Citation Point, 168Jump To The Next Citation Point, 150Jump To The Next Citation Point] (see also [87, 128]), i.e.,

Γ i P (ρ) = κiρ forρi−1 ≤ ρ < ρi(1 ≤ i ≤ n), (117 )
where n is the number of the pieces used to parameterize an EOS and ρi denote boundary densities for which appropriate characteristic values are assigned. Here, ρ0 = 0 and ρn → ∞. κi are the polytropic constants and Γ i the adiabatic indices for each piece.

At each boundary density, ρ = ρi (i = 1,⋅ ⋅⋅,n − 1), the pressure is required to be continuous, i.e., κiρΓ i = κi+1ρΓ i+1 i i. Thus, if one gives κ1, Γ i, and ρi (i = 1,⋅⋅⋅,n), the EOS is totally determined. For the zero-temperature EOS, the first law of thermodynamics (dπœ€ = − Pd (1 βˆ•ρ) or dh = dP βˆ•ρ) holds, and thus, πœ€ and h are also determined except for the choice of the integration constants, which are fixed by the continuity condition of πœ€ (hence equivalently h) at each ρi.

Recently, it has been shown that the piecewise polytropic EOS composed of one piece in the crust region and three pieces in the core region approximately reproduces most of nuclear-theory-based EOS at a high density [167Jump To The Next Citation Point]. Here, three pieces in the core region are required to reproduce a high-mass NS for which inner and outer cores could have different stiffness due to the variation of the properties of the high-density nuclear matter (ρ ≳ 2 × 1014 gβˆ•cm3). Thus, if one focuses on a NS of relatively low mass, a smaller number of pieces is acceptable; see [168Jump To The Next Citation Point, 107Jump To The Next Citation Point] for the two-pieces case. Table 4 lists the parameters employed in [168Jump To The Next Citation Point, 107Jump To The Next Citation Point], and Figure 10View Image shows the relation between the mass and the radius of a spherical NS for these piecewise polytropic EOS. We note that using a single polytrope for the low-density EOS is justified to the extent that the radius and deformability of the NS as well as resulting gravitational waveforms in the merger phase are insensitive to the low-density EOS.


Table 4: The parameters and key ingredients for a piecewise polytropic EOS employed in [107Jump To The Next Citation Point]. Γ 2 is the adiabatic index in the core region and p is the pressure at the fiducial density 14.7 3 ρfidu = 10 gβˆ•cm, which determines the polytropic constant κ2 of the core region and ρ1: the critical rest-mass density separating the crust and core regions. Mmax is the maximum mass of a spherical NS for a given EOS. R135 and π’ž135 are the circumferential radius and the compactness of the NS with MNS = 1.35M βŠ™.
Model Γ 2 log10p [3 gβˆ•cm] ρ1 [14 3 10 gβˆ•cm] Mmax [M βŠ™ ] R135 [km] π’ž135
2H 3.0 13.95 0.7033 2.835 15.23 0.1309
H 3.0 13.55 1.232 2.249 12.27 0.1624
HB 3.0 13.45 1.417 2.122 11.61 0.1718
HBs 2.7 13.45 1.069 1.926 11.57 0.1723
HBss 2.4 13.45 0.6854 1.701 11.45 0.1741
B 3.0 13.35 1.630 2.003 10.96 0.1819
Bs 2.7 13.35 1.269 1.799 10.74 0.1856
Bss 2.4 13.35 0.8547 1.566 10.27 0.1940

View Image

Figure 10: The relation between the mass and the circumferential radius of a spherical NS with piecewise polytropic EOS for which parameters are described in Table 4. For comparison, the curve for the Γ = 2 polytropic EOS with 2 −16 κβˆ•c = 2 × 10 in cgs units is also plotted (dotted curve). The figure is taken from [107Jump To The Next Citation Point]. Note that if the observational result of [51] (if the presence of an ≈ 1.97M βŠ™ NS) is confirmed by 100%, some of the soft EOS displayed here will be ruled out.

In the presence of shocks during the merger phase, the assumption of zero-temperature breaks down. Thus, in addition to the piecewise polytropic part, a correction term has to be added. An often-used minimum prescription is to simply add the following term to the pressure [199, 168, 107Jump To The Next Citation Point]

Pth = (Γ th − 1)ρπœ€th, (118 )
where Γ th is an adiabatic index for the thermal part, and πœ€th = πœ€ − πœ€cold with πœ€cold being determined by the piecewise polytropic EOS. This type of the piecewise polytropic EOS is employed by the KT group for a wide variety of parameters [107Jump To The Next Citation Point, 109Jump To The Next Citation Point].

To self-consistently take into account thermal effects, which play an important role in the evolution of a disk formed after tidal disruption occurs, the best way is probably to employ a finite-temperature EOS, which is derived by a detailed model based on nuclear physics. This type of EOS is usually described in table form, e.g., as

P = P (ρ,T,Ye )and πœ€ = πœ€(ρ,T, Ye), (119 )
where T and Ye are the temperature and the electron fraction [119, 188Jump To The Next Citation Point, 189Jump To The Next Citation Point]. Here, approximately speaking, ρ is determined by the continuity equation (110View Equation), and T by the internal energy, which is essentially determined by the energy equation (112View Equation) in general relativity. Ye is determined by the continuity equation for the electron fraction (114View Equation). To solve this equation correctly, we have to take into account neutrino emission/absorption, by which the electron fraction is varied. Multidimensional numerical-relativity simulations have reached a level of incorporating the effects of neutrino emission/absorption only quite recently [183Jump To The Next Citation Point, 185Jump To The Next Citation Point]. Rather, most multidimensional numerical-relativity simulations have been performed without solving neutrino transfer equations and the equation for Ye [146, 53, 54, 74Jump To The Next Citation Point, 149], but an a priori assumption for Ye as a function of other variables such as ρ and T is imposed. Along this line (under the assumption that Ye is unchanged in the fluid rest frame or Ye is determined by the β-equilibrium), the CCCW group performed the first simulation taking into account the finite-temperature effect [57Jump To The Next Citation Point]. Recently, Sekiguchi has developed a leakage scheme of neutrinos in general relativistic simulation [183, 184, 185] and solves the equation for Ye taking into account neutrino emission for the first time. This technique will be applied to the simulation of BH-NS binaries in the near future.

3.1.6 Adaptive mesh refinement

There are two important length scales in the problem of compact binary coalescence. One is the scale associated with the size of compact objects, which is ∼ GMBH βˆ•c2 for BH and RNS ∼ 5– 8GMNS βˆ•c2 for NS. The other is the wavelength of gravitational waves, λ. For a binary in a circular orbit of angular velocity Ω, the wavelength for the dominant quadrupole mode (l = |m | = 2 mode) is

( −3 )− 1 λ = πc-≈ 105Gm0-- Gc---Ωm0-- , (120 ) Ω c2 0.03
where m0 = MBH + MNS. A longterm numerical-relativity simulation is typically performed with initial angular velocity − 3 Gc Ωm0 ≈ 0.02 for NS-NS binaries, and ∼ 0.03 for BH-NS binaries. Thus, λ is larger than Gm0 βˆ•c2 by a factor of ≳ 100.

To accurately compute the evolution of an orbiting BH and NS resolving the structure of each object, the grid spacing should be smaller than 2 Δ0 ∼ GMBH βˆ•20c and ∼ RNS βˆ•40, respectively. In addition, the outer boundary along each axis of the computation domain should be larger than the wavelength of gravitational waves to accurately take into account gravitational radiation reaction and to accurately extract gravitational waves in the wave zone. Then, if one adopts a uniform grid with grid spacing Δ0, the total grid number in one direction should be larger than

λ [ ( Gc −3Ωm ) −1 --- ≈ max 2100 (1 + Q− 1) ---------0 , Δ0 0.03 ( c2RNS βˆ•GMNS )− 1( Gc− 3Ωm0 ) −1 ] 700(1 + Q ) -------------- ---------- . (121 ) 6 0.03
In three spatial dimensions, the required grid number is more than (2λβˆ•Δ0 )3. Thus, more than 10003 in the grid number is necessary. For such a simulation, a quite high computational cost is required even in the present computational resources. It should also be pointed out that in such studies, a survey for a large parameter space composed of the mass and the spin of the BH, and the mass and EOS of the NS is necessary (see Section 3.2). For this purpose, it is very important to make every effort to save computational costs for each run.

Therefore, to efficiently perform a large number of numerical simulations with the finite-differencing method, adaptive mesh refinement (AMR), such as that proposed by Berger and Oliger [23], is an indispensable technique. In this technique, one prepares several refinement levels of Cartesian boxes (or other geometrical domains) with different grid resolutions; usually in the smaller boxes, the grid resolution is higher. At such a fine refinement level, the BH and the NS are evolved with sufficiently high resolution with grid spacing ≲ Δ0. On the other hand, the propagation of gravitational waves, which have a wavelength much longer than the size of the compact objects, is followed in large-size boxes (in the coarser levels) with a large grid spacing, which is still much smaller than the gravitational wavelength (say, ∼ λβˆ•10).

Mesh refinement techniques in numerical relativity have been developed by several groups (e.g., [93, 182Jump To The Next Citation Point, 35, 228Jump To The Next Citation Point]). The AEI and UIUC groups use the Carpet module of the Cactus code [182] developed primarily by Schnetter, the KT group the SACRA code [228Jump To The Next Citation Point] developed by Yamamoto, and the LBPLI group the HAD code [7]. Recently, two independent codes (Whisky, which employs Carpet, and SACRA) were compared, performing simulations of NS-NS binaries [14], and it was illustrated that the AMR scheme in Carpet and SACRA work well at the same level. The CCCW group employs a two-grid technique, which was originally proposed by the Meudon–Valencia–MPA team [52]. In this approach, geometric-field quantities are evolved with the multi-patch pseudo-spectral method and the hydrodynamic equations are solved using a standard finite-volume scheme on a second grid [58Jump To The Next Citation Point]. In solving the equations for the geometric field, the computational domain encompassing the local wave zone is prepared, while the hydrodynamic equations are solved in the region where the matter presents.


  Go to previous page Go up Go to next page