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:
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.
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.
The notation of “Co” means a corotating NS, and that of “Ir” an irrotational NS.
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).
|Group & Ref.||(1)||(2)||(3)||(4)||(5)||(6)|
|CCCW ||, SH||KS||Ir||Ex|
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 .
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
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  (see also [77, 86, 162, 163]). The AEI, KT, and UIUC groups employ the BSSN formalism, while the CCCW and LBPLI  groups use the GH formalism. These two formalisms are different from the standard 3+1 formalism , 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 , is in a sense a modified version of the 3+1 formalism . 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 . Using the harmonic or generalized harmonic gauge condition,.
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 , 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 . 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 . 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.,.
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[175, 176, 172, 173, 170, 146, 183, 184, 185]).
The hydrodynamic equation is schematically written in the form. 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) , 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 . 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  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, 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.,
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 . 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.
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][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[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 . 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
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
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 , 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  developed primarily by Schnetter, the KT group the SACRA code  developed by Yamamoto, and the LBPLI group the HAD code . Recently, two independent codes (Whisky, which employs Carpet, and SACRA) were compared, performing simulations of NS-NS binaries , 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 . 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 . 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
This work is licensed under a Creative Commons License.