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 “” or “ ” means the single polytropic EOS with an adiabatic index of or .

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

- (2)
- The spatial background metric and the extrinsic curvature
The notation of “CF” means the conformally-flat condition , and that of “M” means the maximal slicing .

The notation of “KS” indicates the Kerr–Schild metric. In this case, the extrinsic curvature 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 and misalignment angle between the spin axis and orbital angular momentum vector.
- (6)
- The mass ratio

Group & Ref. | (1) | (2) | (3) | (4) | (5) | (6) |

KT [202] | CF,M | Co | Pu | |||

KT [203] | CF,M | Co | Pu | |||

KT [197] | CF,M | Ir | Pu | |||

KT [228] | CF,M | Ir | Pu | |||

KT [194] | CF,M | Ir | Pu | |||

KT [107] | PwPoly | CF,M | Ir | Pu | ||

KT [109] | PwPoly | CF,M | Ir | Pu | ||

UIUC [62] | CF,M | Ir | ExPu | |||

UIUC [63] | CF,M | Ir | ExPu | |||

CCCW [58] | KS | Ir | Ex | |||

CCCW [57] | , SH | KS | Ir | Ex | ||

CCCW [74] | KS | Ir | Ex | |||

LBPLI [41] | CF,M | Ir | Ex | |||

AEI [154] | CF,M | Ir | ExPu | |||

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, 61], 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 and orbital angular velocity . (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, and fits to a function of the form

where the parameters , , , , and are all determined by the fitting. The part denotes the smooth inspiral, while the part denotes the unwanted oscillations due to the eccentricity of the orbit. For a nearly circular Newtonian orbit, the eccentricity of the orbit can be written as where . This implies that reducing the orbital eccentricity is equivalent to reducing . (3) The final step is to add the corrections to the approaching velocity and the orbital angular velocity using the parameters, which appear in Equation (105). 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 and the initial orbital angular velocity as 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.There are currently two formalisms for evolving the spacetime metric. One is the Baumgarte–Shapiro–Shibata–Nakamura (BSSN) formalism [140, 195, 18] together with the moving-puncture gauge condition [39, 15, 35, 136], and the other is the generalized harmonic (GH) formalism [129] (see also [77, 86, 162, 163]). The AEI, KT, and UIUC groups employ the BSSN formalism, while the CCCW and LBPLI [6] groups use the GH formalism. These two formalisms are different from the standard 3+1 formalism [230], 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 () and the extrinsic curvature () but also a conformal factor of the three metric ( or or or ), the trace part of the extrinsic curvature (), and the first spatial derivative of the three metric ( or where 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, , and the conformal traceless extrinsic curvature, , 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 [25]. Using the harmonic or generalized harmonic gauge condition,

where is the arbitrary function, Einstein’s equation is written in a hyperbolic form of the spacetime metric . The main modification, analogous to the introduction of or in the BSSN formalism, is treating 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 to make the GH formulation first-order [129].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(, , and or are new constraints. For the GH formalism, is a new constraint.

For handling a BH, the original BSSN formalism was not satisfactory. Re-defining a spatial conformal factor [39], adopting more than fourth-order finite-differencing schemes, and employing an appropriate moving-puncture gauge condition [39, 15, 221, 35] 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.

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 and internal energy (or ), i.e.,

where is the stress-energy tensor. For the perfect fluid, it is written as 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 [73].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

where denotes a lepton number fraction per baryon and 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, 146, 183, 184, 185]).The hydrodynamic equation is schematically written in the form

where represents the set of the evolved variables and are the associated transport fluxes. The third term of Equation (115) 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, 14, 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.

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, 116, 117] and see [51] 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, 168, 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, , 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,

with . In this EOS, the initial condition for the NS is determined by the polytropic EOS, (see Section 2). However, this EOS is known to disagree with typical nuclear-theory-based EOS for a small value of (e.g., [107] and Figure 10).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 [167, 168, 150] (see also [87, 128]), i.e.,

where is the number of the pieces used to parameterize an EOS and denote boundary densities for which appropriate characteristic values are assigned. Here, and . are the polytropic constants and the adiabatic indices for each piece.At each boundary density, , the pressure is required to be continuous, i.e., . Thus, if one gives , , and , the EOS is totally determined. For the zero-temperature EOS, the first law of thermodynamics ( or ) holds, and thus, and are also determined except for the choice of the integration constants, which are fixed by the continuity condition of (hence equivalently ) at each .

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 [167]. 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 (). Thus, if one focuses on a NS of relatively low mass, a smaller number of pieces is acceptable; see [168, 107] for the two-pieces case. Table 4 lists the parameters employed in [168, 107], and Figure 10 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.

Model | [] | [] | [km] | |||

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 |

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, 107]

where is an adiabatic index for the thermal part, and with 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 [107, 109].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

where and are the temperature and the electron fraction [119, 188, 189]. Here, approximately speaking, is determined by the continuity equation (110), and by the internal energy, which is essentially determined by the energy equation (112) in general relativity. is determined by the continuity equation for the electron fraction (114). 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 [183, 185]. Rather, most multidimensional numerical-relativity simulations have been performed without solving neutrino transfer equations and the equation for [146, 53, 54, 74, 149], but an a priori assumption for as a function of other variables such as and is imposed. Along this line (under the assumption that is unchanged in the fluid rest frame or is determined by the -equilibrium), the CCCW group performed the first simulation taking into account the finite-temperature effect [57]. Recently, Sekiguchi has developed a leakage scheme of neutrinos in general relativistic simulation [183, 184, 185] and solves the equation for 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.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 for BH and 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 ( mode) is

where . A longterm numerical-relativity simulation is typically performed with initial angular velocity for NS-NS binaries, and for BH-NS binaries. Thus, is larger than by a factor of .To accurately compute the evolution of an orbiting BH and NS resolving the structure of each object, the grid spacing should be smaller than and , 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 , the total grid number in one direction should be larger than

In three spatial dimensions, the required grid number is more than . Thus, more than 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 . 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, ).

Mesh refinement techniques in numerical relativity have been developed by several groups (e.g., [93, 182, 35, 228]). The AEI and UIUC groups use the Carpet module of the Cactus code [182] developed primarily by Schnetter, the KT group the SACRA code [228] 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 [58]. 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.

Living Rev. Relativity 14, (2011), 6
http://www.livingreviews.org/lrr-2011-6 |
This work is licensed under a Creative Commons License. E-mail us: |