Constructing quasiequilibrium initial data for binary neutron stars with arbitrary spins
Abstract
In general neutron stars in binaries are spinning. Recently, a new quasiequilibrium approximation that includes a rotational velocity piece for each star has been proposed to describe binary neutron stars with arbitrary rotation states in quasicircular orbits. We have implemented this approximation numerically for the first time, to generate initial data for neutron star binaries with spin. If we choose the rotational velocity piece such that it equals the Newtonian rigid rotation law, we obtain stars with fluid 4velocities that have expansion and shear of approximately zero, as one would expect for quasiequilibrium configurations. We also use the new approach to construct and study initial data sequences for irrotational, corotating and fixed rotation binaries.
pacs:
04.20.Ex, 04.30.Db, 97.60.Jd, 97.80.FkI Introduction
Binary neutron stars are at the intersection of two of the most fascinating topics in astrophysics: Gamma ray bursts and gravitational wave astronomy. Binary neutron star mergers (together with black hole neutron star mergers) have been proposed as potential engines for short duration gamma ray bursts Narayan et al. (1992); Piran (2000, 1999); Zhang and Meszaros (2004); Piran (2005); Nakar (2007); Berger (2011); Rezzolla et al. (2011). These are likely generated in the massive accretion disks around the merger remnant: a larger neutron star or a black hole. In addition, binary neutron star systems are one of the most promising sources for gravitational wave detectors such as LIGO Abbott et al. (2009); LIGO (), Virgo Acernese et al. (2008); VIRGO () or GEO GEO 600 (). Several of these detectors have been operating over the last few years, while several others are in the planning or construction phase Schutz (1999). During the inspiral regime, when the two stars are still well separated they can be well approximated by postNewtonian theory. Later, when the stars get close, their matter distributions eventually merge together to form a single differentially rotating object. Depending on the total mass, the two progenitors’ spins, the equation of state and the strength of magnetic fields this object can either promptly collapse to a black hole, or form a hypermassive neutron star. The hypermassive neutron star is supported against collapse by differential rotation. It can survive for many dynamical timescales, while angular momentum is gradually transported from the inner to the outer parts. Eventually the hypermassive neutron star will also collapse to a black hole surrounded by a massive torus, that is more massive than in the prompt collapse case. Such systems could supply the energy required for a short gamma ray burst Narayan et al. (1992); Piran (2000, 1999); Zhang and Meszaros (2004); Piran (2005); Nakar (2007); Rezzolla et al. (2011).
In order to make predictions about the last few orbits and the merger of such systems, fully nonlinear numerical simulations of the Einstein Equations are required. To start such simulations we need initial data that describe the binary a few orbits before merger. The emission of gravitational waves tends to circularize the orbits Peters and Mathews (1963); Peters (1964). Thus, during the inspiral, we expect the two neutron stars to be in quasicircular orbits around each other with a radius that shrinks on a timescale much larger than the orbital timescale. This means that the initial data should have an approximate helical Killing vector . In general these neutron stars will be spinning. In the case of the double pulsar PSR J07373039 Lyne et al. (2004) the spin period of the faster spinning star will be at merger Tichy (2011) and should thus not be neglected. In the quasicircular regime the orbital timescale will be much shorter than the spin precession time scale, thus we can assume that the spins are approximately constant.
To incorporate these ideas and to construct numerical initial data for binary neutron stars with arbitrary spins and masses we will use an approach introduced in Tichy (2011). In this approach the stars are given spin by choosing a rotational velocity for each star. In this way it is possible to construct stars with both rigid or differential rotation. In equilibrium we of course expect the stars to be rigidly rotating such that the expansion and shear of the fluid 4velocities of each star vanish. We find that this can be achieved (to good approximation) by setting the rotational velocity of each star equal to the Newtonian rigid rotation law.
We also construct and discuss initial data sequences for irrotational, corotating and fixed rotation binaries.
Spin will have a noticeable effect on the inspiral and merger of the binary if the spin period is within an order of magnitude of the orbital period. Since the orbital period is on the order of a few milliseconds during the late inspiral, we expect interesting spin effects for spin periods on the order of a few dozen milliseconds or less. To date we have observed only ten binary neutron star systems, thus it is not clear yet how likely such spin periods are. Some population synthesis models Oslowski et al. (2011) suggest that radio observable pulsars in neutron star binaries do have a distribution of spin periods that extends down to about 15ms. Other population synthesis models for pulsars Kiel et al. (2008) come to similar findings. However, since such models involve many parameters that are used to describe sometimes poorly understood physical processes it may be too early to draw definite conclusions about what spin periods can be expected in neutron star binaries. One parameter that may be of particular importance and that illustrates this uncertainty is the magnetic field decay timescale . In the models it is assumed that the magnetic field of each star decays exponentially on this timescale. Since neutron stars spin down due to magnetic dipole radiation, magnetic field decay can have important effects for the spin down rate and thus the expected spin periods before merger. Unfortunately the value of is controversial. In the first model mentioned Oslowski et al. (2011) the magnetic field decay timescale has to be in order to fit observations, while in the other model Kiel et al. (2008) one needs .
Throughout we will use units where . Latin indices such as run from 1 to 3 and denote spatial indices, while Greek indices such as run from 0 to 3 and denote spacetime indices. The paper is organized as follows. Sec. II lists the General Relativistic equations that govern binary spinning neutron stars described by perfect fluids. In Sec. III we briefly describe what algorithm we use to numerically implement these equations. We then present some numerical results in Sec. IV. In particular we describe how one should choose the rotational velocity of each star. We also present particular examples in the form of several initial data sequences. We conclude with a discussion of our method in Sec. V. Some derivations are relegated to the appendices A and B.
Ii Binary neutron stars with arbitrary rotation states
In this section we briefly describe the equations governing binary neutron stars in arbitrary rotation states in General Relativity. These equations were derived in Tichy (2011).
We use the ArnowittDeserMisner (ADM) decomposition of Einstein’s equations and describe the gravitational fields (i.e. the 4metric ) in terms of the 3metric , lapse , shift and the extrinsic curvature . We further assume that the neutron star matter is a perfect fluid with stressenergy tensor
(1) 
Here is the mass density (which is proportional the number density of baryons), is the pressure, is the internal energy density divided by and is the 4velocity of the fluid. Assuming a polytropic equation of state
(2) 
and defining the specific enthalpy
(3) 
, and can all be expressed in terms of . We find
(4) 
where we have used the abbreviation
(5) 
The fluid 4velocity is expressed in terms of the 3velocity
(6) 
which in turn is split into a irrotational piece and a rotational piece
(7) 
where is the derivative operator compatible with the 3metric .
In order to simplify the problem and to obtain elliptic equations we make several assumptions. The first is the existence of an approximate helical Killing vector , such that
(8) 
We also assume similar equations for scalar matter quantities such as . For a spinning star, however, is nonzero. Instead we assume that
(9) 
so that the time derivative of the irrotational piece of the fluid velocity vanishes in corotating coordinates. We also assume that
(10) 
and
(11) 
which describe the fact that the rotational piece of the fluid velocity is constant along the world line of the star center.
These approximations together with the further assumptions of maximal slicing
(12) 
and conformal flatness
(13) 
yield the following coupled equations:
(14) 
(15) 
(16) 
(17) 
and
(18) 
Here , , and we have introduced
(19) 
(20)  
(21) 
where we sum over repeated spatial indices irrespective of whether they are up or down and where is a constant of integration that in general has a different value inside each star. Below we will denote these two values by and . Notice that in an inertial frame the approximate helical Killing vector has the components
(22) 
Here denotes the center of mass position of the system, and is the orbital angular velocity, which we have chosen to lie along the direction.
The elliptic equations (14), (15), (16) and (17) above have to be solved subject to the boundary conditions
(23) 
at spatial infinity, and
(24) 
at each star surface. Note that the rotational piece of the fluid velocity can be freely chosen.
Once the equations (14), (15), (16), (17) and (18) are solved we know (and thus the matter distribution) and the fluid 3velocity via Eq. (7). The 3metric is obtained from Eq. (13) and the extrinsic curvature is given by
(25) 
Our approach, while more general, has certain similarities with the approach in Baumgarte and Shapiro (2009a, b), hereafter BS. For example our Eq. (7) reduces to Eq. (52) of BS if we assume . Furthermore, using the approximations in Eqs. (9), (10) and (11), we are able to find the first integral (18) of the Euler equation. This means that within our approximations the Euler equation has vanishing curl. That is precisely what has to hold if the approach in BS (which converts the Euler equation into an elliptic equation by taking its divergence) is to succeed. The problem with the BS approach is that it simply assumes instead of our Eqs. (9), (10) and (11). This leads to extra terms in the Euler equation, which cause a nonvanishing curl. Thus, the BS approach is not entirely consistent since it requires zero curl of the Euler equation, while at the same time it uses an Euler equation with nonvanishing curl. Furthermore, the boundary condition given by BS for their new elliptic equation has to be imposed at the star surface. Since the location of the star surface is another unknown function one needs an equation or an algorithm to determine this location. No such algorithm is given by BS. The usual iterative approach where we simply search for the surface where in each iteration (see e.g. Sec. III) will not work, because the boundary condition given by BS ensures that occurs at the surface where we impose their boundary condition. So in this kind of iterative procedure the star surface would always remain at the location of the initial guess for the surface.
Iii Numerical method
To construct initial data we have to solve the elliptic equations (14), (15), (16) and (17) together with the algebraic equation (18). This set of equations has a similar structure as for the well known case of irrotational neutron star binaries, which has been solved before Bonazzola et al. (1999); Gourgoulhon et al. (2001); Marronetti et al. (1999); Uryu and Eriguchi (2000); Marronetti and Shapiro (2003); Taniguchi and Gourgoulhon (2002, 2003); Uryu et al. (2006, 2009). In this work we will use the SGRID code Tichy (2006, 2009a, 2009b), which uses pseudospectral methods to accurately compute spatial derivatives. We use the same decomposition into 6 domains as was used in Tichy (2009a) for the case of corotating neutron star binaries. In this approach one star center is put on the positive and the other on the negative axis. Then complicated coordinate transformations are used to transform from Cartesian like coordinates to new coordinates . Here the coordinates and both range from 0 to 1, and is a polar angle measured around the axis. The actual coordinate transformations are different in each domain. This results in two domains that cover the outside of each star (including spatial infinity) for either or . These two domains touch at . Since they contain spatial infinity it is trivial to impose the boundary conditions in Eq. (23). The coordinate transformations contain freely specifiable functions so that one can always make the inner domain boundaries coincide with the star surfaces. The inside of each star is covered by two more domains. One of these stretches from the star surface up to a certain depth inside the star. The other covers the remainder of the star interior. The elliptic equations (14), (15) and (16) need to be solved in all domains, while the matter equations (17) and (18) are solved only inside each star. Two of the domain boundaries always coincide with the neutron star surfaces so that it is straightforward to impose the boundary condition (24) for at each star surface. Notice, however, that Eq. (17) and its boundary condition in Eq. (24) do not uniquely specify a solution . If solves both Eqs. (17) and (24) will be a solution as well. In order to obtain a unique solution we modify the boundary condition by adding the volume integral over the star interior to the left hand side of Eq. (24). Furthermore, in the domains covered by the coordinates we impose the following regularity conditions along the axis:
(28) 
where stands for either , , or , and is the distance from the axis.
In order to solve the elliptic equations (14), (15), (16) and (17) we need a fixed domain decomposition. However, the location of the star surfaces (where ) is not a priorily known, but rather determined by Eq. (18). For this reason we use the following iterative procedure:

We first come up with an initial guess for in each star, in practice we simply choose TolmanOppenheimerVolkoff solutions (see e.g. Chap. 23 in Misner et al. (1973)) for each. For the irrotational velocity potential we choose , where and is the center of the star and the center of mass. We choose the initial orbital angular velocity according to postNewtonian theory.

In order to also solve Eq. (18) we need to know the values of the constants in each star as well as and . We first determine the star centers by finding the maximum of the current along the axis. Since the location of each star center is given by , Eq. (26) [which is equivalent to Eq. (18)] yields
(29) Note that is a function of and . We now update and by solving Eq. (5) for and so that the star centers remain in the same location, when we update according to Eq. (18) or Eq. (26). For this reason Eq. (5) is sometimes referred to as force balance equation. One noteworthy caveat is that we evaluate the derivative of in Eq. (5) for the and before the update. Since we iterate over these steps this approximation does not introduce any errors. However, it is essential for the overall stability of our iterative procedure.

Next, we use Eq. (18) to update in each star, while at the same time adjusting such that the rest mass of each star remains constant. This update is numerically expensive because the domain boundaries need to be adjusted (by changing ) such that they remain at the star surfaces, which change whenever is updated. When we adjust it can be helpful for the stability of the overall iteration to filter out high frequency modes in and to impose . The latter keeps the stars from drifting away from the axis during the iterations.

Finally we go back to step 2.
Iv Numerical results
We have implemented the method described above in the SGRID code Tichy (2006, 2009a, 2009b). In this section we present some numerical results using this code. All our results are presented in units where and we only use in the polytropic equation of state.
iv.1 Choice of rotational velocity piece
As already mentioned we hold the rest masses fixed during the iterations described in Sec. III. Similarly we need to fix the rotational velocity piece that gives rise to the spin in each star. Ideally we would like to choose such that the expansion and shear of the fluid are approximately zero. As we show in Appendices A and B, this is true if the velocity in the corotating frame is of the form
(30) 
where is the location of the star center, which could be defined as the point with the highest rest mass density or as the center of mass of the star. Thus our task is to choose such that Eq. (30) holds. In Tichy (2011) we have speculated that a good choice might be , which can be rewritten in terms of the derivative operator as if we introduce
(31) 
It is clear that
(32) 
satisfies for any function that only depends on the conformal distance from the star’s center. We have tested the simplest case of , and find that the that results from this choice after solving all equations is similar to Eq. (30) but with a position dependent , so that we have differential rotation in the star which leads to nonzero shear. Another simple choice is
(33) 
Numerically, we find that this choice results in a that is very close to the form in Eq. (30).
iv.2 Initial data sequences
In order to test our method we have performed simulations with different rotation states.
In Figs. 2 and 3 we show how the ADM mass and the ADM angular momentum vary as a function of the orbital angular velocity for a binaries with rest masses and . Note that increasing corresponds to decreasing separation. As we can see our numerical results (pluses, stars and crosses) approach the expected postNewtonian results (taken from Schäfer and Wex (1993); Tichy et al. (2003a, b); Tichy and Brügmann (2004)) for nonspinning point particles (solid line) for small . In fact we find that the results for irrotational () stars (marked by stars in Figs. 2 and 3) agree very well with postNewtonian nonspinning point particle results for all we have investigated. On the other hand, for corotating binaries (pluses) we obtain larger and values, especially for larger , which is expected because to maintain corotation the stars have to spin more for higher . Finally, we have also investigated the case of a binary where both stars have the same constant rotational velocity with . This corresponds to a spin period of 14ms [for m/(kg s)]. In Figs. 2 and 3 this case is denoted by crosses. Since here the stars always have rotational velocity we obtain larger and values than in the irrotational case for all . However, since the rotational velocity is always the same we get less and than in the corotating case (pluses) for large , and more and than in the corotating case (pluses) for small .
V Discussion
Realistic neutron stars in binaries are spinning. From observations of millisecond pulsars we know that these spins can be substantial enough to influence the late inspiral and merger dynamics of the binary. The spins might also influence the lifetime of an angular momentum supported hypermassive neutron star that can form after merger. With more angular momentum we expect the hypermassive star to survive for longer. This can have important consequences for both the gravitational waves emitted by the system and also for the likelihood of a gamma ray burst.
The purpose of this paper is to numerically implement and test a new method for the computation of binary neutron star initial data with arbitrary rotation states. This method is derived from the standard matter equations of perfect fluids together with certain quasiequilibrium assumptions. We assume that there is an approximate helical Killing vector and that Lie derivatives of the metric variables with respect to vanish. We also assume that scalar matter variables such as or have Lie derivatives that vanish with respect to . However, since the Lie derivative of the fluid velocity is nonzero for generic spins, we split the fluid velocity into an irrotational piece (derived from a potential ) and a rotational piece , and assume that only the irrotational piece has a vanishing Lie derivative (see Eq. (9)) with respect to . Furthermore we know that the spin of each star remains approximately constant since the viscosity of the stars is insufficient for tidal coupling Bildsten and Cutler (1992). To incorporate this fact, we use Eqs. (10) and (11) which are based on the assumption that is constant along the star’s motion described by the irrotational velocity piece .
From these assumptions we obtain the elliptic equations (14), (15), (16) and (17) together with the algebraic equation (18). The specific enthalpy in Eq. (18) determines the shape of the star surfaces. Since our elliptic solvers work on a fixed domain decomposition where the star surfaces have to coincide with domain boundaries we solve this mixture of elliptic and algebraic equations by iteration. In each iteration we first solve the elliptic equations for a given and then use the algebraic Eq. (18) to update . The stability of this iterative procedure is improved if we do the following. (A) We typically do not take the , and and coming from solving Eqs. (14), (15), (16) and (17) as our new fields. Rather, we take the average of this solution and the values from the previous iteration step as our new fields. In this way , and and change less from one iteration step to the next. (B) We use the force balance condition as given in Eq. (5) to update and .
For each iteration we also need to specify rotational piece of the fluid velocity. We have found that the choice in Eq. (33) works very well in the sense that after numerically solving all equations it leads to a velocity in the corotating frame of the form as in Eq. (30). In Appendices A and B we show that this form of results in a fluid 4velocity with an expansion and shear that are approximately zero, as we would expect for stars in equilibrium. This means that we have found a simple way to generate initial data for neutron star binaries that ensure that the star are spinning and without differential rotation.
We also compare initial data sequences for irrotational, corotating and fixed rotation binaries. We find that our method yields reasonable results. All sequences approach postNewtonian results for large separations. And the binary sequences with fixed rotation yield higher and than corotating configurations for large separations, and lower and than corotating configurations for small separations.
Acknowledgements.
It is a pleasure to thank Sebastiano Bernuzzi for helpful discussions. This work was supported by NSF grants PHY0855315 and PHY1204334.Appendix A Expansion, shear and rotation of the fluid
If denotes the 4velocity of the fluid we can split its derivatives into
(34) 
where
(35) 
and where the expansion, shear, rotation and acceleration are defined as
(36) 
(37) 
(38) 
and
(39) 
If the 4velocity is of the form
(40) 
where is any scalar function, it immediately follows (from ) that
(41) 
(42) 
and
(43) 
For our purposes it is often convenient to write the 4velocity as
(44) 
where in an inertial frame the helical Killing vector is given by Eq. (22) and . Note that leads to . From
(45) 
we see that can be interpreted as the velocity in the corotating frame.
Appendix B Expansion and shear if
Let us assume that is given by Eq. (30). Let us further assume that is small in the sense that
(48) 
with . We know (e.g. from postNewtonian theory) that the shift is also small. For simplicity we assume that
(49) 
Using and Eq. (13) we then find
(50) 
Recall that the Newtonian potential near a mass is given by
(51) 
where is a second mass at a distance in the direction. Since and we find that
(52) 
where we have used the virial theorem and set .
References
 R. Narayan, B. Paczynski, and T. Piran, Astrophys. J. 395, L83 (1992).
 T. Piran, Phys.Rept. 333, 529 (2000), eprint astroph/9907392.
 T. Piran, Phys.Rept. 314, 575 (1999), eprint astroph/9810256.
 B. Zhang and P. Meszaros, Int.J.Mod.Phys. A19, 2385 (2004), eprint astroph/0311321.
 T. Piran, Rev.Mod.Phys. 76, 1143 (2005), eprint astroph/0405503.
 E. Nakar, Phys.Rept. 442, 166 (2007), eprint astroph/0701748.
 E. Berger, New Astron.Rev. 55, 1 (2011), eprint 1005.1068.
 L. Rezzolla, B. Giacomazzo, L. Baiotti, J. Granot, C. Kouveliotou, and M. A. Aloy, Astrophys.J. 732, L6 (2011), eprint 1101.4298.
 B. Abbott et al. (LIGO Scientific Collaboration), Rept. Prog. Phys. 72, 076901 (2009), eprint 0711.3041.
 LIGO, URL http://www.ligo.caltech.edu.
 F. Acernese et al., J. Opt. A: Pure Appl. Opt. 10, 064009 (2008).
 VIRGO, URL http://www.virgo.infn.it.
 GEO 600, URL http://www.geo600.org.
 B. Schutz, Class. Quantum Grav. 16, A131 (1999).
 P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
 P. C. Peters, Phys. Rev. 136, B1224 (1964).
 A. G. Lyne et al., Science 303, 1153 (2004), eprint astroph/0401086.
 W. Tichy, Phys. Rev. D84, 024041 (2011), eprint 1107.1440.
 S. Oslowski, T. Bulik, D. GondekRosinska, and K. Belczynski, Mon. Not. Roy. Astr. Soc. 413, 461 (2011), eprint 0903.3538.
 P. Kiel, J. Hurley, M. Bailes, and J. Murray, Mon. Not. Roy. Astr. Soc. 388, 393 (2008), eprint 0805.0059.
 T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D80, 064009 (2009a), eprint 0909.0952.
 T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D80, 089901(E) (2009b), eprint 0909.0952.
 S. Bonazzola, E. Gourgoulhon, and J.A. Marck, Phys. Rev. Lett. 82, 892 (1999), eprint grqc/9810072.
 E. Gourgoulhon, P. Grandclement, K. Taniguchi, J.A. Marck, and S. Bonazzola, Phys. Rev. D63, 064029 (2001), eprint grqc/0007028.
 P. Marronetti, G. J. Mathews, and J. R. Wilson, Phys. Rev. D60, 087301 (1999), eprint arXiv:grqc/9906088.
 K. Uryu and Y. Eriguchi, Phys. Rev. D61, 124023 (2000), eprint grqc/9908059.
 P. Marronetti and S. L. Shapiro, Phys. Rev. D68, 104024 (2003), eprint grqc/0306075.
 K. Taniguchi and E. Gourgoulhon, Phys. Rev. D66, 104019 (2002), eprint grqc/0207098.
 K. Taniguchi and E. Gourgoulhon, Phys. Rev. D68, 124025 (2003), eprint grqc/0309045.
 K. Uryu, F. Limousin, J. L. Friedman, E. Gourgoulhon, and M. Shibata, Phys. Rev. Lett. 97, 171101 (2006), eprint grqc/0511136.
 K. Uryu, F. Limousin, J. L. Friedman, E. Gourgoulhon, and M. Shibata, Phys. Rev. D80, 124004 (2009), eprint 0908.0579.
 W. Tichy, Phys. Rev. D74, 084005 (2006), eprint grqc/0609087.
 W. Tichy, Class. Quant. Grav. 26, 175018 (2009a), eprint 0908.0620.
 W. Tichy, Phys. Rev. D80, 104034 (2009b), eprint 0911.0973.
 C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation (W. H. Freeman, San Francisco, 1973).
 G. Schäfer and N. Wex, Physics Lett. A 174, 196 (1993).
 W. Tichy, B. Brügmann, M. Campanelli, and P. Diener, Phys. Rev. D67, 064008 (2003a), grqc/0207011.
 W. Tichy, B. Brügmann, and P. Laguna, Phys. Rev. D68, 064008 (2003b), eprint grqc/0306020.
 W. Tichy and B. Brügmann, Phys. Rev. D 69, 024006 (2004), eprint grqc/0307027.
 L. Bildsten and C. Cutler, Astrophys. J. 400, 175 (1992).