Compact star binaries

5.5 Compact star binaries

5.5.1 Formalism

Systems consisting of two compact objects are known to emit gravitational waves. Due to this emission, no closed orbits can exist and the objects follow a spiral-like trajectory. This implies that such systems have no symmetries that can be taken into account and full time evolutions should be made. However, when the objects are relatively far apart, the emission of gravitational waves is small and the inspiral slow. In this regime, one can hope to approximate the real trajectory with a sequence of closed orbits. Moreover, the emission of gravitational waves is known to efficiently circularize eccentric orbits so that only circular orbits are usually considered.

So, a lot of effort has been devoted to the computation of circular orbits in general relativity. This can be done by demanding that the system admit a helical Killing vector $∂t + Ω ∂φ$ , $Ω$ being the orbital angular velocity of the system. Roughly speaking, this means that advancing in time is equivalent to turning the system around its axis. Working in the corotating frame, one is left with a time-independent problem.

Additional approximations must be made in order to avoid any diverging quantities. Indeed, when using helical symmetry, the system has an infinite lifetime and can fill the whole space with gravitational waves, thus causing quantities like the total mass to be infinite. To deal with that, various techniques can be used, the simplest one consisting of preventing the appearance of any gravitational waves. This is usually done by demanding that the spatial metric be conformally flat. This is not a choice of coordinates but a true approximation that has a priori no reason to be verified. Indeed, even for a single rotating black hole, one can not find coordinates in which the spatial three-metric is conformally flat, so that we do not expect it to be the case for binary systems. However, comparisons with post-Newtonian results or non–conformally-flat results tend to indicate that this approximation is relatively good.

Under these assumptions, Einstein’s equations reduce to a set of five elliptic equations for the lapse, the conformal factor and the shift vector. These equations encompass both the Hamiltonian and momentum constraint equations and the trace of the evolution equations. To close the system, one must provide a description of the matter. It is commonly admitted that the fluid is irrotational, the viscosity in neutron stars being too small to synchronize the motion of the fluid with the orbital motion. It follows that the motion of the fluid is described by an additional elliptic equation for the potential of the flow. The matter terms entering the equations via the stress-energy tensor can then be computed, once an equation of state is given. An evolutionary sequence can be obtained by varying the separation between the stars.

5.5.2 Numerical procedure

Up to now, only the Meudon group has solved these equations by means of spectral methods in the case of two neutron stars. Two sets of domains are used, one centered on each star. Each set consists of sphere-like domains that match the surface of the star and extend up to infinity. Functions are expanded in terms of spherical harmonics with respect to the angles $(𝜃,φ )$ and Chebyshev polynomials with respect to the radial coordinates. Each Poisson equation $ΔN = SN$ is split into two parts $ΔN1 = SN1$ and $ΔN2 = SN2$ , such that $SN = SN1 + SN2$ and $N = N1 + N2$ . The splitting is, of course, not unique and only requires that $SNi$ be mainly centered around star $i$ so that it is well described by spherical coordinates around it. The equation labeled $i$ is then solved using domains centered on the appropriate star. The splittings used for the various equations can be found explicitly in Section IV-C of [100 ].

The elliptic equations are solved using the standard approach by the Meudon group found in [109 ]. For each spherical harmonic, the equation is solved using the tau method and the matching between the various domains is made using the homogeneous solutions method (see Section 2.6.4). The whole system of equations is solved by iteration and most of the computational time is spent when quantities are passed from one set of domains to the other by means of a spectral summation (this requires $N 6$ operations, $N$ being the number of collocation points in one dimension). A complete and precise description of the overall procedure can be found in [100 ].

5.5.3 Neutron star binaries

The first sequence of irrotational neutron star binaries computed by spectral means is shown in [40 ]. Both stars are assumed to be polytropes with an index $γ = 2$ . The results are in good agreement with those obtained simultaneously by other groups with other numerical techniques (see, for instance, [24 , 222 ]). One of the important points that has been clarified by [40 ] concerns the evolution of the central density of the stars. Indeed, at the end of the 1990s, there was a claim that both stars could individually collapse to black holes before coalescence, due to the increase of central density as the two objects spiral towards each other. Should that have been true, this would have had a great impact on the emitted gravitational wave signal. However, it turned out that this came from a mistake in the computation of one of the matter terms. The correct behavior, confirmed by various groups and in particular by [40 ], is a decrease in central density as the stars get closer and closer (see Figure I of [40 ]).

If the first sequence computed by spectral methods is shown in [40 ], the complete description and validation of the method are given in [100 ]. Convergence of the results with respect to the number of collocation points is exhibited. Global quantities like the total energy or angular momentum are plotted as a function of the separation and show remarkable agreement with results coming from analytical methods (see Figures 8 – 15 of [100 ]). Relativistic configurations are also shown to converge to their Newtonian counterparts when the compactness of the stars is small (see Figures 16 – 20 of [100 ]).

Newtonian configurations of compact stars with various equations of state are computed for both equal masses [216 ] and various mass ratios [213 ]. One of the main results of the computations concerns the nature of the end point of the sequence. For equal masses, the sequence ends at contact for synchronized binaries and at mass shedding for irrotational configurations. This is to be contrasted with the non–equal-mass case, in which the sequence always ends at the mass shedding limit of the smallest object.

Properties of the sequences in the relativistic regime are discussed in [214 , 215 ]. In [214 ] sequences with $γ = 2$ are computed for various compactness and mass ratios for both synchronized and irrotational binaries. The nature of the end point of the sequences is discussed and similar behavior to the Newtonian regime is observed. The existence of a configuration of maximum binding energy is also discussed. Such existence could have observational implications because it could be an indication of the onset of a dynamic instability. Sequences of polytropes with various indexes ranging from 1.8 to 2.2 are discussed in [215 ]. In particular, the authors are lead to conjecture that, if a configuration of maximum binding energy is observed in the Newtonian regime, it will also be observed in conformal relativity for the same set of parameters.

In [76 ] the authors derive from the sequences computed in [214 ] a method to constrain the compactness of the stars from the observations. Indeed, from the results in [214 ] one can easily determine the energy emitted in gravitational waves per interval of frequency (i.e. the power spectrum of the signal). For large separation, that is, for small frequencies, the curves follow the Newtonian. However, there is a break frequency at the higher end (see Figure 2 of [76 ]). The location of this frequency depends mainly on the compactness of the stars. More precisely, the more compact the stars are, the higher the break frequency is. Should such frequency be observed by the gravitational wave detectors, this could help to put constraints on the compactness of the neutron stars and, thus, on the equation of state of such objects.

5.5.4 Extensions

The framework of [100 ] is applied to more realistic neutron stars in [26 ]. In this work, the equations of state are more realistic than simple polytropes. Indeed, three different equations are considered for the interior, all based on different microscopic models. The crust is also modeled. For all the models, the end point of the evolution seems to be given by the mass shedding limit. However, the frequency at which the shedding occurs depends rather strongly on the EOS. The results are in good agreement with post-Newtonian ones, until hydrodynamic effects begin to be dominant. This occurs at frequencies in the range of 500 – 1000 Hz, depending on the EOS.

Sequences of strange star binaries have also been computed [133 ]. Contrary to the neutron star case, matter density does not vanish at the surface of the stars and one really needs to use surface-fitting domains to avoid any Gibbs phenomenon that would spoil the convergence of the overall procedure. Sequences are computed for both synchronized and irrotational binaries and a configuration of maximum binding energy is attained in both cases. This is not surprising: strange stars are more compact than neutron stars and are less likely to be tidally destroyed before reaching the extremum of energy, making it easier to attain dynamic instability. More detailed results on both neutron star and strange star binaries are discussed in [87 , 91 ].

All the work presented above was done in the conformal flatness approximation. As already stated in Section 5.5.1, this is only an approximation and one expects that the true conformal three-metric will depart from flatness. However, in order to maintain asymptotic flatness of spacetime, one needs to get rid of the gravitational wave content. One such waveless approximation is presented in [196 ] and implemented in [223 ]. Two independent codes are used, one of them being an extension of the work described in [100 ]. The number of equations to be solved is then greater than in conformal flatness (one has to solve for the conformal metric), but the algorithms are essentially the same. It turns out that the deviation from conformal flatness is rather small. The new configurations are slightly further from post-Newtonian results than the conformally-flat ones, which is rather counter-intuitive and might be linked to a difference in the definition of the waveless approximations.