2.1 Formulation

For a solution of BH-NS binaries in quasi-equilibrium, it is at least necessary to solve Einstein’s constraint equations and relativistic hydrostationary equations, supplemented by NS EOS. To derive a better model of quasi-equilibrium states, a part of Einstein’s evolution equation is also solved with an additional assumption. Until now, two approaches have been proposed to constructing BH-NS binaries in quasi-equilibrium. The main difference in those approaches comes from the method for handling the BH singularity. In one approach the BH singularity is eliminated from the computational domain by excising a coordinate sphere and by imposing equilibrium BH boundary conditions [47Jump To The Next Citation Point] (see also [10Jump To The Next Citation Point, 81Jump To The Next Citation Point] for a related isolated-horizon formalism). In the other approach, a BH is handled as a puncture [32Jump To The Next Citation Point]. In this approach, the metric quantities are decomposed into a singular part, which is written analytically and denotes contributions from a BH and a regular part, which is obtained by numerically solving the corresponding basic equations.

These two methods are separately explained below because they give different formulas. We refer to [208Jump To The Next Citation Point, 209Jump To The Next Citation Point, 210Jump To The Next Citation Point] for the excision approach and to [202Jump To The Next Citation Point, 203Jump To The Next Citation Point, 106Jump To The Next Citation Point] for the puncture one. For a more detailed discussion of the decomposition of Einstein’s equation and the formalism, we refer to [29, 231Jump To The Next Citation Point, 44], and for the hydrostatics [80Jump To The Next Citation Point, 214Jump To The Next Citation Point, 219Jump To The Next Citation Point].

2.1.1 Formulation in the excision approach

In the excision approach, Einstein’s equations in the conformal thin-sandwich formalism are solved for constructing quasi-equilibrium configurations [231Jump To The Next Citation Point]. The line element in the 3+1 form is written as

2 μ ν ds = gμνdx dx ( ) = − α2dt2 + γij dxi + βicomdt (dxj + βjcomdt), (22 )
where i βcom is the shift vector seen by a comoving observer. We note that the coordinates in the comoving observer frame are usually chosen for convenience in this field (e.g., [80Jump To The Next Citation Point, 78]).

The spatial metric γij is further decomposed into a conformal factor ψ and a background spatial metric &tidle;γ ij as

γ = ψ4γ&tidle; . (23 ) ij ij
Here, the condition, det(&tidle;γij) = 1, is not always imposed. The shift vector seen by a comoving observer, βi com, can be decomposed as
βi = βi + βi , (24 ) com rot
where i β is the shift vector seen by an inertial observer and i βrot is a rotating shift. The rotating shift is written as βirot = (Ω × R )i, where Ω is the orbital angular velocity vector of the binary system measured at infinity, and R a coordinate vector for which the origin is located at the center of mass of the binary system.

The extrinsic curvature is defined by

1 Kij = − --ℒnγij, (25 ) 2
where ℒn is the Lie derivative along the unit normal to the hypersurface Σt. Following [229, 145], it is split into the trace K and the traceless part Aij as
1 Kij = Aij + -γijK. (26 ) 3
The traceless part is rewritten as
−2 ˆ ij − 10 ˆij Aij = ψ Aij, A = ψ A , (27 )
where ˆAij = &tidle;γikγ&tidle;jlAˆkl. The Hamiltonian constraint then takes the form
5 1 1 5 2 1 − 7 ij Δ&tidle;ψ = − 2πψ ρH + -ψ &tidle;R + --ψ K − --ψ ˆAij ˆA , (28 ) 8 12 8
where &tidle;Δ denotes &tidle;γijD&tidle;i &tidle;Dj, D&tidle;i the covariant derivative with respect to &tidle;γij, and R&tidle; the scalar curvature with respect to &tidle;γij.

Equations (25View Equation), (26View Equation), and (27View Equation) yield

( ) ij ψ6 ij i j j i 2 ij k Aˆ = --- ∂t&tidle;γ + &tidle;D βcom + D&tidle; βcom − -γ&tidle; D&tidle;k βcom . (29 ) 2 α 3
Inserting Equation (29View Equation) into the momentum constraint gives
( ) ( 6 ) &tidle; i 1- ik &tidle; &tidle; j -α- &tidle; ψ-- ij 4 i 4- &tidle;i Δ βcom + 3 &tidle;γ Dk Dj βcom = − ψ6 Dj α ∂tγ&tidle; + 16πα ψ j + 3α D K ( ) ( 6) − D&tidle;i βj + &tidle;Djβi − 2&tidle;γij &tidle;Dk βk -α-&tidle;Dj ψ-- . (30 ) com com 3 com ψ6 α
In addition to the Hamiltonian and momentum constraints, the trace of the evolution equation of the extrinsic curvature is often employed as one of the field equations:
( ) [ 2] ∂tK − ℒβK = − ψ −4 &tidle;D2α + 2 &tidle;Di ln ψD&tidle;i α + α 4 π(ρH + S ) + ψ −12 ˆAij ˆAij + K . (31 ) 3
For a given condition for K and ∂ K t, this equation reduces to an elliptic-type equation for α, which possesses primary information of gravitational potential for an equilibrium or quasi-equilibrium configuration.

The matter terms in the right-hand side of Equations (28View Equation), (30View Equation), and (31View Equation) are derived from the projections of the stress-energy tensor T μν into the spatial hypersurface Σt, defined by

ρH = nμn νTμν, (32 ) ji = − γin T μν, (33 ) μ ν μν Sij = γiμγjνT , (34 ) S = γijS . (35 ) ij
For an ideal fluid for which
T μν = (ρ + ρ 𝜀 + P )uμu ν + P gμν = ρhu μuν + P gμν, (36 )
the following relation holds:
ρH = ρh (αut)2 − P, (37 ) t ji = ρh (αu )ui, (38 ) Sij = ρhuiuj + Pγij. (39 )

The set of Equations (28View Equation) – (31View Equation) has four functions that can be chosen freely; ∂t&tidle;γij, ∂tK, &tidle;γij, and K. For computing quasi-equilibrium states, one usually assumes the presence of a helical Killing vector, μ ξ, and the absence of gravitational waves in the wave zone. Under these assumptions, it is natural to choose the time direction so as to satisfy ξμ = (∂∕∂t)μ (in the comoving frame), and to set ∂t&tidle;γij and ∂tK = 0. For the choice of K and &tidle;γij, two ways have been proposed. The first one is to choose a maximal slicing K = 0 and to adopt a flat metric &tidle;γij = fij, for simplicity. The other choice is to use the Kerr–Schild metric for &tidle;γij and K in the vicinity of the BH or in the whole computational space. One of the advantages of choosing maximal slicing and the flat metric background is that the source term becomes simple and falls off rapidly for r → ∞ enough to obtain accurate results. The disadvantage is that it is not possible to construct the Kerr BH even with a distant orbit, and moreover, the set of Equations (28View Equation) – (31View Equation) may have non-unique solutions [159, 17, 224]. In particular, the spin of a BH has two values for the same spin parameter Ωr (see Equation (49View Equation) for the definition) [132Jump To The Next Citation Point]. The lower branch of the spin as a function of Ωr, which should be physically reasonable because it approaches to the Schwarzschild BH as Ωr → 0, can reach only a ≈ 0.85, much less than the maximum spin of a Kerr BH. The advantage of using the Kerr–Schild metric for &tidle;γ ij and K is that one can calculate a spinning BH with nearly maximum spin, a ≃ 1. The disadvantage is that the source term becomes complicated and falls off slowly for r → ∞. Because of this situation, it is not easy to derive the results as accurately as those with the conformal three metric, if one adopts this background metric in the whole computational space [208Jump To The Next Citation Point]. To recover the accuracy, a modification of the Kerr–Schild metric seems to be necessary; the metric is chosen to be nearly the Kerr–Schild one in the vicinity of the BH, whereas away from the BH, the metric should approach exponentially a conformally-flat metric and a maximal slicing [75Jump To The Next Citation Point, 132Jump To The Next Citation Point].

In the following, we restrict ourselves to the case of maximal slicing and flat spatial background metric for simplicity. Then, Equations (28View Equation), (30View Equation), and (31View Equation) can be written as

5 1- −7 ˆ ˆij Δ-ψ = − 2πψ ρH − 8 ψ AijA , (40 )
Δ βi + 1-∂i(∂ βj) = 16π Φψ3ji + 2Aˆij ∂ (Φ ψ−7) , (41 ) -- 3 j j
7 Δ-Φ = 2πΦ ψ4(ρH + 2S ) + -Φ ψ− 8 ˆAij ˆAij, (42 ) 8
where Δ- and ∂i denote the flat Laplace operator and the partial derivative, and Φ ≡ αψ. Note here that the shift vector βicom was replaced by βi in Equation (41View Equation), because the rotating shift βi rot does not contribute to the equation in the conformally-flat case. Equation (29View Equation) becomes
7 ( ) Aˆij = ψ-- ∂iβj + ∂jβi − 2f ij∂kβk , (43 ) 2Φ 3
where we have replaced i βcom by i β for the same reason as above.

To solve the gravitational field equations (40View Equation), (41View Equation), and (42View Equation), it is necessary to impose appropriate boundary conditions on two different boundaries in the excision approach: outer boundaries at spatial infinity and inner boundaries on the BH horizons. Assuming asymptotic flatness, the boundary conditions at spatial infinity are written as

ψ |r→ ∞ = 1, (44 )
i| β |r→∞ = 0, (45 )
Φ|r→∞ = 1. (46 )
Inner boundary conditions arise from the excision of the BH interior. The assumption that the BH is in equilibrium leads to a set of boundary conditions for the conformal factor and the shift vector [47Jump To The Next Citation Point] (see also [40Jump To The Next Citation Point, 45] as well as the related isolated horizon formalism, e.g., [10, 81]). The boundary condition for the conformal factor is
| ( ) | &tidle;skD&tidle; ln ψ| = − 1- &tidle;hijD&tidle; &tidle;s − ψ2L || , (47 ) k |𝒮 4 i j |𝒮
where i −2 i s ≡ ψ &tidle;s is the outward pointing unit vector normal to the excision surface and hij is the induced metric on the excision surface, hij ≡ ψ4&tidle;hij = γij − sisj. The quantity L is computed from the projection of the extrinsic curvature Kij as L ≡ hijKij. The boundary condition on the normal component of the shift vector is
β⊥ |𝒮 = α|𝒮 . (48 )
Note that the quantity β⊥ is the the normal component of the shift vector in the comoving frame, i βcom. The tangential component must form a conformal Killing vector of the conformal metric &tidle;hij on the excision surface. This can be achieved by choosing them to be Killing vectors of a 2-sphere,
βi|| = Ω &tidle;ξi, (49 ) ∥ 𝒮 r
where Ωr is an arbitrary parameter related to the spin of the BH, and &tidle;ξi is derived by solving the conformal Killing equation on &tidle;h ij. Because we choose a conformally-flat spatial background, Equation (49View Equation) can be rewritten as
i|| i j k β∥ 𝒮 = 𝜖jkΩ rx . (50 )
Here Ωj r is a freely specifiable vector, related to the BH spin, and xk is a Cartesian coordinate centered on the 2-sphere. The quasi-local spin angular momentum associated with an approximate Killing vector of hij is defined by
1 ∮ S(ξ) = --- (Kij − γijK )ξjd2Si, (51 ) 8 π 𝒮
where ξi is the approximate Killing vector that is found by solving the Killing transport equations as described in [55Jump To The Next Citation Point, 40Jump To The Next Citation Point] ([55Jump To The Next Citation Point] described the original formulation and [40Jump To The Next Citation Point] reported a numerical result for BH-BH binaries based on the formulation of [55, 47Jump To The Next Citation Point]), or by a new alternative method for finding approximate Killing vectors of closed 2-spheres [48, 132]. Ωjr is iteratively determined until the quasi-local spin angular momentum of BH relaxes to a required value [40Jump To The Next Citation Point].

According to [47], the boundary condition on the lapse function can be chosen freely. For example, we can choose a Neumann boundary condition

| &tidle;si∂iΦ |𝒮 = 0 (52 )
on the excision surface.

2.1.2 Formulation in the puncture approach

A “puncture” method [32Jump To The Next Citation Point] was proposed by Brandt and Brügmann to describe multiple BH with arbitrary linear momenta and spin angular momenta, extending the original work by Brill and Lindquist [33]. A “moving-puncture approach” [39Jump To The Next Citation Point, 15Jump To The Next Citation Point] was revealed to be quite useful in dynamic simulations. Here, we describe the puncture approach for quasi-equilibrium in the context of BH-NS binaries, which was originally developed by Shibata and UryÅ« [202Jump To The Next Citation Point, 203Jump To The Next Citation Point] and subsequently modified by Kyutoku et al. [106Jump To The Next Citation Point]. The puncture approach employs a mixture of the conformal thin-sandwich decomposition and the conformal transverse-traceless decomposition of Einstein’s equations. The trace part of the extrinsic curvature is set to zero (maximal slicing) and the three metric is assumed to be conformally flat in the work so far.

The basic equations for the gravitational field are Equations. (40View Equation), (41View Equation), and (42View Equation) as in the excision approach. In the puncture approach, the metric quantities, which appear in Equations (40View Equation) and (42View Equation), ψ and Φ, are decomposed into an analytic singular part and a regular part. The former part denotes the contribution from a BH and the latter one is obtained by numerically solving the basic equations. Assuming that the puncture is located at rP = xkP, ψ and Φ are given by

-MP-- ψ = 1 + 2r + ϕ, (53 ) BH Φ = 1 − M-Φ-+ η, (54 ) rBH
where M P and M Φ are positive constants of mass dimension, and r = |xk − xk | BH P is a coordinate distance from the puncture. ϕ and η denote the regular parts of ψ and Φ, respectively, and are determined by solving elliptic equations (see Equations (59View Equation) and (61View Equation) shown below). The quantity MP is an arbitrarily-chosen parameter called the puncture mass, while M Φ is determined by the condition that the ADM mass (MADM) and the Komar mass agree,
∮ ∮ ∂iΦdSi = − ∂iψdSi = 2πMADM. (55 ) r→∞ r→∞
Also, ˆ Aij is decomposed into singular and regular parts as
Aˆ = D&tidle; W + &tidle;D W − 2f &tidle;D W k + KP . (56 ) ij i j j i 3 ij k ij
Here KP ij is the singular part, which denotes a weighted extrinsic curvature associated with the linear momentum, BH Pj, and spin, BH S j, of the BH written by
KP = -3---[lP BH + l PBH − (f − ll )lkPBH ] ij 2r2BH i j j i ij i j k 3 + --3- [li𝜖jkn + lj𝜖ikn]SkBHln. (57 ) rBH
lk = xk∕rBH P and 𝜖ijk is the completely anti-symmetric tensor on Σt.

Wi denotes an auxiliary three-dimensional function and i ij W = f Wj. The equation for Wi is derived by substituting Equation (56View Equation) into the momentum constraint equation. Because the total linear momentum of the binary system should vanish1, the linear momentum of the BH, P BiH, is related to that of the companion NS as

∫ P BH = − j ψ6dV, (58 ) i i
where the right-hand side denotes the (minus) linear momentum of the NS.

The field equations to be solved are summarized as follows:

5 1- −7 ˆ ˆij Δ-ϕ = − 2πψ ρH − 8ψ AijA , (59 )
i 1 i( j) 3 i ij ( −7) Δ-β + --∂ ∂jβ = 16π Φψ j + 2Aˆ ∂j Φ ψ , (60 ) 3
7 Δ-η = 2 πΦψ4 (ρH + 2S ) + -Φ ψ− 8 ˆAij ˆAij, (61 ) 8
ΔW + 1-∂ (∂ W j) = 8π ψ6j . (62 ) i 3 i j i
We note that in the puncture approach Aˆ ij is determined by solving Equation (56View Equation), not Equation (43View Equation), because Aˆij is not straightforwardly defined for the location with α = 0 when adopting Equation (43View Equation). In this approach, the elliptic equation for βi still has to be solved because βi is needed in solving hydrostatic equations.

Equations (59View Equation) – (62View Equation) are elliptic in type, and hence, appropriate boundary conditions have to be imposed at spatial infinity. Because of the asymptotic flatness, the boundary conditions at spatial infinity are written as

| | ϕ | = βi| = η| = W i| = 0. (63 ) r→ ∞ r→ ∞ r→∞ r→∞
Here it is assumed that the equations are solved in the inertial frame.

In contrast to the case in which the excision approach is adopted, the inner boundary conditions do not have to be imposed in the puncture approach. This could be a drawback in this approach, because one cannot impose physical boundary conditions (e.g., Killing horizon boundary conditions) for the BH. However, this could also be an advantage, because we do not have to impose a special condition for the geometric variables, and as a result, flexibility for adjusting a quasi-equilibrium state to a desired state is preserved.

2.1.3 Hydrostatic equations

The hydrostatic equations governing quasi-equilibrium states are the Euler and continuity equations. The matter in the NS interior needs to satisfy those equations. There are several versions for deriving the hydrostatic equations [26, 9, 192, 215, 80Jump To The Next Citation Point, 204Jump To The Next Citation Point]. In this section, we review the version in [204Jump To The Next Citation Point, 219Jump To The Next Citation Point].

The equation of motion is written as

ν ∇ νT μ = 0. (64 )
Assuming a perfect fluid, Equation (36View Equation), the term in the left-hand side of the equation of motion is written as
∇ νT ν = ρu ν∇ν (huμ) + huμ ∇ν (ρuν) + ∇ μP, μ ν ν = ρ [u ∇ν (huμ) + ∇ μh] + huμ∇ ν (ρu ) − ρT ∇ μs, (65 )
where T is the temperature, and s is the entropy per baryon mass. To derive the second line of Equation (65View Equation), we use an equation of local thermodynamic equilibrium,
1 dh = T ds + -dP. (66 ) ρ
Assuming that the entropy per baryon is constant everywhere inside the NS and using Equation (65View Equation), Equation. (64View Equation) yields
u ν∇ (hu ) + ∇ h = 0, (67 ) ν μ μ
where we use the equation of rest-mass conservation,
ν ∇ ν (ρu ) = 0. (68 )

Equation (67View Equation) can be further modified. Defining a spatial velocity μ v in the comoving frame, the 4-velocity is written as

μ t μ μ u = u (ξ + v ), (69 )
where ξμ is a helical Killing vector, and we have a relation, vμnμ = 0. ut denotes the time component of uμ and is calculated through the normalization condition of u μ, uμu μ = − 1. Substituting Equation (69View Equation) into Equation (67View Equation) and using the relation,
ℒ ξ (hu μ) = 0, (70 )
we obtain
( ) h- ℒv (hu μ) + ∇μ ut = 0. (71 )
If the fluid motion in a NS is synchronized with the orbital motion, i.e., the corotating case, we have vμ = 0. Thus, Equation (71View Equation) can be integrated once to yield
h- = constant. (72 ) ut
If the fluid motion in a NS is irrotational, we need to consider the vorticity of the fluid. The relativistic vorticity tensor is defined by
ω μν ≡ ∇ μ(huν) − ∇ ν (hu μ). (73 )
For the irrotational flow, ω μν = 0, there exists a scalar potential Ψ by which the term hu μ is expressed as
hu μ = ∇ μΨ. (74 )
Inserting Equation (74View Equation) into Equation (71View Equation), we have
( ) ∇ μ ℒv Ψ + h- = 0, (75 ) ut
and thus,
-h ℒv Ψ + ut = CE = constant, (76 )
or, in another form,
μ h- v ∇ μΨ + ut = CE. (77 )

It is interesting to note that if we use the Cartan identity,

μ μ ξ ω μν = ℒξ (hu ν) − ∇ ν (huμξ ), (78 )
another form of the integrated Euler equation, e.g., Equation (33) in [80Jump To The Next Citation Point], is obtained. Because of the helical symmetry ℒξ(hu ν) = 0 and the irrotational flow ω μν = 0, the Cartan identity yields
μ ∇ ν (hu μξ ) = 0, (79 )
and thus,
μ hu μξ = − CE. (80 )
This equation is equivalent to Equation (77View Equation).

The next task is to derive an equation for the velocity potential Ψ. The term in the left-hand side of the equation of continuity (68View Equation) is rewritten as

( ---) ∇ ν (ρu ν) = √-1--ℒu ρ√ − g , − g 1 ( t √ --) = -√---ℒv ρu α γ , α γ 1- ( t i) = α Di ραu v . (81 )
Using Equations (69View Equation), (74View Equation), and the expression for the helical Killing vector,
μ μ μ ξ = αn + β com, (82 )
the equation of continuity (68View Equation) is rewritten as
( ) DiD Ψ + (DiΨ − hutβi ) D ln ρ-α − D (hut βi ) = 0, (83 ) i com i h i com
where Equation (81View Equation) is used.

Finally, we comment on the prescription for determining the constant in the right-hand side of Equation (77View Equation). For this task, the central value of the quantities in the left-hand side is usually used. Here the center of the NS is defined as the location of the maximum baryon rest-mass density. Equation (77View Equation) includes one more constant, which should be determined for each quasi-equilibrium figure; the orbital angular velocity as found from Equations (24View Equation), (69View Equation), (74View Equation), and (82View Equation). The method for calculating it will be explained in the next Section 2.1.4.

2.1.4 Orbital angular velocity

The first integral of the Euler equation (77View Equation) includes information of the orbital angular velocity, Ω, through Equations. (24View Equation), (69View Equation), (74View Equation), and (82View Equation). Ω should be determined from the condition that a binary system is in a quasi-equilibrium circular orbit. In the following, we describe two typical methods referring to the rotation axis of the binary system as the Z-axis and to the axis connecting the BH’s and NS’s centers as the X-axis.

In one of two typical methods, a force balance along the X-axis is required. The force balance equation is derived from the condition that the central value of the enthalpy gradient for the NS is zero,

∂ lnh || ------|| = 0,, (84 ) ∂X 𝒪NS
where 𝒪 NS denotes the NS’s center at which the pressure gradient and self-gravitational force of the NS are absent. Thus, Equation (84View Equation) may be regarded as the condition that the gravitational force from the BH at the NS’s center is equal to the centrifugal force associated with the orbital motion, and hence, can be used to determine Ω for a given set of gravitational field variables.

In the other method, Ω is determined by requiring the enthalpy at two points on the NS’s surface along the X-axis to be equal to h = 1 on the surface. At these two points, the pressure is absent. Namely, the sum of the gravitational force from the BH, self-gravitational force from the NS, and the centrifugal force associated with the orbital motion is balanced. These two conditions may be regarded as the conditions that determine C E and Ω, and thus, Ω can be determined for a given set of gravitational field variables. The work by Taniguchi et al. [208Jump To The Next Citation Point, 209Jump To The Next Citation Point, 210Jump To The Next Citation Point, 214Jump To The Next Citation Point] confirmed that in both methods, an accurate numerical result can be computed with a reasonable number of iterations, and that the results by these two methods coincide within the convergence level of the enthalpy. Therefore, both methods work well.

2.1.5 Center of mass of a binary system

Equation (84View Equation) also depends on the location of the center of mass, because the rotating shift i βrot includes R, which is the radial coordinate measured from the center of mass of the binary system. To determine the location of the center of mass of a binary system, in the framework of the excision approach, Taniguchi et al. [208Jump To The Next Citation Point, 209Jump To The Next Citation Point, 210Jump To The Next Citation Point, 214Jump To The Next Citation Point] and Grandclément [82Jump To The Next Citation Point, 83Jump To The Next Citation Point] require that the linear momentum of the system vanishes

∮ i -1- ij P = 8 π K dSj = 0, (85 ) r→∞
where the maximal slicing condition, K = 0, is assumed. The essence of this condition is as follows: for a hypothetical orbital angular velocity, Ω, the total linear momentum of the system depends on the location of the center of mass, and hence, its location is determined by Condition (85View Equation). Once the location of the center of mass is determined in an iteration step, the positions of the BH and NS are moved, keeping the separation, in order for the center of mass of the binary system to locate on the Z-axis.

In the puncture approach, the situation is totally different from the above, because Condition (85View Equation) has already been used to calculate the linear momentum of the BH; see Equation (58View Equation) in Section 2.1.2. In this framework, there is no known natural, physical condition for determining the center of mass of the system. Until now, three methods have been employed to determine the center of mass. In the first method, the dipole part of ψ at spatial infinity is required to be zero [202Jump To The Next Citation Point, 203Jump To The Next Citation Point]. However, it was found that in this condition, the angular momentum derived for a close orbit of Ωm0 ≥ 0.03 is ∼ 2% smaller than that derived by the 3PN approximation for Q = 3. Because the 3PN approximation should be an excellent approximation of general relativity for a fairly distant orbit, as such Ωm0 ≈ 0.03, the obtained initial data deviates from the true quasi-circular state, and hence, the initial orbit would be eccentric.

In the second method, the azimuthal component of the shift vector βϕ at the location of the puncture is required to be equal to − Ω; a corotating gauge condition at the location of the puncture is imposed [197Jump To The Next Citation Point]. This method gives a slightly better result than that of the first method. However, the angular momentum derived for a close orbit of Ωm0 ≥ 0.03 is also ∼ 2% smaller than that derived by the 3PN relation for a larger mass ratio Q ≥ 2. The disagreement is larger for the larger mass ratio. Such initial conditions are likely to deviate from the true quasi-circular state and hence the orbital eccentricity is large as well.

In the last method, the center of mass is determined in a phenomenological manner: One imposes a condition that the total angular momentum of the binary system for a given value of Ωm0 agrees with that derived by the 3PN approximation [106Jump To The Next Citation Point]. This condition can be achieved by appropriately choosing the position of the center of mass. With this method, the drawback in the previous two methods, i.e., the angular momentum becoming smaller than the expected value, is overcome. Recent numerical simulations by the KT group have been performed employing initial conditions obtained by this method, and showed that the binary orbit is not very eccentric with these initial conditions (cf. Section 3.1.1).

2.1.6 Equations of state

A wide variety of EOS has been adopted for the study of quasi-equilibrium states of BH-NS binaries, which are employed as initial conditions of numerical simulations (see Section 3). However, the only EOS used for the study of quasi-equilibrium sequences has been the polytrope,

P = κρ Γ , (86 )
where κ is a polytropic constant and Γ is the adiabatic index. Hereafter, we review only the numerical results in this EOS.

For the polytropic EOS, we have the following natural units, i.e., polytropic units, to normalize the length, mass, and time scales:

R = κ1∕(2Γ −2). (87 ) poly
Because the geometric units with G = c = 1 are adopted, the polytropic units Rpoly normalize all of the length, mass, and time scales.

Even though the EOS used for constructing sequences is only the polytrope, a lot of quasi-equilibria with several EOS have been derived as initial data for the merger simulations. We will summarize those initial data in Section 3.1.1.

2.1.7 Physical quantities

The several key quantities that are necessary in the quantitative analysis of quasi-equilibrium sequences are summarized in this section; the irreducible mass and spin of BH, the baryon rest mass of NS, the ADM mass, the Komar mass, and the total angular momentum of the system. It is reasonable to consider that the irreducible mass and spin of the BH and the baryon rest mass of the NS are conserved during the inspiral of BH-NS binaries. In addition, temperature and entropy of the NS may be assumed to be approximately zero because the thermal effects of not-young NS are negligible to their structure, i.e., it is reasonable to use a fixed cold EOS throughout the sequence. For such a sequence, the ADM mass, the Komar mass, and the total angular momentum of the system vary with the decrease of the orbital separation. These global quantities characterize the quasi-equilibrium sequence.

We classify the study by seven parameters and summarize in Table 1:

The irreducible mass of the BH
 
This is defined by [42Jump To The Next Citation Point]
∘ ----- AEH-- Mirr ≡ 16π , (88)
where AEH is the proper area of the BH event horizon. Because it is not possible to numerically calculate AEH for a study of quasi-equilibrium configuration, people in the numerical relativity community approximate this area with that of the apparent horizon, AAH, which is computed from an integral on the apparent horizon 𝒮,
∫ 4 2 AAH = ψ d x. (89) 𝒮
In the presence of a helical Killing vector, it is reasonable to consider that AAH agrees with AEH (at least approximately).

In the framework of the excision approach, 𝒮 corresponds to the excision surface. On the other hand, it is necessary to determine the apparent horizon in the framework of the puncture approach (although it is not a difficult task).

The spin of the BH
 
The spin is determined by the quasi-local spin angular momentum of the BH, which is given in Equation (51View Equation), i.e.,
∮ S (ξ) = -1- (Kij − γijK )ξjd2Si. (90) 8π 𝒮
The rest mass of the NS
 
This is defined by
∫ t√ --- 3 MB = ρu − gd x, (91)
where g is the determinant of gμν. Assuming the absence of mass loss, this should be conserved.
The ADM mass of the system
 
The ADM mass in isotropic Cartesian coordinates may be computed from
∮ -1- i MADM = − 2π r→ ∞ ∂ ψdSi. (92)
The Komar mass of the system
 
In the conformal thin-sandwich formalism [231] employed so far, the Komar mass is written as
∮ MKomar = 1-- ∂iαdSi, (93) 4π r→∞
where the shift vector is assumed to fall off sufficiently rapidly.
The total angular momentum of the system
 
York [230Jump To The Next Citation Point] defined this quantity by
∮ --1- ( j kl k jl) Ji = 16 π𝜖ijk X K − X K dSl, (94) r→ ∞
where i X is a spatial Cartesian coordinate relative to the center of mass of the binary system.

In addition, the definition of the linear momentum (which is usually set to be zero) is the same as Equation (85View Equation),

∮ P i = 1-- KijdS , (95 ) 8π r→∞ j
where the maximal slicing condition, K = 0, is assumed.

Then, the binding energy of the binary system is often defined by

Eb = MADM − m0, (96 )
where m0 is the ADM mass of the binary system at infinite orbital separation, m0 ≡ MBH + MNS (see Section 1.8 for the definition). For a non-spinning BH, M BH coincides with the irreducible mass Mirr. For a spinning BH, the relation among MBH, Mirr, and a is given by [42]
[ ] 2 1∕2 MBH = Mirr ------------1∕2- . (97 ) 1 + (1 − a2)

In order to measure a global error in the numerical results, the virial error is often defined as the fractional difference between the ADM and Komar masses,

| | |MADM − MKomar | δM ≡ ||----------------||. (98 ) MADM
Here, we note the presence of a theorem, which states that for the helical symmetric spacetime, the ADM mass and Komar mass are equal (e.g., [76, 204Jump To The Next Citation Point]).

2.1.8 A mass-shedding indicator

To determine the orbit at the onset of mass shedding of a NS, Gourgoulhon et al. [80, 211, 212] defined a “sensitive mass-shedding indicator” (in the context of NS-NS binaries) as

-(∂ (lnh-)∕∂r-)eq χ ≡ (∂ (ln h)∕∂r ) . (99 ) pole
Here, the numerator of Equation (99View Equation), (∂(lnh)∕∂r )eq, is the radial derivative of the enthalpy in the equatorial plane at the surface along the direction toward the companion, and the denominator, (∂(ln h)∕∂r) pole, is that at the surface of the pole. The radial coordinate r is measured from the center of the NS. For a spherical NS at infinite separation, χ = 1, while χ = 0 indicates the formation of a cusp, and hence, the onset of mass shedding. Taniguchi et al. [208Jump To The Next Citation Point, 209Jump To The Next Citation Point, 210Jump To The Next Citation Point] analyze this parameter for identifying the mass-shedding limit of BH-NS binaries.
  Go to previous page Go up Go to next page