# Analysis of the Conformally Flat Approximation for Binary Neutron Star Initial Conditions

Authors: In-Saeng Suh; Grant J. Mathews; J. Reese Haywood; N. Q. Lan
Journal: Advances in Astronomy
 (2017)

---

## Abstract

The spatially conformally flat approximation (CFA) is a viable method to deduce initial conditions for the subsequent evolution of binary neutron stars employing the full Einstein equations. Here we analyze the viability of the CFA for the general relativistic hydrodynamic initial conditions of binary neutron stars. We illustrate the stability of the conformally flat condition on the hydrodynamics by numerically evolving ~100 quasicircular orbits. We illustrate the use of this approximation for orbiting neutron stars in the quasicircular orbit approximation to demonstrate the equation of state dependence of these initial conditions and how they might affect the emergent gravitational wave frequency as the stars approach the innermost stable circular orbit.


---

## Body

## 1. Introduction

The epoch of gravitational wave astronomy has now begun with the first detection [1, 2] of the merger of binary black holes by Advanced LIGO [3]. Now that the first ground based gravitational wave detection has been achieved, observations of binary neutron star mergers should soon be forthcoming. This is particularly true as other second generation observatories such as Advanced VIRGO [4] and KAGRA [5] will soon be online. In addition to binary black holes, neutron star binaries are thought to be among the best candidate sources gravitational radiation [6, 7]. The number of such systems detectable by Advanced LIGO is estimated [7–14] to be of order several events per year based upon observed close binary-pulsar systems [15, 16]. There is a difference between neutron star mergers and black hole mergers; however, in that neutron star mergers involve the complex evolution of the matter hydrodynamic equations in addition to the strong gravitational field equations. Hence, one must carefully consider both the hydrodynamic and field evolution of these systems.To date there have been numerous attempts to calculate theoretical templates for gravitational waves from compact binaries based upon numerical and/or analytic approaches (see, e.g., [17–26]). However, most approaches utilize a combination of post-Newtonian (PN) techniques supplemented with quasicircular orbit calculations and then applying full GR for only the last few orbits before disruption. In this paper we analyze the hydrodynamic evolution in the spatially conformally flat metric approximation (CFA) as a means to compute stable initial conditions beyond the range of validity of the PN regime, that is, near the last stable orbits. We establish the numerical stability of this approach based upon many orbit simulations of quasicircular orbits. We show that one must follow the stars for several orbits before a stable quasicircular orbit can be achieved. We also illustrate the equation of state (EoS) dependence of the initial conditions and associated gravitational wave emission.When binary neutron stars are well separated, the post-Newtonian (PN) approximation is sufficiently accurate [27]. In the PN scheme, the stars are often treated as point masses, either with or without spin. At third order, for example, it has been estimated [28–30] that the error due to assuming the stars are point masses is less than one orbital rotation [28] over the ~16,000 cycles that pass through the LIGO detector frequency band [7]. Nevertheless, it has been noted in many works [25, 31–42] that relativistic hydrodynamic effects might be evident even at the separation (~10–100km) relevant to the LIGO window.Indeed, the templates generated by PN approximations, unless carried out to fifth and sixth order [28, 29], may not be accurate unless the finite size and proper fluid motion of the stars are taken into account. In essence, the signal emitted during the last phases of inspiral depends upon the finite size and the equation of state (EoS) through the tidal deformation of the neutron stars and the cut-off frequency when tidal disruption occurs.Numeric and analytic simulations [43–51] of binary neutron stars have explored the approach to the innermost stable circular orbit (ISCO). While these simulations represent some of the most accurate to date, simulations generally follow the evolution for a handful of orbits and are based upon initial conditions of quasicircular orbits obtained in the conformally flat approximation. Accurate templates of gravitational radiation require the ability to stably and reliably calculate the orbit initial conditions. The CFA provides a means to obtain accurate initial conditions near the ISCO.The spatially conformally flat approximation to GR was first developed in detail in [32]. That original formulation, however, contained a mathematical error first pointed out by Flanagan [52] and subsequently corrected in [34]. This error in the solution to the shift vector led to a spurious NS crushing prior to merger. The formalism discussed below is for the corrected equations. Here, we discuss the hydrodynamic solutions as developed in [31–34, 53, 54]. This CFA formalism includes much of the nonlinearity inherent in GR and leads set of coupled, nonlinear, elliptic field equations that can be evolved stably. We also note that an alternative spectral method solution to the CFA configurations was developed by [55, 56], and approaches beyond the CFA have also been proposed [48]. However, our purpose here is to clarify the viability of the hydrodynamic solution without the imposition of a Killing vector or special symmetry. This approach is the most adaptable, for example, to general initial conditions such as that of arbitrarily elliptical orbits and/or arbitrarily spinning neutron stars.Here, we summarize the original CFA approach and associated general relativistic hydrodynamics formalism developed in [32, 34, 53, 54] and illustrate that it can produce stable initial conditions anywhere between the post-Newtonian to ISCO regimes. We quantify how long this method takes to converge to quasiequilibrium and demonstrate the stability by subsequently integrating up to ~100 orbits for a binary neutron star system. We also analyze the EoS dependence of these quasicircular initial orbits and show how these orbits can be used to make preliminary estimates [57] of the gravitational wave signal for the initial conditions.This paper is organized as follows. In Section2 the basic method is summarized and in Section 3 a number of code tests are performed in the quasiequilibrium circular orbit limit to demonstrate the stability of the technique. The EoS dependence of the initial conditions and associated gravitational wave frequency are analyzed in Section 4. Conclusions are presented in Section 5.


## 2. Method

### 2.1. Field Equations

The solution of the field equations and hydrodynamic equations of motion were first solved in three spatial dimensions and explained in detail in the 1990s in [31, 32] and subsequently further reviewed in [53, 58]. Here, we present a brief summary to introduce the variables relevant to the present discussion.One starts with the slicing of space-time into the usual one-parameter family of hypersurfaces separated by differential displacements in a time-like coordinate as defined in the (3+1) ADM formalism [59, 60].In Cartesianx, y, z isotropic coordinates, proper distance is expressed as(1)ds2=-α2-βiβidt2+2βidxidt+ϕ4δijdxidxj,where the lapse function α describes the differential lapse of proper time between two hypersurfaces. The quantity βi is the shift vector denoting the shift in space-like coordinates between hypersurfaces. The curvature of the metric of the 3-geometry is described by a position-dependent conformal factor ϕ4 times a flat-space Kronecker delta (γij=ϕ4δij). This conformally flat condition (together with the maximal slicing gauge, tr⁡Kij=0) requires [60](2)2αKij=Diβj+Djβi-23δijDkβk,where Kij is the extrinsic curvature tensor and Di are 3-space covariant derivatives. This conformally flat condition on the metric provides a numerically valid initial solution to the Einstein equations. The vanishing of the Weyl tensor for a stationary system in three spatial dimensions guarantees that a conformally flat solution to the Einstein equations exists.One consequence of this conformally flat approximation to the three-metric is that the emission of gravitational radiation is not explicitly evolved. Nevertheless, one can extract the gravitational radiation signal and the back reaction via a multipole expansion [32, 61]. An application to the determination of the gravitational wave emission from the quasicircular orbits computed here is given in [57]. The advantage of this approximation is that conformal flatness stabilizes and simplifies the solution to the field equations.As a third gauge condition, one can choose separate coordinate transformations for the shift vector and the hydrodynamic grid velocity to separately minimize the field and matter motion with respect to the coordinates. This set of gauge conditions is key to the present application. It allows one to stably evolve up to hundreds and even thousands of binary orbits without the numerical error associated with the frequent advocating of fluid through the grid.Given a distribution of mass and momentum on some manifold, then one first solves the constraint equations of general relativity at each time for a fixed distribution of matter. One then evolves the hydrodynamic equations to the next time step. Thus, at each time slice a solution to the relativistic field equations and information on the hydrodynamic evolution is obtained.The solutions for the field variablesϕ, α, and βi reduce to simple Poisson-like equations in flat space. The Hamiltonian constraint [60] is used to solve for the conformal factor ϕ [32, 62](3)∇2ϕ=-2πϕ5W2ρ1+ϵ+P-P+116πKijKij.In the Newtonian limit, the RHS is dominated [32] by the proper matter density ρ, but in strong fields and compact neutron stars there are also contributions from the internal energy density ϵ, pressure P, and extrinsic curvature. The source is also significantly enhanced by the generalized curved-space Lorentz factor W(4)W=αUt=1+∑Ui2ϕ41/2,where Ut is the time component of the relativistic four velocity and Ui are the covariant spatial components. This factor, W, becomes important near the last stable orbit as the specific kinetic energy of the stars rapidly increases.In a similar manner [32, 62], the Hamiltonian constraint, together with the maximal slicing condition, provides an equation for the lapse function,(5)∇2αϕ=2παϕ53W2ρ1+ϵ+P-2ρ1+ϵ+3P+716πKijKij.Finally, the momentum constraints yields [60] an elliptic equation for the shift vector [34, 52],(6)∇2βi=∂∂xi13∇·β+4πρ3i,where(7)ρ3i=4αϕ4Si+14π∂ln⁡α/ϕ6∂xj∂βi∂xj+∂βj∂xi-23δij∂βk∂xk.Here Si are the spatial components of the momentum density one-form as defined below.We note that, in early applications of this approach, the source for the shift vector contained a spurious term due to an incorrect transformation between contravariant and covariant forms of the momentum density as was pointed out in [34, 52]. As illustrated in those papers, this was the main reason why early hydrodynamic calculations induced a controversial additional compression on stars causing them to collapse to black holes prior to inspiral [31]. This problem no longer exists in the formulation summarized here.


### 2.2. Relativistic Hydrodynamics

To solve for the fluid motion of the system in curved-space time it is convenient to use an Eulerian fluid description [63]. One begins with the perfect fluid stress-energy tensor in the Eulerian observer rest frame,(8)Tμν=Pgμν+ρ1+ϵ+PUμUν, where Uμ is the relativistic four velocity one-form.By introducing the usual set of Lorentz contracted state variables it is possible to write the relativistic hydrodynamic equations in a form which is reminiscent of their Newtonian counterparts [63]. The hydrodynamic state variables are the coordinate baryon mass density(9)D=Wρ; the coordinate internal energy density(10)E=Wρϵ; the spatial three velocity(11)Vi=αUiϕ4W-βi;and the covariant momentum density(12)Si=D+E+PWUi.In terms of these state variables, the hydrodynamic equations in the CFA are as follows: the equation for the conservation of baryon number takes the form(13)∂D∂t=-6D∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6DVj. The equation for internal energy evolution becomes(14)∂E∂t=-6E+PW∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6EVj-P∂W∂t+1ϕ6∂∂xjϕ6WVj. Momentum conservation takes the form(15)∂Si∂t=-6Si∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6SiVj-α∂P∂xi+2αD+E+PWW-1W∂log⁡ϕ∂xi+Sj∂βj∂xi-WD+E+PW∂α∂xi-αWD+ΓE∂χ∂xi,where the last term in (15) is the contribution from the radiation reaction potential χ as defined in [32, 57]. In the construction of quasistable orbit initial conditions, this term is set to zero. Including this term would allow for a calculation of the orbital evolution via gravitational wave emission in the CFA. However, there is no guarantee that this is a sufficiently accurate solution to the exact Einstein equations. Hence, the CFA is useful for the construction of initial conditions.


### 2.3. Angular Momentum and Orbital Frequency

In the quasicircular orbit approximation (neglecting angular momentum in the radiation field), this system has a Killing vector corresponding to rotation in the orbital plane. Hence, for these calculations the angular momentum is well defined and given by an integral over the space-time components of the stress-energy tensor [64]; that is,(16)Jij=∫Ti0xj-Tj0xidV. Aligning the z-axis with the angular momentum vector then gives(17)J=∫xSy-ySxdV.To find the orbital frequency detected by a distant observer corresponding to a fixed angular momentum we employ a slightly modified derivation of the orbital frequency compared to that of [53]. In asymptotically flat coordinates the angular frequency detected by a distant observer is simply the coordinate angular velocity; that is, one evaluates(18)ω≡dϕdt=UϕU0.In the ADM conformally flat (3+1) curved space, our only task is then to deduce Uϕ from code coordinates. For this we make a simple polar coordinate transformation keeping our conformally flat coordinates, so(19)Uϕ=ΛνϕUν=xUy-yUxx2+y2. Now, the code uses covariant four velocities, Ui=giνUν=βiU0+ϕ4Ui. This gives Ui=Uiβi(W/α)/ϕ4. Finally, one must density weight and volume average ω over the fluid differential volume elements. This gives(20)ω=∫d3xϕ2D+ΓEα/WxUy-yUx-xβy-yβx/x2+y2∫d3xϕ6D+ΓE.This form differs slightly from that of [53] but leads to the similar results.A key additional ingredient, however, is the implementation of a grid three velocityVGi that minimizes the matter motion with respect to Ui and βi. Hence, the total angular frequency to a distant observer ωtot=ω+ωG, and in the limit of rigid corotation, ωtot→ωG, where ωG=xVy+yVx.For the orbit calculations illustrated here we model corotating stars, that is, no spin in the corotating frame. This minimizes matter motion on the grid. However, we note that there is need at the present time of initial conditions for arbitrarily spinning neutron stars and the method described here is equally capable of supplying those initial conditions.As a practical approach the simulation [32] of initial conditions is best run first with viscous damping in the hydrodynamics for sufficiently long time (a few thousand cycles) to relax the stars to a steady state. Subsequently, one can run with no damping. In the present illustration we examine stars at large separation which are in quasiequilibrium circular orbits and stable hydrodynamic configurations. In the initial relaxation phase the orbits are circularized by damping any radial velocity components. During the evolution, the rigorous conservation of angular momentum is imposed by adjusting the orbital angular velocity in (20) such that (17) remains constant. The simulated orbits described in this work span the time from the last several minutes up to orbit inspiral. Here, we illustrate the stability of the multiple orbit hydrodynamic simulation and examine where the initial conditions for the strong field orbit dynamics computed here deviates from the post-Newtonian regime.


## 2.1. Field Equations

The solution of the field equations and hydrodynamic equations of motion were first solved in three spatial dimensions and explained in detail in the 1990s in [31, 32] and subsequently further reviewed in [53, 58]. Here, we present a brief summary to introduce the variables relevant to the present discussion.One starts with the slicing of space-time into the usual one-parameter family of hypersurfaces separated by differential displacements in a time-like coordinate as defined in the (3+1) ADM formalism [59, 60].In Cartesianx, y, z isotropic coordinates, proper distance is expressed as(1)ds2=-α2-βiβidt2+2βidxidt+ϕ4δijdxidxj,where the lapse function α describes the differential lapse of proper time between two hypersurfaces. The quantity βi is the shift vector denoting the shift in space-like coordinates between hypersurfaces. The curvature of the metric of the 3-geometry is described by a position-dependent conformal factor ϕ4 times a flat-space Kronecker delta (γij=ϕ4δij). This conformally flat condition (together with the maximal slicing gauge, tr⁡Kij=0) requires [60](2)2αKij=Diβj+Djβi-23δijDkβk,where Kij is the extrinsic curvature tensor and Di are 3-space covariant derivatives. This conformally flat condition on the metric provides a numerically valid initial solution to the Einstein equations. The vanishing of the Weyl tensor for a stationary system in three spatial dimensions guarantees that a conformally flat solution to the Einstein equations exists.One consequence of this conformally flat approximation to the three-metric is that the emission of gravitational radiation is not explicitly evolved. Nevertheless, one can extract the gravitational radiation signal and the back reaction via a multipole expansion [32, 61]. An application to the determination of the gravitational wave emission from the quasicircular orbits computed here is given in [57]. The advantage of this approximation is that conformal flatness stabilizes and simplifies the solution to the field equations.As a third gauge condition, one can choose separate coordinate transformations for the shift vector and the hydrodynamic grid velocity to separately minimize the field and matter motion with respect to the coordinates. This set of gauge conditions is key to the present application. It allows one to stably evolve up to hundreds and even thousands of binary orbits without the numerical error associated with the frequent advocating of fluid through the grid.Given a distribution of mass and momentum on some manifold, then one first solves the constraint equations of general relativity at each time for a fixed distribution of matter. One then evolves the hydrodynamic equations to the next time step. Thus, at each time slice a solution to the relativistic field equations and information on the hydrodynamic evolution is obtained.The solutions for the field variablesϕ, α, and βi reduce to simple Poisson-like equations in flat space. The Hamiltonian constraint [60] is used to solve for the conformal factor ϕ [32, 62](3)∇2ϕ=-2πϕ5W2ρ1+ϵ+P-P+116πKijKij.In the Newtonian limit, the RHS is dominated [32] by the proper matter density ρ, but in strong fields and compact neutron stars there are also contributions from the internal energy density ϵ, pressure P, and extrinsic curvature. The source is also significantly enhanced by the generalized curved-space Lorentz factor W(4)W=αUt=1+∑Ui2ϕ41/2,where Ut is the time component of the relativistic four velocity and Ui are the covariant spatial components. This factor, W, becomes important near the last stable orbit as the specific kinetic energy of the stars rapidly increases.In a similar manner [32, 62], the Hamiltonian constraint, together with the maximal slicing condition, provides an equation for the lapse function,(5)∇2αϕ=2παϕ53W2ρ1+ϵ+P-2ρ1+ϵ+3P+716πKijKij.Finally, the momentum constraints yields [60] an elliptic equation for the shift vector [34, 52],(6)∇2βi=∂∂xi13∇·β+4πρ3i,where(7)ρ3i=4αϕ4Si+14π∂ln⁡α/ϕ6∂xj∂βi∂xj+∂βj∂xi-23δij∂βk∂xk.Here Si are the spatial components of the momentum density one-form as defined below.We note that, in early applications of this approach, the source for the shift vector contained a spurious term due to an incorrect transformation between contravariant and covariant forms of the momentum density as was pointed out in [34, 52]. As illustrated in those papers, this was the main reason why early hydrodynamic calculations induced a controversial additional compression on stars causing them to collapse to black holes prior to inspiral [31]. This problem no longer exists in the formulation summarized here.


## 2.2. Relativistic Hydrodynamics

To solve for the fluid motion of the system in curved-space time it is convenient to use an Eulerian fluid description [63]. One begins with the perfect fluid stress-energy tensor in the Eulerian observer rest frame,(8)Tμν=Pgμν+ρ1+ϵ+PUμUν, where Uμ is the relativistic four velocity one-form.By introducing the usual set of Lorentz contracted state variables it is possible to write the relativistic hydrodynamic equations in a form which is reminiscent of their Newtonian counterparts [63]. The hydrodynamic state variables are the coordinate baryon mass density(9)D=Wρ; the coordinate internal energy density(10)E=Wρϵ; the spatial three velocity(11)Vi=αUiϕ4W-βi;and the covariant momentum density(12)Si=D+E+PWUi.In terms of these state variables, the hydrodynamic equations in the CFA are as follows: the equation for the conservation of baryon number takes the form(13)∂D∂t=-6D∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6DVj. The equation for internal energy evolution becomes(14)∂E∂t=-6E+PW∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6EVj-P∂W∂t+1ϕ6∂∂xjϕ6WVj. Momentum conservation takes the form(15)∂Si∂t=-6Si∂log⁡ϕ∂t-1ϕ6∂∂xjϕ6SiVj-α∂P∂xi+2αD+E+PWW-1W∂log⁡ϕ∂xi+Sj∂βj∂xi-WD+E+PW∂α∂xi-αWD+ΓE∂χ∂xi,where the last term in (15) is the contribution from the radiation reaction potential χ as defined in [32, 57]. In the construction of quasistable orbit initial conditions, this term is set to zero. Including this term would allow for a calculation of the orbital evolution via gravitational wave emission in the CFA. However, there is no guarantee that this is a sufficiently accurate solution to the exact Einstein equations. Hence, the CFA is useful for the construction of initial conditions.


## 2.3. Angular Momentum and Orbital Frequency

In the quasicircular orbit approximation (neglecting angular momentum in the radiation field), this system has a Killing vector corresponding to rotation in the orbital plane. Hence, for these calculations the angular momentum is well defined and given by an integral over the space-time components of the stress-energy tensor [64]; that is,(16)Jij=∫Ti0xj-Tj0xidV. Aligning the z-axis with the angular momentum vector then gives(17)J=∫xSy-ySxdV.To find the orbital frequency detected by a distant observer corresponding to a fixed angular momentum we employ a slightly modified derivation of the orbital frequency compared to that of [53]. In asymptotically flat coordinates the angular frequency detected by a distant observer is simply the coordinate angular velocity; that is, one evaluates(18)ω≡dϕdt=UϕU0.In the ADM conformally flat (3+1) curved space, our only task is then to deduce Uϕ from code coordinates. For this we make a simple polar coordinate transformation keeping our conformally flat coordinates, so(19)Uϕ=ΛνϕUν=xUy-yUxx2+y2. Now, the code uses covariant four velocities, Ui=giνUν=βiU0+ϕ4Ui. This gives Ui=Uiβi(W/α)/ϕ4. Finally, one must density weight and volume average ω over the fluid differential volume elements. This gives(20)ω=∫d3xϕ2D+ΓEα/WxUy-yUx-xβy-yβx/x2+y2∫d3xϕ6D+ΓE.This form differs slightly from that of [53] but leads to the similar results.A key additional ingredient, however, is the implementation of a grid three velocityVGi that minimizes the matter motion with respect to Ui and βi. Hence, the total angular frequency to a distant observer ωtot=ω+ωG, and in the limit of rigid corotation, ωtot→ωG, where ωG=xVy+yVx.For the orbit calculations illustrated here we model corotating stars, that is, no spin in the corotating frame. This minimizes matter motion on the grid. However, we note that there is need at the present time of initial conditions for arbitrarily spinning neutron stars and the method described here is equally capable of supplying those initial conditions.As a practical approach the simulation [32] of initial conditions is best run first with viscous damping in the hydrodynamics for sufficiently long time (a few thousand cycles) to relax the stars to a steady state. Subsequently, one can run with no damping. In the present illustration we examine stars at large separation which are in quasiequilibrium circular orbits and stable hydrodynamic configurations. In the initial relaxation phase the orbits are circularized by damping any radial velocity components. During the evolution, the rigorous conservation of angular momentum is imposed by adjusting the orbital angular velocity in (20) such that (17) remains constant. The simulated orbits described in this work span the time from the last several minutes up to orbit inspiral. Here, we illustrate the stability of the multiple orbit hydrodynamic simulation and examine where the initial conditions for the strong field orbit dynamics computed here deviates from the post-Newtonian regime.


## 3. Code Validation

### 3.1. Code Tests

To evolve stars at large separation distance it is best [53] to decompose the grid into a high resolution domain with a fine matter grid around the stars and a coarser domain with an extended grid for the fields. Figure 1 shows a schematic of this decomposition from [54].Figure 1
Schematic representation of the field and hydrodynamics grid used in the simulations. The inner blue grid represents the higher resolution matter grid and the outer white grid represents the field grid.As noted in [53] it is best to keep the number of zones across each star between 25 and 40 [54]. This keeps the error in the numerics below 0.5%. It has also been pointed out [53] that an artificial viscosity (AV) shock capturing scheme has an advantage over Riemann solvers because only about half as many zones are required to accurately resolve the stars when an AV scheme is employed compared to a Riemann solver.The time stepsdtn are taken as the minimum of the time step as determined by several conditions. Each condition is also multiplied by a number less than one to accommodate the nonlinear nature of the equations.The first condition is known as the Courant condition, that is, a search over all zonesi for the zone with the minimum sound crossing time:(21)dt1=min⁡dxbiCsi, where Csi is the sound speed in the ith zone.The Newtonian sound speed is given by the variation of pressure with density. In relativity the wave speed is given instead by the adiabatic derivative of the pressure with respect to the relativistic inertial density. In terms of relativistic variables the local sound speed in zonei then becomes [53](22)Cs=ΓiΓi-1EiDi+ΓiEi.The second condition is a search for the zone with minimum time for material to flow across a zone(23)dt2=min⁡dxaiVi. This constraint is introduced to ensure stability and accuracy in the numerical advection calculation.The third condition is introduced to maintain stability of the artificial viscosity algorithm. The viscous equations are analogous to a diffusion equation in four velocity with a diffusion coefficientD≈k1dxi|δUi|, where δUi≡Ui+1-Ui. We then can define a minimum viscous diffusion time across a zone derived from the stability condition for explicit diffusion equations(24)dt3=14min⁡WidxbiδUi.The time stepdt is then assigned to be some fraction (referred to as the Courant parameter) of the minimum of these three conditions(25)dt=kmin⁡dt1,dt2,dt3,Obviously, smaller values for k increase the accuracy of the calculation but also increase the computation time.Figure2 shows a plot of orbital velocity versus time for various values of the Courant parameter. This figure establishes that the routines for the hydrodynamics are stable (e.g., changing the Courant condition has little to effect) as long as k≤0.5.Figure 2
Comparison of the orbital angular velocityω versus time for different values of the Courant parameter k. As can be seen, the simulations with k = 0.25–0.5 result in stable runs that converge to the same value, implying that a smaller k or equivalently a smaller δt is not necessary and would only use extra CPU time. For comparison, we plot a simulation with k=0.6 to show that the stability is lost for k>0.5.Figure3 illustrates the difference between the central density ρc and the central density (ρ52) at the highest resolution of 52 zones for single isolated stars. This is expressed as the fractional error as a function of the number of zones across a star. This plot was calculated using the relatively soft MW EoS, that is, the zero temperature and zero neutrino chemical potential limit of the EoS that has previously been used to model core-collapse supernovae [32, 53, 65].Figure 3
Plot of the error in the central density versus the number of zones across the star. It is clear that there is only a 1% error with ≈15 zones across the star. Increasing the number of zones across the star so that there are >35 zones across the star produces less than a 0.1% error.This figure illustrates that here is only a 1% error in central density with ≈15 zones across the star, while increasing the number of zones across the star to >35 produces less than a 0.1%. In the illustrations below we maintaink=0.5 and ≈25 zones across each star as the best choice for both speed and accuracy needed to compute stable orbital initial conditions.


### 3.2. Orbit Stability

As an illustration of the orbit stability Figure4 shows results from a simulation [54] in which the angular momentum was fixed at J=2.7×1011 cm2 and the Courant parameter set to k=0.5. For this orbital calculation the MW EoS was employed and each star was fixed at a baryon mass of MB=1.54M⊙ and a gravitational mass in isolation of MG=1.40M⊙.Figure 4
Plot of the orbital angular velocity,ω in geometrized units (cm−1) versus cycle number. When ω stops changing the simulation has reached a circular binary orbit solution. This run was extended to over 30,000 cycles, corresponding to ≈20 orbits.Figure4 shows the evolution of the orbital angular velocity ω versus computational cycle for the first 30,000 code cycles corresponding to ≈20 orbits. The stars were initially placed on the grid using a solution of the TOV equation in isotropic coordinates for an isolated star. The stars were initially set to be corotating but were allowed to settle into their binary equilibrium. Notice that it takes ~5,000 cycles, corresponding to ~3 orbits, just to approach the quasiequilibrium binary solution. Indeed, the stars continued to gradually compact and slightly increase in orbital frequency until ~10 orbits; afterward, the stars were completely stable.Figure5 shows the contours of the lapse function α (roughly corresponding to the gravitational potential) and corresponding density profiles at cycle numbers, 0, 5200, and 25800 (≈0, 5, and 19 orbits). Figure 6 shows the contours of central density and the orientation of the binary orbit corresponding to these cycle numbers. One can visibly see from these figures the relaxation of the stars after the first few orbits and the stability of the density profiles after multiple orbits.Figure 5
Contours of the lapse function (left) and central density (right) at cycle numbers 0 (a), 5,200 (b), and 25,800 (c) corresponding to roughly 0, 5, and 19 orbits.
(a) (b) (c)Figure 6
Contours of the central density for the binary system at the approximate number of orbits as labelled.We note, however, that this orbit is on the edge of the ISCO. As such it could be unstable to inspiral even after many orbits. Figures7 and 8 further illustrate this point. In these simulations various angular momenta were computed with a slightly higher neutron star mass (Mb=1.61M⊙ and Mg=1.44M⊙), but the same MW EoS. In this case the binary system was followed for nearly 100 orbits.Figure 7
Plot of the orbital angular velocity,ω, versus cycle. When ω stops changing with time the simulation has reached a circular binary orbit solution. The run (a for J=2.7×1011 cm2) goes over ~10 obits and then becomes unstable to inspiral and merger after ~104 cycles. The stable two runs (b for J=2.8×1011 cm−2 and c for J=2.9×1011 cm−2) were run for 100,000 cycles and ≈100 orbits.Figure 8
Plot of the central density,ρc, versus the number of orbit. The solid lines from top to bottom are for J=3.0,3.2,3.4,3.6,3.8,4.0×1011 cm2.Figure7 illustrates orbital angular frequency versus cycle number for three representative angular momenta bracketing the ISCO. The orbital separation for the lowest angular momentum (J=2.7×1011 cm−2) shown on Figure 7 is just inside the ISCO. Hence, even though it requires about 10 orbits before inspiral, the orbit is eventually unstable. Similarly, Figure 8 shows the central density versus number of orbits for 6 different angular momenta. Here one can see that orbits with J≥3.0×1011 cm−2 are stable. Indeed, for these cases, after about the first 3 orbits the orbits continue with almost no discernible change in orbit frequency or central density.As mentioned previously, the numerical relativistic neutron binary simulations of [43] all start with initial data that are subsequently evolved in a different manner compared to those with which they were created. One conclusion that may be drawn from the above set of simulations, however, is that the initial data must be evolved for ample time (>3 orbit) for the stars to reach a true quasiequilibrium binary configuration. That has not always been done in the literature.


## 3.1. Code Tests

To evolve stars at large separation distance it is best [53] to decompose the grid into a high resolution domain with a fine matter grid around the stars and a coarser domain with an extended grid for the fields. Figure 1 shows a schematic of this decomposition from [54].Figure 1
Schematic representation of the field and hydrodynamics grid used in the simulations. The inner blue grid represents the higher resolution matter grid and the outer white grid represents the field grid.As noted in [53] it is best to keep the number of zones across each star between 25 and 40 [54]. This keeps the error in the numerics below 0.5%. It has also been pointed out [53] that an artificial viscosity (AV) shock capturing scheme has an advantage over Riemann solvers because only about half as many zones are required to accurately resolve the stars when an AV scheme is employed compared to a Riemann solver.The time stepsdtn are taken as the minimum of the time step as determined by several conditions. Each condition is also multiplied by a number less than one to accommodate the nonlinear nature of the equations.The first condition is known as the Courant condition, that is, a search over all zonesi for the zone with the minimum sound crossing time:(21)dt1=min⁡dxbiCsi, where Csi is the sound speed in the ith zone.The Newtonian sound speed is given by the variation of pressure with density. In relativity the wave speed is given instead by the adiabatic derivative of the pressure with respect to the relativistic inertial density. In terms of relativistic variables the local sound speed in zonei then becomes [53](22)Cs=ΓiΓi-1EiDi+ΓiEi.The second condition is a search for the zone with minimum time for material to flow across a zone(23)dt2=min⁡dxaiVi. This constraint is introduced to ensure stability and accuracy in the numerical advection calculation.The third condition is introduced to maintain stability of the artificial viscosity algorithm. The viscous equations are analogous to a diffusion equation in four velocity with a diffusion coefficientD≈k1dxi|δUi|, where δUi≡Ui+1-Ui. We then can define a minimum viscous diffusion time across a zone derived from the stability condition for explicit diffusion equations(24)dt3=14min⁡WidxbiδUi.The time stepdt is then assigned to be some fraction (referred to as the Courant parameter) of the minimum of these three conditions(25)dt=kmin⁡dt1,dt2,dt3,Obviously, smaller values for k increase the accuracy of the calculation but also increase the computation time.Figure2 shows a plot of orbital velocity versus time for various values of the Courant parameter. This figure establishes that the routines for the hydrodynamics are stable (e.g., changing the Courant condition has little to effect) as long as k≤0.5.Figure 2
Comparison of the orbital angular velocityω versus time for different values of the Courant parameter k. As can be seen, the simulations with k = 0.25–0.5 result in stable runs that converge to the same value, implying that a smaller k or equivalently a smaller δt is not necessary and would only use extra CPU time. For comparison, we plot a simulation with k=0.6 to show that the stability is lost for k>0.5.Figure3 illustrates the difference between the central density ρc and the central density (ρ52) at the highest resolution of 52 zones for single isolated stars. This is expressed as the fractional error as a function of the number of zones across a star. This plot was calculated using the relatively soft MW EoS, that is, the zero temperature and zero neutrino chemical potential limit of the EoS that has previously been used to model core-collapse supernovae [32, 53, 65].Figure 3
Plot of the error in the central density versus the number of zones across the star. It is clear that there is only a 1% error with ≈15 zones across the star. Increasing the number of zones across the star so that there are >35 zones across the star produces less than a 0.1% error.This figure illustrates that here is only a 1% error in central density with ≈15 zones across the star, while increasing the number of zones across the star to >35 produces less than a 0.1%. In the illustrations below we maintaink=0.5 and ≈25 zones across each star as the best choice for both speed and accuracy needed to compute stable orbital initial conditions.


## 3.2. Orbit Stability

As an illustration of the orbit stability Figure4 shows results from a simulation [54] in which the angular momentum was fixed at J=2.7×1011 cm2 and the Courant parameter set to k=0.5. For this orbital calculation the MW EoS was employed and each star was fixed at a baryon mass of MB=1.54M⊙ and a gravitational mass in isolation of MG=1.40M⊙.Figure 4
Plot of the orbital angular velocity,ω in geometrized units (cm−1) versus cycle number. When ω stops changing the simulation has reached a circular binary orbit solution. This run was extended to over 30,000 cycles, corresponding to ≈20 orbits.Figure4 shows the evolution of the orbital angular velocity ω versus computational cycle for the first 30,000 code cycles corresponding to ≈20 orbits. The stars were initially placed on the grid using a solution of the TOV equation in isotropic coordinates for an isolated star. The stars were initially set to be corotating but were allowed to settle into their binary equilibrium. Notice that it takes ~5,000 cycles, corresponding to ~3 orbits, just to approach the quasiequilibrium binary solution. Indeed, the stars continued to gradually compact and slightly increase in orbital frequency until ~10 orbits; afterward, the stars were completely stable.Figure5 shows the contours of the lapse function α (roughly corresponding to the gravitational potential) and corresponding density profiles at cycle numbers, 0, 5200, and 25800 (≈0, 5, and 19 orbits). Figure 6 shows the contours of central density and the orientation of the binary orbit corresponding to these cycle numbers. One can visibly see from these figures the relaxation of the stars after the first few orbits and the stability of the density profiles after multiple orbits.Figure 5
Contours of the lapse function (left) and central density (right) at cycle numbers 0 (a), 5,200 (b), and 25,800 (c) corresponding to roughly 0, 5, and 19 orbits.
(a) (b) (c)Figure 6
Contours of the central density for the binary system at the approximate number of orbits as labelled.We note, however, that this orbit is on the edge of the ISCO. As such it could be unstable to inspiral even after many orbits. Figures7 and 8 further illustrate this point. In these simulations various angular momenta were computed with a slightly higher neutron star mass (Mb=1.61M⊙ and Mg=1.44M⊙), but the same MW EoS. In this case the binary system was followed for nearly 100 orbits.Figure 7
Plot of the orbital angular velocity,ω, versus cycle. When ω stops changing with time the simulation has reached a circular binary orbit solution. The run (a for J=2.7×1011 cm2) goes over ~10 obits and then becomes unstable to inspiral and merger after ~104 cycles. The stable two runs (b for J=2.8×1011 cm−2 and c for J=2.9×1011 cm−2) were run for 100,000 cycles and ≈100 orbits.Figure 8
Plot of the central density,ρc, versus the number of orbit. The solid lines from top to bottom are for J=3.0,3.2,3.4,3.6,3.8,4.0×1011 cm2.Figure7 illustrates orbital angular frequency versus cycle number for three representative angular momenta bracketing the ISCO. The orbital separation for the lowest angular momentum (J=2.7×1011 cm−2) shown on Figure 7 is just inside the ISCO. Hence, even though it requires about 10 orbits before inspiral, the orbit is eventually unstable. Similarly, Figure 8 shows the central density versus number of orbits for 6 different angular momenta. Here one can see that orbits with J≥3.0×1011 cm−2 are stable. Indeed, for these cases, after about the first 3 orbits the orbits continue with almost no discernible change in orbit frequency or central density.As mentioned previously, the numerical relativistic neutron binary simulations of [43] all start with initial data that are subsequently evolved in a different manner compared to those with which they were created. One conclusion that may be drawn from the above set of simulations, however, is that the initial data must be evolved for ample time (>3 orbit) for the stars to reach a true quasiequilibrium binary configuration. That has not always been done in the literature.


## 4. Sensitivity of Initial Condition Orbital Parameters to the Equation of State

### 4.1. Equations of State

One hope in the forthcoming detection of gravitational waves is that a sensitivity exists to the neutron star equation of state. For illustration we utilize several representative equations of state often employed in the literature. These span a range from relatively soft to stiff nuclear matter. These are used to illustrate the EoS dependence of the initial conditions. One EoS often employed is that of a polytrope, that is,p=KρΓ, with Γ=2, where in cgs units, K=0.0445(c2/ρn), and ρn=2.3×1014 g cm−3. These parameters, with ρc=4.74×1014 g cm−3, produce an isolated star having a radius = 17.12 km and baryon mass = 1.5M⊙. As noted in previous sections we utilize the zero temperature and zero neutrino chemical potential MW EoS [32, 53, 65]. The third is the equation of state developed by Lattimer and Swesty [66] with two different choices of compressibility, one having compressibility K=220 MeV and the other having K=375 MeV. We denote these as LS 220 and LS 375. The fourth EoS has been developed by Glendenning [67]. This EoS has K=240 MeV, which is close to the experimental value [68]. We denote this EoS as GLN. Table 1 illustrates [54] the properties of isolated neutron stars generated with each EoS. For each case the baryon mass was chosen to obtain a gravitational mass for each star of 1.4M⊙.Table 1
Table presenting central density, baryon mass, and gravitational mass for the five adopted equations of state.
EoS ρ c (×1015 g cm−3) M B (M⊙) M G (M⊙) GLN 1.56 1.54 1.40 MW 1.39 1.54 1.40 LS 220 0.698 1.54 ~1.40 LS 375 0.492 1.54 ~1.40 Γ = 2 polytrope 0.474 1.50 1.40In Figure9 we plot the equation of state index Γ versus density, ρ, for the various EoSs considered here. These are compared to the simple polytropic Γ = 2 EoS often employed in the literature.Figure 9
EoS indexΓ versus central density for various equations of state. Large Γ implies a stiff EoS.


### 4.2. EoS Dependence of the Initial Condition Orbit Parameters

Table2 summarizes the initial condition orbit parameters [54] at various fixed angular momenta for the various equations of state. In the case of orbits unstable to merger, this table lists the orbit parameters just before inspiral. These orbits span a range in specific angular momenta J/M02 of ~5 to 10. The equations of state listed in Tables 1 and 2 are in approximate order of increasing stiffness from the top to the bottom.Table 2
Orbital parameters for each EoS.
EoS J (cm2) ω (rad s−1) d p (km) d c (km) M A D M (M⊙) ρ c (g cm−3) GLN 2.7 × 1 0 11 666.5 71.62 57.67 1.390 1.73 × 1 0 15 2.8 × 1 0 11 592.34 77.82 62.81 1.391 1.69 × 1 0 15 3.0 × 1 0 11 475.05 88.06 73.53 1.394 1.61 × 1 0 15 3.2 × 1 0 11 391.75 100.34 84.31 1.396 1.56 × 1 0 15  MW 2.6 × 1 0 11 780.92 65.22 51.52 1.391 1.67 × 1 0 15 2.7 × 1 0 11 671.85 71.18 57.24 1.393 1.62 × 1 0 15 2.8 × 1 0 11 602.80 76.94 61.86 1.394 1.60 × 1 0 15 3.0 × 1 0 11 482.30 86.91 72.36 1.396 1.55 × 1 0 15 3.5 × 1 0 11 300.46 116.13 100.8 1.399 1.44 × 1 0 15 3.8 × 1 0 11 235.72 136.93 119.74 1.401 1.39 × 1 0 15  LS 220 2.7 × 1 0 11 523.59 90.77 77.34 1.403 7.18 × 1 0 14 2.8 × 1 0 11 472.08 97.53 83.08 1.404 7.14 × 1 0 14 3.0 × 1 0 11 389.96 109.78 94.84 1.405 7.06 × 1 0 14 3.2 × 1 0 11 327.04 122.51 107.10 1.407 6.98 × 1 0 14  LS 375 2.7 × 1 0 11 490.09 97.09 83.92 1.404 5.00 × 1 0 14 2.8 × 1 0 11 442.40 103.95 90.04 1.405 4.98 × 1 0 14 3.0 × 1 0 11 366.67 116.65 102.50 1.406 4.95 × 1 0 14 3.2 × 1 0 11 307.80 130.72 115.60 1.407 4.92 × 1 0 14  Polytrope 1.8 × 1 0 11 804.70 63.30 51.20 1.395 6.78 × 1 0 14 2.1 × 1 0 11 826.03 67.85 55.18 1.396 7.00 × 1 0 14 2.3 × 1 0 11 762.37 74.64 61.72 1.397 6.55 × 1 0 14 2.5 × 1 0 11 624.33 85.87 72.71 1.399 6.24 × 1 0 14 2.6 × 1 0 11 532.83 94.04 80.45 1.400 6.17 × 1 0 14 2.7 × 1 0 11 477.19 101.34 86.95 1.400 6.05 × 1 0 14As expected, the central densities are much higher for the relatively soft MW and GLN equations of state. Also, the orbit angular frequencies are considerably lower for the extended mass distributions of the stiff equations of state than for the more compact soft equations of state. These extended mass distributions induce a sensitivity of the emergent gravitational wave frequencies and amplitude due to the strong dependence of the gravitational wave frequency to the quadrupole moment of the mass distribution.


### 4.3. Gravitational Wave Frequency

The physical processes occurring during the last orbits of a neutron star binary are currently a subject of intense interest. As the stars approach their final orbits it is expected that the coupling of the orbital motion to the hydrodynamic evolution of the stars in the strong relativistic fields could provide insight into various physical properties of the coalescing system [58, 69]. In this regard, careful modeling of the initial conditions is needed which includes both the nonlinear general relativistic and hydrodynamic effects as well as a realistic neutron star equation of state.Figure10 shows the EoS sensitivity of the gravitational wave frequency f=ω/π as a function of proper separation dp between the stars for the various orbits and equations of state summarized in Table 2. These are compared with the circular orbit condition in the (post)5/2-Newtonian approximation, hereafter PN, analysis of reference [70]. In that paper a search was made for the inner most stable circular orbit in the absence of radiation reaction terms in the equations of motion. This is analogous to the calculations performed here which also analyzes orbit stability in the absence of radiation reaction.Figure 10
Computed gravitational wave frequency,f, versus proper separation for each EoS as labelled. The black line corresponds to the (post)5/2-Newtonian estimate. Frequencies obtained from the stiff and polytropic equations of state do not deviate by more than ~10% from the PN prediction until a frequency greater than ~300 Hz. The grey line is an extrapolation of the frequencies obtained using the soft MW and GLN EoSs. These begin to deviate by more than 10% from the PN prediction at a frequency of ~100 Hz.In the (post)5/2-Newtonian equations of motion, a circular orbit is derived by setting time derivatives of the separation, angular frequency, and the radial acceleration to zero. This leads to the circular orbit condition [70] (26)ω02=mA0dh3, where ω0 is the circular orbit frequency and m=2MG0, dh is the separation in harmonic coordinates, and A0 is a relative acceleration parameter which for equal mass stars becomes(27)A0=1-32mdh3-778mdh+ω0dh2+74ω0dh2.Equations (26) and (27) can be solved to find the orbit angular frequency as a function of harmonic separation dh. The gravitational wave frequency is then twice the orbit frequency, f=ω0/π.Although this is a gauge-dependent comparison, for illustration we show in Figure10 the calculated gravitational wave frequency compared to the PN expectation as a function of proper binary separation distance up to 200km. One should keep in mind, however, that there is some uncertainty in associating proper distance with the PN parameter (m/r). Hence, a comparison with the PN results is not particularly meaningful. It is nevertheless instructive to consider the difference in the numerical results as the stiffness of the EoS is varied. The polytropic and stiff EoSs begin to deviate (by >10%) from the softer equations of state (MW and GLN) for a gravitational wave frequency as low as ~100 Hz and more or less continue to deviate as the stars approach the ISCO at higher frequencies.Indeed, a striking feature of Figure10 is that as the stars approach the ISCO, the frequency varies more slowly with diminishing separation distance for the softer equations of state. A gradual change in frequency can mean more orbits in the LIGO window and hence a stronger signal to noise (cf. [57]).Also, for a soft EoS the orbit becomes unstable to inspiral at a larger separation. At least part of the difference between the soft and stiff EoSs can be attributed to the effects of the finite size of the stars which are more compact for the soft equations of state [37].We note that, for comparable angular momenta, our results are consistent with the EoS sensitivity study of [37] based upon a set of equations of state parameterized by a segmented polytropic indices and an overall pressure scale. Their calculations, however, were based upon two independent numerical relativity codes. The similarity of their simulations to the results in Table 2 further confirms the broad validity of the CFA approach when applied to initial conditions. For example, their orbit parameters are summarized in Table II of [37]. Their softest EoS is the Bss221 which corresponds to an adiabatic index of Γ=2.4 for the core and a baryon mass of 1.501M⊙ and an ADM mass of 1.338M⊙ per star for a specific angular momentum of 1.61×1011 cm2 (in our units) with a corresponding gravitational wave frequency of 530 Hz at a proper separation of 46 km. This EoS is comparable to the polytropic, MW, and GLN EoSs shown on Figure 10. For example, our closest orbit with the Γ=2 polytropic EoS corresponds to a specific angular momentum of 1.8×1011 cm2 and an ADM mass of 1.39M⊙ compared to their ADM mass of 1.34M⊙ at J=1.6×1011 cm2 for the same baryon mass of 1.5M⊙. Although, for the softer EoSs, their results are for a closer orbit than the numerical points given on Figure 10, an extension of the grey line fit to the numerical simulations of the soft EoSs would predict a frequency of 540 Hz at the same proper separation of 46 km compared to 530 Hz in the Bss221 simulation of [37].The main parameter characterizing the last stable orbit in the post-Newtonian calculation is the ratio of coordinate separation to total mass (in isolation)dh/m. The analogous quantity in our nonperturbative simulation is proper separation to gravitational mass, dP/m. The separation corresponding to the last stable orbit in the post-Newtonian analysis does not occur until the stars have approached 6.03m. For MG0=1.4M⊙ stars, this would correspond to a separation distance of about 25 km. In the results reported here the last stable orbit occurs somewhere just below 7.7mG0 at a proper separation distance of dP≈30 km for both the polytropic and the MW stars.


## 4.1. Equations of State

One hope in the forthcoming detection of gravitational waves is that a sensitivity exists to the neutron star equation of state. For illustration we utilize several representative equations of state often employed in the literature. These span a range from relatively soft to stiff nuclear matter. These are used to illustrate the EoS dependence of the initial conditions. One EoS often employed is that of a polytrope, that is,p=KρΓ, with Γ=2, where in cgs units, K=0.0445(c2/ρn), and ρn=2.3×1014 g cm−3. These parameters, with ρc=4.74×1014 g cm−3, produce an isolated star having a radius = 17.12 km and baryon mass = 1.5M⊙. As noted in previous sections we utilize the zero temperature and zero neutrino chemical potential MW EoS [32, 53, 65]. The third is the equation of state developed by Lattimer and Swesty [66] with two different choices of compressibility, one having compressibility K=220 MeV and the other having K=375 MeV. We denote these as LS 220 and LS 375. The fourth EoS has been developed by Glendenning [67]. This EoS has K=240 MeV, which is close to the experimental value [68]. We denote this EoS as GLN. Table 1 illustrates [54] the properties of isolated neutron stars generated with each EoS. For each case the baryon mass was chosen to obtain a gravitational mass for each star of 1.4M⊙.Table 1
Table presenting central density, baryon mass, and gravitational mass for the five adopted equations of state.
EoS ρ c (×1015 g cm−3) M B (M⊙) M G (M⊙) GLN 1.56 1.54 1.40 MW 1.39 1.54 1.40 LS 220 0.698 1.54 ~1.40 LS 375 0.492 1.54 ~1.40 Γ = 2 polytrope 0.474 1.50 1.40In Figure9 we plot the equation of state index Γ versus density, ρ, for the various EoSs considered here. These are compared to the simple polytropic Γ = 2 EoS often employed in the literature.Figure 9
EoS indexΓ versus central density for various equations of state. Large Γ implies a stiff EoS.


## 4.2. EoS Dependence of the Initial Condition Orbit Parameters

Table2 summarizes the initial condition orbit parameters [54] at various fixed angular momenta for the various equations of state. In the case of orbits unstable to merger, this table lists the orbit parameters just before inspiral. These orbits span a range in specific angular momenta J/M02 of ~5 to 10. The equations of state listed in Tables 1 and 2 are in approximate order of increasing stiffness from the top to the bottom.Table 2
Orbital parameters for each EoS.
EoS J (cm2) ω (rad s−1) d p (km) d c (km) M A D M (M⊙) ρ c (g cm−3) GLN 2.7 × 1 0 11 666.5 71.62 57.67 1.390 1.73 × 1 0 15 2.8 × 1 0 11 592.34 77.82 62.81 1.391 1.69 × 1 0 15 3.0 × 1 0 11 475.05 88.06 73.53 1.394 1.61 × 1 0 15 3.2 × 1 0 11 391.75 100.34 84.31 1.396 1.56 × 1 0 15  MW 2.6 × 1 0 11 780.92 65.22 51.52 1.391 1.67 × 1 0 15 2.7 × 1 0 11 671.85 71.18 57.24 1.393 1.62 × 1 0 15 2.8 × 1 0 11 602.80 76.94 61.86 1.394 1.60 × 1 0 15 3.0 × 1 0 11 482.30 86.91 72.36 1.396 1.55 × 1 0 15 3.5 × 1 0 11 300.46 116.13 100.8 1.399 1.44 × 1 0 15 3.8 × 1 0 11 235.72 136.93 119.74 1.401 1.39 × 1 0 15  LS 220 2.7 × 1 0 11 523.59 90.77 77.34 1.403 7.18 × 1 0 14 2.8 × 1 0 11 472.08 97.53 83.08 1.404 7.14 × 1 0 14 3.0 × 1 0 11 389.96 109.78 94.84 1.405 7.06 × 1 0 14 3.2 × 1 0 11 327.04 122.51 107.10 1.407 6.98 × 1 0 14  LS 375 2.7 × 1 0 11 490.09 97.09 83.92 1.404 5.00 × 1 0 14 2.8 × 1 0 11 442.40 103.95 90.04 1.405 4.98 × 1 0 14 3.0 × 1 0 11 366.67 116.65 102.50 1.406 4.95 × 1 0 14 3.2 × 1 0 11 307.80 130.72 115.60 1.407 4.92 × 1 0 14  Polytrope 1.8 × 1 0 11 804.70 63.30 51.20 1.395 6.78 × 1 0 14 2.1 × 1 0 11 826.03 67.85 55.18 1.396 7.00 × 1 0 14 2.3 × 1 0 11 762.37 74.64 61.72 1.397 6.55 × 1 0 14 2.5 × 1 0 11 624.33 85.87 72.71 1.399 6.24 × 1 0 14 2.6 × 1 0 11 532.83 94.04 80.45 1.400 6.17 × 1 0 14 2.7 × 1 0 11 477.19 101.34 86.95 1.400 6.05 × 1 0 14As expected, the central densities are much higher for the relatively soft MW and GLN equations of state. Also, the orbit angular frequencies are considerably lower for the extended mass distributions of the stiff equations of state than for the more compact soft equations of state. These extended mass distributions induce a sensitivity of the emergent gravitational wave frequencies and amplitude due to the strong dependence of the gravitational wave frequency to the quadrupole moment of the mass distribution.


## 4.3. Gravitational Wave Frequency

The physical processes occurring during the last orbits of a neutron star binary are currently a subject of intense interest. As the stars approach their final orbits it is expected that the coupling of the orbital motion to the hydrodynamic evolution of the stars in the strong relativistic fields could provide insight into various physical properties of the coalescing system [58, 69]. In this regard, careful modeling of the initial conditions is needed which includes both the nonlinear general relativistic and hydrodynamic effects as well as a realistic neutron star equation of state.Figure10 shows the EoS sensitivity of the gravitational wave frequency f=ω/π as a function of proper separation dp between the stars for the various orbits and equations of state summarized in Table 2. These are compared with the circular orbit condition in the (post)5/2-Newtonian approximation, hereafter PN, analysis of reference [70]. In that paper a search was made for the inner most stable circular orbit in the absence of radiation reaction terms in the equations of motion. This is analogous to the calculations performed here which also analyzes orbit stability in the absence of radiation reaction.Figure 10
Computed gravitational wave frequency,f, versus proper separation for each EoS as labelled. The black line corresponds to the (post)5/2-Newtonian estimate. Frequencies obtained from the stiff and polytropic equations of state do not deviate by more than ~10% from the PN prediction until a frequency greater than ~300 Hz. The grey line is an extrapolation of the frequencies obtained using the soft MW and GLN EoSs. These begin to deviate by more than 10% from the PN prediction at a frequency of ~100 Hz.In the (post)5/2-Newtonian equations of motion, a circular orbit is derived by setting time derivatives of the separation, angular frequency, and the radial acceleration to zero. This leads to the circular orbit condition [70] (26)ω02=mA0dh3, where ω0 is the circular orbit frequency and m=2MG0, dh is the separation in harmonic coordinates, and A0 is a relative acceleration parameter which for equal mass stars becomes(27)A0=1-32mdh3-778mdh+ω0dh2+74ω0dh2.Equations (26) and (27) can be solved to find the orbit angular frequency as a function of harmonic separation dh. The gravitational wave frequency is then twice the orbit frequency, f=ω0/π.Although this is a gauge-dependent comparison, for illustration we show in Figure10 the calculated gravitational wave frequency compared to the PN expectation as a function of proper binary separation distance up to 200km. One should keep in mind, however, that there is some uncertainty in associating proper distance with the PN parameter (m/r). Hence, a comparison with the PN results is not particularly meaningful. It is nevertheless instructive to consider the difference in the numerical results as the stiffness of the EoS is varied. The polytropic and stiff EoSs begin to deviate (by >10%) from the softer equations of state (MW and GLN) for a gravitational wave frequency as low as ~100 Hz and more or less continue to deviate as the stars approach the ISCO at higher frequencies.Indeed, a striking feature of Figure10 is that as the stars approach the ISCO, the frequency varies more slowly with diminishing separation distance for the softer equations of state. A gradual change in frequency can mean more orbits in the LIGO window and hence a stronger signal to noise (cf. [57]).Also, for a soft EoS the orbit becomes unstable to inspiral at a larger separation. At least part of the difference between the soft and stiff EoSs can be attributed to the effects of the finite size of the stars which are more compact for the soft equations of state [37].We note that, for comparable angular momenta, our results are consistent with the EoS sensitivity study of [37] based upon a set of equations of state parameterized by a segmented polytropic indices and an overall pressure scale. Their calculations, however, were based upon two independent numerical relativity codes. The similarity of their simulations to the results in Table 2 further confirms the broad validity of the CFA approach when applied to initial conditions. For example, their orbit parameters are summarized in Table II of [37]. Their softest EoS is the Bss221 which corresponds to an adiabatic index of Γ=2.4 for the core and a baryon mass of 1.501M⊙ and an ADM mass of 1.338M⊙ per star for a specific angular momentum of 1.61×1011 cm2 (in our units) with a corresponding gravitational wave frequency of 530 Hz at a proper separation of 46 km. This EoS is comparable to the polytropic, MW, and GLN EoSs shown on Figure 10. For example, our closest orbit with the Γ=2 polytropic EoS corresponds to a specific angular momentum of 1.8×1011 cm2 and an ADM mass of 1.39M⊙ compared to their ADM mass of 1.34M⊙ at J=1.6×1011 cm2 for the same baryon mass of 1.5M⊙. Although, for the softer EoSs, their results are for a closer orbit than the numerical points given on Figure 10, an extension of the grey line fit to the numerical simulations of the soft EoSs would predict a frequency of 540 Hz at the same proper separation of 46 km compared to 530 Hz in the Bss221 simulation of [37].The main parameter characterizing the last stable orbit in the post-Newtonian calculation is the ratio of coordinate separation to total mass (in isolation)dh/m. The analogous quantity in our nonperturbative simulation is proper separation to gravitational mass, dP/m. The separation corresponding to the last stable orbit in the post-Newtonian analysis does not occur until the stars have approached 6.03m. For MG0=1.4M⊙ stars, this would correspond to a separation distance of about 25 km. In the results reported here the last stable orbit occurs somewhere just below 7.7mG0 at a proper separation distance of dP≈30 km for both the polytropic and the MW stars.


## 5. Conclusions

The relativistic hydrodynamic equilibrium in the CFA remains as a viable approach to calculate the initial conditions for calculations of binary neutron stars. In this paper we have illustrated that one must construct initial conditions that have run for at least several orbits before equilibrium is guaranteed. We have demonstrated that beyond the first several orbits the equations are stable over many orbits (~100). We also have shown that such multiple orbit simulations are valuable as a means to estimate the location of the ISCO prior to a full dynamical calculation. Moreover, we have examined the sensitivity of the initial condition orbit parameters and initial gravitational wave frequency to the equation of state. We have illustrated how the initial condition orbital properties (e.g., central densities, orbital velocities, and binding energies) and location of the ISCO are significantly affected by the stiffness of the EoS.



---

*Source: 6127031-2017-01-09.xml*