Initial data for binary neutron stars with adjustable eccentricity
Niclas Moldenhauer, Charalampos M. Markakis, Nathan K. Johnson-McDaniel, Wolfgang Tichy, Bernd Bruegmann
IInitial data for binary neutron stars with adjustable eccentricity
Niclas Moldenhauer, ∗ Charalampos M. Markakis,
1, 2, † Nathan K.Johnson-McDaniel, ‡ Wolfgang Tichy, § and Bernd Br¨ugmann ¶ Theoretical Physics Institute, University of Jena, Max-Wien-Platz 1, 07743 Jena, Germany Mathematical Sciences, University of Southampton, Southampton, SO17 1BJ, United Kingdom Physics Department, Florida Atlantic University, 777 Glades Road, Boca Raton, FL 33431, USA (Dated: November 6, 2018)Binary neutron stars in circular orbits can be modeled as helically symmetric, i.e., stationary ina rotating frame. This symmetry gives rise to a first integral of the Euler equation, often employedfor constructing equilibrium solutions via iteration. For eccentric orbits, however, the lack of helicalsymmetry has prevented the use of this method, and the numerical relativity community has oftenresorted to constructing initial data by superimposing boosted spherical stars without solving theEuler equation. The spuriously excited neutron star oscillations seen in evolutions of such data arisebecause such configurations lack the appropriate tidal deformations and are stationary in a linearlycomoving—rather than rotating—frame. We consider eccentric configurations at apoapsis that areinstantaneously stationary in a rotating frame. We extend the notion of helical symmetry to ec-centric 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 approxi-mate instantaneous first integrals of the Euler equation throughout each companion. We use theseintegrals as the basis of a self-consistent iteration of the Einstein constraints to construct conformalthin-sandwich initial data for eccentric binaries. We find that the spurious stellar oscillations arereduced by at least an order of magnitude, compared with those found in evolutions of superposedinitial data. The tidally induced oscillations, however, are physical and qualitatively similar to ear-lier evolutions. Finally, we show how to incorporate radial velocity due to radiation reaction in ourinscribed helical symmetry vectors, which would allow one to obtain truly non-eccentric initial datawhen our eccentricity parameter e is set to zero. PACS numbers: 04.20.Ex, 04.25.dk, 04.30.Db, 97.60.Jd
I. INTRODUCTION
The coalescences of binary neutron stars are a promi-nent source for ground-based gravitational wave detec-tors. (See Table 1 in [1] for a review of populationsynthesis predictions and [2] for more recent work.) Inparticular, binary neutron stars are the only compactbinary sources relevant for ground-based detectors thathave been observed to date (via electromagnetic observa-tions of binary pulsars). Neutron star binary coalescencesare also interesting beyond gravitational wave astronomyas potential progenitors of short gamma-ray bursts andsources of the r -process material that enriches the inter-stellar medium with heavy elements [3].All the known binary neutron stars are at least some-what eccentric (eccentricities between 0 .
085 and 0 .
681 forthe systems that will merge within a few hundred mil-lion years; see Tables 2 and 3 in [4]), though these willbe highly circular when they merge, since gravitationalradiation reaction efficiently circularizes the orbit [5];see Table 1 in [6]. However, there are possible (thoughlikely rare) scenarios in which neutron stars can merge ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] with nonnegligible eccentricity, either because they areformed with a high eccentricity and small periapsis dis-tance by dynamical interactions in dense stellar regions,such as globular clusters [7–10] or have their eccentricityexcited by, e.g., the Kozai mechanism in a hierarchicaltriple [11–14] ([15] also treats the latter case, but onlyconsiders the case of binary black holes). Gravitationalwaves from highly eccentric compact binary systems ex-hibit a repeated burst structure, which poses challengesfor gravitational wave astronomy, but also offers poten-tial rewards, as has been explored in a number of works:McWilliams, Pretorius, and collaborators [16, 17] discussstrategies for detection, while Kyutoku and Seto [18] findimprovements in the accuracy of premerger sky local-ization and timing, compared to the quasicircular case.Loutrel, Yunes, and Pretorius [19] consider bursts fromhighly eccentric binaries as a regime for testing generalrelativity, while Tsang [9] has considered the possibilityof obtaining electromagnetic flares from crust crackingduring close encounters.The first full numerical relativity evolutions of highlyeccentric binary neutron stars were carried out by Gold etal. [20], with a further study by East and Pretorius [21].These systems have also been simulated with Newtonianmethods in [8, 22]. In addition, there have been full nu-merical relativity simulations of highly eccentric blackhole-neutron star (BHNS) binaries [23, 24] and binaryblack holes, e.g., [25, 26]. However, all the evolutionsof highly eccentric binaries with neutron stars have used a r X i v : . [ g r- q c ] O c t inconsistent initial data, due to the difficulty of generaliz-ing the standard procedure for quasicircular orbits, whereone uses the binary’s approximate helical Killing vectorto solve the Euler equation via a first integral. In par-ticular, Gold et al. [20] used a superposition of boostedspherical stars, which leads to relatively large constraintviolations in addition to not giving the appropriate mat-ter configuration. The Princeton group [21, 23, 24] solvesthe constraints to obtain their initial data, as describedin [27], but they do not solve the Euler equation. Itwould obviously be desirable to obtain consistent initialdata for these configurations. In particular, it is possiblethat the tidally induced oscillations of the neutron star(s)found in these evolutions are affected by the initial spuri-ous oscillations of the neutron star from imperfect initialdata.Here we present a method to construct binary neutronstar initial data with arbitrary eccentricity, by generaliz-ing the helical Killing vector to a pair of inscribed helicalsymmetry vectors, appropriate for the more general situ-ation of an eccentric orbit at apoapsis. We then provide afirst proof-of-principle numerical implementation of thismethod for equal-mass binaries; our method is applica-ble to arbitrary mass ratios or BHNS binaries as well.This self-consistent method yields constraint-solved ini-tial data (in the Isenberg-Wilson-Mathews approxima-tion [28, 29], i.e., assuming spatial conformal flatness),where both the geometry and matter are momentarilystationary in a rotating frame. We also give an addi-tional generalization to include radial velocity, thoughwe do not implement this numerically in this paper. Asa test of the method, we show that it produces resultswith the expected physical properties both in the qua-sicircular limit, where we compare with the results ofdata calculated with the standard method, as well as fornonzero values of the eccentricity.For simplicity, we have considered irrotational binariesin the current numerical implementation and made theapproximation of a homogeneous velocity field. How-ever, it is possible (and relatively straightforward) todrop these approximations and even add spin to the con-struction, following [30, 31]. In particular, the assump-tion of a homogeneous velocity field is made merely forconvenience, so that we can use a Cartesian multigridelliptic solver without surface-fitted coordinates, but wedemonstrate that it is a reasonable approximation forsufficiently separated binaries. The assumption of irrota-tional stars is standard and is reasonable for a first study:While neutron stars can spin quite quickly (at least up to716 Hz [32]), the spins in known binary neutron stars aremuch more modest, at most 44 Hz for the more massivestar in the double pulsar (see, e.g., Table 2 in [4]), andall of these will decrease further due to spin-down be-fore the stars merge. Moreover, the viscosity of neutronstar matter is far too low for the stars to experience anysignificant tidal spin-up, as established by Kochanek [33]and Bildsten and Cutler [34]. Thus, it has been standardto consider irrotational flow in modeling binary neutron stars (see, e.g., [1]), since the system’s orbital frequencyten to twenty orbits before merger is (cid:38)
100 Hz in thequasicircular case.However, the extent to which the relatively smallspins of the members of observed binary neutron starsis a selection effect remains unclear, particularly becausethe known population is so small (only 9 systems [4]).Moreover, if one did form a binary with even a mod-estly rapidly spinning neutron star, the spin-down beforemerger might not be very significant: Fast-spinning neu-tron stars are thought to have had their spin increased byaccretion (a process known as recycling), which also re-duces their external magnetic field, and thus reduces thestars’ spin-down, as well (see, e.g., [35]). In particular,[30] finds that the more massive star in the double pulsaris expected to have spun down only to 37 Hz at merger,and [36] has found that spins of about this magnitudecan have a sizable effect on the system’s dynamics.The addition of spin may be particularly interestingfor eccentric systems, since one way of forming such bi-naries is through dynamical assembly in dense stellar sys-tems, such as globular clusters, and globular clusters area fertile breeding ground for millisecond pulsars, includ-ing the fastest pulsar known. For instance, all 23 knownpulsars in the rich globular cluster 47 Tucanae have spinfrequencies greater than 125 Hz, and all but three areabove 200 Hz [37].The paper is structured as follows: We first review thenecessary portions of perfect fluid hydrodynamics andthe 3+1 split of the Einstein equations in Secs. II and III,before describing the specific construction we use to ob-tain an approximate first integral to the Euler equation inSec. IV. We then describe the numerical implementationof the method in Sec. V and evaluate its performance inthe quasicircular case in Sec. VI before giving examplesof eccentric binaries in Sec. VII. We discuss and concludein Sec. VIII, and give some ancillary results for nonrela-tivistic incompressible binaries in the Appendix.We use the following notation throughout: We useGreek letters α, β, γ, δ, ... and µ, ν, κ, λ, ... for abstractand concrete spacetime indices, respectively. We also useLatin letters a, b, c, ... and i, j, k, ... for abstract and con-crete spatial indices, respectively. We raise and lowerconcrete spatial indices with the flat conformal metric,while all other indices are raised and lowered with thephysical metric; the summation convention is always inforce. We shall also use index-free notation when con-venient, denoting vectors (spatial or spacetime) usingboldface. We employ units with G = c = M (cid:12) = 1 al-most exclusively, except that we show the appearancesof G explicitly for clarity when making some Newtoniancalculations in Sec. IV C. II. PERFECT FLUID MODELA. Thermodynamic quantities
We consider a spacetime ( M , g αβ ), i.e., a manifold M endowed with a Lorentzian metric g αβ . Furthermore,we assume that this spacetime is globally hyperbolic, soit possesses a Cauchy surface (and, indeed, can be fo-liated by Cauchy surfaces). Part of this spacetime isoccupied by a perfect fluid, characterized by the energy-momentum tensor T αβ = ( (cid:15) + p ) u α u β + pg αβ , (2.1)where u α is the four-velocity, (cid:15) is the proper energy den-sity, and p the fluid pressure. Moreover, we assume thatthe fluid is a simple fluid , i.e., that all the thermodynamicquantities depend only on the entropy density σ and onthe proper baryon number density n . In particular, therelation (cid:15) = (cid:15) ( σ, n ) . (2.2)is called the equation of state (EOS) of the fluid. Thetemperature T and the baryon chemical potential µ arethen defined by T := ∂(cid:15)∂σ and µ := ∂(cid:15)∂n . (2.3)Then, the first law of thermodynamics can be written as d(cid:15) = µ dn + T dσ. (2.4)As a consequence, p is a function of ( σ, n ) entirely deter-mined by (2.2): p = − (cid:15) + T σ + µn. (2.5)Let us introduce the specific enthalpy, h := (cid:15) + pρ = µm b + T s, (2.6)where m b = 1 . × − kg (the atomic mass unit, whichcan be taken as the average nucleon mass) is the baryonrest mass, ρ is the rest-mass density ρ := m b n, (2.7)and s is the specific entropy: s := σρ . (2.8)The second equality in (2.6) is an immediate consequenceof (2.5). From Eqs. (2.4)–(2.8), we obtain the thermody-namic relations [38, 39] d(cid:15) = hdρ + ρT ds, dp = ρ ( dh − T ds ) . (2.9)To describe the matter inside a neutron star, we haveto make a choice for the EOS in order to close the system. For the present study, we restrict attention to polytropes.Specifically, we assume that p = κρ /n , (2.10)where κ is the polytropic constant and n = − is thepolytropic index (Γ is the adiabatic index). Then, onemay express ρ , p , and (cid:15) as functions of the specific en-thalpy hρ = (cid:20) h − κ (1 + n ) (cid:21) n , (2.11a) p = κ (cid:20) h − κ (1 + n ) (cid:21) n , (2.11b) (cid:15) = (cid:20) n ( h − n (cid:21) (cid:20) h − κ (1 + n ) (cid:21) n . (2.11c)For simplicity, we shall only consider the case n = 1 (i.e.,Γ = 2) in the numerical examples of this paper, sincethis is close to the effective polytropic index of realisticnuclear physics equations of state [40]. Treatment of real-istic EOSs is straightforward using a piecewise polytropicapproximation [41]. B. Euler equation
The relativistic Euler equation follows from the con-servation law of energy-momentum: ∇ α T αβ = 0 , (2.12)where ∇ α denotes the covariant derivative compatiblewith the metric g αβ . Using Eqs. (2.4)–(2.9), the diver-gence of the fluid energy-momentum tensor (2.1) can bedecomposed as ∇ α T αβ = π β ∇ α ( ρu α ) + ρ [ u α (d π ) αβ − T ∇ β s ] . (2.13)Here, π α = hu α denotes the canonical momentum 1-form of a fluid element, while its exterior derivative(d π ) αβ = ∇ α π β − ∇ β π α denotes the canonical vorticity2-form. Invoking the baryon number conservation law ∇ α ( ρu α ) = 0 , (2.14)Eqs. (2.12) and (2.13) yield the Euler equation, u α (d π ) αβ = T ∇ β s, (2.15)written here in the Carter-Lichnerowicz form [38, 42–45]. This particular form is quite useful when the fluidconfiguration possesses certain symmetries [46]. In the“dust” limit ( h → T → u β yields u β ∇ β s = 0 , (2.16)implying that specific entropy is constant along flow lines.If the fluid is barotropic , the thermodynamic quantitiesdepend only on the proper baryon number density n (or,equivalently, rest mass density ρ ). Particular examplesare cold ( T = 0) or homentropic ( s = const . ) fluids.This assumption is appropriate for inspiralling neutron-star binaries, as shock heating is absent and the fluidtemperature is much lower than the Fermi temperature[47]. In the remainder of this paper, we shall restrictour attention to barotropic flows, for which the Eulerequation (2.15) simplifies to u α (d π ) αβ = 0 . (2.17)Barotropic fluid streamlines are geodesics of a Rieman-nian manifold with metric h g αβ [48]. Indeed, the fluidelement action S = (cid:90) τ τ L ( x, u ) dτ = − (cid:90) τ τ h ( x ) (cid:113) − g αβ ( x ) u α u β dτ (2.18)can be minimized (most easily via covariant techniques[45]) to obtain the Euler equation (2.17), with canonicalvelocity given by u α = dx α /dτ and canonical momentumgiven by π α = ∂L/∂u α = hu α . C. Noether’s theorem & Bernoulli’s principle
If the ε -family of infinitesimal coordinate transforma-tions x α → x α + ε k α is a continuous symmetry of thefluid element action (2.18), then Noether’s theorem im-plies that the quantity E = − k α π α (2.19)is conserved along streamlines. Indeed, using the equa-tion of motion (2.17) and the constraint u α u α = −
1, onefinds u α ∇ α E = − h u α u β L k ( h g αβ ) = 0 , (2.20)where L k denotes the Lie derivative with respect to thevector k . In geometrical terms, this result follows fromthe fact that k α is a Killing vector of the conformal met-ric h g αβ . The conservation of the quantity (2.19) is ageneralization of Bernoulli’s principle, which is recoveredin the Newtonian limit if the Killing vector generatestime translations that leave the flow unchanged (i.e., ifthe flow is stationary) [38]. Note, however, that a Killingsymmetry only guarantees a weak Bernoulli principle, inthe sense that the quantity (2.19) is conserved only alongstreamlines, but could differ from one streamline to thenext. In order to obtain a strong
Bernoulli principle,i.e., a quantity conserved throughout the fluid, a secondcondition (such as irrotationality or rigidity) is required.This issue will be revisited in Sec. IV.
III. GRAVITATIONAL FIELD EQUATIONSA. Extended conformal thin-sandwich formulation
We consider a spacetime M = R × Σ which is foliatedby a family of spacelike surfaces Σ t . Making the standard3+1 decomposition in a chart { t, x i } , the spacetime met-ric takes the form ds = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) , (3.1)where α is the lapse and β α is the shift vector. These arerelated to the unit normal n α of the three dimensionalspatial hypersurface Σ t and the timelike vector t α via t α = αn α + β α . The shift is purely spatial and satisfies β α n α = 0, while we define the spatial metric γ ab ( t ) byrestricting the projection tensor γ αβ = g αβ + n α n β to Σ t .In the extended conformal thin-sandwich (XCTS) for-mulation [28, 49–54] one decomposes the spatial metricinto a spatial conformal metric ¯ γ ab and a conformal fac-tor ψ defined by γ ab = ψ ¯ γ ab . One also decomposes theextrinsic curvature of the foliation as K ab := − L n γ ab = A ab + 13 γ ab K, (3.2)where K := K aa is the trace of the extrinsic curvatureand A ab = − ψ α (cid:2) ∂ t ¯ γ ab − ( ¯ Lβ ) ab (cid:3) (3.3)is the rescaled traceless part of the extrinsic curvature.Here, ∂ µ denotes a partial derivative with respect to thecoordinate x µ on M , while( ¯ Lβ ) ab = ¯ D a β b + ¯ D b β a −
23 ¯ γ ab ¯ D c β c (3.4)is the traceless part of L β ¯ γ ab and ¯ D a is the covariantderivative compatible with ¯ γ ab .Decomposing Einstein’s equations G αβ = 8 πT αβ (where G αβ is the Einstein tensor), following the XCTSformulation, we take a set of five equations( G αβ − πT αβ ) n α n β = 0 , (3.5a)( G αβ − πT αβ ) γ cα n β = 0 , (3.5b)( G αβ − πT αβ ) (cid:18) γ αβ + 12 n α n β (cid:19) = 0 , (3.5c)and solve them for the five metric coefficients { ψ, α, β a } on the initial slice Σ . The constraint equations (3.5a)and (3.5b), along with (3.5c) can be written in the formof elliptic equations with nonlinear source terms:¯ D ψ = ψ R − ψ (cid:18) A ab A ab − K + 16 πρ H (cid:19) , (3.6a)¯∆ L β a = ( ¯ Lβ ) ab ¯ D b ln( αψ − ) + 16 παψ J a − αψ − ¯ D b ( α − ψ ∂ t ¯ γ ab ) + 43 α ¯ D a K, (3.6b)¯ D ( αψ ) = αψ (cid:20) A ab A ab + 512 K + 2 π ( ρ H + 2 J ) (cid:21) − ψ ( ∂ t − β b ¯ D b ) K + 18 αψ ¯ R, (3.6c)where ¯ D := ¯ D a ¯ D a , ¯ R is the Ricci scalar of the spatialconformal metric ( ¯ R ab is its Ricci tensor) and ¯∆ L β a :=¯ D β a + ¯ D a ( ¯ D b β b ) + ¯ R ab β b .The matter source terms are ρ H , J i , and J , which cor-respond to the energy density, the momentum flux, andthe trace of the stress tensor. They are defined as pro-jections of the stress-energy tensor T αβ and thus can bewritten as ρ H = T αβ n α n β , (3.7a) J c = − T αβ γ cα n β , (3.7b) J = T αβ γ αβ . (3.7c)Using the perfect fluid stress-energy tensor (2.1) and ourassumption of spatial conformal flatness, we obtain ρ H = ρh ( αu t − + (cid:15), (3.8a) J i = ρh α ( u t ) ψ ( β i + u i /u t ) , (3.8b) J = ρh [( αu t ) −
1] + 3 p. (3.8c)If one imposes maximal slicing K = 0 , (3.9a) ∂ t K = 0 , (3.9b)and assumes spatial conformal flatness (Isenberg-Wilson-Mathews [IWM] approximation [28, 29])¯ γ ab = f ab , (3.10a) ∂ t ¯ γ ab = 0 , (3.10b)(where f ab is the metric of flat space) and preserves theseconditions in time (at least for an infinitesimally smalltime interval), then the XCTS equations (3.6) simplifyconsiderably. In Cartesian coordinates ( f ij = δ ij ), theyreduce to [54] ∂ i ∂ i ψ = − ψ ( A ij A ij + 16 πρ H ) , (3.11a) ∂ j ∂ j β i + 13 ∂ i ∂ j β j = 2 ψ A ij ∂ j ( αψ − ) + 16 παψ J i , (3.11b) ∂ i ∂ i ( αψ ) = αψ (cid:20) A ij A ij + 2 π ( ρ H + 2 J ) (cid:21) . (3.11c) Here we raise and lower indices with the flat conformalmetric, and will do so for all other concrete spatial in-dices.Note that, in the literature (as reviewed in, e.g., [54]),authors traditionally invoke a time-like Killing symmetryor quasi-equilibrium to justify the conditions (3.9b) and(3.10b). Typically, the maximal slicing and spatial con-formal flatness conditions (3.9a) and (3.10a) are imposedat a later stage. However, swapping the order of assump-tions makes the Killing symmetry redundant. That is,if one imposes the conditions (3.9a), (3.10a) from thebeginning and preserves these conditions in time, thenEqs. (3.9b), (3.10b) follow without assuming Killing sym-metry or quasi-equilibrium. A notion of stationarity willbe introduced in Sec. IV for the fluid sector, but it is notnecessary for the gravity sector of our system if the IWMapproximation is employed.Note that the IWM approximation, while technicallyconvenient, does not allow for purely outgoing gravita-tional radiation, as would be present in an isolated binaryin nature. In particular, the metric of a non-spinningbinary system is known to no longer be spatially confor-mally flat beyond the first post-Newtonian approxima-tion: See, e.g., the discussion in [55, 56] in the quasicir-cular case, and [57] for some work including eccentricity.The assumption of conformal flatness is thus thought tobe responsible for at least some of the initial spuriousradiation observed at the beginning of all numerical rela-tivity simulations of compact binaries. (One finds reduc-tions in some components of the initial spurious radiationwhen one drops the assumption of conformal flatness inthe binary black hole case, e.g., [58–62].) The waveless approach [63, 64] involves a (constraint-solved) construc-tion of binary neutron star data that does not assumespatial conformal flatness, with some unpublished evolu-tions [65], but in general this aspect has not been stud-ied nearly as well for binary neutron stars as for binaryblack holes. Gravitational waves can naturally be ac-commodated in the fully constrained formulation of gen-eral relativity [66–69]. Nevertheless, the overall physicsof the simulation is not significantly affected by this as-pect of the initial data (see [70] in the binary black holecase), though the high-frequency spurious radiation candecrease the accuracy of the simulation [71, 72]. IV. STATIONARY FLUID APPROXIMATIONA. Circular orbits & helical symmetry
Binaries on circular orbits possess a helical Killing vec-tor which generates time translations in a rotating frame: k α = t α + Ω ϕ α = t α + Ω ( x y α − y x α ) . (4.1)The vectors t = ∂ t , x = ∂ x , and y = ∂ y generate trans-lations in the t , x , and y directions, respectively, while ϕ = ∂ ϕ generates rotations in the ϕ direction (i.e., aboutthe z -axis); Ω is the orbital frequency. For circular or-bits, the system appears stationary in a frame corotatingwith the binary. Hence, by virtue of Noether’s theorem(2.19) and (2.20), the energy in a rotating frame , E = − k α π α (4.2)is conserved along streamlines. The conserved quan-tity E , often called the injection energy [47], is analo-gous to the Jacobi constant of motion of test particlesaround Newtonian circular binaries [45]. For irrotationalor corotating binaries on circular orbits, this quantity isconstant throughout each star and is extremely usefulfor constructing equilibrium models numerically via self-consistent field methods [73–78]. As mentioned earlier,this conservation law is equivalent to (a relativistic gen-eralization of) the strong Bernoulli principle and followsquickly from the Cartan identity [38, 45] applied to anyvector field k α that Lie-derives the flow: L k π α = k β (d π ) βα + ∇ α ( k β π β ) = 0 . (4.3)For flows where u α is parallel to k α (such as rigid rotationor rigid translation), the first term in the above equationvanishes by virtue of the Euler equation (2.17). In themore relevant case where the flow is irrotational , i.e., thecanonical vorticity vanishes, we have π α = ∇ α Ψ ⇔ (d π ) αβ = ∇ α π β − ∇ β π α = 0 , (4.4)for some velocity potential Ψ. The Euler equation (2.17)is thus automatically satisfied and the first term in theidentity (4.3) again vanishes, implying that the injectionenergy (4.2) is constant throughout the star: ∇ α E = 0 . (4.5)For many binary neutron stars, the stars’ spin frequencyis much smaller than the orbital frequency for the lastten to twenty orbits before merger, making irrotationalflow a natural approximation, as is discussed in Sec. I. B. Constant three-velocity approximation
For any barotropic flow, Kelvin’s circulation theoremguarantees that a flow that is initially irrotational willremain irrotational [47, 79]. This result is exact for time-dependent spacetimes without symmetries. Hence, inwhat follows, we shall allow for eccentricity but we will re-tain the assumption of irrotational flow . To construct ir-rotational initial data, one typically substitutes Eq. (4.4)into the continuity equation (2.14) and numerically solvesthe resulting equation ∇ α (cid:16) ρh ∇ α Ψ (cid:17) = 0 (4.6)for the velocity potential. This equation requires bound-ary conditions on the star surface and thus surface-fitted coordinates are typically used. For elliptic solvers basedon Cartesian multigrid methods, this is technically diffi-cult.Nevertheless, multigrid methods are widespread, asthey provide a simple way to test a new method and allowone to solve on the grid used for evolution, avoiding theneed to interpolate as well as the need for surface-fittedcoordinates (which are required for spectral solvers, suchas lorene [80] and sgrid [81]). In particular, there is amultigrid method implemented in the bam code [82, 83]that we also use for evolutions. For such multigrid im-plementations, instead of solving for the velocity poten-tial, it can be convenient to approximate the fluid three-velocity measured by coordinate observers as homoge-neous , i.e., constant throughout the fluid. In particular,if the neutron stars are initially at apoapsis, with thecenter of each star located on the x -axis, then we assertthat each companion initially moves rigidly along the y -direction with instantaneous four-velocity field approxi-mated by u α = u t ( t α + v y y α ) . (4.7)Here, the parameter v y := u y /u t = dy/dt denotes the in-stantaneous three-velocity of a fluid element measured bya coordinate observer, which we approximate as constantthroughout the star. Note that, in general, this parame-ter has a different sign and magnitude for each star. Forirrotational incompressible flows, this quantity is exactly constant throughout the star in the Newtonian limit (cf.the discussion in the Appendix). For irrotational, rela-tivistic, compressible flows, this assumption is valid inan approximate sense, as illustrated in Fig. 1. A directcomparison with exact irrotational initial data (obtainedby solving for the velocity potential) shows that the pa-rameter v y is approximately constant to an accuracy of ∼
1% when the neutron stars are 67 (cid:39)
100 km apart(measured by the coordinate separation of their centers;recall that we take M (cid:12) = 1). This accuracy degrades to ∼
10% when the separation decreases to 31 . (cid:39)
47 km.The Lorentz factor in the above equation is determinedfrom the normalization condition g αβ u α u β = −
1, whichyields u t = [ − g tt − g ty v y − g yy ( v y ) ] − / . (4.8)The Lorentz factor u t is not assumed to be constant.We note that the constant three-velocity approxima-tion (4.7) is optional . It merely provides a way of easilyconstructing initial data for approximately irrotationalstars, on circular or eccentric orbits, using multigrid el-liptic solvers. One can still opt to solve Eq. (4.6) for theexact velocity potential with a solver that uses surface-fitted coordinates, using a method outlined at the end ofthe next subsection. x . . . . . . . . . . v y BAM , d = 31 . BAM , d = 67 . SGRID , d = 31 . SGRID , d = 67 . FIG. 1. Comparing equilibrium data constructed with theconstant fluid velocity approximation (using the implementa-tion of our method in bam ) to data constructed by solvingfor the exact velocity potential with sgrid . We show thefluid velocity along the x -axis for two different values of d ,the coordinate separation of the stars’ centers. One can see adeviation of ∼
10% for d = 31 .
2, but this deviation decreasesto ∼
1% for d = 67. We have cut the data at the surface ofthe star, denoted by the vertical thin dotted lines, since thevelocity is only well-defined in the star’s interior. C. Eccentric orbits & instantaneous helicalsymmetry
Binaries on eccentric orbits are not stationary in a ro-tating frame and thus lack helical symmetry. Hence, noself-consistent method for constructing initial data ex-isted for eccentric binaries to date, and numerical rela-tivity groups have resorted to using a superposition ofboosted spherical stars as initial data, possibly solvingthe constraints, but not solving the Euler equation. Useof such initial data entails not one, but three distinctphysical approximations regarding the fluid configura-tion: (i) that the stars are spherical and thus lack tidaldeformation, (ii) that the fluid is stationary in a linearlycomoving frame, and (iii) that the initial velocity fieldis homogeneous, i.e., constant throughout the fluid. Theabove approximations are only valid at infinitely largeseparation and are violated as the stars inspiral towardseach other. In addition, the metric construction as a lin-ear superposition of boosted spherical star data entails(iv) violations of the Einstein constraint equations.Thus, it remains an open question whether the f -modeoscillations observed in simulations of highly eccentric bi-nary neutron stars [20, 21] and black hole-neutron starbinaries [23, 24] are spuriously excited due to inconsis-tent initial data. In what follows, it will be demonstratedthat the above approximations are distinct and may infact be relaxed one by one. This allows one to examinewhether removing certain approximations removes spuri-ous oscillations in the early part of inspiral simulations.In particular, it will be shown that dropping assumptions(i), (iv) and retaining assumptions (ii), (iii), leads to os- cillations of the same magnitude as for boosted sphericalstars. On the other hand, the oscillations are nearly elim-inated if assumption (iii) is retained but assumption (ii)is dropped. That is, stationarity in a linearly comovingframe is the main source of error in boosted sphericalstar initial data. In the circular limit, the data shouldbe stationary in a rotating frame . One can quantify thediscrepancy between these two symmetries by testing theformer symmetry against exact circular initial data sta-tionary in a rotating frame; the discrepancy (10–20%)is comparable in magnitude to the central density oscil-lations in simulations. Thus, for eccentric binaries, themain difficulty is to define a suitable notion of stationar-ity, i.e., to obtain a generalization of the vector field (4.1)for which the energy (4.2) is approximately constant.We shall assume that the neutron stars are initially atapoapsis, with the center (i.e., point of maximum den-sity) of the star of mass m , m located respectively atposition x = a (1 + e ) m m + m + x cm , (4.9a) x = − a (1 + e ) m m + m + x cm (4.9b)on the x -axis, where x cm denotes the initial positionof the center of mass (on the x -axis). Here, a , b , and e = (1 − b /a ) / denote the semimajor axis, semiminoraxis, and eccentricity of the ellipse traced by the vectorjoining the two star centers and the binary’s center ofmass. We have selected the apoapsis for our construc-tion of initial data because this represents a moment of time symmetry of the radial motion, in the sense thatthe radial velocity vanishes. In addition, maximizing thedistance between the two stars happens to maximize theaccuracy of our approximations.We seek a vector field k α that approximately Lie-derives the flow. If such a vector exists, and the flowis irrotational, then the Cartan identity (4.3) will givean injection energy of the form (4.2) that is approxi-mately constant throughout the fluid. In light of Eq. (4.7)and the fact that the difference between a circular orbitand an eccentric orbit lies in the magnitude of the y -component of the initial 3-velocity at the apoapsis, our ansatz consists of generalizing the helical vector field(4.1) by adding a boost along the y -direction. Thisamounts to a change in the center of rotation of the heli-cal vector, and we shall use both viewpoints interchange-ably in what follows. These considerations lead us tointroduce what we term an instantaneously inscribed he-lical vector (or, more informally, a “helliptical” vector ), k α = t α + ω ϕ α + λ y α = t α + ω [( x − x c ) y α − y x α ] , (4.10) The parameter e estimated via this formula is only used as inputto monotonically control the eccentricity of the orbit obtainedupon evolving the initial data. Due to finite-size and relativis-tic effects, the actual orbits deviate from closed ellipses and theorbital eccentricity deviates from the Newtonian point particlelimit, cf. Secs. IV D and VII C. with a rotation frequency ω that now differs from the or-bital frequency. The displacement x c , or boost parameter λ = − ω x c , will be determined on physical grounds (with λ and x c different for each star). The choice ω = 0 and λ = v y would lead to k α = t α + v y y α , implying station-arity in a linearly comoving frame. As mentioned earlier,constructing initial data with this assumption yields spu-rious oscillations of similar magnitude to boosted spheri-cal stars and thus explains why assumption (ii) describedabove must be dropped. Instead, for eccentric binaries,the parameter ω should be nonzero and λ (or x c ) shouldvanish in the circular limit.The ω parameter can be determined by requiring thatEq. (4.3) holds exactly for incompressible binaries onNewtonian eccentric orbits (cf. the discussion in the Ap-pendix). This yields v y , = (1 − e ) ω ( x , − x cm ) . (4.11)The positions of the star centers, x , , are given by (4.9)and ω = (1 + e ) − (1 − e ) − / ¯Ω , (4.12)where ¯Ω = 2 π/T is the mean motion and T is the orbitalperiod. The λ or x c parameter can then be determinedby requiring k α to be initially parallel to u α at the starcenter. Substituting (4.9) into (4.10) and comparing to(4.7) yields λ , = − ω x c , = − v y , e − e − ω x cm , (4.13)with v y = v y , given by (4.11). As expected, for e = 0the inscribed helical symmetry vectors (4.10) yield initialdata for circular orbits, while e = 1 corresponds to zerotangential velocity v y , = 0 for a fixed ω , giving a head-on collision starting from rest. [Note that ω blows up as e (cid:37)
1. This can be seen in Eq. (4.11), which our methodreproduces in the Newtonian limit, as discussed below,even though we do not use that equation to determine ω .However, one can simply set ω = 0 and v y , = 0 in ourmethod to obtain an exactly head-on collision.]The parameters ω and λ in the inscribed helical sym-metry vectors in Eq. (4.10) have been determined such In particular, this choice would mean that ( ∂ t + v y ∂ y ) h vanishes.For circular orbits, one can check this against the exact enthalpyprofile, which satisfies ( ∂ t + Ω ∂ ϕ ) h = 0. Using this exact re-lation to eliminate the time derivative, the previous expressionbecomes ( ∂ t + v y , ∂ y ) h = Ω[( x , − x ) ∂ y + y∂ x ] h for star 1 , x = x , , y = 0. In addition, the force equation [ofthe form (4.20)] one obtains from the injection energy (4.2) thatis constructed from k α = t α + v y y α violates true force balanceat the center of the star. The violation stems from the absenceof centrifugal forces, which leave gravitational forces unbalanced.This results in density oscillations of order 10 −
20% in simula-tions, as illustrated in Fig. 4. xy a b r c x x c x cm FIG. 2. Illustration of the approximation of the orbits byusing circles inscribed into the orbital ellipse in a way thattheir curvature is the same as the one of the ellipse. Weshow the scaled semimajor axis a = a m M , semiminor axis b = b m M , and the radius r c = b /a and center x c of theinscribed circle as well as the center x of one star. The centerof mass is denoted by x cm . that we get the correct answer for irrotational incom-pressible stars moving on ellipses (see the Appendix).There is, however, a completely different way to obtainthis result. All we need are the following three reasonableassumptions: (i) There exists a vector k α that approxi-mately Lie-derives the flow. (ii) k α is along the motionof the star center. (iii) Each star center moves along asegment of an elliptic orbit at apoapsis.Assumption (ii) is absolutely necessary, otherwise k a can never be an approximate Killing vector. Assumption(iii) specifies what orbit we want. It seems reasonablethat at least approximately we should have Newtonianand thus elliptic orbits. Since we only need a small seg-ment of an orbit near apoapsis, we will approximate thissegment by the circle inscribed into the elliptical orbitthere (see Fig. 2), i.e., the circle that that has the samecurvature radius R c as the ellipse at apoapsis. Fromelementary considerations it is clear that R c has to be R c = (1 − e ) A for an ellipse with semimajor axis A andeccentricity e . In Newtonian theory it is well known thattwo particles of masses m and m that orbit around eachother, move on ellipses with semimajor axes a = d e and a = d e , where d and d are the distances of theparticles from the center of mass at apoapsis. Thus theradii of the inscribed circles have to be r c , = (1 − e ) d , . (4.14)Obviously these two inscribed circles are not centeredon the center of mass, but on the points x c , = x , ∓ r c , = x cm + e ( x , − x cm ) , (4.15)where we have used d , = | x , − x cm | and assumed thatapoapis occurs on the x -axis. (The upper and lower signscorrespond to the subscripts 1 and 2, respectively.) As-sumption (ii) then tells us that the approximate Killingvector must have the form k α , = t α + ω , [( x − x c , ) y α − y x α ] (4.16)near each star. While the expressions for k α and k α inEq. (4.16) look different, we will show next that ω = ω ,so that far from the stars (where x (cid:29) x c , ) there is onlyone approximate Killing vector. From the definition ofthe center of mass, we obtain m d = m d and thususing Eq. (4.14) we find m r c = m r c . (4.17)If we assume that the center of mass is at rest, momen-tum conservation demands that m ω r c = m ω r c . To-gether with Eq. (4.17), this implies that ω = ω =: ω. (4.18)Finally, in order to satisfy assumption (ii), v y in Eq. (4.7)must be chosen to be v y , = ± ωr c , = (1 − e ) ω ( x , − x cm ) . (4.19)The value of ω is usually derived from a “force balance”equation, which has the form of Eq. (4.5) applied at thecenter of each star. In the Newtonian limit, this equationreads ¯ D i E = 0 , (4.20)where the injection energy in star 1 is given by E = 12 v + Φ + h − ω [( x − x c ) v y − y v x ] (4.21)and Φ is the Newtonian gravitational potential. Here¯ D a is the covariant derivative compatible with the Eu-clidian 3-metric f ab in E . (We have used the samenotation for this covariant derivative as for the covari-ant derivative compatible with the conformal 3-metric,since we took the conformal 3-metric to be flat.) If thestars are far apart and thus almost spherical, then theorbits are almost elliptic and we can approximate thepotential due to star 2 as that of a point mass, given byΦ = − Gm / (cid:107) x − x (cid:107) , and neglect the gradient of star1’s potential at its center. (We show factors of G explic-itly here and during the rest of this discussion to makethe distinction between the gravitational and centripetalforces immediately apparent.) If we use this expressionin Eq. (4.20) at the star center x = x where the enthalpyis maximum, ∂ x h = 0, we find G m r − ωv y = 0 , (4.22)where r is the separation of the stars. Using v y = v y from Eq. (4.19) we find G m m r = m ω r c , (4.23) which implies that ω is such that the centripetal forceneeded to keep star 1 on the inscribed circular orbit isprovided by the gravitational force due to star 2.Thus the ansatz (4.10) admits a beautifully simple geo-metrical interpretation. The vector k α can be interpretedas an instantaneously inscribed helical vector field , thatgenerates time translations in a frame rotating about apoint ( x c , , ,
0) given by Eq. (4.15). When projectedonto a spatial slice Σ t , the integral curves of this vec-tor field are circles inscribed into the elliptic trajectoryat the apoapsis , with their center at ( x c , , ,
0) and a ra-dius r c , equal to the radius of curvature of the ellipse atthe apoapsis. Contrary to the circular case, the integralcurves of this vector field do not coincide with the (ec-centric) orbit of the star center (though they are tangentto each other at t = 0, which suffices for constructing ini-tial data), and the energy (4.2) is not conserved by theflow. Instead, as one may see by operating with k α ∂ α on(A16), this energy is constant spatially [Eq. (4.20) is sat-isfied throughout the fluid] at the initial time t = 0 butnot for t >
0. Because the spatial derivatives of the in-jection energy E approximately vanish initially, we inferthat its time derivative also vanishes instantaneously (i.e. ∂ t E = − k a , ∂ a E = 0 at t = 0). Eq. (2.20) is then satisfiedand, in this sense, k α , may be considered approximateKilling vectors for the region of spacetime occupied byeach star.In principle, if one wishes to construct initial data atany time (not necessarily at apoapsis), one could use amore general ansatz, k α = t α + ω [( x − x c ) y α − ( y − y c ) x α ],with the parameters ω , x c and y c determined by assum-ing that the orbit is instantaneously circular , i.e., by con-structing a circle inscribed at the elliptic trajectory at thepoint of interest. We have not tested this more generalconstruction, as nonvanishing radial derivatives and closeseparation are expected to degrade accuracy. Instead, forreasons detailed above, we opt to construct initial dataat apoapsis and set y c = 0. Note that the stars will likelyobtain a small initial radial velocity component when thedata are evolved, from the “kick” due to the initial spu-rious radiation. However, at the relatively large separa-tions we are considering, this radial velocity componentwill be small, so the stars’ orbit will still begin very closeto periapsis, thus retaining all of the favorable featureswhich led us to choose that point.With the ansatz (4.10), the first integral (4.2) to therelativistic Euler equation (2.17) becomes E = − h { u t + ω [( x − x c ) u y − y u x ] } , (4.24)where u µ = g µν u ν = u t ( g µt + v y g µy ) (4.25)if the constant 3-velocity approximation (4.7) is used, or u µ = h − ∇ µ Ψ (4.26)if one solves Eq. (4.6) for the exact velocity poten-tial. In this case, time derivatives are eliminated in fa-vor of spatial derivatives using the replacements ∂ t Ψ → −E − k a D a Ψ, ∂ t ρ → − k a D a ρ [with k a denoting thespatial part of the vector field (4.10) and D a the covari-ant derivative compatible with the physical 3-metric γ ab ],resulting in the equation [47] D a D a Ψ = ( β a + k a ) D a Λ − [ D a Ψ − Λ( β a + k a )] D a ln αρh , (4.27)where Λ := α − [ E +( β a + k a ) D a Ψ]. For fixed h and ρ ( h ),this equation is elliptic [the principal part is γ ab + (cid:96) a (cid:96) b , (cid:96) a := α − ( β a + k a ), which is clearly positive definite,since the 3-metric is positive definite]. It can be solvediteratively for Ψ using a numerical Poisson solver, withthe right hand side (including the additional terms oneobtains if one does not assume spatial conformal flatness)treated as a fixed source in each iteration, as in [64]. Theboundary condition[ D a Ψ − Λ( β a + k a )] D a ρ | S = 0 (4.28)is imposed on the surface of the star [47].Taking the Newtonian limit of the energy (4.24) (cf.the discussion in the Appendix) and applying Eq. (4.20)at the center of one star [given by (4.9) and defined asthe point of maximum specific enthalpy, ∂ x h | x = x , = 0]yields a force balance equation ∂ x E| x = x , = 0. It is re-assuring and straightforward to check that, for inversesquare gravitational forces, this equation amounts to Ke-pler’s third law for eccentric binaries¯Ω = G ( m + m ) a , (4.29)where a is the semimajor axis of the ellipse traced by thevector joining the two star centers. D. Radiation reaction & radial velocity
When the eccentricity parameter e is set to zero,the ansatz (4.10) reduces to the helically symmetricansatz (4.1). However, because this expression neglectsthe radial velocity due to radiation reaction, responsi-ble for binary inspiral, evolutions of helically symmetricdata exhibit residual eccentricity, e.g., the orbital sepa-ration acquires an oscillatory contribution and does notdecrease monotonically in time [84]. In light of the pre-ceding discussion, a way to incorporate radial velocityis to include a (constant as measured by coordinate ob-servers) three-velocity v x along the radial (initially x -)direction in u α and k α , thereby replacing Eqs. (4.7) and(4.10) by the ansatz u α = u t ( t α + v x x α + v y y α ) (4.30)and k α = t α + ω ( x − x c ) y α − ( ωy − v x ) x α . (4.31)The relation between ωx c and v y is again given byEq. (4.11), which guarantees that the condition u α = u t k α is satisfied at the center x = x , of each star. The radial velocity parameter v x can be determined,for example, by post-Newtonian or effective-one-bodytheory, with or without tidal (finite-size) corrections forthe given equation of state. This has been sufficient tosignificantly decrease the eccentricity in simulations inthe past (for example, the simulations performed in [85]used a Lorentz boost of the initial data based on point-particle post-Newtonian values for the radial velocity andled to significantly lower eccentricity than the simulationsin [86] that did not incorporate radial velocity). Alter-natively, if one wishes to obtain truly quasicircular data,then one may set e = 0 and an initial value for v x in theabove equations, evolve the resulting data for a periodof time adequate to determine the orbital eccentricity,adjust v x to reduce the eccentricity, and iterate this pro-cess until the resulting eccentricity is sufficiently small,similarly to what has been done for binary black holesimulations [87–91]. As stated earlier, the approximation (4.30) is meant tobe used when elliptic solvers without surface-fitted coor-dinates are employed. If surface-fitted solvers are avail-able, one may instead solve the elliptic equation (4.27)for the velocity potential, with k α given by (4.31). Ineither case, the approximate first integral to the Eulerequation is given by (4.2), and provides the basis for aself-consistent iteration that will be outlined below.If the stars are represented by compact monopolesources that inspiral towards their center of mass with atime-dependent radial velocity, then one can straightfor-wardly show that the vector (4.31) Lie-derives the Newto-nian gravitational field in a zone near each monopole. Inthis sense, the vector field is approximately Killing neareach compact star. Note that this is not true globallyand that the vector field is different for each star; thisdoes not limit our formulation since k α is merely usedfor the hydrodynamics inside each star and not for thegravitational field equations. V. NUMERICAL METHODA. Elliptic solver
To construct binary neutron star initial data (i.e., ex-pressions for the five elliptic quantities ψ , α , and β i andthe matter density profile, in our constant 3-velocity ap-proximation), we solve the five elliptic equations (3.11)together with the first integral of the Euler equation(4.24), where the latter has to be fulfilled throughoutthe stars. We solve these equations by iteration usinga self-consistent field method. In each step we approx-imate derivatives by standard second-order finite differ-encing operators in a full approximation storage (FAS) While the present manuscript was near completion, we were in-formed that recently this procedure was also implemented in [84]for binary neutron star simulations. olliptic code [93], and therefore yield higher con-vergence orders. Another way to increase the accuracyand obtain spectral convergence would be to implementthe method in sgrid .The iteration process is highly sensitive to the initialguess and cannot be started with an arbitrary set of val-ues. Solutions to Einstein’s equations for isolated non-rotating neutron stars are well known and will serve as aninitial guess for our iteration. We proceed as in [20]. Wefirst construct two single relativistic non-rotating spher-ical stars (solutions to the Tolman-Oppenheimer-Volkoff[TOV] equations [94, 95]) with the same baryonic massor central density (depending on what we fix during theiteration) as that desired for the stars in the initial data.We then boost these stars with a Lorentz transformationin the ± y direction to give the appropriate orbital motionand superpose the resulting 4-metrics by g (sup) µν = g (1) µν + g (2) µν − η µν , (5.1)where g ( A ) µν denotes the metric of star A ∈ { , } (in-cluding the boost) and η µν is the Minkowski metric. Weextract initial values for the elliptic quantities from thesuperposed 4-metric and also initialize the matter en-thalpy profile using the TOV solution. Since the spatialmetric is not conformally flat, due to the boost, we sim-ply take ψ to be the xx component of the 3-metric. Notethat one can usually take the boost to be zero withoutaffecting the convergence if the stars are widely separated(as in all the eccentric runs we show in Sec. VII). How-ever, a nonzero boost is necessary to obtain convergencewhen the stars are close (as in the d = 31 . ∼ ∼
10. We set Dirichlet bound-ary conditions for the elliptic variables ( ψ , α , and β i ) atour outer boundary by using values from the superposedTOV metric there, which we find produces better results(e.g., better agreement with sgrid at the boundary inthe quasicircular case) than just using the values thesevariables would have at infinity (i.e., α = ψ = 1 and β i = 0). B. Iteration scheme
In each iteration step we have to compute four con-stants from the integrated Euler equation (4.24) and its derivative with respect to x , which—evaluated at thestar’s center—yields the force-balance equation0 = { ∂ x u t + ω [( x − x c ) ∂ x u y + u y − y ∂ x u x ] }| x = x , . (5.2)The constants of interest are the orbital frequency ω , thecenter of mass x cm , and the injection energy of each star E , , which is given as the constant of integration. Weare free to make arbitrary choices for the central densityand the separation of the stars in advance and fix themthroughout the iteration. Note that it is also possibleto fix the total rest-mass instead of the central density,which we do when computing sequences. We evaluateboth the first integral to the Euler equation (4.24) andthe force balance-equation (5.2) at the centers of the twostars, which are located at fixed positions x , . For un-equal mass stars, we have to use a root finder inside ouroverall iterative scheme to obtain these constants, but forequal mass stars this system is degenerate and x cm cansimply be set to zero, allowing us to use algebraic solu-tions without an additional call to a root finder (thoughwe still solve by iteration overall).Using (4.25), the first integral can be written as E , = − hu t {− α + ψ [ β i β i + v , β y − ω y β x + ω ( x − e ( x , − x cm ) − x cm )( β y + v , )] } (5.3)and the force-balance equation can be rewritten in thesame way. The latter can be solved for ω algebraicallyif we set h (cid:48) ( x ) = 0, which ensures that the maximumdensity stays at the center of the star. (While the force-balance equation contains u t , which depends on ω , wesolve it with u t fixed, and then update u t later in theiteration.) The frequency ω can now be substituted into(4.24), along with value of the enthalpy h obtained fromthe fixed central density, to obtain the injection energy E .Using these values for ω and E , we are able to computethe enthalpy density profile by solving (5.3) for h . Wethen use the equation of state to obtain the mass density ρ ( h ). We finally update u t using (4.8) [and (4.25) tolower the index] and the new values of the constants andsolve the force-balance equation for ω again with the new u t , iterating over these steps until the change in ω fallsbelow numerical accuracy, which usually happens after afew iterations.Having gathered all the necessary constants, we cancompute the source terms using (3.8) and then use themultigrid scheme to solve the elliptic equations. In thisstep, we employ softening, i.e., instead of taking the fullvalue of the updated variable, we use a weighted aver-age of the old and new variables. Specifically, all ellipticvariables X are set using X = 0 . X new + 0 . X old , sim-ilar to the softening used in [31]. We then return to thecomputation of the constants, which are no longer valid,since the elliptic variables have now changed. This pro-cess is iteratively repeated until the change in the ellipticvariables falls below a prescribed threshold.2 VI. QUASICIRCULAR ORBITS
Here we give examples of quasicircular initial data con-structed with our method, and verify the code’s resultsagainst those obtained using sgrid . We summarize theproperties of the initial data sets we consider in this sec-tion and the following one in Table I, which also givesthe labels we use for the different sets.We compare the e = 0 limit of our data with heli-cal Killing vector initial data constructed using the spec-tral code sgrid [81], which solves for the velocity poten-tial using surface-fitted coordinates and also compacti-fies the grid to include spatial infinity. For the purposesof our comparison, we use n A = n B = 24, n φ = 8, n c = 16 points in sgrid (see [81] for details about sgrid ’s grid structure), which is sufficiently accurate,due to the code’s spectral convergence.We perform evolutions of the data using the bam code, which is a finite difference adaptive mesh refine-ment code for evolving the Einstein equations [82], andincludes a high resolution shock-capturing module tosolve the equations of relativistic hydrodynamics [83].Specifically, we use the same evolution setup as in [96],with the following differences: We use second-order spa-tial finite differencing for the geometry, consistent withthe order of the multigrid algorithm, and fourth-orderRunge-Kutta integration in time, along with fourth-orderKreiss-Oliger dissipation (with a factor of 0 . / ecc0 dataset. Therefore, we considerthree different resolutions, with a finest grid spacing of0 . ∼
500 with 5 levels of mesh refinement(where each level doubles the resolution), which is suf-ficient for our purposes. One can use more refinementlevels for highly accurate data intended for evolution andgravitational-wave extraction. Fig. 3 shows the expectedsecond-order convergence in a one dimensional compari-son of the momentum constraint. Here we plot the largestof the components, which is the y -component, since thestars are initially moving in the ± y -direction. The con-vergence behavior for the other components or the Hamil-tonian constraint is similar. At the surface of the stars wecan see some deviation from perfect convergence, includ-ing spikes in the constraint violations at the surface itself(shown in Sec. VII A in the eccentric case), which are cutoff here to show the central behavior in detail. These x − . − . − . − . . D y × D y . x D y . x × D y . x × FIG. 3. The y -component of the momentum constraint ( D y )for the ecc0 setup. We plot D y along the x -axis (which passesthrough the centers of both stars) for three grid spacings of∆ x = 0 . x , and ∆ x/ features are undesirable, but they are not easy to removein an implementation without surface-fitted coordinates.The convergence in the eccentric case is further detailedin Sec. VII A, which includes a discussion of the spikesand a comparison of the convergence of the Hamiltonianand momentum constraints.As mentioned earlier, using initial data inconsistentwith the hydrodynamic properties of the system can leadto spurious oscillations of the neutron stars, which wouldcontaminate the gravitational wave signal. As shown inFig. 4, the oscillations obtained in the evolution of the sgrid data (which solves for the velocity potential in ad-dition to solving the constraints) are negligible comparedto those obtained when evolving superimposed sphericalTOV stars, which exceed ∼ sgrid and subsequently evaluating the orbits and tun-ing the boost parameter to iteratively lower the orbitaleccentricity (via bisection). Note that it is possible toreduce the spurious oscillations even for simple superim-posed TOV data by changing the stars’ shapes, as foundby Tsatsin and Marronetti [99], who adjusted the coor-dinates, matter density, and velocity in an ad hoc buttunable way. This allowed them to reduce the oscilla-tions by an order of magnitude, even without explicitlysolving the hydrodynamic and constraint equations.On the other hand, solving the hydrodynamic and con-straint equations, but assuming stationarity in a linearlycomoving frame does not significantly reduce the spuri-3 TABLE I. Parameters for the initial data sets considered in this paper. Here m , denotes the baryonic mass of one of the stars(recall that we are only considering the equal-mass case in this paper), d denotes the initial coordinate separation of the stars’centers, e is the eccentricity parameter set in the initial data, and λ denotes the boost parameter used in the inscribed helicalsymmetry vectors. Additionally, κ is the scale parameter in the polytropic EOS, (∆ x ) min denotes the finest grid spacing, and“points” denotes the number of points used in each direction on each of l max refinement levels (as well as the fundamental gridlevel l = 0, giving l max + 1 levels total); the levels with l ≥ l mv are moving. (We do not give l mv for the sequence data sets seq0 – seq0.9 , which we do not evolve.) We name the sets using their eccentricity, with markers for the cases with a differentchoice of the boost parameter (“v”), data used for a sequence (“seq”), and an additional high-resolution case (“high”).Name m , d e λ κ (∆ x ) min points l max l mv ecc0 1.620 31 . − ωx c .
65 0.09375, 0.1875, 194, 98, 5 10.375 50ecc0v 1.620 31 . v y e . , ..., .
44] 0 , . , . , . − ωx c · · · ecc e . . , . , . , . , − ωx c
100 0.25 130 6 30 . , . , . . . − ωx c
100 0.125 258 6 3 t . . . . . . ρ m a x / ρ m a x , unsolved superpos SGRID ω = 0 , λ = − ωx c ω = 0 , λ = v FIG. 4. Comparison of the oscillations of the star, measuredusing the maximum density at each timestep, ρ max , normal-ized by the maximum density at t = 0, ρ max, 0 . We showevolutions of sgrid data (solid black) and the correspondingdataset constructed with our method assuming stationarity ina rotating frame ( ecc0 , blue dashed). One can clearly see theimprovement over the strongly oscillating curves of superim-posed boosted spherical stars (red dot dashed) and data com-puted using stationarity in a linearly comoving frame ( ecc0v ,green dotted). Note that the latter data set was only evolvedfor a short time, since we were only interested in the spuriousoscillations. ous oscillations. As a noteworthy caveat, one should keepin mind that the configurations that assume stationar-ity in a linearly comoving frame do not converge easily;one must use significant softening (over-relaxation) andcarefully adjust the order in which the equations are it-erated. Although the errors become smaller and the so-lution seems to converge after a few iterations, the errorseems to saturate and the solution tends to diverge af-ter a large number of iterations if the error tolerance issmall. Since in this case, as discussed in Sec. IV C, trueforce balance is lacking (a major source of instability,cf. footnote 2), we did not pursue this approach further. By employing assumptions consistent with the hydro-dynamic properties of the system, i.e., assuming sta-tionarity in a rotating frame at apoapsis (discussed inSec. IV C), the density oscillations exhibited in simula-tions were reduced by an order of magnitude, i.e., to ∼ VII. ECCENTRIC ORBITSA. Convergence
In this section, we perform evolutions of initial datasets constructed using the pair of inscribed helical sym-metry vectors (4.10) for nonzero eccentricity. We firstconsider initial data for the ecc0.915(high) cases (whoseparameters are given in Table I). The grid setup is arealistic one that could be used for production-qualityevolutions, with 6 levels of mesh refinement (beyond thecoarsest grid level) and an outer boundary at ∼ . Weconsider two resolutions: The lower resolution has 130 points in each refinement level and a finest resolution of0 .
25 (with 64 points across the star). The higher resolu-tion has 258 points in each refinement level and a finestresolution of 0 .
125 (with 128 points across the star). Ittook 0 . . ∼
800 simula-4 − . − . − . . . . H × H . H . ×
25 30 35 40 45 50 55 x − . − . − . . . . . D y × D y . D y . × FIG. 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 Hamil-tonian constraint ( H ) and the y -component of the momentumconstraint ( D y ) for two resolutions, with grid spacings of 0 . .
125 in the finest box and demonstrate second-order con-vergence of the constraints computed with sixth-order finitedifferencing stencils, the same order used in the evolution ofthis dataset we show later. (Recall that the initial data codeis only second order accurate.) tion time) was ∼
25 and ∼
50 times longer than thetime it took to solve for the initial data for the low- andhigh-resolution runs, respectively. For longer evolutionsthis ratio is even larger [e.g., the total wall clock time tomerger for the (low-resolution) ecc0.45 case is ∼ ecc0.915 case, while the initial data solvetook about the same time for both eccentricites].In Fig. 5, we show convergence of the Hamiltonianand momentum constraints in the finest box (surround-ing one star) along the axis passing through both stars’centers. (We chose this dataset with a relatively largevalue of e = 0 .
915 as a representative for all other ec-centricities. The Hamiltonian and momentum constraintfor other values of e are almost the same as for the caseat hand.) The convergence is clearly of second order, asexpected, apart from some spikes at the surface of thestar, which are also to be expected, since the density ofthe n = 1 polytropic stars we consider has a cusp atthe surface (i.e., it is not differentiable there). Moreover,we have chosen to compute the constraints here usingsixth-order finite differencing, the order we use in theevolutions (since the use of higher-order finite differenc-ing improves their accuracy, e.g., reducing the constraintviolations during the evolution), and one would expectthe use of higher-order finite differencing to amplify anysuch features. t . . . . . . ρ m a x / ρ m a x , ρ . ρ . × ρ . FIG. 6. The initial oscillations of a star in a binary with e = 0 . .
25 and 0 . ecc0.915 and ecc0.915high )to demonstrate how the density oscillations decrease with in-creasing resolution at approximately second order. The no-tation ρ . × denotes that the oscillations (not the totalmaximum density) are multiplied by 2 , as discussed in thetext. This increase of the finite difference order changes theshape of the constraint violations (cf. Fig. 3 and the bot-tom panel of Fig. 5—there is not much difference dueto the value of the eccentricity), in addition to changingtheir magnitude (increasing the maximum magnitude atthe center of the star for the Hamiltonian constraint, butdecreasing it for the momentum constraint). This differ-ence in shape is likely to be expected, since the remain-ders from second-order finite differencing are relativelylarge here. The use of higher-order finite differencing alsocreates spikes in the Hamiltonian constraint at the star’ssurface—there are only some slight wiggles present at thesurface when the Hamiltonian constraint is computed us-ing second-order finite differencing. However, increasingthe order of finite differencing decreases the size of thespikes at the surface of the star in the momentum con-straint.We also show the improvement of the spurious den-sity oscillations with increased resolution in Fig. 6. Weevolved the ecc0.915(high) initial data described above(now using sixth-order spatial finite differencing, as dis-cussed further in Sec. VII C) and monitored the max-imum density. We see a clear improvement in the os-cillations when doubling the resolution and the con-vergence can be estimated by multiplying the oscilla-tions by the appropriate scaling factor, i.e., consider-ing 4∆ ρ max + ρ max , , where ∆ ρ max = ρ max − ρ max , forsecond-order convergence with a factor of 2 differencein the grid spacing. Apart from some smaller superim-posed features, which can be seen, e.g., around t = 50 or t = 180, this scaling shows that the oscillations decreasewith increasing resolution with almost second-order con-vergence. Of course, we do not expect the oscillations to5completely converge away in the continuum limit, sincewe have still assumed spatial conformal flatness and haveneglected the radial component of the velocity from ra-diation reaction. However, we might expect that in thecontinuum limit these oscillations would be at the samesmall level seen for the sgrid quasicircular data in Fig. 4,in which case it makes sense to compute the convergenceorder assuming that the oscillations are zero in the con-tinuum limit. Indeed, this expectation is borne out bythe results shown in the figure and preliminary resultsfrom our implementation of the method in sgrid . B. Eccentric Sequences
As a check of our results, we compute constant-rest-mass sequences for equal mass stars of a fixed baryonicmass m b = 1 .
625 for varying eccentricity e (the seq0 – seq0.9 data sets in Table I); an isolated star with thatbaryonic mass has a gravitational mass of M i = 1 . E b = M ADM − M , where M = 2 M i and M ADM denotes the Arnowitt-Deser-Misner (ADM) mass, anasymptotic quantity that gives a measure of the totalmass of the spacetime. The ADM mass is defined via anintegral at spatial infinity (see, e.g., [100]), and thus isgenerally obtained by extrapolation in numerical codesthat do not use compactified coordinates to include spa-tial infinity on the grid (see, e.g., [82]). In our currentsituation, we found that the resolution of the outer gridswas insufficient to allow us to obtain accurate results fromextrapolation. We thus chose to obtain the ADM massfrom a single sufficiently large extraction radius (thoughnot too large, to avoid errors due to low resolution).Here we can use a simplified formula for the ADM massapplicable to our spatially conformally flat case, givenin Eq. (16) of [100], which gives significantly better re-sults with no extrapolation than the standard expression(given in, e.g., Eq. (7) of [100]). Specifically, in empiricaltests with sgrid data, we found that changing the ex-traction radius from r = 150 to r = 500 leads to a ∼ M ADM when using the standard expression,while this deviation is less than 0 .
01% when using the ex-pression that takes advantage of conformal flatness. Thisis to be expected, since (as discussed around Eq. (16)in [100]) the simplified expression can be evaluated atany radius in a region where the conformal factor satis-fies the Laplace equation, and the conformal factor in ourdata satisfies the Laplace equation to a good approxima-tion in the region in question [see Eq. (3.11a)], since thematter source is zero there, and the A ij A ij term will besmall (it falls off asymptotically as the shift squared, andthe shift goes to zero at infinity). The standard measurefor the ADM angular momentum is sufficiently accuratefor finite radii, so that we do not use extrapolation here,as well (the deviation caused by the change of extractionradius considered above is around 0 .
30 35 40 45 50 coordinate separation d . . . . . . ω e = 0 . e = 0 . e = 0 . e = 0 . .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . M ¯Ω0 . . . . . J / M .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . M ¯Ω − . − . − . − . − . − . . . E b / M FIG. 7. Sequences for equal-mass binary neutron starswith varying eccentricity e (data sets seq0 – seq0.9 in Ta-ble I). These sequences are computed with fixed baryonicmasses, yielding isolated stars with gravitational masses of M i = 1 . M = 2 M i . From top to bottom thequantities shown are the rotation ω as a function of the coor-dinate separation d of the two stars’ centers, and the ADM an-gular momentum J ADM and binding energy E b = M ADM − M as functions of the normalized mean motion M ¯Ω. The angu-lar momentum has been normalized by M , while the bind-ing energy is normalized by M . In addition, in the lower twoplots, the expected Newtonian behavior is plotted in dashedlines (black line instead of color scheme for the bottom plot,since the Newtonian prediction is independent of e ) and thepost-Newtonian (3PN) results in dotted lines for e = 0. monly supplied by plotting the dependence of the binding6energy E b and ADM angular momentum J ADM on theorbital frequency Ω. In the eccentric case, we will usethe mean motion ¯Ω, rather than the frequency ω whichappears within the two inscribed helical symmetry vec-tors (4.10). The reason is that it is ¯Ω, not ω , that satisfiesKepler’s third law (4.29) in the Newtonian limit. Specif-ically, increasing either ω or ¯Ω corresponds to decreasingthe binary’s initial coordinate separation, but only ¯Ω de-creases with increasing eccentricity e , as expected fromthe Newtonian limit, while ω ∝ / √ − e in the Newto-nian limit. The top panel of Figure 7 gives a plot of ω ,as it is the quantity used in the code, which illustratesthis increase with increasing e .Figure 7 shows these sequences and compares them tothose expected from Newtonian theory. Considering thequasicircular limit e = 0, we can see a qualitative agree-ment of the results. We also find the expected improve-ment in this agreement when we compare with the thirdorder post-Newtonian (PN) curve (errors of ∼ (See, e.g., Fig. 3 in [64] for comparison of a quasicircularsequence with 3PN predictions using a code with surface-fitted coordinates that solves for the velocity potential.)One sees that the angular momentum decreases as theeccentricity increases in the middle panel of Fig. 7. Thisbehavior is predicted by the Newtonian limit, as is the in-dependence of the binding energy E b on the eccentricityseen (approximately) in the bottom panel of that figure.Specifically, the Newtonian expressions in terms of ¯Ω are E b M = − η M ¯Ω) / , (7.1a) JM = η (1 − e ) / ( M ¯Ω) − / , (7.1b)where η := m m /M denotes the symmetric mass ratio,with M the total mass, so η = 1 / a = [(1 − e ) / (1 + e ) ] / ¯Ω into Eqs. (2.2)in Mora and Will [101]. The expression for the angularvelocity at aphelion, Ω a , in terms of ¯Ω comes from thestandard Newtonian expressions above Eqs. (2.2) in Moraand Will, noting that ¯Ω = M/a , where a is the binary’ssemimajor axis, by Kepler’s third law.It would be obvious to compare these results to thePN calculations from Mora and Will, too. However, While there are now 4PN results available for the energy and an-gular momentum in the quasicircular case [102], which are givenexplicitly in Appendix A of [84], we had initially intended to carryout the comparison with PN results in the eccentric case, as well,for which the 3PN computations in Mora and Will [101] seemedthe obvious choice. Moreover, given the effects of our constantthree-velocity approximation and finite boxes, comparisons withhigher-order PN would not necessarily be too illuminating. the PN results obtained in that manner are not well-behaved for head-on collisions. Since we approach suchconfigurations for increasing eccentricities, while Moraand Will approach an unbound parabolic orbit (with anonzero angular momentum), the comparison is not ap-propriate. While one could also consider the alterna-tive quasi-Keplerian parametrization of eccentric orbitsreviewed in Sec. 10 of [102], which remains well-behavedeven for vanishing angular momentum through 1PN, itis not clear how to relate the quantities used to describethe orbit in this parametrization to our ¯Ω variable usingonly data at apoapsis. Moreover, as noted by Sperhake etal. [103], the zero angular momentum limit is special evenin the Newtonian limit, since all three types of orbits(elliptic, parabolic, and hyperbolic) degenerate to head-on collisions when the angular momentum vanishes, soit is perhaps not surprising that post-Newtonian resultsbehave strangely there. (Sperhake et al. also give addi-tional caveats about comparing eccentric post-Newtonianresults with numerical relativity simulations.)
C. Trajectories and waveforms
To verify that the initial data obtained with ourmethod actually show the desired features for nonzeroeccentricities, it is useful to consider the star’s trajecto-ries. We define the trajectory of a star to be the co-ordinate position of the local minimum of the lapse ateach time step. In Fig. 8 we show that the eccentric-ity parameter e has the expected dramatic influence onthe orbits of the stars (as measured by their trajecto-ries). We choose a series of initial data sets with fixedcentral enthalpy h = 0 . . d = 80, but vary theeccentricity parameter. For the evolutions we again usethe bam code with the same settings given in Sec. VI,except that we now use sixth-order spatial finite differ-encing to increase accuracy, as was done for binary blackhole evolutions in Husa et al. [104]. Here we use eighth-order dissipation (with the same factor of 0 . et al. ’s choice of fourth-order dissipation (made forreasons of speed).Since increasing e yields a smaller tangential velocityof the star at apoapsis, the orbits become less circular,as the stars fall faster towards each other. In general,the number of orbits the stars perform before mergerwill decrease for eccentric orbits, since the configurationapproaches a head-on collision as one increases e . How-ever, it is possible to find interesting configurations wherethe stars undergo one or more encounters before merger,as found in [20]. We illustrate this with an evolutionwith two encounters before merger (the ecc0.45 case) inFig. 8, showing the gravitational waves for this case inFig. 9. For the present illustration, we have chosen to7present the waves extracted at a finite radius, though wehave checked that extrapolation to infinity (using radiifrom 300 to 900) only produces significant differences inthe high-frequency part of the merger signal, where thegrid spacing at the outermost extraction radii is likelytoo large to accurately transport waves with these highfrequencies.While it would seemingly be desirable to compare oneor more of the trajectories shown in Fig. 8 with a tra-jectory for the same case computed with superposeddata (cf. the trajectories shown in Fig. 1 of Gold etal. [20]), we do not do so, since the comparison might ac-tually be more confusing than illuminating. In particular,these are coordinate trajectories, and thus are a gauge-dependent quantity, and while the gauge conditions forboth evolutions are the same, the initial gauge is not,since we initialize the lapse and shift in both cases usingthe values given by the initial data. Indeed, the initialportions of the trajectories for the two cases look ratherdifferent, so while the qualitative zoom-whirl features inthe evolution are the same as found for the superposeddata, the initial portion of the trajectory for a run whichhas the same qualitative behavior of the trajectories andwaveform looks much more eccentric with the superposeddata than with the new data. (The Gold et al. [20] re-sults also show more eccentric tracks for situations forwhich the trajectory and waveform have the same quali-tative behavior, but here the primary difference is likelythat the initial coordinate separation of the stars is al-most twice as large as the one we consider here, while themasses are the same.)The gravitational waveform shown in Fig. 9 revealsthe same key features found in Gold et al. [20], specifi-cally the high frequency signals between the bursts frompre-merger encounters: We have also checked that thefrequency of these oscillations agrees with the f -modefrequency of an isolated nonrotating star with the samebaryonic mass and equation of state, as found by Gold etal. [We estimated this frequency using the fits givenin [105] and the values of 1 .
399 and 9 .
586 for the iso-lated star’s gravitational mass and areal radius; note thatEq. (6) in [105] contains a typographical error, which iscorrected in Eq. (14) of [106].]Since our method is only exact in the limit of New-tonian point particles, the eccentricity e evol of the orbitobtained when evolving the initial data is not expectedto be the eccentricity parameter e ID used in the con-struction of the data (i.e., the parameter we have beenreferring to as e so far). Therefore it would be usefulto have a way to determine the eccentricity obtained inthe evolution and relate this to the input eccentricity,allowing one to obtain a specific eccentricity, if desired.First, we will consider the different methods available todetermine e evol . While there is no known definition ofeccentricity in the comparable mass case in full generalrelativity, we can give a quantitative measure of the ec-centricity of a given evolution by fitting an ellipse to ashort section of the trajectory near the beginning of the − −
10 0 10 20 30 40 x − y e = 0 . e = 0 . e = 0 . e = 0 . FIG. 8. The trajectories of one star for different valuesof the eccentricity parameter e . The evolutions are basedon the initial data sets ecc0.45 , ecc0.5 , ecc0.73 , and ecc0.915 (see Table I) which are identical, except for the value of theeccentricity parameter e . While larger values of e lead toa rapid merger, for smaller values of the eccentricity (suchas e = 0 . et al. [20]. t − r − . − . − . . . . . R e ( r Ψ , ) e = 0 . FIG. 9. The gravitational waveform (in the form of the l = m = 2 mode of the Newman-Penrose scalar Ψ ) for the ecc0.45 case (see Table I). We extracted the waveform at adistance r = 500 from the binary’s center-of-mass, and shiftthe time axis by r = 500 to account (approximately) for thewaves’ travel time. The two small bursts (at t − r (cid:39) f -mode oscillations of the stars are visibleinbetween the bursts. orbit, as soon after the initial relaxation has completed aspossible. We expect that the trajectory will be most ap-proximately elliptical there, since radiation reaction andother strong-gravity effects will not have had much timeto affect the orbit. Besides fitting the track to an ellipsedirectly, it is also possible to measure the position angle ofthe trajectory and the proper distance between the starsand then fit an ellipse to the coordinates of the “properdistance trajectory” one obtains in this manner. Thisalternative method gives some indication of the extentto which the determination is contaminated by gauge ef-8 . . . . . . . e ID . . . . . . . e e v o l d coord d prop FIG. 10. The relation of the eccentricity parameter e ID in-put into the code (i.e., e in previous plots) and the “output”eccentricity of the evolutions e evol measured using two differ-ent fits to the initial portion of the trajectory, computed forthe ecc0.45 – ecc0.96 data sets in Table I. The fit to an ellipseusing the coordinate distance d coord gives quite similar resultsto the fit using the proper distance d prop and the trajectory’sposition angle. The black dot-dashed line shows the “ideal”relation e ID = e evol . The d coord and d prop curves agree wellenough that one can use them to obtain an estimate of theeccentricity of the system being evolved. fects (though the proper distance is not gauge invariant,since it is not computed along a geodesic). There arevarious other methods for eccentricity estimation in theliterature (see the references cited in [90], and also [84]for more recent work), but most of them are only appli-cable to small eccentricities and also often require severalorbits, so we do not consider them here.Note that while the Newtonian definition of eccentric-ity we use, e = (1 − b /a ) / , seemingly requires know-ing both the semimajor and semiminor axes of the orbit, a and b , this is not actually the case. The method ofperforming a fit to a segment of the orbit that we usehere does not require us to define either of quantitiesfor the orbit as a whole—for which they would likelybe undefined, or at least difficult to define, due to or-bital precession—but only to obtain them for the ellipsethat is the best fit to the segment of the orbit we con-sider, where they are given directly by the fit. Obvi-ously, the resulting eccentricity estimate should merelybe taken as a reasonable way of measuring the eccentric-ity which gives a qualitative feel for how eccentric theorbit is, rather than anything fundamental. In particu-lar, since this estimate uses the trajectories, which aregauge-dependent, one might obtain very different resultsfor this estimate if one evolved the same data using acode which used substantially different gauge conditions,e.g., the Princeton group’s code, which uses generalizedharmonic coordinates, and has been used to study highlyeccentric systems [21, 23, 24].We give the relation between the eccentricity in theevolution we measure using these methods ( e evol ) and theeccentricity parameter input to the code ( e ID ) in Fig. 10. From experimentation with the interval over which wefit the ellipse, we can estimate the error introduced byusing different intervals for the fit (varying both the lo-cation and size of the interval). Specifically, we vary theinterval length between 200 and 3000 points, which corre-sponds to coordinate displacements from ∼ . ∼ ∼
30 to ∼ ∼ .
4. Thisagreement is sufficient for the purpose of constructing ec-centric orbits, since it merely serves as a rough estimateof the expected eccentricity obtained in the evolution ofthe data. Note also that we have computed Fig. 10 forthe same choice of stellar masses as in Fig. 8 (i.e., usingthe ecc0.45 – ecc0.96 data sets in Table I), but the resultsone obtains for different stellar masses are very similar.In particular, a reduction of ∼
10% in the mass only re-sulted in a change of ∼ .
5% in e evol , independent of themethod used to compute it.Additionally, note that the method we have used todetermine the eccentricity is only applicable for e (cid:38) . e , one obtains inaccurate results due to radia-tion reaction. In particular one sees both a larger devia-tion of the eccentricity measured with the coordinate andproper distances for small e , as well as a noticeable off-set (of ∼ .
3) for quasicircular data (though this offset issmaller than the value of ∼ . e = 0). Such a large offset is not seen when one uses anyof the (previously mentioned) eccentricity determinationmethods that are specialized to small eccentricities.Altogether, Fig. 10 shows that the eccentricities weobtain in evolutions behave as expected (i.e., increasemonotonically as e ID increases). VIII. SUMMARY AND OUTLOOK
There are certain scenarios in which binary neutronstars can merge without having shed all of their eccen-tricity, e.g., due to dynamical capture, and simulations ingeneral relativity are the only way to model such mergersaccurately (necessary to study, e.g., their gravitationalwaves, ejecta, and merger remnants). In this paper wehave given the first method capable of providing consis-tent initial data for such systems (i.e., initial data thatsolves both the constraint equations of general relativityand the Euler equation). Our method proceeds by gen-eralizing the approximate helical Killing vector that isused to solve the Euler equation via its first integral inthe quasicircular case to a pair of inscribed helical sym-metry vectors (one for each star), which allows us to pro-vide initial data for binary neutron stars with arbitraryeccentricity. We find that the initial spurious oscillationsfound in evolutions of inconsistent data are reduced by9an order of magnitude or more (with higher resolution)using our consistent data, which assumes stationarity ina rotating frame, while the spurious oscillations remainif one assumes stationarity in a linearly comoving frame.We also find that the oscillations induced by the tidaldeformation at each close encounter are indeed physicaland not qualitatively altered compared to those observedin earlier evolutions of initial data that did not solve theEuler equation [20, 21].Considering the foundations of the method, we givetwo motivations for the inscribed helical symmetry vec-tors we introduce. In one derivation, we add a boostto the standard approximate helical Killing vector usedfor quasicircular initial data to adjust the binary’s ve-locity at apoapsis, allowing one to control the binary’seccentricity. In the other, more geometrical derivation,we show how the same vector arises from approximatingan elliptical orbit at apoapsis using an inscribed circle.We also show that the fact that there are two differentinscribed helical symmetry vectors (one for each star)does not spoil the derivation of the extended conformalthin-sandwich equations, by showing that one can obtainthese equations without assuming Killing symmetry byimposing spatial conformal flatness and maximal slicingfrom the outset. Additionally, we give a further exten-sion of the method that allows one to add radial velocityto the stars, so that (at least in principle) one can obtainconsistent binary neutron star initial data with arbitraryinitial tangential and radial velocities.For our first implementation of the method we intro-duced, we chose to use a second-order Cartesian multigridsolver (i.e., without surface-fitted coordinates), for sim-plicity and comparative speed. Without surface-fittedcoordinates, one cannot easily solve for the velocity po-tential, so we have taken the 3-velocity to be constant,which we show is a good approximation if the stars arenot too close. However, this is not a requirement of themethod, and one can easily solve the equation for thevelocity potential, as well, if one is using a code thatemploys surface-fitted coordinates (e.g., sgrid [81]).There are many potential extensions of the method, asdiscussed in the paper, and we are already in the pro-cess of implementing and testing a number of them. Afirst and straightforward step will be to construct binarieswith unequal mass stars. In addition, we are currentlyextending the sgrid code to use our inscribed helicalsymmetry vectors, in order to take advantage of sgrid ’sspectral accuracy and surface-fitted coordinates, whichwould allow us to solve for the velocity potential easily.The sgrid implementation will also naturally allow usto use more realistic EOSs, modeled as piecewise poly-tropes, which sgrid has recently been extended to han-dle. Furthermore, this implementation should make itpossible to add arbitrary spin to the stars, which mightbe especially pertinent for eccentric systems, as discussedin Sec. I.We have also implemented the generalized vector thatincludes radial velocity and now need to investigate the properties of the data we obtain from it. In particular,since the generalized vector allows one to modify both theradial and tangential velocity components of the stars, itshould allow us to obtain low-eccentricity initial data,similar to the work done for black hole-neutron star bi-naries in [77] and very recently for binary neutron starsin [84], or the various well-established methods for ec-centricity reduction for binary black holes [87–91]. Low-eccentricity binary neutron star initial data are particu-larly important from a gravitational wave data analysispoint of view: The residual eccentricity in current simu-lations is large enough to bias determination of the tidaldeformation [86], which would provide a valuable con-straint on the poorly-known equation of state of cold,dense nuclear matter.It may also be interesting to consider PN corrections tothe Newtonian expressions for the orbital motion used inderiving the inscribed helical symmetry vectors in orderto obtain sequences that yield better agreement with PNpredictions. This would also facilitate comparisons withanalytic techniques, such as PN or EOB formulationswith tidal corrections and eccentricity. We note, however,that neglecting PN corrections to the orbital motion inour construction mostly affects the relation between thevalue of the eccentricity parameter e used in the codeand the eccentricity obtained in the simulation, and notthe accuracy of the initial data or simulation: One canalways iterate over e to obtain any desired eccentricity inthe simulation.Even without the possible extensions of the method,one may already make certain useful investigations withthe current initial data. In particular, in the near future,we intend to revisit and extend the studies of Gold etal. [20] and East and Pretorius [21] (e.g., concerning prop-erties of the merger remnant and ejecta) in order to de-termine qualitative and quantitative changes upon usingimproved initial data, before going on to study more gen-eral scenarios (considering eccentricity reduction, addingspin, etc.). The ability to construct self-consistent initialdata for eccentric binary neutron stars opens the doorto studying many interesting physical situations, in boththe high- and low-eccentricity regimes, without signifi-cant limitations in accuracy. ACKNOWLEDGMENTS
We are grateful to K¯oji Ury¯u for his guidance andfor pointing out that initial data can be constructedwithout assuming equilibrium if the IWM approxima-tion and maximal slicing are employed. We also thankSarp Akcay, Leor Barack, Thomas Baumgarte, Sebas-tiano Bernuzzi, John L. Friedman, Eric Gourgoulhon,Alexandre Le Tiec, Marcus Thierfelder, and CliffordWill for fruitful discussions and suggestions. Addition-ally, we thank the anonymous referee for useful com-ments. We gratefully acknowledge support from theDFG SFB/Transregio 7 “Gravitational Wave Astron-0omy”, DFG Research Training Group 1523/1 “Quantumand Gravitational Fields”, STFC grant PP/E001025/1,NSF grant PHY-1305387, and the supercomputing grantfrom the John von Neumann Institute for Computing(NIC) provided on JUROPA at the J¨ulich Supercomput-ing Centre (JSC).
Appendix: Injection energy and velocity potentialfor nonrelativistic incompressible binaries
The nonrelativistic Euler equation for barotropic flowscan be written in the Crocco form [38] ∂ t v a + v b ( ¯ D b v a − ¯ D a v b ) = − ¯ D a H, (A.1)where H = v + h N + Φ is the Hamiltonian of a fluidelement with specific enthalpy h N , v a = dx a /dt is its ve-locity, Φ is the gravitational potential, and ¯ D a is the co-variant derivative compatible with the Euclidian 3-metric f ab in E . (As before, we have used the same notationfor this covariant derivative as for the covariant deriva-tive compatible with the conformal 3-metric, since wetook the conformal 3-metric to be flat.) If the flow is ir-rotational, v a = ¯ D a Ψ, then Eq. (A.1) has a first integral ∂ t Ψ = − H. (A.2)For incompressible flows, the specific enthalpy is givenby h N = (cid:82) dp/ρ = p/ρ , where p is the pressure, and themass density ρ is constant. Then, the continuity equation ∂ t ρ + ¯ D a ( ρv a ) = 0 (A.3)simplifies to a Laplace equation for the velocity potential,¯ D a ¯ D a Ψ = 0 . (A.4)In what follows, we construct analytic solutions toEqs. (A.2) and (A.4) for binaries on circular and eccentricorbits.
1. Circular orbits
For irrotational incompressible binaries on circular or-bits, all fluid elements move on circles with different cen-ters, but with the same radius R , and with the samespeed v = Ω R . Then, Eqs. (A.2) and (A.4) have theexact solutionΨ( t, r ) = −E t + v ( t ) · r = −E t − Ω R ( x sin Ω t − y cos Ω t )(A.5)where E is the injection energy (which is constant in bothspace and time) [47], v ( t ) = ˙ R = − Ω R (sin Ω t ˆ x − cos Ω t ˆ y ) = ∇ Ψ (A.6)is the fluid velocity and R ( t ) = R cos Ω t ˆ x + R sin Ω t ˆ y (A.7) is the position of the star’s center relative to the cen-ter of mass (assumed here to coincide with the origin).Here we use ∇ for the index-free version of ¯ D a . It isstraightforward to check that( ∂ t + Ω ∂ ϕ )Ψ = [ ∂ t + Ω( x ∂ y − y ∂ x )]Ψ = −E . (A.8)This equation can be interpreted as Eq. (A.2) trans-formed to a rotating frame. Alternatively, it may in-terpreted as a first integral of the equation( ∂ t + Ω L ϕ ) v a = 0 , (A.9)which follows from helical symmetry, i.e., stationarity ina rotating frame. The conserved injection energy followsfrom Eqs. (A.2) and (A.8) and reads E = H − Ω( xv y − yv x ) . (A.10)Taking the gradient of this equation and evaluating at thecenter of the star ( R, ,
0) at t = 0 gives a force balanceequation ∂ x E| x = R = ∂ x Φ | x = R − Ω R = 0 , (A.11)which yields Kepler’s third law for inverse square forces.
2. Eccentric orbits
We wish to generalize this derivation to eccentric bina-ries. In this case, the position of the stellar center relativeto the center of mass (assumed again to coincide with theorigin) is given by R ( t ) = [ a cos ζ ( t ) + ae ] ˆ x + b sin ζ ( t ) ˆ y , (A.12)where ζ ( t ) is the eccentric anomaly, related to the meananomaly ¯Ω t via the Kepler equation¯Ω t = ζ ( t ) + e sin ζ ( t ) . (A.13)Here, a , b , e and ¯Ω are the semi-major axis, semi-minoraxis, eccentricity, and mean motion of the orbit of onestar, respectively. For simplicity, we have chosen to studyan effectively one body problem by assuming an extrememass ratio, so the other, massive star (and thus the centerof mass) is at the origin, which is chosen to be the leftfocus of the ellipse. It is straightforward to relax theextreme mass ratio assumption and recover the two bodyequations, but we defer this until the end of this section.We have also assumed that the values ζ = 0 and ζ = π correspond to apoapsis and periapsis, respectively. TheKepler equation (A.13) has a series solution ζ ( t ) = ¯Ω te + 1 + ( ¯Ω t ) e e + 1) + O ( t ) . (A.14)The fluid velocity is homogeneous and given by v ( t ) = ˙ R = − ˙ ζ ( t )[ a sin ζ ( t ) ˆ x − b cos ζ ( t ) ˆ y ] = ∇ Ψ(A.15)1and the velocity potential is given byΨ( t, r ) = −E t + v ( t ) · r = −E t − ˙ ζ ( t )[ a x sin ζ ( t ) − b y cos ζ ( t )] . (A.16)If we operate on the above expression with ∂ t + k i ∂ i = ∂ t + ω∂ ϕ + λ∂ y , where k i is the spatial part of our in-scribed helical symmetry vector (4.10), and demand thatthe resulting expression be constant throughout the starat t = 0, i.e., ∇ E = 0, we obtain v = (1 − e ) ω a (A.17) with ω given by Eq. (4.12). If, in addition, we demandthat k i = vy i at the star center x = a (1 + e ) at t = 0, weobtain Eq. (4.13). It is straightforward to check that theforce balance equation ∂ x E| x = a (1+ e ) = 0 (A.18)applied to the star center for inverse square forces yieldsKepler’s third law for eccentric binaries. To recover thetwo-body equations, it suffices to rescale the ellipse bya factor depending on the mass of each companion, asindicated by Eq. (4.9). Then, Eq. (A.17) is replaced byEqs. (4.11), while e , ω , and ¯Ω remain unchanged. [1] J. A. Faber and F. A. Rasio, Living Rev. Relativ-ity , 8 (2012), arXiv:1204.3858 [gr-qc], [2] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel,K. Belczynski, C. Fryer, D. Holz, T. Bulik, and F. Pan-narale(2014), arXiv:1405.7016 [astro-ph.HE][3] D. Eichler, M. Livio, T. Piran, and D. N. Schramm,Nature , 126 (1989)[4] K. Postnov and L. Yungelson, Living Rev. Relativity , 3 (2014), arXiv:1403.4754 [astro-ph.HE][5] P. C. Peters, Phys. Rev. , B1224 (1964)[6] I. Kowalska, T. Bulik, K. Belczynski, M. Dominik,and D. Gondek-Rosinska, Astron. Astrophys. , A70(Mar. 2011), arXiv:1010.0511 [astro-ph.CO][7] R. M. O’Leary, B. Kocsis, and A. Loeb, Mon. Not. R.Astron. Soc. , 2127 (2009), arXiv:0807.2638 [astro-ph][8] W. H. Lee, E. Ramirez-Ruiz, and G. van de Ven, Astro-phys. J. , 953 (2010), arXiv:0909.2884 [astro-ph.HE][9] D. Tsang, Astrophys. J. , 103 (2013),arXiv:1307.3554 [astro-ph.HE][10] J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astro-phys. J. , 71 (2014), arXiv:1308.2964 [astro-ph.HE][11] F. Antonini and H. B. Perets, Astrophys. J. , 27(Sep. 2012), arXiv:1203.2938 [astro-ph.GA][12] S. Naoz, B. Kocsis, A. Loeb, and N. Yunes, Astrophys.J. , 187 (2013), arXiv:1206.4316 [astro-ph.SR][13] N. Seto, Phys. Rev. Lett. , 061106 (2013),arXiv:1304.5151 [astro-ph.CO][14] J. M. Antognini, B. J. Shappee, T. A. Thompson, andP. Amaro-Seoane, Mon. Not. R. Astron. Soc. , 1079(Mar. 2014), arXiv:1308.5682 [astro-ph.HE][15] F. Antonini, N. Murray, and S. Mikkola, Astrophys. J. , 45 (Jan. 2014), arXiv:1308.3674 [astro-ph.HE][16] W. E. East, S. T. McWilliams, J. Levin, and F. Preto-rius, Phys. Rev. D , 043004 (2013), arXiv:1212.0837[gr-qc][17] K. S. Tai, S. T. McWilliams, and F. Pretorius(2014),arXiv:1403.7754 [gr-qc][18] K. Kyutoku and N. Seto, Mon. Not. R. Astron. Soc. , 1934 (Jul. 2014), arXiv:1312.2953 [astro-ph.HE][19] N. Loutrel, N. Yunes, and F. Pretorius(2014),arXiv:1404.0092 [gr-qc][20] R. Gold, S. Bernuzzi, M. Thierfelder, B. Br¨ugmann,and F. Pretorius, Phys. Rev. D , 121501(R) (2012),arXiv:1109.5128 [gr-qc] [21] W. E. East and F. Pretorius, Astrophys. J. Lett. ,L4 (2012), arXiv:1208.5279 [astro-ph.HE][22] S. Rosswog, T. Piran, and E. Nakar, Mon. Not. R.Astron. Soc. , 2585 (2013), arXiv:1204.6240 [astro-ph.HE][23] B. C. Stephens, W. E. East, and F. Pretorius, Astro-phys. J. Lett. , L5 (2011), arXiv:1105.3175 [astro-ph.HE][24] W. E. East, F. Pretorius, and B. C. Stephens, Phys.Rev. D , 124009 (2012), arXiv:1111.3055 [astro-ph.HE][25] R. Gold and B. Br¨ugmann, Classical Quantum Gravity , 084035 (2010), arXiv:0911.3862 [gr-qc][26] R. Gold and B. Br¨ugmann, Phys. Rev. D , 064051(2013), arXiv:1209.4085 [gr-qc][27] W. E. East, F. M. Ramazano˘glu, and F. Pretorius, Phys.Rev. D , 104053 (2012), arXiv:1208.3473 [gr-qc][28] J. A. Isenberg, Int. J. Mod. Phys. D , 265 (2008),arXiv:gr-qc/0702113[29] J. R. Wilson and G. J. Mathews, Frontiers in Numeri-cal Relativity (Cambridge University Press, Cambridge,1989) ISBN 0521366666[30] W. Tichy, Phys. Rev. D , 024041 (2011),arXiv:1107.1440 [gr-qc][31] W. Tichy, Phys. Rev. D , 064024 (2012),arXiv:1209.5336 [gr-qc][32] J. W. T. Hessels, S. M. Ransom, I. H. Stairs, P. C. C.Freire, V. M. Kaspi, and F. Camilo, Science , 1901(2006), arXiv:astro-ph/0601337 [astro-ph][33] C. Kochanek, Astrophys. J. , 234 (1992)[34] L. Bildsten and C. Cutler, Astrophys. J. , 175 (1992)[35] T. M. Tauris, in Evolution of Compact Binaries , Astro-nomical Society of the Pacific Conference Series, Vol.447, edited by L. Schmidtobreick, M. R. Schreiber,and C. Tappert (2011) p. 285, arXiv:1106.0897 [astro-ph.HE][36] S. Bernuzzi, T. Dietrich, W. Tichy, and B. Br¨ugmann,Phys. Rev. D , 104021 (2014), arXiv:1311.4443 [gr-qc][37] [38] E. Gourgoulhon, EAS Publ. Ser. , 43 (2006),arXiv:gr-qc/0603009[39] N. Andersson and G. Comer, Living Rev. Relativity ,1 (2007), arXiv:gr-qc/0605010[40] J. M. Lattimer and M. Prakash, Astrophys. J. , 426(2001) [41] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Fried-man, Phys. Rev. D , 124032 (2009), arXiv:0812.2163[astro-ph][42] J. L. Synge, Proc. Lond. Math. Soc. , 376 (1937),reprinted in Gen. Relativ. Gravit. , 2177 (2002)[43] A. Lichnerowicz, Ann. Scient. ´Ecole Norm. Sup. , 285(1941)[44] A. H. Taub, Arch. Ration. Mech. Anal. , 312 (Jan.1959)[45] B. Carter, in General Relativity: An Einstein CentenarySurvey , edited by S. Hawking and W. Israel (CambridgeUniversity Press, Cambridge, England, 1979)[46] E. Gourgoulhon, C. Markakis, K. Ury¯u, andY. Eriguchi, Phys. Rev. D , 104007 (2011),arXiv:1101.3497 [gr-qc][47] J. L. Friedman and N. Stergioulas, Rotating Relativis-tic Stars (Cambridge University Press, New York, NY,2013) ISBN 978-0-521-87254-6[48] A. Lichnerowicz,
Relativistic hydrodynamics andmagnetohydrodynamics: Lectures on the ex-istence of solutions (W. A. Benjamin, NewYork, NY, 1967) [49] T. W. Baumgarte, G. B. Cook, M. A. Scheel, S. L.Shapiro, and S. A. Teukolsky, Phys. Rev. D , 7299(1998), arXiv:gr-qc/9709026[50] J. W. York, Phys. Rev. Lett. , 1350 (1999), arXiv:gr-qc/9810051[51] K. Ury¯u and Y. Eriguchi, Phys. Rev. D , 124023(2000), arXiv:gr-qc/9908059[52] H. P. Pfeiffer and J. W. York, Phys. Rev. D , 044022(2003), arXiv:gr-qc/0207095[53] X. Huang, C. Markakis, N. Sugiyama, and K. Ury¯u,Phys. Rev. D , 124023 (2008), arXiv:0809.0673 [astro-ph][54] T. W. Baumgarte and S. L. Shapiro, Numerical Rel-ativity: Solving Einstein’s Equations on the Computer (Cambridge University Press, Cambridge, 2010)[55] N. K. Johnson-McDaniel, N. Yunes, W. Tichy, and B. J.Owen, Phys. Rev. D , 124039 (2009), arXiv:0907.0891[gr-qc][56] G. Sch¨afer, in Mass Motion Gen. Relativ. Proc. CNRSSch. Orleans/France , edited by L. Blanchet, A. Spal-licci, and B. Whiting (2009) p. 167, arXiv:0910.2857[57] R. Rieth and G. Sch¨afer(1996), arXiv:gr-qc/9603043[58] G. Lovelace, Classical Quantum Gravity , 114002(2009), arXiv:0812.3132 [gr-qc][59] B. C. Mundim, B. J. Kelly, Y. Zlochower, H. Nakano,and M. Campanelli, Classical Quantum Gravity ,134003 (2011), arXiv:1012.0886 [gr-qc][60] G. Reifenberger and W. Tichy, Phys. Rev. D , 064003(2012), arXiv:1205.5502 [gr-qc][61] T. Chu, Phys. Rev. D , 064062 (2014),arXiv:1310.7900 [gr-qc][62] W. Throwe et al. (unpublished)[63] K. Ury¯u, F. Limousin, J. L. Friedman, E. Gourgoulhon,and M. Shibata, Phys. Rev. Lett. , 171101 (2006),arXiv:gr-qc/0511136[64] K. Ury¯u, F. Limousin, J. L. Friedman, E. Gourgoul-hon, and M. Shibata, Phys. Rev. D , 124004 (2009),arXiv:0908.0579 [gr-qc][65] M. Shibata (unpublished) [66] S. Bonazzola, E. Gourgoulhon, P. Grandcl´ement, andJ. Novak, Phys. Rev. D , 104007 (2004), arXiv:gr-qc/0307082[67] I. Cordero-Carri´on, P. Cerd´a-Dur´an, H. Dimmelmeier,J. L. Jaramillo, J. Novak, and E. Gourgoulhon, Phys.Rev. D , 024017 (Jan 2009), arXiv:0809.2325 [gr-qc][68] I. Cordero-Carri´on, P. Cerd´a-Dur´an, and J. M. Ib´a˜nez,J. Phys. Conf. Ser. , 012078 (2011)[69] I. Cordero-Carri´on, P. Cerd´a-Dur´an, and J. M. Ib´a˜nez,Phys. Rev. D , 044023 (2012), arXiv:1108.0571 [gr-qc][70] B. Garcia, G. Lovelace, L. E. Kidder, M. Boyle, S. A.Teukolsky, M. A. Scheel, and B. Szilagyi, Phys. Rev. D , 084054 (2012), arXiv:1206.2943 [gr-qc][71] Y. Zlochower, M. Ponce, and C. O. Lousto, Phys. Rev.D , 104056 (2012), arXiv:1208.5494 [gr-qc][72] F. Zhang and B. Szil´agyi, Phys. Rev. D , 084033(2013), arXiv:1309.1141 [gr-qc][73] R. H. Price, C. Markakis, and J. L. Friedman, J. Math.Phys. , 073505 (2009), arXiv:0903.3074 [astro-ph.SR][74] S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Phys.Rev. Lett. , 892 (1999), arXiv:gr-qc/9810072[75] S. Teukolsky, Astrophys. J. , 442 (1998), arXiv:gr-qc/9803082[76] E. Gourgoulhon(1998), arXiv:gr-qc/9804054[77] F. Foucart, L. E. Kidder, H. P. Pfeiffer, andS. A. Teukolsky, Phys. Rev. D , 124051 (2008),arXiv:0804.3787 [gr-qc][78] M. Shibata, Phys. Rev. D , 024012 (1998), arXiv:gr-qc/9803085[79] C. Markakis, Rotating and binary relativistic stars withmagnetic field , Ph.D. thesis, University of Wisconsin-Milwaukee (2011)[80] E. Gourgoulhon, P. Grandcl´ement, J.-A. Marck, J. No-vak, and K. Taniguchi, [81] W. Tichy, Classical Quantum Gravity , 175018(2009), arXiv:0908.0620 [gr-qc][82] B. Br¨ugmann, J. A. Gonz´alez, M. Hannam, S. Husa,U. Sperhake, and W. Tichy, Phys. Rev. D , 024027(2008), arXiv:gr-qc/0610128[83] M. Thierfelder, S. Bernuzzi, and B. Br¨ugmann, Phys.Rev. D , 044012 (2011), arXiv:1104.4751 [gr-qc][84] K. Kyutoku, M. Shibata, and K. Taniguchi(2014),arXiv:1405.6207 [gr-qc][85] J. S. Read, C. Markakis, M. Shibata, K. Ury¯u, J. D. E.Creighton, and J. L. Friedman, Phys. Rev. D , 124033(2009), arXiv:0901.3258 [gr-qc][86] J. S. Read, L. Baiotti, J. D. E. Creighton, J. L. Fried-man, B. Giacomazzo, K. Kyutoku, C. Markakis, L. Rez-zolla, M. Shibata, and K. Taniguchi, Phys. Rev. D ,044042 (2013), arXiv:1306.4065 [gr-qc][87] H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom,G. Lovelace, and M. Scheel, Classical Quantum Gravity , S59 (2007), arXiv:gr-qc/0702106[88] S. Husa, M. Hannam, J. A. Gonz´alez, U. Sperhake,and B. Br¨ugmann, Phys. Rev. D , 044037 (2008),arXiv:0706.0904 [gr-qc][89] W. Tichy and P. Marronetti, Phys. Rev. D , 024012(2010), arXiv:1010.2936 [gr-qc][90] A. H. Mrou´e and H. P. Pfeiffer(2012), arXiv:1210.2958[gr-qc][91] M. P¨urrer, S. Husa, and M. Hannam, Phys. Rev. D ,124051 (2012), arXiv:1203.4258 [gr-qc][92] W. H. Press, S. A. Teukolsky, W. T. Vetterling, andB. P. Flannery, Numerical Recipes: The Art of Scientific Computing , 3rd ed. (Cambridge University Press, NewYork, NY, 2007)[93] P. Galaviz, B. Br¨ugmann, and Z. Cao, Phys. Rev. D ,024005 (2010), arXiv:1004.1353 [gr-qc][94] R. C. Tolman, Phys. Rev. , 364 (1939)[95] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. ,374 (1939)[96] S. Bernuzzi, M. Thierfelder, and B. Br¨ugmann, Phys.Rev. D , 104030 (2012), arXiv:1109.3611 [gr-qc][97] S. Bernuzzi, A. Nagar, M. Thierfelder, andB. Br¨ugmann, Phys. Rev. D , 044030 (2012),arXiv:1205.3403 [gr-qc][98] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao,W. Tichy, and B. Br¨ugmann, Phys. Rev. D , 084057(2013), arXiv:1212.2901 [gr-qc][99] P. Tsatsin and P. Marronetti, Phys. Rev. D , 064060(2013), arXiv:1303.6692 [gr-qc] [100] N. ´O Murchadha and J. W. York, Phys. Rev. D , 2345(1974)[101] T. Mora and C. M. Will, Phys. Rev. D , 104021(2004), arXiv:gr-qc/0312082[102] L. Blanchet, Living Rev. Relativity , 2 (2014),arXiv:1310.1528 [gr-qc][103] U. Sperhake, E. Berti, V. Cardoso, J. A. Gonz´alez,B. Br¨ugmann, and M. Ansorg, Phys. Rev. D , 064069(2008), arXiv:0710.3823 [gr-qc][104] S. Husa, J. A. Gonz´alez, M. Hannam, B. Br¨ugmann, andU. Sperhake, Classical Quantum Gravity , 105006(2008), arXiv:0706.0740 [gr-qc][105] H. K. Lau, P. T. Leung, and L. M. Lin, Astrophys. J. , 1234 (2010), arXiv:0911.0131 [gr-qc][106] Y.-H. Sham, L.-M. Lin, and P. T. Leung, Astrophys. J.781