# EVOLUTIONARY SEQUENCES OF IRROTATIONAL BINARY NEUTRON STARS

## 1 Introduction

Inspiraling neutron star binaries are expected to be among the strongest sources of gravitational radiation that could be detected by the interferometric detectors currently under construction (GEO600, LIGO, TAMA and Virgo). These binary systems are therefore subject to numerous theoretical studies. Among them are (i) Post-Newtonian (PN) analytical treatments (e.g. , , ) and (ii) fully relativistic hydrodynamical treatments, pioneered by the works of Oohara and Nakamura (see e.g. ) and Wilson et al. . The most recent numerical calculations, those of Baumgarte et al.  and Marronetti et al. , rely on the approximations of (i) a quasiequilibrium state and (ii) of synchronized binaries. Whereas the first approximation is well justified before the innermost stable orbit, the second one does not correspond to physical situations, since it has been shown that the gravitational-radiation driven evolution is too rapid for the viscous forces to synchronize the spin of each neutron star with the orbit  as they do for ordinary stellar binaries. Rather, the viscosity is negligible and the fluid velocity circulation (with respect to some inertial frame) is conserved in these systems. Provided that the initial spins are not in the millisecond regime, this means that close binary configurations are better approximated by zero vorticity (i.e. irrotational) states than by synchronized states.

Dynamical calculations by Wilson et al.  indicate that the neutron stars may individually collapse into a black hole prior to merger. This unexpected result has been called into question by a number of authors (see Ref.  for a summary of all the criticisms and their answers). Recently Flanagan  has found an error in the analytical formulation used by Wilson et al. . This error may be responsible of the observed radial instability. As argued by Mathews et al. , one way to settle this crucial point is to perform computations of relativistic irrotational configurations. We have performed recently such computations . They show no compression of the stars, although the central density decreases much less than in the corotating case. In the present report, we give more details about the results presented in Ref.  and extend them to the compactification ratio $`M/R=0.17`$ (the results of Ref.  have been obtained for a compactification ratio $`M/R=0.14`$).

## 2 Analytical formulation of the problem

### 2.1 Basic assumptions

We have proposed a relativistic formulation for quasiequilibrium irrotational binaries  as a generalization of the Newtonian formulation presented in Ref. . The method was based on one aspect of irrotational motion, namely the counter-rotation (as measured in the co-orbiting frame) of the fluid with respect to the orbital motion (see also Ref. ). Since then, Teukolsky  and Shibata  gave two formulations based on the definition of irrotationality, which implies that the specific enthalpy times the fluid 4-velocity is the gradient of some scalar field  (potential flow). The three formulations are equivalent; however the one given by Teukolsky and by Shibata greatly simplifies the problem. Consequently we used it in the present work.

The irrotational hypothesis amounts to say that the co-momentum density is the gradient of a potential

$$h𝐮=\mathrm{\Psi },$$
(1)   
where $`h`$ and $`𝐮`$ are respectively the fluid specific enthalpy and fluid 4-velocity.

Beside the physical assumption (1), two simplifying approximations are introduced:

1. The spacetime is supposed to have a helicoidal symmetry , which means that the orbits are exactly circular and that the gravitational radiation content of spacetime is neglected.
2. The spatial 3-metric is assumed to be conformally flat (Wilson-Mathews approximation ), so that the full spacetime metric reads

$$ds^2=(N^2B_iB^i)dt^22B_idtdx^i+A^2f_{ij}dx^idx^j,$$
(2)   
where $`f_{ij}`$ is the flat space metric.

The Killing vector corresponding to hypothesis (i) is denoted by $`𝐥`$ (cf. Fig. 1).

Approximation (i) is physically well motivated, at least up to the innermost stable orbit. Regarding the second approximation, it is to be noticed that (i) the 1-Post Newtonian (PN) approximation to Einstein equations fits this form, (ii) it is exact for arbitrary relativistic spherical configurations and (iii) it is very accurate for axisymmetric rotating neutron stars . An interesting discussion about some justifications of the Wilson-Mathews approximation may be found in . A stronger justification may be obtained by considering the 2.5-PN metric obtained by Blanchet et al.  for point mass binaries. Using Eq. (7.6) of Ref. , the deviation from a conformally flat 3-metric (which occurs at the 2-PN order) can be computed at the location of one point mass (i.e. where it is expected to be maximal), the 3-metric $`h_{ij}`$ being written as

$$h_{ij}=A^2f_{ij}+h_{ij}^{2\mathrm{P}\mathrm{N}}+h_{ij}^{2.5\mathrm{PN}}.$$
(3)   
The result is shown in Table 1 for two stars of $`1.4M_{}`$ each. It appears that at a separation as close as $`30\mathrm{km}`$, where the two stars certainly almost touch each other, the relative deviation from a conformally flat 3-metric is below $`2\%`$.

### 2.2 Equations to be solved

We refer to Ref.  for the presentation of the partial differential equations (PDEs) which result from the assumptions presented above. In the present report, let us simply mention a point which seems to have been missed by various authors: the existence of the first integral of motion

$$h𝐥𝐮=\mathrm{const}.$$
(4)   
does not result solely from the existence of the helicoidal Killing vector $`𝐥`$. Indeed, Eq. (4) is not merely the relativistic generalization of the Bernoulli theorem which states that $`h𝐥𝐮`$ is constant along each single streamline and which results directly from the existence of a Killing vector without any hypothesis on the flow. In order for the constant to be uniform over the streamlines (i.e. to be a constant over spacetime), as in Eq. (4), some additional property of the flow must be required. One well known possibility is rigidity (i.e. $`𝐮`$ colinear to $`𝐥`$) . The alternative property with which we are concerned here is irrotationality \[Eq. (1)\]. This was first pointed out by Carter .

## 3 Numerical procedure

### 3.1 Description

The numerical procedure to solve the PDE system is based on the multi-domain spectral method presented in Ref. . We simply recall here some basic features of the method:

* Spherical-type coordinates $`(\xi ,\theta ^{},\phi ^{})`$ centered on each star are used: this ensures a much better description of the stars than with Cartesian coordinates.
* These spherical-type coordinates are surface-fitted coordinates: i.e. the surface of each star lies at a constant value of the coordinate $`\xi `$ thanks to a mapping $`(\xi ,\theta ^{},\phi ^{})(r,\theta ,\phi )`$ (see  for details about this mapping). This ensures that the spectral method applied in each domain is free from any Gibbs phenomenon.
* The outermost domain extends up to spatial infinity, thanks to the mapping $`1/r=(1\xi )/(2R_0)`$. This enables to put exact boundary conditions on the elliptic equations for the metric coefficients: spatial infinity is the only location where the metric is known in advance (Minkowski metric).
* Thanks to the use of a spectral method  in each domain, the numerical error is evanescent, i.e. it decreases exponentially with the number of coefficients (or equivalently grid points) used in the spectral expansions, as shown in Fig. 2.

The PDE system to be solved being non-linear, we use an iterative procedure. This procedure is sketched in Fig. 3. The iteration is stopped when the relative difference in the enthalpy field between two successive steps goes below a certain threshold, typically $`10^7`$ (cf. Fig. 4).

The numerical code is written in the LORENE language , which is a C++ based language for numerical relativity. A typical run make use of $`N_r=33`$, $`N_\theta =21`$, and $`N_\phi =20`$ coefficients (= number of collocation points, which may be seen as number of grid points) in each of the domains on the multi-domain spectral method. 8 domains are used : 3 for each star and 2 domains centered on the intersection between the rotation axis and the orbital plane. The corresponding memory requirement is 155 MB. A computation involves $`250`$ steps (cf. Fig. 4), which takes 9 h 30 min on one CPU of a SGI Origin200 computer (MIPS R10000 processor at 180 MHz). Note that due to the rather small memory requirement, runs can be performed in parallel on a multi-processor platform. This especially useful to compute sequences of configurations.

### 3.2 Tests passed by the code

In the Newtonian and incompressible limit, the analytical solution constituted by a Roche ellipsoid is recovered with a relative accuracy of $`10^9`$, as shown in Fig. 2. For compressible and irrotational Newtonian binaries, no analytical solution is available, but the virial theorem can be used to get an estimation of the numerical error: we found that the virial theorem is satisfied with a relative accuracy of $`10^7`$. A detailed comparison with the irrotational Newtonian configurations recently computed by Uryu & Eriguchi  will be presented elsewhere. Regarding the relativistic case, we have checked our procedure of resolution of the gravitational field equations by comparison with the results of Baumgarte et al.  which deal with corotating binaries \[our code can compute corotating configurations by setting to zero the velocity field of the fluid with respect to the co-orbiting observer\]. We have performed the comparison with the configuration $`z_A=0.20`$ in Table V of Ref. . We used the same equation of state (EOS) (polytrope with $`\gamma =2`$), same value of the separation $`r_C`$ and same value of the maximum density parameter $`q^{\mathrm{max}}`$. We found a relative discrepancy of $`1.1\%`$ on $`\mathrm{\Omega }`$, $`1.4\%`$ on $`M_0`$, $`1.1\%`$ on $`M`$, $`2.3\%`$ on $`J`$, $`0.8\%`$ on $`z_A`$, $`0.4\%`$ on $`r_A`$ and $`0.07\%`$ on $`r_B`$ (using the notations of Ref. ).

## 4 Numerical results

### 4.1 Equation of state and compactification ratio

As a model for nuclear matter, we consider a polytropic equation of state (EOS) with an adiabatic index $`\gamma =2`$:

$$p=\kappa (m_\mathrm{B}n)^\gamma ,e=m_\mathrm{B}n+p/(\gamma 1),$$
(5)   
where $`p`$, $`n`$, $`e`$ are respectively the fluid pressure, baryon density and proper energy density, and $`m_\mathrm{B}=1.66\times 10^{27}\mathrm{kg}`$, $`\kappa =1.8\times 10^2\mathrm{J}\mathrm{m}^3\mathrm{kg}^2`$. This EOS is the same as that used by Mathews, Marronetti and Wilson (Sect. IV A of Ref ).

The mass – central density curve of static configurations in isolation constructed upon this EOS is represented in Fig. 5. The three points on this curve corresponds to three configurations studied by our group and that of Marronetti, Mathews and Wilson:

* The configuration of baryon mass $`M_\mathrm{B}=1.625M_{}`$ and compactification ratio $`M/R=0.14`$ is that considered in the dynamical study of Mathews, Marronetti and Wilson  and in the quasiequilibrium study of our group (Ref.  and this paper).
* The configuration of baryon mass $`M_\mathrm{B}=1.85M_{}`$ and compactification ratio $`M/R=0.17`$ is studied in the present paper.
* The configuration of baryon mass $`M_\mathrm{B}=1.95M_{}`$ and compactification ratio $`M/R=0.19`$ has been studied recently by Marronetti, Mathews and Wilson <sup>2</sup><sup>2</sup>2Marronetti et al.  use a different value for the EOS constant $`\kappa `$: their baryon mass $`M_\mathrm{B}=1.55M_{}`$ must be rescaled to our value of $`\kappa `$ in order to get $`M_\mathrm{B}=1.95M_{}`$. Apart from this scaling, this is the same configuration, i.e. it has the same compactification ratio $`M/R=0.19`$ and its relative distance with respect to the maximum mass configuration, as shown in Fig. 5, is the same. by means of a new code for quasiequilibrium irrotational configurations.

### 4.2 Irrotational sequence with $`M/R=0.14`$

In this section, we give some details about the irrotational sequence $`M_\mathrm{B}=1.625M_{}`$ presented in Ref. . This sequence starts at the coordinate separation $`d=110\mathrm{km}`$ (orbital frequency $`f=82\mathrm{Hz}`$), where the two stars are almost spherical, and ends at $`d=41\mathrm{km}`$ ($`f=332\mathrm{Hz}`$), where a cusp appears on the surface of the stars, which means that the stars start to break. The shape of the surface at this last point is shown in Fig. 6.

The velocity field with respect to the co-orbiting observer, as defined by Eq. (52) of Ref. , is shown in Fig. 7. Note that this field is tangent to the surface of the star, as it must be.

The lapse function $`N`$ (cf. Eq. 2) is represented in Fig. 8. The coordinate system $`(x,y,z)`$ is centered on the intersection between the rotation axis and the orbital plane. The $`x`$ axis joins the two stellar centers, and the orbital is the $`z=0`$ plane. The value of $`N`$ at the center of each star is $`N_\mathrm{c}=0.64`$.

The conformal factor $`A^2`$ of the 3-metric \[cf. Eq. (2)\] is represented in Fig. 9. Its value at the center of each star is $`A_\mathrm{c}^2=2.20`$.

The shift vector of nonrotating coordinates, $`𝐍`$, (defined by Eq. (9) of Ref. ) is shown in Fig. 10. Its maximum value is $`0.10c`$.

The $`K_{xy}`$ component of the extrinsic curvature tensor of the hypersurfaces $`t=\mathrm{const}`$ is shown in Fig. 11. We chose to represent the $`K_{xy}`$ component because it is the only component of $`K_{ij}`$ for which none of the sections in the three planes $`x=0`$, $`y=0`$ and $`z=0`$ vanishes.

The variation of the central density along the $`M_\mathrm{B}=1.625M_{}`$ sequence is shown in Fig. 12. We have also computed a corotating sequence for comparison (dashed line in Fig. 12). In the corotating case, the central density decreases quite substantially as the two stars approach each other. This is in agreement with the results of Baumgarte et al. . In the irrotational case (solid line in Fig. 12), the central density remains rather constant (with a slight increase, below $`0.1\%`$) before decreasing. We can thus conclude that no tendency to individual gravitational collapse is found in this case. This contrasts with results of dynamical calculations by Mathews et al.  which show a central density increase of $`14\%`$ for the same compactification ratio $`M/R=0.14`$.

### 4.3 Irrotational sequence with $`M/R=0.17`$

In order to investigate how the above result depends on the compactness of the stars, we have computed an irrotational sequence with a baryon mass $`M_\mathrm{B}=1.85M_{}`$, which corresponds to a compactification ratio $`M/R=0.17`$ for stars at infinite separation (second heavy dot in Fig. 5). The result is compared with that of $`M/R=0.14`$ in Fig. 13. A very small density increase (at most $`0.3\%`$) is observed before the decrease. Note that this density increase remains within the expected error ($`2\%`$, cf. Sect. 2.1) induced by the conformally flat approximation for the 3-metric, so that it cannot be asserted that this effect would remain in a complete calculation.

Marronetti, Mathews and Wilson  have recently computed quasiequilibrium irrotational configurations by means of a new code. They use a higher compactification ratio, $`M/R=0.19`$ (third heavy dot in Fig. 5). They found a central density increase as the orbit shrinks much pronounced than that we found for the compactification ratio $`M/R=0.17`$: $`3.5\%`$ against $`0.3\%`$. We will present irrotational sequences with the compactification ratio $`M/R=0.19`$ and compare with the results by Marronetti et al.  in a future article.

## 5 Innermost stable circular orbit

An important parameter for the detection of a gravitational wave signal from coalescing binaries is the location of the innermost stable circular orbit (ISCO), if any. In Table 2, we recall what is known about the existence of an ISCO for extended fluid bodies. The case of two point masses is discussed in details in Ref. .

For Newtonian binaries, it has been shown  that the ISCO is located at a minimum of the total energy, as well as total angular momentum, along a sequence at constant baryon number and constant circulation (irrotational sequences are such sequences). The instability found in this way is dynamical. For corotating sequences, it is secular instead . This turning point method also holds for locating ISCO in relativistic corotating binaries . For relativistic irrotational configurations, no rigorous theorem has been proven yet about the localization of the ISCO by a turning point method. All what can be said is that no turning point is present in the irrotational sequences considered in Sect. 4: Fig. 14 shows the variation as the orbit shrinks of the ADM mass of the spatial hypersurface $`t=\mathrm{const}`$ (which is a measure of the total energy, or equivalently of the binding energy, of the system) for the $`M_\mathrm{B}=1.625M_{}`$ sequence. Clearly, the ADM mass decreases monoticaly, without showing any turning point. Figure 15 shows the evolution of the total angular momentum along the same sequence. Again there is no turning point. The same behaviour holds for the $`M_\mathrm{B}=1.85M_{}`$ sequence, as shown in Figs. 16 and 17.

## 6 Conclusion and perspectives

We have computed evolutionary sequences of quasiequilibrium irrotational configurations of binary stars in general relativity. The evolution of the central density of each star have been monitored as the orbit shrinks. For a compactification ratio $`M/R=0.14`$, the central density remains rather constant (with a slight increase, below $`0.1\%`$) before decreasing. For a higher compactification ratio $`M/R=0.17`$ (i.e. stars closer to the maximum mass configuration), a very small density increase (at most $`0.3\%`$) is observed before the density decrease. It can be thus concluded that no substantial compression of the stars is found, which means that no tendency to individually collapse to black hole prior to merger is observed. Moreover, the observed density increase remains within the expected error ($`2\%`$, cf. Sect. 2.1) induced by the conformally flat approximation for the 3-metric, so that it cannot be asserted that this effect would remain in a complete calculation.

No turning point has been found in the binding energy or angular momentum along evolutionary sequences, which may indicate that these systems do not have any innermost stable circular orbit (ISCO).

All these results have been obtained for a polytropic EOS with the adiabatic index $`\gamma =2`$. We plan to extend them to other values of $`\gamma `$ in the near future. We also plan to abandon the conformally flat approximation for the 3-metric and use the full Einstein equations, keeping the helicoidal symmetry in a first stage.

## Acknowledgments

We would like to thank Jean-Pierre Lasota for his constant support and Brandon Carter for illuminating discussions. The numerical calculations have been performed on computers purchased thanks to a special grant from the SPM and SDU departments of CNRS.

## References