Binary Neutron Stars in Quasi-equilibrium

Binary Neutron Stars in Quasi-equilibrium

Abstract

Quasi-equilibrium sequences of binary neutron stars are constructed for a variety of equations of state in general relativity. Einstein’s constraint equations in the Isenberg-Wilson-Mathews approximation are solved together with the relativistic equations of hydrostationary equilibrium under the assumption of irrotational flow. We focus on unequal-mass sequences as well as equal-mass sequences, and compare those results. We investigate the behavior of the binding energy and total angular momentum along a quasi-equilibrium sequence, the endpoint of sequences, and the orbital angular velocity as a function of time, changing the mass ratio, the total mass of the binary system, and the equation of state of a neutron star. It is found that the orbital angular velocity at the mass-shedding limit can be determined by an empirical formula derived from an analytic estimation. We also provide tables for 160 sequences which will be useful as a guideline of numerical simulations for the inspiral and merger performed in the near future.

binaries: close – equation of state – stars: neutron

1 Introduction

Coalescing binary neutron stars are among the most promising sources of gravitational waves for ground-based laser-interferometric gravitational-wave detectors such as LIGO (Brown et al., 2004), GEO600 (Lück et al., 2006), TAMA300 (Ando et al., 2005), and VIRGO (Acernese et al., 2007). Merger of binary neutron stars, together with that of black hole-neutron star binaries, is also considered to be one of the candidates for the central engines of short-hard gamma-ray bursts (Narayan et al., 1992). These facts motivate us to study coalescing binary neutron stars.

Binary neutron stars evolve as a result of gravitational radiation reaction and eventually merge. This evolutionary sequence is usually divided into three stages, depending on the characteristic timescales associated with orbital period and gravitational radiation reaction, as well as on the tidal effects for each neutron star. The first stage is the adiabatic, inspiral phase. In this phase, the timescale of orbital shrink due to the emission of gravitational waves is much longer than the orbital period, and thus, the binary system evolves adiabatically. In addition, each neutron star can be treated as a point mass, because the orbital separation is much larger than the neutron star radius and hence the tidal deformation of the neutron star is negligible. In this phase, a post-Newtonian approximation together with the point particle approximation is a robust tool for determining the orbital evolution and for computing gravitational waveforms (see e.g. Blanchet (2006) and references therein).

The second stage is called the intermediate phase or the quasi-equilibrium phase. In this phase, the binary system is considered to be still in the adiabatic, inspiral phase, but we need to take into account tidal effects on each neutron star, i.e., hydrodynamic effects in neutron stars, as well as full effects of general relativity, because the orbital separation between two neutron stars is only a few times of the neutron-star radius and thus they are in a strong two-body gravitational field. One of the aims of the present paper is to contribute to the understanding of this phase. We will explain more details about the purpose of the present paper later.

The last stage is the merger phase, for which the timescale of orbital shrink becomes shorter than the orbital period and thus the evolution of the system proceeds in a dynamical manner. Furthermore, the system becomes highly general relativistic, because the compactness of the system, defined by the ratio of the gravitational radius to the radius of the system, becomes larger than . To clarify the merger phase, numerical relativity is the unique approach. Since the first fully general relativistic merger simulation was performed by Shibata (1999), huge effort has been devoted in this research field (Shibata & Uryū, 2000, 2002; Shibata et al., 2003, 2005; Shibata & Taniguchi, 2006; Duez et al., 2003; Miller et al., 2004; Anderson et al., 2008a, b; Yamamoto et al., 2008; Liu et al., 2008; Baiotti et al., 2008; Giacomazzo et al., 2009; Kiuchi et al., 2009, 2010). (See e.g. Oechslin & Janka (2007) and Oechslin et al. (2007) for works focusing on the micro-physics in a neutron star but not in fully general relativistic framework.)

Now we return to the intermediate phase and explain the purposes for studying the quasi-equilibrium phase of binary neutron stars in general relativity in detail. There are two roles for this study. One is to clarify the physical conditions in this phase, for example, how large the tidal deformation of a neutron star is, when the mass-shedding from the neutron star occurs, and what the orbital angular velocity at the mass-shedding limit is. The other is to provide initial data for studying the merger phase by numerical relativity simulations. Numerical relativity, in which Einstein’s evolution equations are solved, requires initial data that satisfy Einstein’s constraint equations and also that should be as physical as possible. Obviously, it is important to derive accurate initial data for a scientific study. Constructing such initial data is just obtaining relativistic binary neutron stars in quasi-equilibrium.

The first effort on this issue was devoted to constructing corotating binary neutron stars in general relativity because implementing a numerical code for computing such solutions is relatively easy (Baumgarte et al., 1997, 1998; Marronetti et al., 1998; Usui et al., 2000; Usui & Eriguchi, 2002; Taniguchi & Gourgoulhon, 2002b, 2003; Tichy, 2009). Kochanek (1992) and Bildsten & Cutler (1992), however, found that the timescale of coalescence driven by the gravitational-radiation reaction is much shorter than that of synchronization due to the viscosity in a neutron star. This implies that if the spin angular velocity of neutron stars is much smaller than the orbital angular velocity in a late inspiral phase, we can regard the rotation state of a neutron star to be approximately irrotational for the subsequent phase until the merger sets in. Additionally, any neutron star spins down due to electromagnetic radiation during its life from birth to the coalescence. The spin down timescale of a neutron star in a known binary is at longest as short as the coalescing timescale (larger than yrs; Lorimer, 2008). Moreover, the spin period of neutron stars in a known binary is always larger than 20 ms which is times larger than the orbital period in the late inspiral phase just prior to merger, 2–3 ms. Therefore, we can conclude that the irrotational flow is physically more realistic.2

Under the assumption of irrotation, formulation for solving relativistic hydrostatic equations was derived (Bonazzola et al., 1997; Asada, 1998; Shibata, 1998; Teukolsky, 1998). Soon after the formulation was derived, quasi-equilibrium sequences of irrotational binary neutron stars were calculated (Bonazzola et al., 1999; Marronetti et al., 1999; Uryū & Eriguchi, 2000; Uryū et al., 2000; Gourgoulhon et al., 2001; Taniguchi & Gourgoulhon, 2002b, 2003; Bejger et al., 2005), based on the Isenberg-Wilson-Mathews (IWM) approximation to general relativity (Isenberg, 1978, 2008; Wilson & Mathews, 1989). (See Shibata et al. (2004) for an advanced formulation and Uryū et al. (2006, 2009) for the results.) Even though a lot of sequences have been, so far, calculated, systematic survey has not yet been done. In particular, unequal-mass, irrotational binary neutron stars with an equation of state other than single polytrope has not been studied in detail. In Taniguchi & Gourgoulhon (2002b, 2003), quasi-equilibrium sequences of unequal-mass binaries were calculated, but a polytropic equation of state (EOS) was used. In Bejger et al. (2005) and Uryū et al. (2009), non-polytropic EOSs were used, but quasi-equilibrium sequences of non-equal-mass binaries were not computed. Actually, unequal-mass, irrotational binary neutron stars with realistic EOSs in quasi-circular orbits were constructed and used initial data for merger simulations in Shibata et al. (2005), Shibata & Taniguchi (2006), and Kiuchi et al. (2009, 2010).3 We, however, have not constructed sequences, and rather computed only some initial data sets for each neutron star mass. The purpose of the present paper is to complete the issue and to provide a database of the sequences.

To compute unequal-mass binary systems of arbitrary mass ratio, we develop a new code for the present research, because the numerical code that we developed for the previous works (Taniguchi & Gourgoulhon, 2002b, 2003) had a problem with calculating binary systems composed of significantly different-mass neutron stars in general relativity even though the problem was not in Newtonian computation (Taniguchi & Gourgoulhon, 2002a). As we will explain in Section 2, the method to determine the center of mass of unequal-mass binary systems, i.e., the position of the rotation axis, caused the problem in the previous code, but we have overcome it by employing a new method.

In addition, we employ a wide variety of EOSs; piecewise polytropic EOSs (Read et al., 2009a, b), tabulated realistic EOSs derived from nuclear physics, and fitted EOSs to the tabulated realistic EOSs. Some of the first and second EOSs were, respectively, used in Uryū et al. (2009) and in Bejger et al. (2005), but we adopt a wider set of EOSs in this paper. Furthermore, we systematically study the unequal-mass binaries, whereas Bejger et al. (2005) and Uryū et al. (2009) focused only on the equal-mass case.

This paper is organized as follows. We briefly review the basic equations and explain the improvement of the numerical code in Section 2. In Section 3, the results for the code test are shown. In Section 4, we show numerical results and discuss the effects of EOS on each sequence. Section 5 is devoted to a summary. Throughout this paper we adopt geometrized units with , where denotes the gravitational constant and the speed of light. Latin and Greek indices denote purely spatial and spacetime components, respectively.

2 Formulation

In this section, we briefly describe the basic equations to be solved and the method to calculate a quasi-equilibrium configuration in circular orbits. For more details, we would like to recommend readers to refer to Cook (2000), Baumgarte & Shapiro (2003), and Gourgoulhon (2007) for the part of gravitational field equations, and Gourgoulhon et al. (2001) for that of hydrostatic equations.

2.1 Gravitational field equations

The line element in 3+1 form is written as

 ds2 = gμνdxμdxν, = −α2dt2+γij(dxi+βidt)(dxj+βjdt),

where is the spacetime metric, is the lapse function, is the shift vector, and is the spatial metric induced on a spatial hypersurface . Note here that the direction of the shift vector is the normally used one which has a sign opposite to that used in Gourgoulhon et al. (2001) and Taniguchi & Gourgoulhon (2002b, 2003). The spatial metric is further decomposed into the conformal factor and a background metric , and is written as

 γij=ψ4~γij. (2)

The extrinsic curvature is defined by

 Kij=−12Lnγij (3)

where is the Lie derivative along the unit normal to the hypersurface . We split it into the trace and the traceless part as

 Kij=Aij+13γijK. (4)

We further rewrite the traceless part as

 Aij = ψ−2~Aij, (5a) Aij = ψ−10~Aij, (5b)

where . The Hamiltonian constraint, then, takes the form

 ~∇2ψ=−2πψ5ρH+18ψ~R+112ψ5K2−18ψ−7~Aij~Aij, (6)

where denotes , the covariant derivative with respect to , and the scalar curvature with respect to . The matter term, , in Equation (6) is calculated by taking a projection of the stress-energy tensor. (See Equation (11) below.)

Equations (3), (4), and (5b) yield

 ~Aij=ψ62α(∂t~γij+~∇iβj+~∇jβi−23~γij~∇kβk). (7)

Using Equation (7) for the traceless part of the extrinsic curvature, the momentum constraint is written as

 ~γjk~∇j~∇kβi+13~γik~∇k(~∇jβj) = −αψ6~∇j(ψ6α∂t~γij)+16παψ4ji+43α~∇iK −(~∇iβj+~∇jβi−23~γij~∇kβk)αψ6~∇j(ψ6α),

where the matter term, , is calculated by taking a projection of the stress-energy tensor. (See Equation (12) below.) In addition to the Hamiltonian and momentum constraints, we solve the trace of the evolution equation of the extrinsic curvature

 ∂tK−LβK=−ψ−4(~∇i~∇iα+2~∇ilnψ~∇iα) +α[4π(ρH+S)+ψ−12~Aij~Aij+K23], (9)

where is a matter term defined by Equation (14) below.

As mentioned, the matter terms in the above equations are calculated by taking the projections of the stress-energy tensor into the spatial hypersurface . In the present paper, we assume an ideal fluid and adopt the form of as

 Tμν=(ρ+ρϵ+P)uμuν+Pgμν, (10)

where is the fluid 4-velocity, is the baryon rest-mass density, is the specific internal energy, and is the pressure. Defining the future-directed unit normal to as , the projections of can be written as

 ρH = nμnνTμν, (11) ji = −γiμnνTμν, (12) Sij = γiμγjνTμν, (13) S = γijSij. (14)

The set of equations, (6)–(9), has four functions that we can choose freely; , , , and . For simplicity, we chose a maximal slicing and adopt a flat metric for the spatial background metric in the present paper. We then assume the presence of a helical Killing vector, , and the absence of gravitational waves in the wavezone. Under these assumptions, it is natural to choose the time direction so as to satisfy , and to set and . Then, the basic equations are written as

 Δ––ψ=−2πψ5ρH−18ψ−7~Aij~Aij, (15) Δ––βi+13ηik∂k(∂jβj)=16πΦψ3ji +2~Aij∂j(Φψ−7), (16) Δ––Φ=2πΦψ4(ρH+2S)+78Φψ−8~Aij~Aij, (17) ~Aij=ψ72Φ(ηik∂kβj+ηjk∂kβi−23ηij∂kβk),

where denotes the flat Laplacian, the flat partial derivative, and . The derived formulation is called the IWM approximation (Isenberg, 1978, 2008; Wilson & Mathews, 1989). Note here that a similar but more general formulation was recently proposed and started to be called the extended conformal thin-sandwich decomposition (XCTS; Pfeiffer & York, 2003), mainly in the field of vacuum spacetime, i.e., black hole spacetime. In this new formulation, the conformal flatness is not assumed. (See Gourgoulhon (2007) for a more detailed explanation.)

It may be worthy to note that there is an effort to construct quasi-equilibrium binary neutron stars beyond the IWM approximation. Shibata et al. (2004) proposed a formalism in which all the Einstein equations are solved under the assumption of a helical Killing vector and artificially imposing asymptotic flatness. This formalism was used to solve quasi-equilibrium binary neutron stars in Uryū et al. (2006, 2009).

When solving the set of equations, (15)–(LABEL:eq:taij), we need to impose the outer boundary conditions for , , and . Because our numerical grids include spatial infinity by compactifying the outermost domain, we can impose the exact outer boundary conditions as

 ψ = 1, (19) βi = (\boldmathΩ×\boldmathR)i, (20) Φ = 1, (21)

where is the orbital angular velocity and is the radial coordinate measured from the center of mass of the binary system. Although the shift vector seen by a co-orbiting observer, , diverges at infinity, this diverging term,

 βirot≡(\boldmathΩ×\boldmathR)i, (22)

does not affect the set of equations we solve. If we define the shift vector seen by an inertial observer as , the shift vector seen by a co-orbiting observer can be written as

 βi=βiiner+βirot. (23)

Because we have the relations,

 Δ––βirot=0, (24) ∂jβjrot=0, (25)

and

 ηik∂kβjrot+ηjk∂kβirot=0, (26)

substituting Equation (23) into Equations (16) and (LABEL:eq:taij) yields

 Δ––βiiner+13ηik∂k(∂jβjiner)=16πΦψ3ji +2~Aij∂j(Φψ−7), (27) ~Aij=ψ72Φ(ηik∂kβjiner+ηjk∂kβiiner −23ηij∂kβkiner). (28)

The outer boundary condition for is then at spatial infinity. We actually solve for .

2.2 Hydrostatic equations

The hydrostatic equations governing the quasi-equilibrium state are the Euler and continuity equations. For both irrotational and synchronized motions, the Euler equation can be integrated once to give

 hαγγ0=constant, (29)

where is the fluid specific enthalpy, is the Lorentz factor between the co-orbiting and Eulerian observers, and is the Lorentz factor between the fluid and co-orbiting observers. If we define the 4-velocity of the co-orbiting observer by , the Lorentz factors are written as

 γ0 = −nμvμ=(1−γijUi0Uj0)−1/2, (30) γ = −vμuμ (31) = γ0(1−γijUiUj0)(1−γijUiUj)−1/2,

where is the orbital 3-velocity with respect to the Eulerian observer,

 Ui0=βiα, (32)

and denotes the fluid 3-velocity with respect to the Eulerian observer,

 Ui=ψ−4αuth~∇iΨ, (33)

for irrotational binary systems. Here is the time component of the fluid 4-velocity and is the velocity potential which is calculated by solving the equation of continuity written as

 (34)

where is the covariant derivative with respect to . Note that the fluid 3-velocity corresponds to the orbital 3-velocity for synchronized binary systems.

For the determination of the constant on the right-hand side of Equation (29), we use the central value of the quantities on its left-hand side. The center of a neutron star is defined as the location of the maximum baryon rest-mass density in the present paper. Equation (29) includes one more constant which should be determined for each quasi-equilibrium figure; the constant is the orbital angular velocity as we find from Equations (32), (23), and (22). The method for calculating it will be explained in the next section.

2.3 Orbital angular velocity and the center of mass of a binary system

The method for determining the orbital angular velocity is as follows: we first set the rotation axis of the binary system to be the -axis, and the line connecting the centers of mass of each neutron star to be the -axis. Requiring a force balance along the -axis, we impose a condition of quasi-circular orbit for the binary system. The force balance equation is obtained by setting the central values of the gradient of enthalpy to be zero for each star,

 ∂lnh∂X∣∣Oa=0, (35)

where () denotes the center of each neutron star. Because Equation (29) includes through Equation (22), the force balance Equation (35) may be regarded as the equation for determining the orbital angular velocity. Equation (35) also depends on the location of the center of mass, because Equation (22) includes which is the radial coordinate measured from the center of mass of the binary system. For equal-mass binaries, the force balance equations for each star degenerate and the location of the center of mass becomes trivial, i.e., we can set it to the middle between two stars. For unequal-mass binaries, on the other hand, we have a couple of equations, Equation (35) for and 2, for two parameters of the orbital angular velocity and the location of the center of mass. In the previous papers (Taniguchi & Gourgoulhon, 2002b, 2003), those parameters were determined by solving the couple of equations as stated in Section II B of Taniguchi & Gourgoulhon (2002a). This method works for Newtonian binary systems (Taniguchi & Gourgoulhon, 2002a) and also in the case that the difference in mass of the neutron stars is small for general relativistic binary systems. However, if the difference in mass of the neutron stars is significantly large, the coupled equations, Equation (35) for and 2, would not be simultaneously satisfied at earlier steps of computational iteration because the state of the binary neutron stars is far from equilibrium, and as a result, the computation would fail to achieve the convergence to a solution.

To avoid such crush of computation for small mass ratios, where denotes the Arnowitt-Deser-Misner (ADM) mass for a spherical star in isolation, we adopt the same method as used for black hole-neutron star binaries, described in Taniguchi et al. (2006, 2007, 2008), to determine the location of the center of mass; we require that the linear momentum of the system vanishes

 Pi=18π∮∞KijdSj=0. (36)

Here we have assumed maximal slicing condition, . Once the location of the center of mass is determined in an iteration step, we move the position of each star, keeping the separation, in order for the center of mass of the binary system to locate on the -axis.

For the computation of the orbital angular velocity, we keep the method that Equation (35) is satisfied. As the readers may realize, Equation (35) gives two values of the orbital angular velocity for the unequal-mass case because there are two equations in Equation (35). Even though those values of the orbital angular velocity are very close (the relative difference is within the convergence level), they are slightly different because of numerical error. In the present paper, we just take an average of the two values.

It may be worthy to comment on another method for determining the orbital angular velocity. While our method for calculating is to require the force balance Equation (35), it is possible to obtain by requiring the enthalpy at two points on the neutron star’s surface to be equal, i.e., on the surface. We confirm that the results of those two methods coincide within the convergence level of the enthalpy.

2.4 Global quantities and a mass-shedding indicator

A sequence of binary neutron stars should be constructed for a fixed baryon rest mass of each star,

 M(a)B=∫star aρut√−gd3x,a=1,2, (37)

as the orbital separation decreases. This is because we regard the baryon rest mass as conserved as the orbital separation decreases due to the emission of gravitational waves. Along such a constant-baryon-rest-mass sequence, we then monitor three global quantities: the ADM mass, the Komar mass, and the total angular momentum, as well as a sensitive mass-shedding indicator of a star (see Equation (49)).

The ADM mass in isotropic Cartesian coordinates is written as

If we use Equation (15) and Gauss’ theorem, the ADM mass can be written in terms of volume integral as

Both of Equations (38) and (39) give the same results relative to the convergence level of the computation.

The Komar mass is written as

 MKomar=14π∮∞∂iαdSi, (40)

where we use the fact that the shift vector falls off rapidly enough to be neglected from Equation (40). Using the boundary conditions that and at infinity and the definition , we can rewrite Equation (40) as

 MKomar=14π∮∞(∂iΦ−∂iψ)dSi. (41)

Using Equations (15) and (17), the Komar mass can be also written in terms of volume integral as

 MKomar = 14π∫V[2πψ4(Φ+ψ)ρH+4πΦψ4S (42) +18ψ−7(7Φψ−1+1)~Aij~Aij]dV.

The total angular momentum of the binary system is calculated to give

 Ji=116πϵijk∮∞(XjKkl−XkKjl)dSl, (43)

where is a spatial Cartesian coordinate relative to the center of mass of the binary system. This equation can be also rewritten in the form of volume integral as

 Ji=ϵijk∫Vψ10XjjkdV, (44)

where we use the momentum constraint equation. Similarly, the linear momentum (36) is written as

 Pi=∫Vψ10jidV. (45)

As we mentioned, both of the global quantities calculated by surface integral at infinity and by volume integral give the same results within the convergence level of the computation. In the present paper, we show the results by the volume integral.

The binding energy of the binary system is defined as

where is the ADM mass of the binary system at infinite orbital separation, as defined by the sum of the ADM mass of two isolated neutron stars with the same baryon rest mass,

To measure a global error in the numerical results, we define the error in the virial theorem as the fractional difference between the ADM and Komar masses,

We refer to as the virial error, and use it to measure the magnitude of numerical error in the ADM mass.

Finally, a sensitive mass-shedding indicator is defined as

 χ≡(∂(lnh)/∂r)eq(∂(lnh)/∂r)pole. (49)

Here, the numerator of Equation (49), , is the radial derivative of the enthalpy in equatorial plane at the surface along the direction toward the companion star, and the denominator, , is that at the surface of the pole. The radial coordinate is measured from the center of the corresponding neutron star. For spherical stars at infinite separation, the indicator takes , while indicates the formation of a cusp, and hence the onset of mass-shedding. Note here that because our numerical code is based on a spectral method, it is impossible to construct cusp-like configurations, because it is accompanied by the Gibbs phenomena. This is also the case for the configuration with smaller values of . Thus, we stop constructing a sequence when reaches .

2.5 Equations of state

In this section, we summarize four types of EOSs employed in this paper.

Polytropic EOSs

The first EOS is a polytrope

 P=κρΓ, (50)

where is a polytropic constant and is the adiabatic index. This EOS has been often used in simply modeling binary neutron stars in quasi-equilibrium. We use this EOS only in testing our new code. The specific internal energy for the polytrope is written as

 ϵ=κρΓ−1Γ−1. (51)

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

 Rpoly=κ1/(2Γ−2). (52)

Because geometrized units with are adopted, the polytropic units normalize all of the length, mass, and time scales.

Piecewise polytropic EOSs

The second EOS is a piecewise polytrope introduced by Read et al. (2009a, b). In the present paper, we set the number of polytrope segments to two. Then, the EOS is written as

 P = κ0ρΓ0,(0≤ρ<ρ0) (53) P = κ1ρΓ1,(ρ0≤ρ) (54)

where the dividing density is close to the nuclear density of order (see below). The adiabatic index of the crust is set to a fixed value as , whereas we choose three values for the adiabatic index of the core, , 2.7, and 3.0. The polytropic constant of the crust, , is set to in cgs units. The polytropic constant of the core, , is calculated by requiring the continuity of the pressure at the dividing density as

 κ1=κ0ρΓ0−Γ10. (55)

The dividing density is calculated by setting the fiducial density, , and the pressure at the fiducial density, . We take the fiducial density as where is in cgs units. Because is larger than the dividing density , the EOS has the form of at the fiducial density. Using this equation and Equation (55), the dividing density is obtained as

 log10ρ0=log10P1−14.7×Γ1−log10κ0Γ0−Γ1. (56)

The specific internal energy in the crust and that in the core are, respectively, written as

 ϵ0 = κ0ρΓ0−1Γ0−1 (57) ϵ1 = a1+κ1ρΓ1−1Γ1−1, (58)

where is a constant which is calculated by requiring the continuity of the enthalpy at the dividing density .

We summarize the adiabatic index of the core , the logarithm of the pressure at the fiducial density , the dividing density , and the constant in Table 1; we employ nine parameter sets in the present work. Figure 1 plots the mass-radius relation of spherical stars for the chosen EOSs and indicates that a wide variety of EOSs are modeled.

Tabulated realistic EOSs

We use tabulated EOSs for zero-temperature nuclear matter which are derived by using various theories of dense nuclear matter and different solution methods of the many-body problem in nuclear physics. As described in Bejger et al. (2005), the EOS of Baym et al. (1971) is used for , that of Haensel & Pichon (1994) for , and that obtained by Douchin & Haensel (2001) for in the neutron star crust, where is the density at the crust-core interface. (See Chamel & Haensel (2008) for a review of the EOS in the neutron star crust.)

For the EOS of the neutron star core, we select six EOSs: APR (Akmal et al., 1998), BBB2 (Baldo et al., 1997), BPAL12 (Zuo et al., 1999), FPS (Friedman & Pandharipande, 1981), GNH3 (Glendenning, 1985), and SLy4 (Douchin & Haensel, 2001). Those tabulated EOSs are interpolated by using the Hermite interpolation which is basically the same method as in Swesty (1996). Figure 2 shows the mass-radius relation for spherical stars in these EOSs.

Fitted EOSs to the tabulated realistic EOSs

The method of an interpolation for the tabulated realistic EOSs is not unique. Haensel & Potekhin (2005) introduced a fitting formula using an analytic function. We modified their method (Shibata et al., 2005) and derived a new fitting formula for FPS and SLy4 to satisfy the first law of thermodynamics. In Shibata & Taniguchi (2006), we also derived a fitting formula for APR. In the present paper, we construct quasi-equilibrium sequences for those EOSs, fitAPR, fitFPS, and fitSLy4, and compare the results with those by the tabulated EOSs.

Figure 3 compares the mass-radius relation for spherical stars with those by the tabulated EOSs, APR, FPS, and SLy4. A good agreement is found between two corresponding results for each EOS.

3 Code Tests

We implemented a new numerical code based on the spectral method library, LORENE, developed by the Meudon relativity group. (See the LORENE Web site, http://www.lorene.obspm.fr/, for more detailed explanations of this code library.) Our new numerical code was tested to check its ability for accurately computing unequal-mass binary systems composed of significantly different-mass neutron stars as well as equal-mass ones with similar accuracy to those obtained by the old code used in Gourgoulhon et al. (2001), Taniguchi & Gourgoulhon (2002b, 2003), and Bejger et al. (2005).

3.1 Convergence tests

The first test is to check if the convergence of the ADM mass of an equal-mass binary neutron star is achieved with increasing the resolution. We choose the polytrope and set the baryon rest mass of each star to in polytropic units. Two different orbital separations, and 3.045 in polytropic units, are chosen. Those coordinate separations are, respectively, about 8.965 and 3.736 times larger than the coordinate radius of an isolated, spherical neutron star with the same baryon rest mass. We choose these separations because the former one, , is larger than the farthest separation in the equal-mass sequences we show in the present paper () and the latter one, , is an intermediate separation between the farthest and the closest ones.

Figures 4 and 5 compare the ADM mass at and , respectively, for different resolutions (for different number of collocation points in the terminology of the spectral method). We choose five resolutions, , , , , and , where , , and denote the number of collocation points for the radial, polar, and azimuthal directions, respectively. The horizontal axis of these figures denotes the cube root of the total number of collocation points, . The error bar is drawn for an estimated error size derived from the virial error.

It is found from Figure 4 that the ADM mass for () is in approximately convergent level because its value is approximately identical to that for (). The central value of the ADM mass for () is in the error bar of that for . This implies that the number of collocation points, , is sufficiently large when computing the binary with smaller separations than within the fractional error of for the ADM mass, which satisfies the required accuracy in our present computation.

Figure 5 shows that the ADM mass for is in approximately convergent level because its value is approximately identical to that for at the separation of (). The central value for () is in the error bar of that for . This implies that it is safe to decrease the number of collocation points for smaller separations than within the accuracy required in the present computation.

From these convergence tests, we decide to use the number of collocation points of for larger separations and for closer ones, keeping the number of collocation points for the radial direction.

3.2 An unequal-mass binary system composed of significantly different-mass stars

To illustrate that our new code can compute a binary system composed of significantly different-mass stars, we show in Figure 6 a quasi-equilibrium configuration of binary neutron stars with the polytrope whose baryon rest masses are and 0.15, respectively. The compactness of those stars when they have a spherical shape is, respectively, and 0.1452, where is the circumferential radius. Because the maximum baryon rest mass of the polytrope is , the model of is an extremely light neutron star. We compute this model just to demonstrate the ability of our code. In the figure, the star on the left-hand side has and that on the right-hand side does 0.15. The thick solid circles are the location of the neutron star’s surface. The orbital angular velocity, the binding energy, the total angular momentum, the virial error, and the linear momentum of this figure are , , , , and , respectively. The relative difference of the binding energy from that obtained by the third post-Newtonian (3PN) approximation at the same orbital angular velocity is about . Here we define

 δEb≡(Eb)num(Eb)3PN−1, (59)

where and denote the binding energy obtained by the numerical computation and that by the 3PN approximation, respectively. The relative difference of the total angular momentum from that obtained by the 3PN approximation at the same orbital angular velocity is , defining similar to Equation (59). The relative error in baryon rest mass of the less massive star is and that of the more massive star is . We do not find any problem to accurately construct a binary system composed of two significantly different-mass neutron stars.

3.3 Comparison with the results obtained by the old code

We compare the linear momentum of an unequal-mass binary neutron star along a quasi-equilibrium sequence for a model of piecewise polytrope (PwPoly30-1345) with masses of and in Figure 7. This figure compares the linear momentum of the -component, where we assume that the centers of mass of the neutron stars are located on the -axis (see Figure 6 about the location of each star). This shows that the present results give by more than two orders smaller values than those calculated by the old code through the sequence. This improvement results from the change in the solution method of the center of mass for achieving a convergence. Note that the linear momentum of the - and -directions is zero within the machine precision because of the imposed symmetries.

We also compare the relative difference of the binding energy from that obtained by the 3PN approximation in Figure 8. The definition of the relative difference is shown in Equation (59). The (black) solid curve with the filled circles denotes the results calculated by the new code, and the (red) dashed curve with the open squares is those computed by the old code. The error bar is drawn for an estimated error size derived from the virial error. It is found from Figure 8 that the results by the new code have a factor of 2–5 smaller error bar than those by the old code except for very close separations (). Additionally, when we used the old code, we could not compute sufficiently converged figures for larger separations (), because the code fails to correctly determine the center of mass during the computational iterations. For smaller mass ratio than the model shown in Figure 8, , the old code also fails to achieve the convergence, because the code crushes at earlier steps of computational iteration.

We conclude that our new code gives results as accurate as those obtained by the old code. In addition, it can compute models with smaller mass ratios that the old code cannot.

4 Numerical Results

Quasi-equilibrium sequences for 18 EOSs are computed, choosing three total masses, , , and , and three mass ratios for each total mass. The computation is performed with the collocation points of for larger separations and for closer ones. The number of domains which cover the computational region around each star is six for larger separations and five for closer ones. The results are summarized in Appendix B.

4.1 Contours of baryon rest-mass density

Figures 912 show contours of baryon rest-mass density for 4 EOSs in the equatorial plane. Masses of the stars on the left- and right-hand sides are and , respectively. All the figures are drawn for the closest separation that we can choose. The - and -axes are measured by the coordinate length in km units. Figure 9 is for the tabulated EOS of APR and Figure 10 is for that of GNH3. The former EOS gives relatively compact stars, while the latter one produces rather less compact stars. Figures 11 and 12 are selected among the EOSs of piecewise polytrope, PwPoly30-1345 and PwPoly24-1345. Both of the piecewise polytropes produce stars of a similar size for a given mass, , but the former model, PwPoly30-1345, has a stiff EOS for the core () while the latter one, PwPoly24-1345, has a soft EOS (). This difference in is reflected strongly in the structure of more massive star for which the central density is high and the effect of is appreciated. By contrast, a striking difference is not seen for the less massive star.

4.2 Binding energy and total angular momentum

Figure 13 shows the binding energy of quasi-equilibrium sequences for five piecewise polytropic EOSs. Those sequences are calculated for binary neutron stars composed of equal-mass stars of . All of the piecewise polytropes we select here have , but the value of varies from 13.95 to 13.15. The thick (red) short-dashed, thick (green) long-dashed, thick (blue) dot-dashed, thick (violet) dot-dot-dashed, and thick (magenta) dot-dash-dashed curves denote, respectively, the results for (PwPoly30-1395), 13.55 (PwPoly30-1355), 13.45 (PwPoly30-1345), 13.35 (PwPoly30-1335), and 13.15 (PwPoly30-1315). The thin (black) solid curve denotes the results in the 3PN approximation. The total angular momentum for the same EOSs is also plotted in Figure 14. The sequences are terminated just before the stars reach the mass-shedding limit, because the spectral method we use has a problem in handling a cusp-like figure. We will discuss the endpoint in detail in Section 4.3.

Figure 13 shows that the orbital angular velocity at the closest separation increases from to as the value of decreases from 13.95 to 13.15. Because the model with (PwPoly30-1395) has the largest radius for a spherical neutron star among the models in Figure 13, the effect of tidal deformation is the largest. We understand this behavior with the help of a Newtonian analytic estimation as follows: by equating the gravity of a neutron star (NS1) attracting on a test mass on the neutron star’s surface with the tidal force of the companion neutron star (NS2) attracting on the test mass, we obtain the separation at which the mass-shedding will occur for NS1. The separation is written in the form

 (60)

where , , , and denote a constant of order unity, the NS1’s radius, the NS1’s mass, and the NS2’s mass, respectively. Eliminating the separation at the mass-shedding from the Keplerian angular velocity , we have

 (61)

where is the total mass. Taking the ratio of the orbital angular velocity at the mass-shedding for PwPoly30-1395 to that for PwPoly30-1315, we obtain

 (MtotΩ)1395(MtotΩ)1315=(r1315NS1r1395NS1)3/2. (62)

Noting that the masses are and is the circumferential radius, we have the ratio as 0.51 from Table 2. This value agrees approximately with the obtained ratio of the orbital angular velocity at the closest separation,

 (M0Ω)1395(M0Ω)1315≃0.0290.056≃0.52. (63)

As a result of the smaller orbital angular velocity at the closest separation, the model with (PwPoly30-1395) has the nondimensional angular momentum as large as , as shown in Figure 14. By contrast, the most compact model with (PwPoly30-1315) has , much smaller than 0.96, because the binary system can come closer than that composed of less compact stars. These differences suggest that the merger process will depend strongly on the EOS: for the stiffer EOS, the nondimensional angular momentum at the onset of merger may be too large to form a black hole soon, whereas for the softer EOS, it may be small enough that the merged object collapses to a black hole in a dynamical timescale ms.

Figure 15 is drawn for the binding energy of equal-mass binary neutron stars with . The chosen EOSs are the piecewise polytrope with a fixed value of , but the value of is varied from 3.0 to 2.4. The thick (red) short-dashed, thick (green) long-dashed, and thick (blue) dot-dashed curves denote the results for (PwPoly30-1345), 2.7 (PwPoly27-1345), and 2.4 (PwPoly24-1345), respectively. For these three models, the radius of a spherical star for has almost the same value but very slightly more compact for than for (see Figure 1 and Table 2). It is found that the orbital angular velocity at the closest separation increases from to 0.046 as the adiabatic index of the neutron star core is decreased from to 2.4. This variation of the orbital angular velocity is smaller than that by the change of as seen in Figure 13. This implies that for neutron stars of mass , the variation in , which affects the stellar radius, results in a larger effect on the determination of the orbital angular velocity at the closest separation than the effect of . The reason for this is that the maximum density for the neutron star mass of is not so large that their structure depends weakly on , which determines the stiffness of the core EOS.

Figure 16 compares the binding energy of an equal-mass binary neutron star with that of an unequal-mass one. The EOS we choose here is the piecewise polytrope with and (PwPoly30-1345). The thick (red) dashed curve denotes the result for an equal-mass binary with , while the thick (blue) dot-dashed one is for an unequal-mass binary composed of and stars. The thin (black) solid and thin (green) dotted curves are the results of the 3PN approximation for the equal-mass and unequal-mass binaries, respectively. When the mass ratio decreases to the value as small as , the fractional binding energy and orbital angular velocity at the closest separation decrease by about 10%. The reason why the orbital angular velocity at the closest separation decreases for the unequal-mass case is that the less massive star is tidally deformed by the companion more massive star and starts shedding mass at a larger separation than that for the equal-mass case. The binding energy decreases for the unequal-mass case because it is proportional to the reduced mass according to the results of the 3PN approximation (the ratio of the unequal-mass case to the equal-mass one is ) and the orbital angular velocity at the termination point of the sequence should be smaller.

In Figure 17, the total angular momentum for the same models as in Figure 16 is shown along both equal-mass and unequal-mass sequences. The sequence of the total angular momentum for the unequal-mass case is located below that of the equal-mass case at the same orbital angular velocity. This is also because the total angular momentum is proportional to the reduced mass, according to the results of the 3PN approximation. However, because the orbital angular velocity at the termination point of the sequence is smaller for the unequal-mass case, its total angular momentum is approximately the same value as that for the equal-mass case coincidentally.

Finally, we show the effect of the total mass of binary neutron stars on the binding energy in Figure 18. The EOS we choose for this figure is one of the tabulated realistic EOSs, the APR EOS. Figure 18 shows the results for three total masses, , , and . As the neutron star mass increases, the star becomes more compact and less subject to tidal disruption. For a binary system composed of more massive neutron stars, the two stars are necessary to come closer each other to reach their mass-shedding limit. This makes the orbital angular velocity at the mass-shedding limit increase for binary systems with massive neutron stars. Because of this effect, the orbital angular velocity at the closest separation is about 0.038 for whereas it is about 0.054 for . We can understand this behavior by using Equation (61). By substituting each total mass and circumferential radius, we obtain the ratio of the orbital angular velocity as

 (MtotΩ)2.4M⊙(MtotΩ)3.0M⊙≃(2.4M⊙/11.37 km3.0M⊙/11.31 km)3/2≃0.71. (64)

The actual ratio of the orbital angular velocity at the closest separation we obtained is

 (M0Ω)2.4M⊙(M0Ω)3.0M⊙≃0.0380.054≃0.70. (65)

This value agrees again with the estimated ratio given in Equation (64).

4.3 Endpoint of sequences

We construct quasi-equilibrium sequences for three cases of total mass, three mass ratios for each total mass, and 18 EOSs. For all the sequences, we do not find the turning point of the binding energy and total angular momentum, which represents the innermost stable circular orbit of the binary system. Instead, the sequences terminate at the mass-shedding limit of a less massive star. (For equal-mass binaries, two neutron stars reach the mass-shedding limit at the same time.) As explained before, we stop constructing the sequences at because our spectral method code cannot handle a cusp-like figure which appears at the mass-shedding limit. This implies that the closest separation with the largest value of we can calculate is not the actual endpoint of the sequence.

To determine the orbital angular velocity at the mass-shedding limit, we extrapolate a curve of the mass-shedding indicator as a function of . The procedure is as follows: we first make a fitting polynomial equation of as a function of for sequences. Then, we extrapolate the fitting polynomial equation to and determine the value of at . The extrapolated value of at is defined as . Figure 19 shows such extrapolated curves for the piecewise polytropic EOS with and (PwPoly30-1345). The thick (black) solid, thick (red) dashed, and thick (blue) dot-dashed curves denote the computed sequences for equal-mass binary neutron stars with masses of , , and , respectively. The thin (black) dotted, thin (red) dashed, and thin (blue) dot-dashed curves denote the extrapolated curves.

For unequal-mass binaries, we apply the same method of extrapolation to less massive stars which will be tidally disrupted by their companion more massive stars. Figure 20 shows the sequences and extrapolated curves for unequal-mass binary neutron stars for the piecewise polytropic EOS with and (PwPoly30-1345). The (black) solid, (red) dashed, and (blue) dot-dashed curves denote the computed sequences with masses of versus versus , versus , and versus