Initial data for binary neutron stars with adjustable eccentricity

Binary neutron stars in circular orbits can be modeled as helically symmetric, i.e., stationary in a rotating frame. This symmetry gives rise to a first integral of the Euler equation, often employed for constructing equilibrium solutions via iteration. For eccentric orbits, however, the lack of helical symmetry has prevented the use of this method, and the numerical relativity community has often resorted to constructing initial data by superimposing boosted spherical stars without solving the Euler equation. The spuriously excited neutron star oscillations seen in evolutions of such data arise because such configurations lack the appropriate tidal deformations and are stationary in a linearly comoving—rather than rotating—frame. We consider eccentric configurations at apoapsis that are instantaneously stationary in a rotating frame. We extend the notion of helical symmetry to eccentric orbits, by approximating the elliptical orbit of each companion as instantaneously circular, using the ellipse’s inscribed circle. The two inscribed helical symmetry vectors give rise to approximate instantaneous first integrals of the Euler equation throughout each companion. We use these integrals as the basis of a self-consistent iteration of the Einstein constraints to construct conformal thin-sandwich initial data for eccentric binaries. We find that the spurious stellar oscillations are reduced by at least an order of magnitude, compared with those found in evolutions of superposed initial data. The tidally induced oscillations, however, are physical and qualitatively similar to earlier evolutions. Finally, we show how to incorporate radial velocity due to radiation reaction in our inscribed helical symmetry vectors, which would allow one to obtain truly noneccentric initial data when our eccentricity parameter $e$ is set to zero.

Figure 1

Comparing equilibrium data constructed with the constant fluid velocity approximation (using the implementation of our method in bam) to data constructed by solving for the exact velocity potential with sgrid. We show the fluid velocity along the $x$ axis for two different values of $d$, the coordinate separation of the stars’ centers. One can see a deviation of $\sim 10\%$ for $d=31.2$, but this deviation decreases to $\sim 1\%$ for $d=67$. We have cut the data at the surface of the star, denoted by the vertical thin dotted lines, since the velocity is only well defined in the star’s interior.

Figure 2

Illustration of the approximation of the orbits by using circles inscribed into the orbital ellipse in a way that their curvature is the same as the one of the ellipse. We show the scaled semimajor axis $a_1=a\frac{m_2}{M}$, semiminor axis $b_1=b\frac{m_2}{M}$, and the radius $r_{{\mathrm{c}}_1}=b_1^2/a_1$ and center $x_{{\mathrm{c}}_1}$ of the inscribed circle as well as the center $x_1$ of one star. The center of mass is denoted by $x_{\mathrm{cm}}$.

Figure 3

The $y$ component of the momentum constraint (${\mathcal{D}}^y$) for the ecc0 setup. We plot ${\mathcal{D}}^y$ along the $x$ axis (which passes through the centers of both stars) for three grid spacings of $\mathrm{\Delta }x=0.1875$, $2\mathrm{\Delta }x$, and $\mathrm{\Delta }x/2$ in the finest box. The constraint violations were computed with second-order finite differencing, which represents the accuracy of the multigrid algorithm, so we scale the two finer resolutions as appropriate for second-order convergence. Note that the feature on the right side of the plot is due to inaccuracies at the surface of the star, leading to the spikes that we can also see in the eccentric case (see Sec. 7a) but do not show here to focus on the convergence in the strong-field interior of the star.

Figure 4

Comparison of the oscillations of the star, measured using the maximum density at each time step, ${\rho }_{\max }$, normalized by the maximum density at $t=0$, ${\rho }_{\max ,0}$. We show evolutions of sgrid data (solid black line) and the corresponding data set constructed with our method assuming stationarity in a rotating frame (ecc0, blue dashed line). One can clearly see the improvement over the strongly oscillating curves of superimposed boosted spherical stars (red dot-dashed line) and data computed using stationarity in a linearly comoving frame (ecc0v, green dotted line). Note that the latter data set was only evolved for a short time, since we were only interested in the spurious oscillations.

Figure 5

The initial constraint violations for the ecc0.915 data in the finest box (surrounding one star) along the $x$ axis, which passes through both stars’ centers. We show the Hamiltonian constraint ($\mathcal{H}$) and the $y$ component of the momentum constraint (${\mathcal{D}}^y$) for two resolutions, with grid spacings of 0.25 and 0.125 in the finest box and demonstrate second-order convergence of the constraints computed with sixth-order finite differencing stencils, the same order used in the evolution of this data set we show later. (Recall that the initial data code is only second order accurate.)

Figure 6

The initial oscillations of a star in a binary with $e=0.915$, illustrated by considering the maximum density at each time step, normalized by the initial maximum density. We show this for two resolutions with finest grid spacings of 0.25 and 0.125, respectively (i.e., ecc0.915 and ecc0.915high) to demonstrate how the density oscillations decrease with increasing resolution at approximately second order. The notation ${\rho }_{0.125}\times 2^2$ denotes that the oscillations (not the total maximum density) are multiplied by $2^2$, as discussed in the text.

Figure 7

Sequences for equal-mass binary neutron stars with varying eccentricity $e$ (data sets seq0–seq0.9 in Table 1). These sequences are computed with fixed baryonic masses, yielding isolated stars with gravitational masses of $M_i=1.5149$; we define $M=2M_i$. From top to bottom the quantities shown are the rotation $\omega$ as a function of the coordinate separation $d$ of the two stars’ centers, and the ADM angular momentum $J_{\mathrm{ADM}}$ and binding energy $E_b=M_{\mathrm{ADM}}-M$ as functions of the normalized mean motion $M\overline{\mathrm{\Omega }}$. The angular momentum has been normalized by $M^2$, while the binding energy is normalized by $M$. In addition, in the lower two plots, the expected Newtonian behavior is plotted in dashed lines (black line instead of color scheme for the bottom plot, since the Newtonian prediction is independent of $e$) and the post-Newtonian (3PN) results in dotted lines for $e=0$.

Figure 8

The trajectories of one star for different values of the eccentricity parameter $e$. The evolutions are based on the initial data sets ecc0.45, ecc0.5, ecc0.73, and ecc0.915 (see Table 1) which are identical, except for the value of the eccentricity parameter $e$. While larger values of $e$ lead to a rapid merger, for smaller values of the eccentricity (such as $e=0.45$), one can obtain one or more encounters before merger, as found by Gold et al. [20].

Figure 9

The gravitational waveform (in the form of the $l=m=2$ mode of the Newman-Penrose scalar ${\mathrm{\Psi }}_4$) for the ecc0.45 case (see Table 1). We extracted the waveform at a distance $r=500$ from the binary’s center of mass, and shift the time axis by $r=500$ to account (approximately) for the waves’ travel time. The two small bursts (at $t-r\simeq 900$ and 2100) correspond to close encounters before merger; the tidally induced $f$-mode oscillations of the stars are visible between the bursts.

Figure 10

The relation of the eccentricity parameter $e_{\mathrm{ID}}$ input into the code (i.e., $e$ in previous plots) and the “output” eccentricity of the evolutions $e_{\text{evol}}$ measured using two different fits to the initial portion of the trajectory, computed for the ecc0.45–ecc0.96 data sets in Table 1. The fit to an ellipse using the coordinate distance $d_{\text{coord}}$ gives quite similar results to the fit using the proper distance $d_{\text{prop}}$ and the trajectory’s position angle. The black dot-dashed line shows the “ideal” relation $e_{\mathrm{ID}}=e_{\text{evol}}$. The $d_{\text{coord}}$ and $d_{\text{prop}}$ curves agree well enough that one can use them to obtain an estimate of the eccentricity of the system being evolved.