A&A 431, 297-306 (2005)
DOI: 10.1051/0004-6361:20041441

Impact of the nuclear equation of state on the last orbits of binary neutron stars

M. Bejger¹ - D. Gondek-Rosinska^1,2 - E. Gourgoulhon² - P. Haensel^1,2 - K. Taniguchi^3, - J. L. Zdunik¹

1 - N. Copernicus Astronomical Center, Polish Academy of Sciences, Bartycka 18, 00-716 Warszawa, Poland
2 - Laboratoire de l'Univers et de ses Théories, UMR 8102 du CNRS, Observatoire de Paris, Université Paris 7, 92195 Meudon Cedex, France
3 - Department of Earth Science and Astronomy, Graduate School of Arts and Sciences, University of Tokyo, Komaba, Meguro, Tokyo 153-8902, Japan

Received 9 June 2004 / Accepted 13 October 2004

Abstract
We present calculations of quasi-equilibrium sequences of irrotational binary neutron stars based on realistic equations of state (EOS) for the whole neutron star interior. Three realistic nuclear EOSs of various softness and based on different microscopic models have been joined with a recent realistic EOS of the crust, giving in this way three different EOSs of a neutron-star interior. Computations of quasi-equilibrium sequences are performed within the Isenberg-Wilson-Mathews approximation to general relativity. For all evolutionary sequences, the innermost stable circular orbit (ISCO) is found to be given by mass-shedding limit (Roche lobe overflow). The EOS dependence on the last orbits is found to be quite important: for two $1.35~M_{\odot}$ neutron stars, the gravitational wave frequency at the ISCO (which marks the end of the inspiral phase) ranges from 800 Hz to 1200 Hz, depending upon the EOS. Detailed comparisons with 3rd order post-Newtonian results for point-mass binaries reveals a very good agreement until hydrodynamical effects (dominated by high-order functions of frequency) become important, which occurs at a frequency ranging from 500 Hz to 1050 Hz, depending upon the EOS.

Key words: dense matter - equation of state - binaries: close - stars: neutron - gravitational waves

   
1 Introduction

Detection of gravitational waves (GW) with the use of ground-based laser interferometers, such as VIRGO (Acernese et al. 2004), LIGO (Abbot et al. 2004), GEO600 (Hewitson et al. 2003), or TAMA300 (Ando et al. 2001) will provide crucial information about various astrophysical objects. Among them, coalescing binary neutron stars (NS) are particularly interesting for the use of GWs to probe their interiors. Matched filtering techniques involving high-order post-Newtonian effects will provide the individual masses and spins of the NSs. The last few orbits before the final merger are dominated by the strong tidal forces acting between the components. The tidal deformation of NSs shape and of the matter distribution in the stellar interior are expected to depend rather strongly on the poorly known equation of state (EOS) of dense matter. The GW signal carries therefore some imprint of the EOS. In particular, the final frequency of the last stable circular orbit is strongly correlated with the compactness parameter M/R (Faber et al. 2002; Taniguchi & Gourgoulhon 2003), and thus can provide direct constraints on the theory of dense hadronic matter above the nuclear saturation density. A related event potentially rich in information about the dense matter EOS is the merger of a neutron star and a black hole (Prakash & Lattimer 2004; Prakash et al. 2004). In addition to locating the transition from inspiral to the merger phase, computations of the last stable orbits of binary NSs are also required for providing initial data to compute the dynamical merger phase (Shibata et al. 2003, and references therein).

The last orbits of inspiraling binary neutron stars have been studied by a number of authors in the quasi-equilibrium approximation, and in the framework of Isenberg-Wilson-Mathews (IWM) approximation of general relativity (see Baumgarte & Shapiro 2003, for a review). The quasi-equilibrium assumption approximates the evolution of the system by a sequence of exactly circular orbits (as the time evolution of the orbit is much larger than the orbital period). The IWM approximation amounts to using a conformally flat spatial metric (the full spacetime metric remaining non-conformally flat), which reduces the problem to solving only five of the ten Einstein equations. Within these two approximations, the most realistic studies have considered irrotational binaries, as opposed to synchronized ones, as the viscosity of neutron star matter is far too low to ensure synchronization during the late stage of the inspiral (Bildsten & Cutler 1992; Kochaneck 1992).

The quasi-totality of the existing quasi-equilibrium IWM studies (Bonazzola et al. 1999; Taniguchi & Gourgoulhon 2002b; Marronetti et al. 1999; Gourgoulhon et al. 2001; Uryu et al. 2000; Taniguchi & Gourgoulhon 2003; Uryu & Eriguchi 2000) employ a polytropic EOS to model the neutron star interior. The only exception is the recent work of Oechslin et al. (2004) who have used two EOSs: (i) a pure nuclear matter EOS, based on a relativistic mean field model; and (ii) a "hybrid'' EOS with a phase transition to quark matter at high density. At 2 $\times$ $10^{14}~{\rm g~ cm}^{-3}$ (i.e. in the vicinity of nuclear density), these two EOSs are matched with a polytropic one with adiabatic index $\gamma=2.86$. This last assumption of Oechslin et al. is somewhat ad hoc, because the EOS of the neutron star crust is very different from a polytrope, and its local adiabatic index is much smaller. Indeed, the crust polytropic index within the inner crust (which contains some 99.9% of the total crust mass) varies from $\gamma\simeq 0.5$ near the neutron drip point to $\gamma\simeq 1.6$ in the bottom layers near the crust-core interface (Douchin & Haensel 2001). The crustal EOS plays a crucial role in defining the mass-shedding limit which marks the end point of quasi-equilibrium binary configurations.

In the present article, we study the last orbits of irrotational binary neutron star systems by using a set of three dense matter EOSs which are representative of the contemporary many-body theory of dense matter. Contrary to Oechslin et al. (2004), we describe the neutron star crust by means of a realistic EOS obtained in the many-body calculations. As in the works mentioned above, we use the quasi-equilibrium and IWM approximations.

The paper is arranged in the following way: in Sect. 2 the differences between EOSs used in the computation are briefly described. Section 3 is devoted to a description of the numerical methods used to obtain the quasi-equilibrium orbital sequences. In Sect. 4 the results are presented, Sect. 5 contains discussion and final remarks.

   
2 Description of the equations of state

   
2.1 The EOS of the crust

The outer layer of a neutron star contains neutron-rich nuclei, which due to Coulomb repulsion form a crystal lattice if the matter temperature is below the corresponding melting temperature. This layer is called the neutron star crust, and extends up to the density at the crust-core interface, $\rho_{\rm cc} \sim 10^{14}~{\rm g~cm^{-3}}$. The precise value of  $\rho_{\rm cc}$ is model dependent, and varies within (0.6-1.4) $\times$ $10^{14}~{\rm g~ cm}^{-3}$. The EOS of the crust is rather well established (for a review, see, e.g., Haensel 2003). In the present paper the EOS of the crust is composed of three segments. For densities smaller than  $10^{8}~{\rm g~cm^{-3}}$ we used the EOS of Baym et al. (1971). For $10^8~{\rm g~cm^{-3}}<\rho<10^{11}~{\rm g~cm^{-3}}$ we applied the EOS of Haensel & Pichon (1994), which makes maximal use of the experimental masses of neutron-rich nuclei. Finally, for densities $10^{11}~{\rm g~cm^{-3}}<\rho<\rho_{\rm cc}$we used the EOS of Douchin & Haensel (2001). For neutron stars of $M=1.35~{M}_\odot$, the crust contains less than 2% of the stellar mass. However, it is the region most easily deformed under the action of the tidal forces resulting from the gravitational field produced by the companion star. Below the melting temperature, elastic shear terms in the stress tensor are nonzero, but they are two orders of magnitude smaller than the main diagonal pressure term (Haensel 2001). Dissipative and thermal effects accompanying matter flow inside neutron stars will be briefly discussed in Sect. 5; they are expected to be small at the quasi-equilibrium evolution stage. Therefore, to a very good approximation, the crust layer behaves in the tidal force field as an ideal degenerate fluid, described by a zero temperature EOS.

An important quantity which actually determines the response of the crust layer to the compression or decompression of matter is the adiabatic index $\gamma=(n/p){\rm d}p/{\rm d}n$, where p is the local pressure and n the corresponding baryon number density. In the outer crust, the pressure is determined by the ultra-relativistic electron gas, so that $\gamma=4/3$. However, the outer crust contains only 10^-5 of the NS mass. Some thousand times more, i.e., about  $0.03~M_\odot$, is contained in the inner crust, composed of a lattice of heavy neutron-rich nuclei immersed in a neutron gas and an electron gas. Its adiabatic index depends strongly on the density, varying from $\gamma\simeq 0.5$ near the neutron drip point (extreme softening of the EOS), up to $\gamma\simeq 1.6$ close to  $\rho_{\rm cc}$ (Douchin & Haensel 2001).

   
2.2 The EOS of the core

The neutron star core contains matter of a density significantly larger than the normal nuclear density, equal to the saturation density of infinite symmetric nuclear matter, $\rho_0=2.8$ $\times$ $10^{14}~{\rm g~ cm}^{-3}$, and corresponding to the baryon number density $n_0=0.16~{\rm fm^{-3}}$. For $\rho>\rho_0$ the EOS of the core is poorly known, and this uncertainty grows rapidly with increasing density. The theoretical EOSs, derived using different theories of dense matter and different methods of solution of the many-body problem, differ significantly at  $10^{15}~{\rm g~cm^{-3}}$, characteristic of the central cores of neutron stars under study. In the present paper, we will limit ourselves to neutron star cores consisting of nucleons and hyperons, i.e., baryons whose properties are known from terrestrial experiments. The speculative case of exotic phases of dense matter (kaon and pion condensation, quark deconfinement), whose existence is - as of this publication - not substantiated either by experiments or by astronomical observations, will be considered in future publications. The most extreme case of hypothetical strange stars built of self-bound strange quark matter will also be studied separately.

Three EOSs of the core were considered. Two of them may be considered as soft and stiff extremes of the EOSs of matter composed of nucleons, electrons and muons. The first, BPAL12, is of phenomenological type and can be considered as a soft extreme of the nucleonic EOS of NS matter (Prakash et al. 1997; Bombaci 1995). The second, APR, is based on precise variational calculations and includes realistic two-body (Argonne A18) and three-body (Urbana UIX) nucleon interactions (Akmal et al. 1998). We considered also one EOS in which hyperons are present at $\rho>\rho_{\rm H}$, where the threshold for the hyperon appearance $\rho_{\rm H}\simeq 2\rho_0$ (model 3 of Glendenning 1985); it will be referred to as the GNH3 EOS. This EOS was obtained using the Relativistic Mean Field model of baryonic matter. For $\rho<\rho_{\rm H}$ (nucleons only) this EOS is very stiff but causal ( $v_{\rm sound}<c$). The appearance of hyperons strongly softens the EOS as compared to the pure nucleon case. The hyperons soften the EOS because they are more massive than nucleons, and when they start to fill their Fermi sea they are slow and replace the highest energy nucleons.

	Figure 1: The pressure p against the energy density $\rho$ for the EOSs used in the paper: APR (dashed line), GNH3 (dotted line), and BPAL12 (solid line). The empty circles correspond to the central density of a non-rotating stellar model with a gravitational mass equal to  $1.35~M_{\odot}$ (Table 1).
Open with DEXTER

All three EOSs are displayed in Fig. 1. As we mentioned they are very different because they assume very different strong interaction models at supra-nuclear densities.

   
2.3 Mass-radius relation for static NSs

The differences between the three EOSs are reflected in the M(R) curves of statically stable models shown in Fig. 2 (M being the gravitational mass and R the area radius).

	Figure 2: Gravitational mass of static isolated stars against area radius for the APR (dashed line), GNH3 (dotted line) and BPAL12 (solid line) EOSs.
Open with DEXTER

The BPAL12 EOS is overall quite soft, and therefore produces NSs with R decreasing rapidly with increasing M. The BPAL12 EOS yields the smallest R(M) for $M>1.2~M_\odot$. The maximum allowable mass for this EOS is marginally consistent with observations, $M_{\rm max}=1.46~{M}_\odot$. The APR EOS is stiff, with its radius staying nearly constant at $R\simeq 11.5~$km for $M=(0.6{-}1.6)~{M}_\odot$, and yields high $M_{\rm max} =2.2~{M}_\odot$. The BPAL12 and APR EOSs are soft and stiff extremes within our set of the EOSs as far as the values of  $M_{\rm max}$ are concerned. However, their M(R) intersect at $M=1.2~M_\odot$. Therefore, these EOSs produce NSs with similar values of the compaction parameter M/R for $M\simeq 1.2~M_\odot$.

The GNH3 EOS is different to the two previous ones, because its M(R) curve is composed of two different segments. The lower-mass segment ( $M<1.3~{M}_\odot$) consists of stars with no or with only small hyperon cores. The radius stays nearly constant at $R\simeq 14$ km for $M=(0.5{-}1.3)~{M}_\odot$. This segment is connected via a "knee'' with a high-mass segment consisting of neutron stars with increasingly larger soft hyperon cores. Along this high-mass (hyperon-softened) segment, R decreases rapidly with increasing M, reaching a very flat maximum at $M\simeq 1.9~{M}_\odot$.

Table 1: Properties of isolated static neutron stars of gravitational mass $M=1.35~{M}_\odot$ for the three EOSs used in our computations. $M/R\equiv GM/Rc^2$ is the compaction parameter, R is the areal radius, $M_{\rm B}$ is the baryon mass and  $\rho _{\rm c}$ is the central energy density, respectively.

Different aspects of the EOSs show up if we consider the $1.35~M_{\odot}$ static stars. Their parameters, calculated using three EOSs, are given in Table 1. A NS configuration is a functional of the EOS at the densities $\rho$ smaller than the central density  $\rho _{\rm c}$. We see that for the GNH3 EOS the central density is barely higher than the hyperon threshold  $\rho_{\rm H}$, and therefore the EOS inside the NS is very stiff, actually the stiffest of all three. The highest stiffness of the GNH3 EOS in the interior of a  $1.35~M_{\odot}$ NS can be clearly recognized by looking at Figs. 1 and 3. Looking at Fig. 1, we see that this EOS has the highest P at any $\rho<\rho_{\rm c}$. Additional information is contained in Fig. 3: this EOS has particularly large $\gamma$ for $\rho\sim \rho_0$, i.e., in the outer layers of the NS core, which are therefore quite "inflated'' in comparison with those in the other two NS models. These features are responsible for the particularly large R for the GNH3 EOS. The GNH3 configuration has the largest R and therefore the smallest compaction parameter M/R. The BPAL12 EOS is clearly the softest, with R smaller by 27% than for the GNH3 EOS. The APR EOS, which is moderately stiff in this mass range, yields R which is only 8% larger than for the BPAL12 model.

The differences in stiffness reflect the characteristics of the nuclear model underlying each of the EOSs of the NS core. Using the density-dependent adiabatic index $\gamma$ it is particularly easy to visualize these differences, see Fig. 3. The strong drop in $\gamma$ above  $\rho_{\rm H}$ reflects the hyperon softening in the GNH3 EOS. The values of $\gamma$ increasing up to the maximum at $\rho=\rho_{\rm c}$ tell us about the stiffening of the APR EOS at $\rho\sim 2\rho_0$, in contrast to the behavior of the GNH3 EOS which softens close to this density. Finally, the BPAL12 EOS remains close to a polytrope with $\gamma\simeq 2.2$, except for a small region around $\rho_0$.

	Figure 3: The adiabatic index $\gamma$ against the energy density $\rho$ for the EOSs used in our computations: APR (dashed line), GNH3 (dotted line), and BPAL12 (solid line). The empty circles correspond to the central density of a non-rotating stellar model with a gravitational mass equal to  $1.35~M_{\odot}$.
Open with DEXTER

   
3 Numerical method

3.1 Numerical code for close binary systems

The present computations of close binary neutron star systems rely on the assumption of a quasi-equilibrium state (helical Killing vector approximation), with irrotational flow of the fluid, and a conformally flat spatial 3-m (Isenberg-Wilson-Mathews approximation). We construct the quasi-equilibrium sequences of binary neutron stars described by the realistic EOSs using a numerical code which solves the five coupled, nonlinear, elliptic equations for the gravitational field, supplemented by an elliptic equation for the velocity potential of irrotational flows. The code has been used successfully for calculating the final phase of the inspiral of binary neutron stars described by the polytropic equation of state (Taniguchi & Gourgoulhon 2002b,a,2003; Taniguchi et al. 2001). This code is built upon the C++ library L ORENE (http://www.lorene.obspm.fr) and can be downloaded freely from the L ORENE CVS repository, as Lorene/Codes/Bin_star/coal.C. The complete description of the resulting general relativistic equations, the whole algorithm, as well numerous tests of the code can be found in Gourgoulhon et al. (2001). Additional tests have been presented in Sect. 3 of Taniguchi & Gourgoulhon (2003).

The numerical technique relies on a multi-domain spectral method with surface-fitted coordinates. We have used one domain for each star and 3 (resp. 4) domains for the space around them for a large (resp. small) separation. In each domain, the number of collocation points of the spectral method is chosen to be N_r $\times$ $N_{\theta}$ $\times$ $N_{\varphi}$ = 25 $\times$ 17 $\times$ 16, where N_r, $N_{\theta}$, and  $N_{\varphi}$ denote the number of points in the radial, polar and azimuthal directions respectively. The accuracy of the computed relativistic models was estimated using a relativistic generalization of the the Virial Theorem (Friedman et al. 2002; see also Sect. 3.A of Taniguchi & Gourgoulhon 2003). The virial relative error was a few times 10^-5.

3.2 The velocity potential of irrotational flows

Let us briefly discuss a technical difference between the computations of the realistic EOS irrotational binary NS and the polytropic case considered by Bonazzola et al. (1999), Marronetti et al. (1999), Uryu & Eriguchi (2000), Uryu et al. (2000), Gourgoulhon et al. (2001) and Taniguchi & Gourgoulhon (2002b,2003). The difference stems from the fact that the realistic EOS presented in Sect. 2 are given in tabulated form, and a certain thermodynamic coefficient $\zeta$, (required in the computation of the velocity potential of the irrotational fluid flow, cf. Sect. 2 of Gourgoulhon et al. 2001) is not given explicitly by the tabulated EOS.

As already discussed in Sect. 2, we adopt the approximation of the perfect fluid for the form of the stress-energy tensor, and we represent the NS matter by a zero-temperature EOS. In view of this, we can use the Gibbs-Duhem identity:

$$% \frac{{\rm d} p}{e+p}=\frac{1}{h}{\rm d} h$$ (1)

where $e=\rho c^2$ is the proper energy density of the fluid, p is the pressure, $h=(e+p)/(m_{\rm b}n)$ is the specific enthalpy, with ndenoting the fluid baryon number density and $m_{\rm b}$ the mean baryon mass.

Following Gourgoulhon et al. (2001), we write the equation for the coefficient $\zeta$ as

$$% \zeta = \frac{{\rm d(ln}~H)}{{\rm d(ln}~n)} = \frac{{\rm d(... ...rm d(ln}~n)} = \frac{{\rm d(ln}~H)}{{\rm d(ln}~p)}\cdot\gamma,$$ (2)

where $\gamma$ is the adiabatic index. The quantity $H:=\ln h$appears in the first integral of fluid motion and denotes the log-enthalpy of the fluid, and the value of ${\rm d(ln}H)/{\rm d(ln}p)$ can easily be evaluated by means of Eq. (1) which, as a strict thermodynamic relation, is exact.

In practice, we used two methods to get the adiabatic index $\gamma=(n/p){\rm d}p/{\rm d}n$. In the first method we used an analytical formula for the adiabatic index being an approximation of tabulated values of $\gamma$. In the second method the adiabatic index was obtained directly from a tabulated EOS by taking the derivative of the second order polynomial which goes through three consecutive (p, n) EOS points. The value of the ${\rm d}p/{\rm d}n$ multiplied by the n/p ratio is evaluated at the middle point, and the resulting discrete values of $\gamma$ are ready to be interpolated, similar to the other quantities from the tabulated EOS.

The second method proved to be very robust, and acts as a consistency check; it was tested with EOS for which the precise values of the adiabatic index were obtained from the microscopic considerations, namely for the SLy EOS (Douchin & Haensel 2000).

   
4 Numerical results

   
4.1 Quasi-equilibrium sequences of neutron stars described by realistic EOS

In this section we present the numerical results for evolutionary sequences of close neutron star binaries described by the three realistic EOS introduced in Sect. 2. By evolutionary sequence, we mean a sequence of quasi-equilibrium configurations of decreasing separation d and with constant baryon mass $M_{\rm B}$. Such a sequence is expected to approximate the true evolution of binary neutron stars, which is entirely driven by the reaction to gravitational radiation and hence occur at fixed baryon number and fluid circulation (zero in the irrotational case considered here). We calculated evolutionary sequences of binary systems composed of two identical neutron stars (equal mass system), with the baryon mass $M_{\rm B}$corresponding to the gravitational mass $M:=M_1=M_2=1.35~M_{\odot}$for a static isolated star (see Table 1 for the values of $M_{\rm B}$). We have selected the gravitational mass $M=1.35~{M}_\odot$ for two reasons: (i) it agrees with the "average NS mass'' obtained from observations of radio pulsars in binary system; and (ii) it allows us to compare our results with calculations made by other authors for polytropic models (Faber et al. 2002).

For the present discussion, let us recall the definition of the ADM mass of the system, which measures the total energy content of a slice $t={\rm const}$of spacetime (hypersurface $\Sigma_t$):

$$% M_{\rm ADM} := \frac{1}{16\pi} \oint_\infty \left[ {\cal D... ...al D}_i \left( f^{kl} \gamma_{kl} \right) \right] {\rm d}S^i,$$ (3)

where $\gamma_{ij}$ is the metric induced by the spacetime metric on $\Sigma_t$, f_ij is a flat metric on $\Sigma_t$ (the condition of asymptotic flatness being $\gamma_{ij} \rightarrow f_{ij}$), ${\cal D}_i$ is the covariant derivative associated with f_ij, and the integral is to be taken on a sphere at spatial infinity. In the case considered here of a conformally flat spatial metric, $\gamma_{ij} = \Psi^4 f_{ij}$, the surface integral (3) can be converted into a volume integral over the whole hypersurface $\Sigma_t$:

$$% M_{\rm ADM} := \int_{\Sigma_t} \Psi^5 \left( E +\frac{1}{16\pi} K_{ij} K^{ij} \right) {\rm d}^3 x,$$ (4)

where E is the total energy density of matter with respect to the observer whose 4-velocity is normal to $\Sigma_t$ and K_ij denotes the extrinsic curvature tensor of $\Sigma_t$.

	Figure 4: Baryon number density isocontours in the coordinate plane y=0 (containing the rotation axis) ( left panels) and in the coordinate plane z=0 (orbital plane) ( right panels), for configurations close to the mass-shedding limit. The upper (resp. middle, lower) panels correspond to the GNH3 (resp. APR, BPAL12) EOS. The thick solid lines denote stellar surfaces.
Open with DEXTER

At infinite separation, the ADM mass of the system, Eqs. (3) and (4), tends toward the sum of the gravitational masses of isolated static stars, and will be denoted by $M_\infty$:

$$% \lim_{d\rightarrow\infty} M_{\rm ADM} = M_\infty := M_1+M_2 = 2.7~M_{\odot}.$$ (5)

We then define the orbital binding energy by

$$% E_{\rm bind} := M_{\rm ADM} - M_\infty.$$ (6)

The variation of $E_{\rm bind}$ along an evolutionary sequence corresponds to the loss of energy via gravitational radiation. Gravitational waves are emitted mostly at twice the orbital frequency: $f_{\rm GW}=2f$.

4.2 Innermost stable circular orbit

Each evolutionary sequence terminates by a mass-shedding point, which marks the end of existence of quasi-equilibrium configurations. The shape of the stars close to this limit is presented in Fig. 4. The mass-shedding is revealed by the formation of a cusp at the stellar surface in the direction of the companion (Roche lobe overflow). This cusp is marginally visible in Fig. 4.

A turning point of $E_{\rm bind}$ along an evolutionary sequence would indicate an orbital instability (Friedman et al. 2002). This instability originates both from relativistic effects (the well-known r=6M last stable orbit of Schwarzschild metric) and hydrodynamical effects (for instance, such an instability exists for sufficiently stiff EOS in the Newtonian regime, see e.g. Taniguchi et al. 2001 and references therein). It is secular for synchronized systems and dynamical for irrotational ones, such as those considered here. Thus the quasi-equilibrium inspiral of binary neutron stars can terminate by either the orbital instability (turning point of  $E_{\rm bind}$) or the mass-shedding limit. In the latter case, there might exist equilibrium configurations beyond the mass-shedding limit, i.e. dumb-bell like configurations (see e.g. Eriguchi & Hachisu 1985). However dynamical calculations for the $\gamma =2$ polytrope have shown that the time to coalescence was shorter than one orbital period for configurations at the mass-shedding limit (Marronetti et al. 2004; Shibata & Uryu 2001). Therefore we may safely define the end of quasi-equilibrium inspiral by the mass-shedding limit in the case where no turning point of  $E_{\rm bind}$ is found along the sequences (which is actually the present case, as we shall discuss below).

The variation of the orbital binding energy along evolutionary sequences is presented in Fig. 5, where the points correspond to the equilibrium configurations of binary system calculated using a numerical method (see Sect. 3) and the lines present our best fit described in detail in Sect. 4.3. In the scale of Fig. 5 there is no visible difference between our numerical results and the fit, i.e. fitting curves pass through the points. We plotted also the binding energy obtained in the framework of Newtonian theory and the 3PN post-Newtonian approximation for point masses (Blanchet 2002). Figure 5 shows clearly that no turning point of  $E_{\rm bind}$ occurs along the evolutionary sequences. Hence there is no orbital instability prior to the mass-shedding limit. We conclude that the innermost stable circular orbit (ISCO) of the computed configurations are given by the mass-shedding limit.

Figure 5: Orbital binding energy $E_{\rm bind} = M_{\rm ADM} - M_\infty$ of the binary system versus frequency of gravitational waves (twice the orbital frequency) along three irrotational quasi-equilibrium sequences. The lines (dotted - GNH3, dashed - APR, solid - BPAL12) were plotted using the fit described in the text, whereas the points represent actual data. Thin solid line touching the bottom-right corner shows the 3rd post-Newtonian approximation for point masses derived by Blanchet (2002). The lower thin dotted curve corresponds to the Newtonian limit for point masses.

Open with DEXTER

Figure 6: Orbital binding energy of the binary system minus the (point mass) Newtonian term $-k_{\rm N} f^{2/3}$ versus frequency of gravitational waves (twice the orbital frequency) along three irrotational quasi-equilibrium sequences. Thin solid line shows the 3rd post-Newtonian approximation for point masses by Blanchet (2002). Slightly below there is the 2PN line (thin, dashed) and the 1PN one (thin, long dash-dotted). For the comparison post-Newtonian approximation by Damour et al. (2000) - 3PN (dashed-dotted thin line) and 2PN (thin, long dashed) are also plotted.

Open with DEXTER

To compare different approaches to the relativistic description of binary systems we present in Fig. 6 the orbital binding energy after subtraction of the Newtonian term $\propto$ $f_{\rm GW}^{3/2}$. This figure shows the effects of general relativity and of finite sizes (hydrodynamics). One can see that binary neutron star systems are quite far from Newtonian and 1PN systems. On the contrary, 2PN and 3PN results by Blanchet (2002) and 3PN results by Damour et al. (2000) are very close to our results for a wide range of frequencies. Note that the difference between 2PN and 3PN approximations are much larger in the Damour et al. (2000) treatment (Effective One Body method) than in that of Blanchet (2002) (standard post-Newtonian expansion).

   
4.3 Analytical fits to the numerical results

Following Faber et al. (2002), we have performed some polynomial fits to each of the computed sequences. This is a very important step in order to obtain the first derivative of the functions required for the energy spectrum of gravitational waves (see below). We used two different approaches to approximate our numerical results. The first one, similar to that presented by Faber et al. (2002), is based on quadratic approximation of numerical results. We decided to make a fit only to the difference between exact results obtained for the ADM mass  $M_{\rm ADM}$ (cf. Eq. (3)) and the prediction of the Newtonian theory - i.e. we made a fit to the function  $E_{\rm bind} + k_{\rm N}f_{\rm GW}^{2/3}$, where the second term corresponds to the Newtonian point-like mass behaviour with $k_{\rm N}=(G\pi/4)^{2/3}M^{5/3} = 4.06$ $\times$ $10^{-4}~M_{\odot}{\rm Hz}^{-2/3}$. We found it sufficient to fit the numerical results by a second order polynomial without any linear term:

$$% E_{\rm bind} = -k_{\rm N}f_{\rm GW}^{2/3}+k_{\rm 2}f_{\rm GW}^2,$$ (7)

i.e. contrary to Faber et al. (2002) we assume k₁=0. The best-fit coefficients $k_{\rm 2}$ are collected in Table 2. Our aim was to obtain accurate formula for gravitational radiation in the region where it is effectively emitted and therefore we have performed the fitting procedure for frequencies higher than 500-600 Hz. This approximation works quite well for high frequencies (i.e. small distance between stars) and there is no need to introduce the linear term, as done by Faber et al. (2002). The advantage of our quadratic formula is its simplicity and good accuracy in the region corresponding to the effective emission of gravitational waves.

However it is possible to find a much better approximation of the numerical results if one takes into account a high order post-Newtonian approximation for the binding energy of point-mass systems. The 3PN formula as obtained by Blanchet (2002) from the standard post-Newtonian expansion reads

 

$$% \frac{E_{\rm bind}^{\rm 3PN}}{M_\infty} = -{1\over 8} ~ {\Omega... ...} + {5\over3072} \left({41} \pi^2 -{285473\over 864}\right) ~ {\Omega_*}^{8/3},$$ (8)

where $\Omega_*$ is the orbital angular frequency expressed in geometrized units:

$$\begin{eqnarray*}\Omega_* :=2\pi {GM_{\infty}\over c^3} f = 2\pi {GM\over c^3} f_{\rm GW}. \end{eqnarray*}$$

The terms in ${\Omega_*}^{2/3}$, ${\Omega_*}^{4/3}$, ${\Omega_*}^{2}$ and  ${\Omega_*}^{8/3}$ in Eq. (8) are respectively the Newtonian, 1PN, 2PN and 3PN term.

Table 2: Parameters of polynomial fits (7) and (9).

	Figure 7: Difference $E_{\rm bind} - E_{\rm bind}^{\rm 3PN}$ between the binding energy of irrotational binary neutron stars built upon a realistic equation of state and the binding energy of binary point masses in the 3PN approximation of Blanchet (2002). The dots correspond to numerical results and the lines to polynomial fits to them (see text for details).
Open with DEXTER

In Fig. 7 we present the difference between our numerical results and the 3PN approximation given by Eq. (8). Formula (8) approximates very well the behavior of a binary system of realistic neutron stars for a very large range of binary periods (notice the scale of the y-axis of Fig. 7!). From Fig. 7, we can define the frequencies  $f_{\rm npm}$ as those frequencies at which the deviation from point-mass behavior becomes important. The values of these frequencies for each of the three considered EOSs are given in Table 3. Of course, they are not precise numbers and should be treated as having a $\sim$ $\pm 50~{\rm Hz}$ uncertainty. We can approximate the results presented in Fig. 7 by the power-law dependence on frequency  $f_{\rm GW}$:

$$% E_{\rm bind} - E_{\rm bind}^{\rm 3PN} = a~\left({f_{\rm GW}\over 1000~{\rm Hz}}\right)^n\cdot$$ (9)

Because of the steep character of the function $E_{\rm bind} - E_{\rm bind}^{\rm 3PN}$seen in Fig. 7, the power n is quite large. The values are listed in Table 2. We assume the integer number of the power n, although this is in principle not obvious from a theoretical point of view. Of course more careful treatment of this approximation is possible (e.g. more terms in the expansion) but from Fig. 7 it is clear that there is no need to do so and the leading term of high power is sufficient.

One can draw an important conclusion from the presented results and their comparison with relativistic approximations for point masses in a binary system. We can expect that taking into account the next orders in a post-Newtonian approximation does not change the energy by an amount larger than the difference between 2PN and 3PN models. As a consequence the large deviation of our numerical results from the 3PN approximation is caused by the effect of a finite size of the star (e.g. tidal forces). The very high power n in relation (9) indicates that, even for small departures from a point mass approximation, high-order tidal effects are very important, and dominate the relation  $E_{\rm bind} (f_{\rm GW})$. Indeed the lowest order tidal term is known to be n=4 (Lai et al. 1994) and the values obtained here are well above this.

4.4 Energy spectrum of gravitational waves for different realistic EOS

We computed the energy spectrum of gravitational waves obtained as the first derivatives of the fitted functions (Eq. (9)). The relation between  ${\rm d}E_{\rm bind}/{\rm d}f$ and the GW frequency  $f_{\rm GW}$ is presented in Fig. 8. In this figure we draw straight lines corresponding to the Newtonian case $\sim$ $f_{\rm GW}^{2/3}$ to find the break frequencies at which the energy spectrum has dropped by 20%, 25%, 40%. These values are important from the point of view of future detections: they show the difference between the amplitude of the real signal and the Newtonian template which allows to calculate the real wave form amplitude from the detector noise. We also compare our results to the 25% case from Faber et al. (2002). There is no visible difference between our models for different EOSs at the break frequency level of 10% (the case considered by Faber et al. 2002) and the situation is then very precisely described by the 3PN formula.

Table 3: Gravitational wave frequencies (in Hz) computed from the calculated data for GNH3, APR and BPAL12 EOSs: the $f_{\rm npm}$ denotes the frequency at which the non-point-mass effects start to be important, $f_{\rm 20}, f_{\rm 25}$ and  $f_{\rm 40}$ are the so-called break-frequencies (see text), whereas  $f_{\rm end}$ is the GW frequency on the final orbit.

   
4.5 Comparison with polytropic EOS

Up to now, all calculations (except those of Oechslin et al. 2004) of the hydrodynamical inspiral and merger phases have been done for the simplified equation of state of dense matter, for the polytropic EOS, where the dependence between pressure and baryon density is given by $p=\kappa {n}^{\gamma}$. It has been shown that the results obtained for these polytropic EOSs depend mainly on the compaction parameter M/R. It is therefore interesting to check if the properties of inspiraling neutron stars described by a realistic EOS can be predicted, in a good approximation, by studying binaries with assumed polytropic EOSs. In order to make such a comparison we construct sequences of binary NSs described by polytropic EOS parametrized by the compaction parameter M/R given in Table 1 for each of the three realistic equations of state. For a given $\gamma$, we calculate the value of the $\kappa$ coefficient which yields the same R at $M=1.35~{M}_\odot$ as that predicted by a selected realistic EOS used in the present paper. The values of $\gamma$, the pressure coefficient $\kappa$, and the baryon masses of polytropic static isolated NSs of $M=1.35~{M}_\odot$ are presented in Table 4.

	Figure 8: Energy spectrum of GW waves emitted by the binary neutron star system versus frequency of gravitational waves (twice the orbital frequency) along three irrotational quasi-equilibrium sequences. The straight lines correspond to the Newtonian dependence of energy multiplied by 1, 0.8, 0.75 and 0.6.
Open with DEXTER

Table 4: Pressure coefficient $\kappa$, adiabatic index $\gamma$ and baryon mass $M_{\rm B}$ for polytropic NSs having the same compactness parameter and mass ( $M=1.35~{M}_\odot$) than NSs described by realistic EOSs (Table 1) ( $\hat{\rho}:=1.66$ $\times$ $10^{14}~{\rm g/cm^3}$ and $\hat{n}:=0.1~{\rm fm^{-3}}$).

Figure 9: Gravitational mass versus stellar radius for sequences of static neutron stars described by the Akmal et al. (1998) EOS (solid line) and polytropic EOS with $\gamma =2$ (dotted line) and $\gamma =2.5$ (dashed line). Configurations with gravitational mass $1.35~M_{\odot}$ (marked by a circle) described by polytropic EOSs have the same compaction parameter, M/R = 0.176, as a neutron star configuration based on the APR EOS.

Open with DEXTER

	Figure 10: Binding energy of a binary neutron star system as a function of gravitational wave frequency for the APR EOS (thick solid line) and two corresponding polytropic EOS with $\gamma =2.5$ (dashed line) and $\gamma =2$ (dotted line). Static isolated NSs have M/R=0.176 and $M=1.35~{M}_\odot$. The thin solid line shows the 3PN approximation for point masses by Blanchet (2002), Eq. (8).
Open with DEXTER

In Fig. 9 we show the dependence of gravitational mass versus radius for isolated non-rotating stars based on the APR EOS, and two different polytropes having M/R = 0.176 and $M=1.35~{M}_\odot$at infinite separation. Finally in Fig. 10 we present the binding energy versus the frequency of gravitational waves at the last orbits of the inspiral for the APR EOS and the two corresponding polytropes. The differences between quasi-equilibrium sequences described by realistic and polytropic sequences are small. For the three different EOSs (APR and two polytropes) the frequency of gravitational waves at the last calculated orbit ( $f_{\rm end}$) is $\sim$ $1140~{\rm Hz}$. Also the binding energy of the system at the mass-shedding point is close to each other, between -0.0372 and  $-0.0366~M_\odot$. However one can see the differences in the frequencies  $f_{\rm npm}$ at which non-point-mass effects start to be important. The  $f_{\rm npm}$ has the smallest value for the polytrope with $\gamma =2.5$ and the highest for $\gamma =2$. Although the matter in the stellar interior is stiffer for the APR EOS ( $\gamma \simeq 3$, see Fig. 3), which is also clear from Fig. 9, the APR curve lies between the $\gamma =2.5$ and $\gamma =2$ polytropes in Fig. 10. This is due to the relatively soft equation of state for the crust for realistic NS models, which makes the response of the crust to tidal forces different from that of a polytrope with $\gamma =2.5$ or $\gamma =2$.

We obtained similar results by comparing two other sequences of realistic EOS - BPAL12 and GNH3 with corresponding polytropic cases. The irrotational flow is weakly affected by the changes in the EOS of the core, but it is expected that the differences should be seen in the merger phase. The outer layers of the star (those with sub-nuclear densities, i.e., the crust), which are properly treated in the present paper have an influence on the properties of the binary system at the last stages of inspiral. However the crucial parameter which determines the energy-frequency spectrum (energy per frequency) of emitted gravitational energy is M/R.

   
5 Concluding remarks

We have presented a set of evolutionary sequences of binary neutron stars based on three selected realistic EOSs. These EOSs are based on modern many-body calculations. Three baryonic EOSs of a neutron-star core have been joined with a recent EOS of a neutron-star crust, and in this way we obtained three different models of the neutron star interior, from the surface to the stellar center. We restricted our analysis to models of neutron star cores without exotic phases (meson condensates, quark matter). In this way, the differences between the core EOSs reflect the uncertainties in the existing theories of the interactions in nuclear matter.

In the present paper we considered only those constituents of dense matter, that have been studied in the laboratory. We did not consider phase transitions to hypothetical exotic phases of dense matter, which were proposed by many authors, but which still remain speculative. Results obtained for the NS-NS binaries with exotic-phase neutron-star cores and realistic envelopes will be considered in future publications. Similarly, the case of a binary involving strange quark stars built of self-bound strange quark matter will also be presented in a separate paper.

We have computed quasi-equilibrium sequences of an irrotational NS-NS binary by keeping the baryon mass constant to a value that corresponds to individual gravitational masses of  $1.35~M_{\odot}$ at infinite separation. For a long time of evolution of the binary system its binding energy is very accurately given by the 3PN post-Newtonian formula for a point-mass system. However the departure from this 3PN model at low binary periods has a quite abrupt character, presumably due to high order tidal effects. The sequences end at the onset of the mass transfer between the stars (i.e. when a cusp forms at the surface of the stars). This point defines the ISCO since no turning point of the binding energy has been found along the sequences, which would have revealed some orbital instability. The gravitational wave frequency at the ISCO is 800 Hz, 1130 Hz and 1230 Hz, for the GNH3, APR and BPAL12 EOSs respectively.

In a recent work based on the numerical integration of the full set of time-dependent Einstein equations, Marronetti et al. (2004) have located the dynamical ISCO by comparing the time evolution of quasi-equilibrium initial data at various separations (see also Fig. 14 of Shibata & Uryu 2001). This defines the true ISCO, as opposed to the "quasi-equilibrium'' ISCO obtained here. For a polytropic EOS with $\gamma =2$ and a compactness parameter M/R=0.14, they obtain the ISCO at an orbital frequency which is 15% lower than the mass-shedding frequency of the quasi-equilibrium sequence of Taniguchi & Gourgoulhon (2002b). This makes us confident that the values of the gravitational wave frequencies given above are quite close to those of the end of the true inspiral.

When comparing our results with those of the recent work of Oechslin et al. (2004), one should stress the basic difference in the EOS of matter at sub-nuclear densities. Oechslin et al. represented the EOS of the crust by a polytrope with the adiabatic index $\gamma=2.86$. In this way, they made their EOS quasi-continuous at the crust-core interface, because their nuclear core EOS is very stiff. However, such an EOS for the crust is very unrealistic, because the real $\gamma$ can be as low as 0.5 near the neutron drip point. As the crust EOS is crucial for the last stable orbits of the NS binary, this implies differences between their and our results, even if for a $M=1.35~{M}_\odot$ model, the compactness M/R resulting from the GNH3 EOS and their EOS are very close. In particular, Oechslin et al. (2004) have found a turning point in the binding energy along their sequences, resulting in a quasi-equilibrium ISCO. This difference is certainly due to the stiffness of the polytropic EOS used by them in the outer layers of the star. For polytropic irrotational binaries a turning point ISCO exists only if $\gamma > 2.5$ (Uryu et al. 2000; Taniguchi & Gourgoulhon 2003).

In our calculations we treated the neutron-star matter as an ideal fluid. In other words, we neglected the elastic shear (in the crust - if not molten) and viscous (in the crust and in the core) terms in the matter stress tensor. There, the terms are believed to be small, but they lead to specific physical phenomena. In particular, the matter flow in NS interior will break the beta equilibrium between baryons and leptons, and this will imply a neutrino burst at the last stage of inspiral (see Haensel 1992). Moreover, dissipative processes will heat the matter. Both effects will be studied in future publications.

All numerical results presented here, including the full binary configurations, are available for downloading from the LORENE main server http://www.lorene.obspm.fr/data/

Acknowledgements

We are grateful to Luc Blanchet for providing us with the formulae for the binding energy at various post-Newtonian orders. We also warmly thank Koji Uryu for useful discussions and careful reading of the manuscript. This work has been funded by the KBN grants 5P03D.017.21, 2P03D.019.24, 1P03D.008.27 and PBZ-KBN-054/P03/2001 and has been partially supported by the Associated European Laboratory Astro-PF (Astrophysics Poland-France) and by the "Bourses de recherche 2004 de la Ville de Paris''. K.T. also acknowledges a Grant-in-Aid for Scientific Research (No. 14-6898) of the Japanese Ministry of Education, Culture, Sports, Science and Technology.

References

• Abbott, B., Abbott, R., Adhikari, R., et al. 2004, Nucl. Instrum. Meth., 517, 154 [NASA ADS] (In the text)
• Acernese, F., Amico, P., Arnaud, N., et al. 2004, Class. Quantum Grav., 21, S385 [CrossRef] (In the text)
• Akmal, A., Pandharipande, V. R., & Ravenhall, D. G. 1998, Phys. Rev. C, 58, 1804 [NASA ADS] [CrossRef] (In the text)
• Ando, M., Arai, K., Takahashi, R., et al. 2001, Phys. Rev. Lett., 86, 3950 [NASA ADS] [CrossRef] (In the text)
• Baumgarte, T. W., & Shapiro, S. L. 2003, Phys. Rep., 376, 41 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Baym, G., Pethick, C., & Sutherland, P. 1971, ApJ, 170, 299 [NASA ADS] [CrossRef] (In the text)
• Bildsten, L., & Cutler, C. 1992, ApJ, 400, 175 [NASA ADS] [CrossRef]
• Blanchet, L. 2002, Phys. Rev. D, 65, 124009 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Bombaci I. 1995, in Perspectives on Theoretical Nuclear Physics, ed. I. Bombaci, A. Bonaccorso, A. Fabrocini, et al., 223
• Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 1999, Phys. Rev. Lett., 82, 892 [NASA ADS] [CrossRef]
• Damour, T., Jaranowski, P., & Schäfer, G. 2000, Phys. Rev. D, 62, 084011 [NASA ADS] [CrossRef] (In the text)
• Douchin, F., & Haensel, P. 2000, Phys. Lett. B, 485, 107 [NASA ADS] [CrossRef] (In the text)
• Douchin, F., & Haensel, P. 2001, A&A, 380, 151 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Eriguchi, Y., & Hachisu, I. 1985, A&A, 142, 256 [NASA ADS] (In the text)
• Faber, J. A., Grandcl $\acute{\rm e}$ment, P., Rasio, F. A., & Taniguchi, K. 2002, Phys. Rev. Lett., 89, 1102
• Friedman, J. L., Uryu, K., & Shibata, S. 2002, Phys. Rev. D, 65, 064035 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
• Glendenning, N. K. 1985, ApJ, 293, 470 [NASA ADS] [CrossRef] (In the text)
• Gourgoulhon, E., Grandcl $\acute{\rm e}$ment, P., Taniguchi, K., Marck, J.-A., & Bonazzola, S. 2001, Phys. Rev. D, 63, 064029 [NASA ADS] [CrossRef]
• Haensel, P. 1992, A&A, 262, 131 [NASA ADS] (In the text)
• Haensel, P. 2001, in Physics of Neutron Star Interiors, ed. D. Blaschke, N. K. Glendenning, & A. Sedrakian, Lecture Notes in Physics (Springer), 127 (In the text)
• Haensel, P. 2003, in Final Stages of Stellar Evolution, ed. J.-M. Hameury, & C. Motch, EAS Publ. Ser., 7, 249 (EDP Sciences) (In the text)
• Haensel, P., & Pichon, B. 1994, A&A, 283, 313 [NASA ADS] (In the text)
• Hewitson, M., Aufmuth, P., Aulbert, C., et al. 2003, Class. Quantum Grav., 20, S581 [CrossRef] (In the text)
• Kochaneck, C. S. 1992, ApJ, 398, 234 [NASA ADS] [CrossRef]
• Lai, D., Rasio, F., & Shapiro, S. L. 1994, ApJ, 420, 811 [NASA ADS] [CrossRef] (In the text)
• Marronetti, P., Mathews, G. J., & Wilson, J. R. 1999, Phys. Rev. D, 60, 087301 [NASA ADS] [CrossRef]
• Marronetti, P., Duez, M. D., Shapiro, S. L., & Baumgarte, T. W. 2004, Phys. Rev. Lett., 92, 141101 [NASA ADS] [CrossRef]
• Oechslin, R., Uryu, K., Poghosyan, G., & Thielemann, F. K. 2004, MNRAS, 349, 1469 [NASA ADS] [CrossRef] (In the text)
• Prakash, M., Bombaci, I., Prakash, M., et al. 1997, Phys. Rep., 280, 1 [NASA ADS] [CrossRef]
• Prakash, M., & Lattimer, J. M. 2004, J. Phys. G: Nucl. Part. Phys., 30, S451 [CrossRef]
• Prakash, M., Ratkovic, S., & Lattimer, J. M. 2004, J. Phys. G: Nucl. Part. Phys., 30, S1279 [NASA ADS] [CrossRef]
• Shibata, M., & Uryu, K. 2001, Phys. Rev. D, 64, 104017 [NASA ADS] [CrossRef]
• Shibata, M., Taniguchi, K., & Uryu, K. 2003, Phys. Rev. D, 68, 084020 [NASA ADS] [CrossRef] (In the text)
• Taniguchi, K., Gourgoulhon, E., & Bonazzola, S. 2001, Phys. Rev. D, 64, 064012 [NASA ADS] [CrossRef]
• Taniguchi, K., & Gourgoulhon, E. 2002a, Phys. Rev. D, 65, 044027 [NASA ADS] [CrossRef]
• Taniguchi, K., & Gourgoulhon, E. 2002b, Phys. Rev. D, 66, 104019 [NASA ADS] [CrossRef]
• Taniguchi, K., & Gourgoulhon, E. 2003, Phys. Rev. D, 68, 124025 [NASA ADS] [CrossRef]
• Uryu, K., & Eriguchi, Y. 2000, Phys. Rev. D, 61, 124023 [NASA ADS] [CrossRef]
• Uryu, K., Shibata, M., & Eriguchi, Y. 2000, Phys. Rev. D, 62, 104015 [NASA ADS] [CrossRef]