Constructing quasi-equilibrium initial data for binary neutron stars with arbitrary spins
In general neutron stars in binaries are spinning. Recently, a new quasi-equilibrium approximation that includes a rotational velocity piece for each star has been proposed to describe binary neutron stars with arbitrary rotation states in quasi-circular 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 4-velocities that have expansion and shear of approximately zero, as one would expect for quasi-equilibrium 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.Fk
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 post-Newtonian 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 non-linear 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 quasi-circular 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 J0737-3039 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 quasi-circular 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 4-velocities 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 Arnowitt-Deser-Misner (ADM) decomposition of Einstein’s equations and describe the gravitational fields (i.e. the 4-metric ) in terms of the 3-metric , lapse , shift and the extrinsic curvature . We further assume that the neutron star matter is a perfect fluid with stress-energy tensor
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 4-velocity of the fluid. Assuming a polytropic equation of state
and defining the specific enthalpy
, and can all be expressed in terms of . We find
where we have used the abbreviation
The fluid 4-velocity is expressed in terms of the 3-velocity
which in turn is split into a irrotational piece and a rotational piece
where is the derivative operator compatible with the 3-metric .
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
We also assume similar equations for scalar matter quantities such as . For a spinning star, however, is non-zero. Instead we assume that
so that the time derivative of the irrotational piece of the fluid velocity vanishes in corotating coordinates. We also assume that
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
and conformal flatness
yield the following coupled equations:
Here , , and we have introduced
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
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.
at spatial infinity, and
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 3-velocity via Eq. (7). The 3-metric is obtained from Eq. (13) and the extrinsic curvature is given by
Notice that Eq. (18) can also be written as
where we have introduced
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 non-vanishing 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 non-vanishing 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:
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 Tolman-Oppenheimer-Volkoff 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 post-Newtonian 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
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
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
It is clear that
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 non-zero shear. Another simple choice is
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 post-Newtonian results (taken from Schäfer and Wex (1993); Tichy et al. (2003a, b); Tichy and Brügmann (2004)) for non-spinning 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 post-Newtonian non-spinning 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 .
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 quasi-equilibrium 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 non-zero 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 4-velocity 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 post-Newtonian 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 PHY-0855315 and PHY-1204334.
Appendix A Expansion, shear and rotation of the fluid
If denotes the 4-velocity of the fluid we can split its derivatives into
and where the expansion, shear, rotation and acceleration are defined as
If the 4-velocity is of the form
where is any scalar function, it immediately follows (from ) that
For our purposes it is often convenient to write the 4-velocity as
where in an inertial frame the helical Killing vector is given by Eq. (22) and . Note that leads to . From
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
with . We know (e.g. from post-Newtonian theory) that the shift is also small. For simplicity we assume that
Using and Eq. (13) we then find
Recall that the Newtonian potential near a mass is given by
where is a second mass at a distance in the -direction. Since and we find that
where we have used the virial theorem and set .
- R. Narayan, B. Paczynski, and T. Piran, Astrophys. J. 395, L83 (1992).
- T. Piran, Phys.Rept. 333, 529 (2000), eprint astro-ph/9907392.
- T. Piran, Phys.Rept. 314, 575 (1999), eprint astro-ph/9810256.
- B. Zhang and P. Meszaros, Int.J.Mod.Phys. A19, 2385 (2004), eprint astro-ph/0311321.
- T. Piran, Rev.Mod.Phys. 76, 1143 (2005), eprint astro-ph/0405503.
- E. Nakar, Phys.Rept. 442, 166 (2007), eprint astro-ph/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 astro-ph/0401086.
- W. Tichy, Phys. Rev. D84, 024041 (2011), eprint 1107.1440.
- S. Oslowski, T. Bulik, D. Gondek-Rosinska, 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 gr-qc/9810072.
- E. Gourgoulhon, P. Grandclement, K. Taniguchi, J.-A. Marck, and S. Bonazzola, Phys. Rev. D63, 064029 (2001), eprint gr-qc/0007028.
- P. Marronetti, G. J. Mathews, and J. R. Wilson, Phys. Rev. D60, 087301 (1999), eprint arXiv:gr-qc/9906088.
- K. Uryu and Y. Eriguchi, Phys. Rev. D61, 124023 (2000), eprint gr-qc/9908059.
- P. Marronetti and S. L. Shapiro, Phys. Rev. D68, 104024 (2003), eprint gr-qc/0306075.
- K. Taniguchi and E. Gourgoulhon, Phys. Rev. D66, 104019 (2002), eprint gr-qc/0207098.
- K. Taniguchi and E. Gourgoulhon, Phys. Rev. D68, 124025 (2003), eprint gr-qc/0309045.
- K. Uryu, F. Limousin, J. L. Friedman, E. Gourgoulhon, and M. Shibata, Phys. Rev. Lett. 97, 171101 (2006), eprint gr-qc/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 gr-qc/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), gr-qc/0207011.
- W. Tichy, B. Brügmann, and P. Laguna, Phys. Rev. D68, 064008 (2003b), eprint gr-qc/0306020.
- W. Tichy and B. Brügmann, Phys. Rev. D 69, 024006 (2004), eprint gr-qc/0307027.
- L. Bildsten and C. Cutler, Astrophys. J. 400, 175 (1992).