# Numerical models of irrotational binary neutron stars in general relativity

## Abstract

We report on general relativistic calculations of quasiequilibrium configurations of binary neutron stars in circular orbits with zero vorticity. These configurations are expected to represent realistic situations as opposed to corotating configurations. The Einstein equations are solved under the assumption of a conformally flat spatial 3-metric (Wilson-Mathews approximation). The velocity field inside the stars is computed by solving an elliptical equation for the velocity scalar potential. Results are presented for sequences of constant baryon number (evolutionary sequences). Although the central density decreases much less with the binary separation than in the corotating case, it still decreases. Thus, no tendency is found for the stars to individually collapse to black hole prior to merger.

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 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 well approximated by zero vorticity (i.e. irrotational) states.

Moreover, 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). As argued by Mathews et al. , one way to settle this crucial point is to perform computations of relativistic irrotational configurations. We present here the first quasiequilibrium irrotational relativistic binary neutron stars models computations.

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.

The hydrodynamic equations may be derived as follows. For a perfect fluid at zero temperature, the momentum-energy conservation equation $`𝐓=0`$ is equivalent to the uniformly canonical equation of motion,

$$𝐮(𝐰)=0$$
(1)   
and

$$(n𝐮)=0,$$
(2)   
where $`𝐮`$ is the fluid 4-velocity and $`𝐰`$ the momentum density 1-form $`𝐰=h𝐮`$, $`h`$ being the fluid specific enthalpy $`h=(e+p)/(m_\mathrm{B}n)`$ ($`𝐰`$ denotes the exterior derivative of w). In the above equations, $`n`$, $`e`$ and $`p`$ denote respectively the fluid proper baryon density, proper total energy density and pressure. It is clear that a potential flow

$$𝐰=\mathrm{\Psi }$$
(3)   
is a solution of the equation of motion (1). Moreover, this particular solution is the relativistic generalization of the classical irrotational flow and, as stated above, corresponds to the physical situation reached by a binary system of neutron stars.

As a first approximation of the relativistic treatment of the problem, we will assume that there exists a helicoidal symmetry . Let us denote by l the associated Killing vector. It is to be noticed that this symmetry is exact at the Newtonian limit. The helicoidal symmetry implies $`_𝐥𝐰=0`$. From Cartan’s identity $`_𝐥𝐰=𝐥𝐰+(𝐥𝐰)`$, the potential form (3) leads immediately to the following first integral of motion

$$𝐥𝐰=\mathrm{Const}.$$
(4)   
This was first pointed out by Carter . Note that this result is not merely the relativistic generalization of the Bernoulli theorem which states that $`𝐥𝐰`$ 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.

The fluid motion is now completely determined by the scalar potential $`\mathrm{\Psi }`$. The equation for $`\mathrm{\Psi }`$ can be derived from the baryon number conservation (2). One obtains:

$$\frac{n}{h}\mathrm{\Psi }+\mathrm{\Psi }\left(\frac{n}{h}\right)=0.$$
(5)   
Within the 3+1 formalism and taking into account the helicoidal symmetry, this last equation is nothing but a Poisson-like equation which reads<sup>*</sup><sup>*</sup>*Latin indices run in the range 1,2,3 and geometrized units ($`G=1`$ and $`c=1`$) are used.

$$nD_iD^i\mathrm{\Psi }+D^inD_i\mathrm{\Psi }=\frac{h\mathrm{\Gamma }_\mathrm{n}}{N}B^iD_in+n\left\{\left(D^i\mathrm{\Psi }+\frac{h\mathrm{\Gamma }_\mathrm{n}}{N}B^i\right)D_i\mathrm{ln}hD^i\mathrm{\Psi }D_i\mathrm{ln}N\frac{B^i}{N}D_i(h\mathrm{\Gamma }_\mathrm{n})\right\}+Knh\mathrm{\Gamma }_\mathrm{n},$$
(6)   
where we have introduced the covariant derivative $`D`$ with respect to the spatial 3-metric, the trace $`K`$ of the extrinsic curvature tensor, the shift vector $`𝐁`$ and the lapse $`N`$ defined by the 3+1 orthogonal decomposition of the helicoidal Killing vector $`𝐥=N𝐧𝐁`$ ($`𝐧`$ being the unit future directed normal vector to the hypersurface $`t=\mathrm{Const}`$) and the Lorentz factor $`\mathrm{\Gamma }_\mathrm{n}=𝐧𝐮`$ of the fluid with respect to the Eulerian observer whose 4-velocity is $`𝐧`$. This equation has been recently derived by Teukolsky  and independently by Shibata . Note that Eq. (6) is independent of the gravitational field equations.

As a first milestone in our project of studying coalescing binary systems, we will adopt the Wilson-Mathews approximation for the form of the metric . This approximation consists in taking a conformally flat 3-metric, so that the full spacetime metric reads

$$ds^2=(N^2B_iB^i)dt^22B_idtdx^i+A^2\eta _{ij}dx^idx^j,$$
(7)   
where $`𝜼`$ is the flat space metric. The field equations reduce now to the Wilson-Mathews equations  for $`N`$, $`𝐁`$ and $`A`$. Right now, it is not obvious whether the Wilson-Mathews approximation is valid in the case of coalescing compact binaries. We presently use this particular form of the metric in order to simplify the problem. However, it is to be noticed that (i) the 1-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 . Finally we chose maximal slicing: $`K=0`$.

Following , we introduce $`\mathrm{\Omega }`$ such that the helicoidal Killing vector l satisfies

$$𝐥=\frac{}{\tau }+\mathrm{\Omega }\frac{}{\varphi },$$
(8)   
where $`\tau `$ and $`\varphi `$ are respectively the time and azimuthal coordinate associated with the asymptotic inertial observer at rest with respect to the binary system (i.e. such that the ADM 3-momentum vanishes on the slices $`\tau =\mathrm{const}`$). Besides, we introduce the non-rotating shift vector $`𝐍`$ defined by

$$𝐁=𝐍\mathrm{\Omega }\frac{}{\varphi }.$$
(9)   
The gravitational field equations are derived within the 3+1 formalism from the Hamiltonian constraint, momentum constraint and trace of the spatial part of the Einstein equation . Introducing $`\nu =\mathrm{ln}N`$ and $`\beta =\mathrm{ln}(AN)`$, they can be written as

$$\underset{¯}{\mathrm{\Delta }}\beta =4\pi A^2S+\frac{3}{4}A^2K_{ij}K^{ij}\frac{1}{2}\left(\overline{}_i\nu \overline{}^i\nu +\overline{}_i\beta \overline{}^i\beta \right)$$
(10)   
$$\underset{¯}{\mathrm{\Delta }}\nu =4\pi A^2(E+S)+A^2K_{ij}K^{ij}\overline{}_i\nu \overline{}^i\beta $$
(11)   
$$\underset{¯}{\mathrm{\Delta }}N^i+\frac{1}{3}\overline{}^i\left(\overline{}_jN^j\right)=16\pi NA^2(E+p)U^i+2NA^2K^{ij}\overline{}_j(3\beta 4\nu ),$$
(12)   
where $`\overline{}`$ is the covariant derivative associated with the flat 3-metric $`𝜼`$ and $`\underset{¯}{\mathrm{\Delta }}=\overline{}^i\overline{}_i`$ is the corresponding LaplacianThroughout the article, we use the notation $`\overline{}^i=\eta ^{ij}\overline{}_j`$.. $`E=\mathrm{\Gamma }_\mathrm{n}^2(e+p)p`$, $`S=3p+(E+p)U_iU^i`$ and $`U^i=D^i\mathrm{\Psi }/(h\mathrm{\Gamma }_\mathrm{n})`$ are respectively the fluid energy density, trace of the stress tensor and fluid 3-velocity, the three of them measured by the Eulerian observer. $`\mathrm{\Gamma }_\mathrm{n}`$ can be computed according to

$$\mathrm{\Gamma }_\mathrm{n}=\left[1+\frac{1}{A^2h^2}\overline{}_i\mathrm{\Psi }\overline{}^i\mathrm{\Psi }\right]^{1/2}=\frac{1}{\sqrt{1U_iU^i}}$$
(13)   
and the extrinsic curvature tensor is computed by means of the identity

$$K^{ij}=\frac{1}{2A^2N}\left\{\overline{}^iN^j+\overline{}^jN^i\frac{2}{3}\eta ^{ij}\overline{}_kN^k\right\},$$
(14)   
which results from the Killing equation for $`𝐥`$.

The matter distribution is determined by the first integral of motion (4). Taking its logarithm leads to

$$H+\nu +\frac{1}{2}\mathrm{ln}\left(1A^2\eta _{ij}\frac{B^iB^j}{N^2}\right)+\mathrm{ln}\mathrm{\Gamma }=\mathrm{Const},$$
(15)   
where $`H:=\mathrm{ln}h`$ and $`\mathrm{\Gamma }`$ is the fluid Lorentz factor with respect to the co-orbiting observer (i.e. the observer whose 4-velocity is collinear to $`𝐥`$):

$$\mathrm{\Gamma }=\mathrm{\Gamma }_\mathrm{n}\frac{1+A^2\eta _{ij}{\displaystyle \frac{B^i}{N}}U^j}{\sqrt{1A^2\eta _{ij}{\displaystyle \frac{B^iB^j}{N^2}}}}.$$
(16)   
Note that Eq. (15) corresponds to Eq. (66) in Ref.  and that the constant which appears at its r.h.s. is nothing but the logarithm of the constant $`C`$ introduced by Teukolsky .

Now, introducing $`\zeta =d\mathrm{ln}H/d\mathrm{ln}n`$, the equation for the fluid velocity potential (6) becomes

$$\zeta H\underset{¯}{\mathrm{\Delta }}\mathrm{\Psi }+\overline{}^iH\overline{}_i\mathrm{\Psi }=A^2h\mathrm{\Gamma }_\mathrm{n}\frac{B^i}{N}\overline{}_iH+\zeta H\left\{\left(\overline{}^i\mathrm{\Psi }+A^2h\mathrm{\Gamma }_\mathrm{n}\frac{B^i}{N}\right)\overline{}_iH\overline{}^i\mathrm{\Psi }\overline{}_i\beta A^2\frac{B^i}{N}\overline{}_i(h\mathrm{\Gamma }_\mathrm{n})\right\}.$$
(17)   
The equations to be solved (10) (11) (12) (17) constitute a system of non-linear Poisson-like equations. Because of the elliptical nature of these equations, we have exhibited the common flat Laplacian operator $`\underset{¯}{\mathrm{\Delta }}`$ which can be solved by means of usual spectral methods (cf. e.g.  or ). Because of the non-linearities, we use an iterative procedure based on a multi-domain spectral method  to get the solution. The numerical code will be described in details in a forthcoming paper . Let us simply mention here some 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`$ (cf. Fig. 6 of Ref. ). 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 in . 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 simply setting $`\mathrm{ln}\mathrm{\Gamma }=0`$ in Eq. (15) and using $`U^i=B^i/N`$ for the fluid 3-velocity\]. 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 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. ).

These tests being passed, we turned towards calculations of irrotational relativistic binaries. We chose to investigate the instability issue raised by Wilson and Mathews  by computing an evolutionary sequence (i.e. a sequence at fixed baryon number) made of irrotational quasiequilibrium models. We took the same configuration than that presented by Mathews, Marronetti and Wilson (Sect. IV A of Ref ), namely two identical stars obeying a $`\gamma =2`$ polytropic EOS \[$`p=\kappa (m_\mathrm{B}n)^\gamma `$, $`e=m_\mathrm{B}n+p/(\gamma 1)`$\] with $`\kappa =1.8\times 10^2\mathrm{J}\mathrm{m}^3\mathrm{kg}^2`$ and each having a baryon mass $`M_\mathrm{B}=1.625M_{}`$. For such parameters, we found that the gravitational mass of a single star in isolation is $`M=1.515M_{}`$ (in agreement with Ref. ), with a central energy density $`e_\mathrm{c}^{\mathrm{inf}}=4.005\rho _{\mathrm{nuc}}c^2`$ ($`\rho _{\mathrm{nuc}}:=1.66\times 10^{17}\mathrm{kg}\mathrm{m}^3`$) and a compactification ratio $`M/R=0.140`$ ($`e_\mathrm{c}`$ and $`M/R`$ are slightly different from that quoted in Ref.  probably due to different values for the constants $`G`$, $`c`$, $`M_{}`$ and $`m_\mathrm{B}`$; we use $`G=6.6726\times 10^{11}\mathrm{m}^3\mathrm{kg}^1\mathrm{s}^2`$, $`c=2.99792458\times 10^8\mathrm{m}\mathrm{s}^1`$, $`M_{}=1.989\times 10^{30}\mathrm{kg}`$ and $`m_\mathrm{B}=1.66\times 10^{27}\mathrm{kg}`$).

We define the coordinate separation $`d`$ as the coordinate distance between the two density maxima. Using the same value as the one considered by Mathews et al. , namely $`d=100\mathrm{km}`$, we found that a $`M_\mathrm{B}=1.625M_{}`$ irrotational configuration in quasiequilibrium at this separation has a total angular momentum of $`J/(2M_\mathrm{B})^2=1.13`$ which is quite similar to the value $`1.09`$ found by Mathews et al. . However, we did not observed any substantial increase of the central (i.e. maximum) density with respect to static stars in isolation, as they did (they report a central density increase of $`14\%`$).

In order to investigate the evolution of a coalescing binary system, we have computed a full sequence at fixed $`M_\mathrm{B}=1.625M_{}`$, starting at the separation $`d=110\mathrm{km}`$ ($`\mathrm{\Omega }/(2\pi )=82\mathrm{Hz}`$) and ending at $`d=41\mathrm{km}`$ ($`\mathrm{\Omega }/(2\pi )=332\mathrm{Hz}`$). We have considered both corotating and irrotational cases. The evolution of the central density along these sequences is shown in Fig. 1. In the corotating case, we found that the central density decreases substantially when the stars approach each other, as expected from previous independent calculations . In the irrotational case, we found that the central density still decreases with the separation but much less than in the corotating case.

All the calculations have been performed with spectral methods using the library LORENE (Langage Objet pour la RElativité NumériquE ). We warmly thank Jean-Pierre Lasota for his support. The calculations presented in this paper have been made possible thanks to a special grant from the SPM and SDU departments of CNRS.