Hydrodynamics in full general relativity with conservative AMR
aa r X i v : . [ g r- q c ] M a y Hydrodynamics in full general relativity with conservative AMR
William E. East , Frans Pretorius , and Branson C. Stephens Department of Physics, Princeton University, Princeton, New Jersey 08544, USA. Center for Gravitation and Cosmology, University of Wisconsin-Milwaukee, Milwaukee, Wisconsin 53211, USA.
There is great interest in numerical relativity simulations involving matter due to the likelihoodthat binary compact objects involving neutron stars will be detected by gravitational wave obser-vatories in the coming years, as well as to the possibility that binary compact object mergers couldexplain short-duration gamma-ray bursts. We present a code designed for simulations of hydrody-namics coupled to the Einstein field equations targeted toward such applications. This code hasrecently been used to study eccentric mergers of black hole-neutron star binaries. We evolve the fluidconservatively using high-resolution shock-capturing methods, while the field equations are solvedin the generalized-harmonic formulation with finite differences. In order to resolve the various scalesthat may arise, we use adaptive mesh refinement (AMR) with grid hierarchies based on truncationerror estimates. A noteworthy feature of this code is the implementation of the flux correction al-gorithm of Berger and Colella to ensure that the conservative nature of fluid advection is respectedacross AMR boundaries. We present various tests to compare the performance of different limitersand flux calculation methods, as well as to demonstrate the utility of AMR flux corrections.
I. INTRODUCTION
The interface between strong field gravity and mat-ter dynamics promises to be one of the important fron-tiers in the coming years. A new generation of gravita-tional wave detectors (LIGO [1], GEO [2],TAMA [3], andVIRGO [4]) are now operational, and within the next fewyears are expected to reach sensitivities that will allowobservations of the Universe in gravitational radiation forthe first time. The prime targets of these observationsare compact object (CO) binaries composed of combi-nations of black holes (BHs) and neutron stars (NSs).Modeling of such sources is a crucial ingredient to real-ize the promise of gravitational wave astronomy. Even ifan event is detected with a high signal-to-noise ratio, re-constructing the dynamics of the system that producedthe signal cannot be done directly but instead will re-quire template banks of theoretical waveforms informedby numerical simulations.Compact object mergers involving NSs are expectedto be significant sources of not only gravitational radia-tion, but also possible progenitors for short-gamma-raybursts (SGRBs) [5–7] and other electromagnetic and neu-trino counterparts [8]. Efforts are already underway touse potential gravitational wave sources as triggers forsearches for electromagnetic transients [9, 10]. Observa-tions would help constrain evolutionary models for theprogenitor stars and their environments. Perhaps mostintriguingly, the observations would give clues to theequation of state (EOS) of matter at nuclear densities (asin NS interiors), which cannot be probed in laboratorieson Earth and is not fully understood at the theoreticallevel (for a broad discussion see for example [11]). Thereason that the gravitational wave signature could con-tain information about the matter EOS (and other detailsabout the internal structure of neutron stars) is that theEOS in general has a significant effect on the bulk motionof matter, and it is this bulk motion that is the mecha-nism by which gravitational waves are produced. Several studies to date have looked into this issue, suggesting theimprint of the EOS on the gravitational waves will bestrong enough to detect [12–23] (though, in some cases,the expected frequencies are higher than the range towhich the current generation of ground-based detectorsare most sensitive, thus limiting the information whichcan be extracted). While CO binaries containing NSsare a particularly interesting class of sources involvinggeneral relativistic (GR) hydrodynamics, they are by nomeans the only such systems. Examples of additionalsystems that have already been considered include BHaccretion tori [24–27] and NS-white dwarf mergers [28].Thoroughly modeling systems like those describedabove would require evolution of the spacetime, the pho-ton and neutrino radiation fields, and the magnetized,relativistic fluid. Even a minimalistic treatment, withthe Einstein equations coupled to the equations of rela-tivistic hydrodynamics, represents a complex, nonlinearsystem of partial differential equations. Numerical simu-lations are thus essential for exploring such strong field,dynamical systems. There is a long history of adaptingsuccessful techniques for simulating Newtonian hydrody-namics to relativistic and general relativistic fluids whichwe will not attempt to summarize (see [29] for a reviewof general relativistic hydrodynamics). Instead, we willbriefly attempt to place the code described in the presentpaper in the context of other recent codes developed forfluids on evolving spacetimes. Several of these codes [32–36] solve the field equationsin the BSSN formulation [37, 38]. The remainder [39, 40]use the generalized-harmonic formulation [41, 42] which Note that our focus is restricted to codes which handle dynam-ically evolving gravitational fields. Such codes, however, fre-quently owe much to earlier, fixed-background evolution codes(see [29]). In addition, advancements such as GR-hydro withmultipatch grids [30] and with GPUs [31] have recently beenmade with fixed-background codes. we also employ; unlike our code, however, these groupsconvert to a fully first-order formulation [43]. Mostgroups use finite-difference methods for the metric evolu-tion and a conservative, high-resolution shock-capturing(HRSC) scheme for the hydro evolution; these unigridalgorithms are then interfaced with some sort of adap-tive mesh refinement (AMR). A notable exception forthe metric evolution is [39], which employs pseudospec-tral methods for the metric and then interpolates to afinite-volume grid for the fluid.Some groups have implemented the MHD equationsin full GR; since these codes all make use of conserva-tive HRSC methods, they may be principally differenti-ated by how they meet the challenge of preserving the ∇ · B = 0 constraint. (A straightforward finite-differenceevolution of the magnetic field would generically lead tomagnetic monopoles and, hence, unphysical behavior.) WhiskyMHD employs constrained transport [34] for thispurpose, which preserves the constraint to machine ac-curacy, whereas the code of [44] uses hyperbolic diver-gence cleaning. Constrained transport, however, requiresspecial interpolation at refinement level boundaries in or-der to preserve the constraint. The Illinois group foundthat a vector-potential formulation of the MHD equa-tions works well when coupled to AMR [45]. This is be-cause the constraint is preserved by construction with thevector-potential, even with the restriction and prolonga-tion operations of AMR (see also [46] for a thorough ex-amination of the electromagnetic gauge condition). Stud-ies indicate that magnetic fields do not significantly affectthe gravitational dynamics of CO mergers (see e.g. [44]),but they could be critical for understanding EM counter-parts including the possible formation of a SGRB engine.A new method to treat the MHD equations was recentlypresented in [47], where ideal MHD is used in high mat-ter density regions ( e.g. inside a NS), while the force-freeapproximation is used elsewhere ( e.g. the magnetosphereof a NS). The authors applied the method to study thecollapse of magnetized hypermassive NSs (which couldbe formed via binary NS mergers) and suggested thatintense EM outbursts could accompany such events.Besides MHD, the other major advances in the phys-ical model for numerical relativity codes have been inthe arena of microphysics. While the Γ = 2 EOS wasthe community standard for quite some time, most codesnow allow for a nuclear theory-based EOS [48, 49] and/oruse various parametrized, piecewise polytropic EOSs in-spired by the range of plausible nuclear EOSs [50, 51].These advances in EOS description primarily affect thecold NS structure, but the group developing the
SACRA code has also begun to account for neutrino transportvia a simplified leakage scheme [12, 52]. The same grouphas also made available a formulation for a more accuratetruncated moment scheme with a variable Eddington fac-tor closure [53], which shows much promise for numericalrelativity simulations with neutrino physics beyond theleakage approximation.Another category of GR hydrodynamics codes employs the conformal-flatness approximation, which is particu-larly useful when supernova simulations are the targetapplication. An example is CoCoNuT/VERTEX, whichincorporates relativistic hydrodynamics, conformally flatgravity, and ray-by-ray neutrino transport [54]. The codeof [55] employs a similar scheme for hydrodynamics andgravity but adds a test magnetic field; this code has beenused to study the magnetorotational instability in super-novae.Newtonian (and semi-Newtonian) [56, 57], conformallyflat [58, 59], and fixed-background [60] SPH codes rep-resent an important, orthogonal approach to studyingCO interactions. SPH has an advantage over Eulerianschemes when a large range of spatial scales is involved.Such a situation may arise in CO mergers when materialis stripped from a star in a tidal interaction and formsan extended tail. On the other hand, Eulerian codes arethe standard approach when strong shocks are present,as would arise in binary NS mergers or disk circulariza-tion. (Recent progress has been made, however, in ap-plying SPH to situations with relativistic shocks [61].) Inaddition, SPH has not (to our knowledge) yet been cou-pled to a code which evolves the full Einstein equations.Nonetheless, comparisons between Eulerian and SPH re-sults could prove very useful on a problem-by-problembasis to characterize the errors in both methods.Though current efforts in GR simulations involvingmatter tend to focus on increasingly complex physicalmodels, there remain many unanswered questions in theastrophysics of compact objects that can be addressedwith a code which solves the Einstein equations coupledto perfect fluid hydrodynamics. We have thus focused ourcode development on hydrodynamics in full GR, whilemaintaining a flexible infrastructure to accommodate ad-ditional physics modules in the future. We evolve thefield equations in the generalized-harmonic formulationusing finite differences. The fluid is evolved conserva-tively using one of several different shock-capturing tech-niques we test here. We have also implemented the hy-drodynamical equations in a manner that is independentof EOS. We make use of AMR by dynamically adaptingthe mesh refinement hierarchy based on truncation er-ror estimates of a select number of the evolved variables.We also utilize Berger and Colella [62] style flux correc-tions (also known as “refluxing”) in order to make theuse of AMR compatible with the conservative nature ofthe hydrodynamic equations. Though AMR flux correc-tions have been implemented in other astrophysical hy-drodynamics codes (such as Athena [63], CASTRO [64],
Enzo [65], and FLASH [66]), to our knowledge this algo-rithm has not been used previously for hydrodynamicssimulations in full general relativity. A further notewor-thy feature of our implementation is that we store cor- Note that “flux correction” here refers to the enforcement of con-servation at AMR boundaries, not the recalculation of fluxes witha more dissipative scheme to preserve stability as in Athena [67]. rections to the corresponding fluid quantity integrated inthe volume of a given cell instead of the flux, allowingfor easy implementation within a computational infras-tructure that supports cell-centered but not face-centereddistributed data structures. The code described here hasrecently been applied to studying BH-NS mergers witheccentricity as may arise in dense stellar systems such asgalactic nuclear clusters and globular clusters [68, 69].In the remainder of this paper we outline our compu-tational methodology for simulating hydrodynamics cou-pled to the Einstein field equations and describe tests ofthis methodology. In Sec. II we review the generalized-harmonic approach to solving the field equations andpresent our methods for conservatively evolving a per-fect fluid coupled to gravity, including our method forinverting the conserved quantities to obtain the primi-tive fluid variables and the implementation of flux correc-tions to enforce the conservation of fluid quantities acrossAMR boundaries. In Sec. III we present simulation re-sults which test these methods, highlight the strengthsand weaknesses of various shock -capturing techniques,and demonstrate the utility of the flux correction algo-rithm.
II. COMPUTATIONAL METHODOLOGY
In this section we begin by explaining the basic equa-tions and variables we use to numerically evolve the Ein-stein equations in Sec. II A and then discuss the conserva-tive formulation of the hydrodynamics equations that weuse in Sec. II B. The evolution of conserved fluid variablesnecessitates an algorithm for inverting these quantities toobtain the primitive fluid variables which we present inSec. II C. Finally in Sec. II D we present the details ofour algorithm for AMR with flux corrections.
A. Solution of the Einstein equations
We solve the field equations in the generalized-harmonic formulation [41, 42] where we fix the coordi-nate degrees of freedom by specifying the evolution ofthe source functions H a := (cid:3) x a . In this form the evo-lution equation for the metric, g ab , becomes manifestlyhyperbolic: g cd ∂ c ∂ d g ab + ∂ b g cd ∂ c g ad + ∂ a g cd ∂ c g bd +2 H ( a,b ) − H d Γ dab + 2Γ cdb Γ dca = − π (2 T ab − g ab T ) (1)where Γ abc is the Christoffel symbol, T ab is the stress-energy tensor, and T is its trace. We evolve the metric,the source functions, and their respective time derivativesusing fourth-order Runge-Kutta where the spatial deriva-tives are calculated using fourth-order accurate finite-difference techniques. In other words, we have reducedthe evolution equations to first order in time so that there are 28 “fundamental” variables { g ab , H a , ∂ t g ab , ∂ t H a } ,but we directly discretize all first and second spatial gra-dients without the introduction of additional auxiliaryvariables.Analytically one can show [70] that if one begins withinitial data that satisfies the Hamiltonian and momen-tum constraints, initially set H a = (cid:3) x a , and then evolvethe metric according to (1) and the source functions ac-cording to some specified differential equations, then theconstraint equation H a − (cid:3) x a = 0 will be satisfied forall time. Numerically this statement will only be true towithin truncation error, which can grow exponentially inblack hole space times; to prevent this we add constraintdamping terms as in [71, 72]. In practice, ensuring that H a − (cid:3) x a is converging to zero for a given numerical sim-ulation run at different resolutions provides an excellentcheck that the numerical solution is indeed converging toa solution of the field equations.As described in [42], the computational grid we use iscompactified so as to include spatial infinity. This waywe can impose boundary conditions on the metric simplyby requiring that it be Minkowski. However we evolvethe metric of the uncompactified coordinates since thecompactified metric is singular at spatial infinity. B. Conservative Hydrodynamics
Coupled to gravity we consider a perfect fluid withstress-energy tensor T ab = ρhu a u b + g ab P , (2)where h := 1 + P/ρ + ǫ is the specific enthalpy and u a isthe four-velocity of the fluid element. The intrinsic fluidquantities ρ , the rest-mass density; P , the pressure; and ǫ , the specific energy are defined in the comoving frameof the fluid element. The equations of hydrodynamics arethen written in conservative form as follows [73]: ∂ t D + ∂ i ( Dv i ) = 0 (3) ∂ t S a + ∂ i (cid:0) √− gT ia (cid:1) = 12 √− gT bc ∂ a g bc (4)where v i is the coordinate velocity, g is the determinant ofthe metric, and the index i runs over spatial coordinatesonly. Note that (4) explicitly contains the time derivativeof the metric for index a = t . The conserved variables D and S a are defined as follows: D := √− gρu t (5) S a := √− gT ta (6)where D is simply the time component of the matter 4-current . In some implementations of the GR (magneto)hydrodynamicequations, see for e.g. [74], the analog of S t in (6) that is evolved In some situations we wish to perform axisymmetricsimulations where we use the symmetry to reduce thecomputational domain to two dimensions. We do thisusing a modification of the Cartoon method [75] as de-scribed in [42], where we take the x -axis as the axis ofsymmetry, and only evolve the z = 0 slice of the space-time. For the hydrodynamics this means that effectivelyeach fluid cell becomes a cylindrical shell, and we use thefact that the Lie derivative of the fluid fields with respectto the axisymmetric killing vector are zero to rewrite thecoordinate divergences in the above equations as ∂ i ( Dv i ) = ∂ x ( Dv x ) + 2 ∂ y ( yDv y ) (7)and similarly for ∂ i ( √− gT ia ) for the t and x components.For the y component there is an additional source term ∂ i ( √− gT iy ) = ∂ x ( √− gT xy ) + 2 ∂ y ( y √− gT yy ) − ( S z v z + √− gP/y ) . (8)By writing the y flux contribution in terms of ∂ y weensure that when we discretize our evolution will be con-servative with respect to the cylindrical shell volume el-ement. We choose a special form for the equation for S z : ∂ t S z + ∂ x ( √− gT xz ) + 2 y ∂ y ( y √− gT yz ) = 0 , (9)since in axisymmetry the quantity yS z is exactly con-served (that is, it has no source term).The conservative evolution system is solved numeri-cally using HRSC schemes. We briefly summarize the dif-ferent methods we have implemented and test in this pa-per, though the references should be consulted for morecomplete details. For calculating intercell fluxes we haveimplemented HLL [76], the Roe solver [77], and the Mar-quina flux [78] method. The HLL method is straight-forward to implement since it does not require the spec-tral decomposition of the flux Jacobian and is based onestimates for the largest and smallest signal velocities.The Roe solver works by solving the linearized Riemannproblem obtained using the flux Jacobian at each cell in-terface (using the so-called Roe average of the left andright states). The Marquina flux method is an extensionof this idea that avoids the artificial intermediate stateand switches to a more viscous local Lax-Friedrich-typemethod from [79] when the characteristic speeds change has the rest-mass density D subtracted off. This could provideimproved results in situations where the rest-mass density is or-ders of magnitude larger than the internal or magnetic energy,and accuracy in these latter quantities is important. Thoughwe have not explored this alternative, in the scenarios studiedhere (in particular since we are not looking at the behavior ofmagnetic fields) the added effect of a small amount of internalrelative to rest energy on the dynamics of the fluid or metricwill be negligible, and we expect either definition of S t to givecomparable accuracy results here. sign across the interface. Since the latter two methodsrequire the spectral decomposition of the flux Jacobian,we give it for our particular choice of conserved variablesin the Appendix. For reconstructing fluid primitive vari-ables at cell faces we have implemented MC and min-mod [80], PPM [81] , and WENO5 [83] , all of whichmay be used interchangeably with any flux method. MCand minmod are both slope limiter methods that reduceto linear reconstruction for smooth flows. Minmod isthe more diffusive of the two. In comparison, PPM andWENO5 are higher-order reconstruction methods. PPMis based on parabolic reconstruction with modificationsto handle contact discontinuities, avoid spurious oscilla-tions from shocks by reducing order, and impose mono-tonicity. WENO5 combines three different three-pointstencils with weights that are determined by a measureof the smoothness of the quantity being reconstructed.The specific fluid quantities that we reconstruct on thecell faces are ρ , u , and W U i , where u := ρǫ , W is theLorentz factor between the local fluid element and an ob-server normal to the constant t hypersurfaces, and U i isthe Eulerian velocity (the explicit form of which is givenin the following section). We choose to reconstruct W U i instead of simply U i since any finite value of this quantitycorresponds to a subluminal velocity.The fluid is evolved in time using second-order Runge-Kutta. Since the fluid is evolved in tandem with themetric, the first and second substeps of the fluid Runge-Kutta time step are chosen to coincide with the first andthird substeps of the metric time step. Since the spa-tial discretization of the fluid equations that we use isonly second-order we choose to use second-order timestepping for the hydrodynamics and we have not yetexperimented with higher-order methods. We still usefourth-order Runge-Kutta for nonvacuum metric evolu-tion (even though for evolutions with matter the overallconvergence rate will be no greater than second-order)both for convenience and because in vacuum dominatedregions we may expect some improvement in accuracy.For general relativistic hydrodynamics we evolve the fluidon a finite subset (though the majority) of the total grid(which as mentioned extends to spatial infinity throughour use of compactified coordinates), and at the outerboundary for the fluid we impose an outflow condition.Finally, as is common practice for this method of sim-ulating hydrodynamics, we require that the fluid densitynever drop below a certain threshold, adding a so-callednumerical atmosphere. We give this numerical atmo-sphere a spatial dependence that makes it less dense ap-proaching the boundaries and choose a maximum value In particular we use the reconstruction parameters presentedin [82]. Specifically, we perform reconstruction with the stencils andweights presented in Section A2 of [84]. Specifically we let ρ atm ( x c , y c , z v ) = ¯ ρ cos ( x c ) cos ( y c ) cos ( z c )where ¯ ρ is a constant, and ( x c , y c , z c ) are the compactified coor-dinates which range from -1 to 1. that makes it dynamically negligible (typically at least 10orders of magnitude below the maximum density). Theatmosphere is initialized using a cold equation of state( e.g. a polytropic equation of state). C. Primitive Inversion
The set of hydrodynamical equations is closed by anEOS of the form P = P ( ρ, ǫ ). While the conserved vari-ables S a and D are simply expressed in terms of fluidprimitive variables ( ρ , P , ǫ , and v i ) and the metric, thereverse is not true. This necessitates a numerical inver-sion to obtain the primitive variables following each up-date of the conserved variables. The method we use issimilar to the one used in [85] for spherical symmetry.First, we decompose the 4-dimensional metric into theusual ADM space plus time form ds = g ab dx a dx b = − α dt + γ ij ( dx i + β i dt )( dx j + β j dt ) (10)where γ ij is the spatial metric, α the lapse function and β i the shift vector. Then, from the metric and conservedvariables we construct two quantities, S := γ ij S i S j = γH W ( W −
1) (11) E := β i S i − S t = √− g ( HW − P ) , (12)where H := ρh and γ is the determinant of the spatialmetric. We reduce the problem of calculating the prim-itive fluid variables from the metric and conserved vari-ables to a one-dimensional root problem, where we beginwith a guess for H and iteratively converge to the correctvalue such that f ( H ) = 0 for some function. From (12)we can choose f ( H ) = E/ √− g − HW + P. (13)Note that given the metric and conserved variables, f ( H )is only a function of H , and can be computed as follows.First, calculate W = (1 + √ / S γH = W ( W − . (14)Then compute ρ and ǫ from ρ = D/ ( √ γW ) , (15)and ǫ = − H ( W − /ρ + W E/ ( Dα ) − , (16)respectively. Once ρ and ǫ are known, P can be obtainedfrom the equation of state, and then f ( H ) above. Aniterative procedure for solving f ( H ) = 0, where f ( H )is calculated as just described, thus gives the primitivevariables ρ , P , and ǫ . The three-velocity can then becomputed from U i = γ ij S j √ γHW , (17) where the Eulerian velocity U i is related to the gridthree-velocity through U i = ( v i + β i ) /α . This inversionscheme is implemented so as to allow any EOS of theform P = P ( ρ, ǫ ); thus, Γ-law, piecewise polytrope, andtabular equations of state such as the finite-temperatureEOS of Shen et al. [86, 87] (for a given electron fraction Y e ) are all supported.In practice we solve for f ( H ) = 0 numerically usingBrent’s method [88], which does not require knowledgeof derivatives and is guaranteed to converge for any con-tinuous equation of state as long as one begins with abracket around the correct solution. This can be usefulwhen dealing with equations of state interpolated fromtabulated values. One can avoid losing accuracy in the ul-trarelativistic and nonrelativistic limit by Taylor expand-ing the above inversion formulas (see [85]), for example,in 1 / Λ and Λ, respectively. We have implemented suchexpansions in our primitive inversion algorithm, thoughwe have not yet made any significant study of the inver-sion calculation in these regimes.In some cases the conserved variables will, due to nu-merical inaccuracies, evolve to a state that does not cor-respond to any physical values for the primitive variables.This causes the inversion procedure to fail. This can hap-pen in very low density regions that are not dynamicallyimportant but still must be addressed. We handle suchsituations using a method similar to that of [73] by ig-noring the value of S t and instead requiring the fluid tosatisfy a cold equation of state. D. AMR with flux corrections
Many of the problems we are interested in applyingthis code to involve a range of length scales, and in manycases we expect the small length scale features not tobe volume filling, for example the individual compactobjects in binary mergers. Such scenarios can be effi-ciently resolved with Berger and Oliger style adaptivemesh refinement [89]. A description of the variant of thealgorithm we use can be found in [90]; here we mentionsome particulars to this implementation, and give a de-tailed description of the extension to ensure conservationacross refinement boundaries.The computational domain is decomposed into a hi-erarchy of uniform meshes, where finer (child) meshesare entirely contained within coarser (parent) meshes.The hierarchy is constructed using (primarily) trunca-tion error (TE) estimates, which are computed within The initial bracket for the root finding is chosen by first checkingif [ H / (1 + δ ) , H (1 + δ )], where H is the value of H computedfor the primitive variables at the previous time step and δ > f ( H ). If it is not, as a failsafe we try successively largerbrackets with [ H / (1 + δ ) n , H (1 + δ ) n ] for n ≥ the Berger and Oliger time subcycling procedure by com-paring the solution obtained on adjacent levels of refine-ment before the coarser levels are overwritten with thesolution from the finer level. Typically we only use theTE of the metric variables, since fluid variables in generaldevelop discontinuities as well as turbulent features thatdo not follow strict convergence. The layout of the AMRhierarchy is then periodically adjusted in order to keepthe TE below some global threshold. In some situationswe also require that a region where the fluid density isabove a certain threshold always be covered by a mini-mum amount of resolution. This can be used to ensure,for example, that the resolution around a NS does nottemporarily drop below some level even if the TE of themetric variables in the neighborhood of the star becomessmall.When setting the values of the metric variables on theAMR boundary of a given child level we interpolate fromthe parent level using third-order interpolation in timeand fourth-order in space. For the cell-centered vari-ables, the outer two cells in each spatial direction (fora refinement ratio of 2) on a child level are initially setusing second order interpolation in time and space fromthe parent level. Following the evolution of the childlevel and flux correction applied to the parent level whenthey are in sync as described below, but before the cell-centered values on the child level are injected into theparent level, the values in the child boundary cells arereset using first-order conservative (spatial) interpolationfrom the parent level ( i.e. the value in the child cell isset to be the same as that of the parent cell in which thechild cell is contained). This ensures that the boundarycells on the child level are consistent with the correspond-ing flux-corrected cells on the parent level but does notaffect the order of convergence of the scheme since thesevalues are not used in the evolution step. During a regridwhen adding cells to the domain of a refined level we alsouse first-order conservative interpolation from the over-lapping parent level to initialize the values of the fluidvariables at new cells (fourth-order interpolation is usedfor the metric variables). Note that the actual domainthat is refined is larger than the volume where the TEestimate is above threshold by a given buffer in any direc-tion. The buffer size and regridding interval are chosen sothat if change in the region of high TE is associated withbulk motion of the solution (e.g. the NS moving throughthe domain), this region will never move by more thatthe size of the buffer between regrids. This ensures thatnew cells (for this kind of flow) are always interpolatedfrom regions of the parent that are below the maximumTE threshold. Thus, though the interpolation operationto initialize new cells is first-order, we find the error itintroduces is negligible ( i.e. , below the maximum desiredTE).AMR boundaries require special treatment in conser-vative hydrodynamics codes however, since the fluxesacross the boundary of a fine-grid region will not exactlymatch the corresponding flux calculated on the coarse- grid due to differing truncation errors. To enforce con-servation, we correct the adjacent coarse grid cells usingthe fine-grid fluxes according to the method of Berger andColella [62]. In the remainder of this section we reviewthe algorithm and outline our specific implementation.We will concentrate on the evolution of D on a 3-dimensional spatial grid, though the remaining conservedfluid quantities are treated the same way, and modifica-tion to different numbers of spatial dimensions is trivial.Equation (3) is evolved numerically at a given resolutionas D n +1 i,j,k = D ni,j,k − ∆ t [( F xi +1 / ,j,k − F xi − / ,j,k ) / ∆ x +( F yi,j +1 / ,k − F yi,j − / ,k ) / ∆ y +( F zi,j,k +1 / − F zi,j,k − / ) / ∆ z ] (18)where D ni,j,k is the volume average of D over the ( i, j, k )cell at time t = n ∆ t ; F xi +1 / ,j,k is the flux F x = Dv x through the ( i + 1 / , j, k ) cell face; ∆ x is the x length ofeach cell and so on for the y and z direction. In prac-tice the flux values will be calculated with some HRSCtechnique combined with Runge-Kutta, but the specificsare not relevant here. Now consider a situation with twosequential levels of refinement, L and L + 1, where level L +1 has a higher resolution with spatial refinement ratioof r in each direction, and its domain is a subset of level L . (In practice, we always take r = 2.) Here we focusthe discussion on a left boundary in the x direction, asillustrated in Fig. 1; boundaries along the right face andother coordinate directions are treated in a like manner.When evolving according to the Berger-Oliger algo-rithm, after each time step of length ∆ t is taken on level L , r steps of length ∆ t/r are taken on level L + 1. Thenthe results obtained on L + 1 are injected into level L where the levels overlap i.e. , the restriction operation isperformed conservatively by setting the value in the par-ent cell to the (coordinate) volume-weighted average ofthe child cells that make up the parent cell. Now on level L , the change in D due to flux going through the cell face( i L + 1 / , j L , k L ) on a timestep will be δD L = − ∆ t ∆ x F xi L +1 / ,j L ,k L ( t n ) . (19)On level L + 1, the change in D in one fine-level time stepdue to flux passing through one of the r cell faces thatmake up this same interface is δD L +1 ,j,k,m = − ∆ t/r ∆ x/r × F xi L +1 +1 / ,j L +1 + j,k L +1 + k ( t n + m ∆ t/r ) . (20)for j , k , and m ∈ { , , . . . , r − } . Now becauseof truncation error, in general the change in the net“mass” δM L := δD L V L within the coarse-level cell For the conserved fluid variable D which we focus on for speci- FIG. 1. A visualization of a refinement level boundary and its treatment in the flux correction algorithm. The top showscells in the x direction on refinement level L while the bottom shows equivalent cells for the L + 1 refinement level (here therefinement ratio is 2). Fluxes are symbolized by arrows. On the bottom level the blue cells (“type B” in the discussion in thetext) and those to the left on level L + 1 are boundary cells and will have their values set by interpolation from level L followingan evolution step on level L . Because of truncation error, subsequent evolution on level L + 1 will give a flux differing from thatcomputed on the parent level L . Consequently, when the new fine grid solution is injected back to the parent level (in cells tothe right of the red/dotted pattern cell), the solution about the boundary will no longer be consistent with the flux previouslycomputed there. The correct this, the fluid quantity in the red/dotted pattern cell is adjusted to exactly compensate for thedifference in flux computed between the coarse and fine levels. at ( i L , j L , k L ) computed with the coarse-level fluxeswill not equal the corresponding quantity δM L +1 := P j,k,m δD L +1 ,j,k,m V L +1 ,j,k,m computed with the fine-level fluxes, where V L is the coordinate volume of thecell ( i L , j L , k L ) and V L +1 ,j,k,m is the coordinate volumeof the cell ( i L +1 , j L +1 + j, k L +1 + k ). Thus, after the val-ues of D on level L + 1 are injected into level L (in cells( i L +1 , j L , k L ) and to the right in this example), the solu-tion on level L will suffer a violation of mass conservationproportional to δM L − δM L +1 . To restore the conserva-tive nature of the algorithm, the idea, described in detailbelow, is to adjust the conservative variable D in thecell ( i L , j L , k L ) post-injection by an amount to exactlycompensate for this truncation error induced difference.The scheme originally proposed in [62] is to definean array that keeps track of a correction to the fluxesthrough cell faces on level L that make up the boundaryof the evolved cells on level L + 1. Consider the casewhere ( i L + 1 / , j L , k L ) is such a face. This face-centered ficity, the value of the quantity integrated within the volume ofa cell in fact represents the rest mass in that cell. Through-out this section we will therefore use the term ‘mass’ to refer tothe value of a conserved fluid variable volume integrated withina cell, though for other conserved fluid variables this will notcorrespond to a physical mass. flux correction array, δF , is initialized with the inverseof the flux in (19), δF = − F xi L +1 / ,j L + j,k L + k , and thenduring the course of taking the r time steps on level L + 1receives corrections from the terms in (20) δF → δF + 1 r X j,k,m F xi L +1 +1 / ,j L +1 + j,k L +1 + k ( t n + m ∆ t/r ) . (21)After the cell values on level L are overwritten by theinjected values on level L + 1 where they overlap, thecells on level L that abut level L + 1 though are notthemselves covered by level L + 1 cells are corrected withthe flux stored in δF .The way we implement the flux correction algorithmis slightly different from this. In particular we wish toavoid the added computational complexity of implement-ing face-centered grid functions, and therefore we keeptrack of a cell-centered correction. The correction is thusalso more naturally represented as a correction to thefluid quantity integrated within the volume of the cell( e.g. for D the rest-mass) rather than a flux. Againreferring to Fig. 1, we define the first few cells at theboundary of level L + 1 as buffer cells since the calcula-tion of flux requires knowledge of the state on both sidesof the interface. These cells will have their values setby interpolation from those in level L . The innermostbuffer cells for the boundary on level L + 1 we call typeB cells (blue cells in the lower half of the figure). Theseare the cells where the level L + 1 contribution to themass correction will be stored. The cell on level L whichcontains the type B cell we will refer to as a type A cell(red, dotted-pattern cell). Type A cells are the ones thatreceive mass corrections in this algorithm. For each cellon each refinement level we use a bitmask grid functionthat indicates whether the cell is one of the above types(A or B), and if so which of the six possible faces (+ x , − x , + y , − y , + z , − z ) abut the boundary. For simplicityin the implementation we do not allow grid hierarchieswhere a cell would be both type A and type B .In the following we outline the additional tasks rel-ative to the basic Berger-Oliger algorithm that need tobe performed with our implementation of Berger-Colella.Following the spirit of these algorithms, we break downthe tasks into those the AMR “driver” code implements,which do not require knowledge of the specific equationsbeing evolved or what physical quantities the variablesrepresent, and conversely the “application” steps thatwould need to be implemented by a unigrid applicationcode plugging into the driver to become AMR-capable.The driver tasks include the following:(i) For the conserved fluid density D , allocate a storagegrid function to keep track of the associated masscorrection δM , i.e. the total correction to D withinthe volume of a given cell.(ii) Upon initialization set all correction arrays δM tozero, and compute the bitmask for the current re-finement hierarchy.(iii) After any regrid, recompute the bitmask array forthe new hierarchy.(iv) During the stage when buffer cells are set for vari-able D at interior boundaries on level L + 1 viainterpolation from level L , also interpolate the cor-rection variable δM , where the latter’s interpola-tion operator simply sets δM in a child cell to be1 /r that of the parent cell (for a three-dimensionalspatial grid).(v) Following injection of arrays D and δM from level L +1 to level L , where the injection operator for δM is an algebraic sum over child cells (a) zero all typeB cells in δM on level L +1, (b) call the applicationroutine (first item in the next list) to apply the masscorrections to D stored in the injected δM to typeA cells on level L , (c) zero all type A cells in δM on level L . In other words, an inner (non-physical) boundary on level L mustbe at least one cell away from any inner boundary on level L − The following are new tasks that the unigrid applicationcode needs to implement:(i) A routine that will add the mass corrections storedin δM to D for all type A cells on a given grid (i.e.,set D L → D L + δM/V L ) .(ii) When taking a single time step on a grid, for anycell marked type A, set δM to minus the changein mass of the cell from fluxes through cell facesindicated by the bitmask. For example, with thecase illustrated in Fig.1 and discussed above aroundEqs. (19) and (20), set δM L = − V L δD L .(iii) When taking a single time step on a grid, forany cell marked type B, add to δM the changein mass of the cell from fluxes through cell facesindicated by the bitmask. For example, with thesame example above, set δM L +1 ,j,k → δM L +1 ,j,k + V L +1 ,j,k,m δD L +1 ,j,k,m .For the GR-hydro equations we have five conservedfluid variables, D and S a . Though the latter do havenonzero source terms — since gravity can be a source(or sink) of energy-momentum — the above algorithmensures there will be no artificial loss/gain in the pres-ence of AMR boundaries due to truncation error fromthe advection terms. III. TESTS
In this section we present a number of tests of themethods presented above. We begin by demonstratingthe fourth-order convergence of the evolution of the Ein-stein equations for vacuum spacetimes before moving onto a number of flat space, relativistic hydrodynamics teststhat probe the treatment of fluid discontinuities. Weconclude with several tests of hydrodynamics in curvedspacetimes.
A. Vacuum evolution
In [42, 91] several tests of convergence of an earlierversion of the code (without hydrodynamics) were pre-sented. However, since then we have updated the evo-lution of the Einstein equations to fourth-order spatial Since we consider D a density and δM a mass, this requiresnormalization by the volume element V L , which the applicationknows. Note that in our code even though we have includedthe uncompactified metric volume element √− g in the definitionof the conservative variables and fluxes, compactification (andin axisymmetry, the cylindrical shell volume element) effectivelymakes the grid non-uniform and so the volume scaling is non-trivial. An alternative implementation could move this correc-tion step to the driver list of tasks, though then the applicationwould need to supply the driver with the array of local volumeelements. differencing and fourth-order Runge-Kutta time differ-encing, so we first show two vacuum tests: a Brill waveevolution [92, 93] and a boosted BH evolution.
1. Brill wave
For the Brill wave test we begin with initial data wherethe spatial line element is given by ds = ψ (cid:16) e B dx + e B y + z r dy +( e B − yzr ( dydz + dzdy ) + e B z + y r dz (cid:17) (22)where r = p y + z , B = 2 Ar exp( − ( r/σ r ) − ( x/σ x ) ),and the value of the conformal factor Ψ is determined bysolving the Hamiltonian constraint. We choose A = 40, σ r = 0 .
16, and σ x = 0 .
12. The initial data is chosen to betime symmetric ( ˙ γ ij = 0) and maximally sliced ( K = 0)with the conformal lapse ˜ α := Ψ − α = 1. The remainingmetric components are chosen to satisfy the harmonicgauge ( (cid:3) x a = 0). This describes a gravitational wavethat initially collapses inward before dispersing. In Fig. 2we show results from convergence tests in axisymmetry atthree resolutions where the medium and high runs had,respectively, 1.5 and 2 times the resolution of the lowrun. The constraint equations ( H a − (cid:3) x a = 0) as well asthe metric components show the expected fourth-orderconvergence.
2. Boosted BH evolution
As an additional vacuum spacetime test we evolved aboosted BH in three dimensions. We began with initialdata describing a BH in harmonic coordinates [94] withboost parameter v = 0 .
25. As described in [42], duringthe evolution we avoid the BH singularity by searchingfor an apparent horizon and excising a region within.To demonstrate convergence we performed this simula-tion at three resolutions, the lowest of which has approx-imately 30 points covering the diameter of the BH. Themedium (high) resolution has 1.5 (2.0) times the numberof points in each dimension, respectively. For all res-olutions we used the same AMR hierarchy, determinedbased on truncation error estimates at the lowest reso-lution, with six levels of 2:1 refinement. In Fig. 3 wedemonstrate that the constraint equations are convergingto zero at fourth-order. When hydrodynamics is includedthe theoretical limiting convergence rate of our code willdrop to second-order (in the absence of shocks). Howeverin vacuum dominated regions, for example the gravita-tional wave zone, one can expect that for the finite res-olutions we can practically achieve the convergence willbe somewhere between second- and fourth-order. l n ( || C a || M ) Low res.Med. res. × (1.5) High res. × (2) −1.2−1−0.8 g tt −4 time (M) ∆ g tt Med−Low res.(High−Med.res) × FIG. 2.
Top : The natural log of the L norm of the con-straint violation, C a := H a − (cid:3) x a , for a Brill wave evolution( i.e. natural log of qR | C a | d x/ R d x ). The three reso-lutions shown are scaled assuming fourth-order convergence.Time is shown in units of M , the total ADM mass of the space-time, and the constraints are multiplied by M to make themdimensionless. The lowest resolution has a grid spacing of h = 1 . M . Middle/Bottom : The value of the metric com-ponent g tt evaluated at ( x, y, z ) = (0 , M,
0) (middle) andthe difference in this quantity between low and medium reso-lution, and medium and high resolution (bottom), the latterscaled so that the two curves should coincide for fourth-orderconvergence.
B. Relativistic hydrodynamic tests in flatspacetime
We have performed a number of standard tests for rel-ativistic, inviscid hydrodynamics that probe how well agiven numerical scheme handles the various discontinu-ities that arise. The best combination of reconstructionand flux calculation methods depends on the problem0 −5 time (M BH ) || C a || M B H Low res.Med. res. × High res. × FIG. 3. The L -norm of the constraint violation ( C a := H a − (cid:3) x a ) in the equatorial plane for a boosted BH simula-tion with v = 0 .
25. The three resolutions shown are scaledassuming fourth order convergence. Time is shown in unitsof M BH , the ADM mass of the BH in its rest frame, and thenorm of the constraints is multiplied by M BH so as to makeit dimensionless. under consideration. We have thus implemented severaloptions and maintained a modular code infrastructureso that they are readily interchangeable and upgradable.While strong shocks such as the ones considered here arenot expected to play an important dynamical role in bi-nary BH-NS mergers, they might be important in otherpotential applications of interest (such as NS-NS grazingimpacts, or understanding EM emission from collisions).Thus, the ability to tailor the reconstruction and fluxmethods to the problem at hand may prove importantin the future. In this section, we closely follow the se-quence of tests used in the development of the RAM codeof Zhang and MacFadyen [95], so that our results maybe compared with theirs. Though they focus on more so-phisticated flux-reconstruction algorithms, their simplermethods (labeled U-PPM and U-PLM, denoting recon-struction of the unknowns with piecewise parabolic andlinear reconstruction, respectively) are comparable to theones we employ.
1. 1D Riemann problems
We first present a series of four relativistic, one-dimensional (1D) Riemann problem tests for which theexact solution is known (see Sections 4.1-4.4 of [95]). Inall cases, the domain is x ∈ [0 ,
1] and there are initiallytwo fluid states, a left and a right, initially separated byan imaginary partition at x = 0 .
5. At t = 0, the parti-tion is removed and the fluid evolves to some new state. A Γ-law EOS is used for all the tests. In Table I wesummarize the initial states and adiabatic indices usedfor the four tests, which we label as RT1 (Riemann Test1), RT2, RT3, and TVT (Transverse Velocity Test). Wecompare the performance of the various combinations ofreconstruction schemes and flux methods to the exact so-lution and summarize the errors and convergence rates inTable II. Exact solutions to these four tests were gener-ated using a solver provided by B. Giacomazzo, which isdescribed in [96]. Taking HLL as our basic flux method,we performed this series of Riemann problem tests withfour reconstruction methods: MC, minmod, WENO5,and PPM. For MC and WENO5, we also explored theeffect of the flux method by running the tests with theRoe solver and Marquina’s method. Most cases have aCourant-Friedrichs-Lewy (CFL) factor of 0.5. However,the Roe solver, when combined with WENO5, does notseem to work well for problems with very strong shocks,such as RT2 and TVT. For a CFL factor of 0.5, we ob-tain acceptable results with Roe only by using a morediffusive limiter (like MC). For RT2 and TVT, we thususe Roe combined with WENO5 with a CFL factor of0.1.All of the methods we considered perform well on RT1,which is a fairly easy test. The lowest overall error occursfor WENO5 reconstruction (though the density profilebetween the shock and the contact discontinuity seemsnot to be as flat as in the other cases). The overall suc-cess of WENO5 may be due to the fact that the shock isrelatively mild and there is an extended rarefaction thatbenefits from the high-order reconstruction. In Fig. 4 wecompare the density profile obtained using HLL and var-ious reconstruction methods to the exact solution. Wenote that the tests which used the Roe or Marquina fluxcalculation with WENO5 do not have the oscillation vis-ible in the plot around x = 0 . Test Γ a P L ρ L v L P R ρ R v R RT1 5/3 13.33 10.0 0.0 10 − − . , . b TABLE I. The initial left and right states for the 1D Riemann problems. a Adiabatic index of EOS b In this case v x =0 but the transverse velocity v y = 0 .
99 is nonzero.Reconstruction Flux method RT1 RT2 RT3 TVTError a Convergence b Error Convergence Error Convergence Error ConvergenceMC HLL 0.034 0.82 0.110 0.59 0.062 0.77 0.238 0.72Roe 0.032 0.82 0.110 0.60 0.052 0.80 0.233 0.72Marquina 0.036 0.82 0.127 0.59 0.056 0.79 0.227 0.76Minmod HLL 0.061 0.86 0.169 0.42 0.054 0.71 0.395 0.76WENO5 HLL 0.033 0.84 0.093 0.76 0.039 0.61 0.191 0.83Roe 0.032 0.85 0.096 0.79 0.039 0.60 0.198 0.81Marquina 0.036 0.85 0.093 0.76 0.038 0.66 0.183 0.82PPM HLL 0.041 0.88 0.133 0.67 0.024 1.01 0.248 0.78TABLE II. 1D Riemann test results. a The L1 norm of the error for resolution N = 400. b The average convergence rate between runs with N = 200 , weak solution of the equations when discontinuities arepresent.)For the transverse velocity test, the initial data areset up as in RT2, except that there is a transverse ve-locity v y = 0 .
99 on the right side of the partition. Thestrong shock propagates into the boosted fluid, and thestructure of the shock is altered, since the velocities inall directions are coupled through the Lorentz factor [98].Again, the reconstruction technique influences the resultmore than the flux calculation. For WENO5 reconstruc-tion, the errors for HLL, Roe, and Marquina are all veryclose in magnitude. WENO5 and PPM yield the bestresults overall. In Fig. 7 we show the density profile atdifferent resolutions for HLL combined with WENO5.
2. 1D Shock-heating problem
We next consider a one-dimensional shock-heatingproblem as in [95], which tests a code’s conservation ofenergy as well as the ability to handle strong shocks.For this problem, the computational domain is x ∈ [0 , x = 1. The fluid movestoward this boundary with an ultrarelativistic initial ve-locity of v = 1 − − . The fluid has an initial density of ρ = 1 . ǫ = 0 . /
3. When the relativistic fluid slams into the wall, its kineticenergy is converted into internal energy behind a shockwhich propagates to the left. Because the fluid is initiallycold, essentially all of the heat is generated through thisconversion. As explained in [95], the shock speed andthe compression ratio of the shock (or equivalently, thepost-shock density) is known analytically. We evaluateour errors by calculating the L1 norm of the density er-rors on the entire computational domain. The averagerate of convergence is also calculated using this measureof error.We performed this test using HLL with five differentreconstruction methods at four different resolutions (200,400, 800, and 1600 zones). Results are shown in Ta-ble III. We find that, due to the extremely strong shock,there is a tendency for post-shock oscillations to formwith less diffusive reconstruction schemes (see Fig. 8).The WENO5 solution is afflicted with severe post-shockoscillations and exhibits poor convergence to the exactcompression ratio. Very diffusive reconstruction schemes(zero slope and minmod) are comparatively quite success-ful and converge rapidly to the exact compression ratio.PPM, with its flattening step, gives the best convergencerate overall.2 x d e n s it y minmodMCWENO-5PPMexact Riemann Test 1
FIG. 4. Density at t = 0 . N = 400. The inset shows the shock and contactdiscontinuity. The exact solution was generated using thecode of [97]. x d e n s it y Riemann Test 2
FIG. 5. Density at t = 0 .
3. Emery step problem
Next we consider the two-dimensional (2D) Emery stepproblem [81, 99], with the setup as in [95]. In this sce-nario, a fluid flows through a wind tunnel at relativisticspeed and hits a step, which is represented by a reflect-ing boundary condition. The computational domain is( x, y ) ∈ [0 , × [0 , − [0 . , × [0 , .
2] where the sub-tracted region represents the step. At the left boundary, x d e n s it y minmodMCWENO-5PPMexact Riemann Test 3
FIG. 6. Density at t = 0 . N = 400. The post-shock oscillations arelargest in the MC case and smallest for PPM. x d e n s it y exact2004008001600 Transverse Test
FIG. 7. Density at t = 0 . inflow conditions are enforced (as in the initial data),while at the right, outflow conditions are enforced. Allremaining boundaries are reflecting. The fluid is initial-ized with density ρ = 1 .
4, velocity v x = 0 . . Reconstruction Error a Convergence b ZERO 950 0.94Minmod 801 0.92MC 1500 0.85PPM 824 0.96WENO5 1670 0.53TABLE III. Shock -heating test results. For this test we alsocompare to zero slope reconstruction, labeled “ZERO.” a The L1 norm of the error for resolution N = 400. b The average convergence rate between runs with N = 200 , d e n s it y ZEROMINMODMCPPMWENO5
Shock Heating Test
FIG. 8. Density at t = 2 . the fluid reflects off the step is distorted by large am-plitude post-shock oscillations. These propagate down-stream, rolling up the boundaries between the differentsolution regions. The higher resolution runs with MCalso have these features, but at shorter wavelengths andlower amplitude. PPM and WENO5 reconstruction per-forms much better, and these results are shown at tworesolutions in Fig. 9. (This figure can be compared tothose of [95, 100].) The PPM results appear slightly bet-ter than WENO5 at a given resolution, likely because ofthe deliberate oscillation suppression in the PPM algo-rithm.
4. 2D shock tube problem
As an additional test of these algorithms’ ability topropagate strong, multidimensional shocks we consider a2D shock tube test. The computational domain ( x, y ) ∈ [0 , × [0 ,
1] is divided into four equal quadrants. The ini- tial fluid states in the lower/upper, left/right quadrantsare ( ρ, P, v x , v y ) LL = (0 . , , , ρ, P, v x , v y ) LR = (0 . , , , . ρ, P, v x , v y ) UL = (0 . , , . , ρ, P, v x , v y ) UR = (0 . , . , , . In this simulation the lower-right and upper-left quad-rants converge on the upper-right quadrant creating apair of curved shocks. We use a Γ = 5 / ×
400 and a CFL factor of 0.5 and are com-parable to [95] and the references therein. Though themain shock features are captured by all of the methodswe considered, oscillations arising from the curved shockfronts are evident in varying degrees. The fourth panel issimilar to the first but contains a refined mesh in the cen-ter that has the same resolution as the other three panels,while the remainder of the domain has half the resolu-tion. Though the majority of the flow is thus effectivelyderefined, the principal features remain the same. Thisis despite the fact that the shocks must travel throughor along refinement boundaries, and the numerical shockspeeds differ slightly on either side of such boundariesdue to the different truncation errors.
C. Hydrodynamic tests in curved spacetime
1. Bondi accretion
As a first test of our code’s ability to simulate rela-tivistic hydrodynamics in the strong field regime, we con-sider Bondi flow. We set up our initial conditions witha stationary solution to spherical accretion onto a blackhole [101]. We use Kerr-Schild coordinates for the blackhole metric. In order to test our code’s ability to convergeto the correct solution we measured how the conserveddensity D differed from the exact solution as a functionof time for three resolutions. The lowest resolution hasa grid spacing of h = 0 . M BH , while the medium andhigh resolutions have twice and 4 times the resolution re-spectively. As shown in Fig. 11, || D − D exact || convergesto zero at second-order. For this test we tried both theMC and WENO5 limiters (with HLL for the flux calcula-tion). Though both had similar levels of error and showedthe expected convergence, WENO5 had larger errors atlow and medium resolutions. This is probably because,at lower resolutions, the larger WENO5 stencil extendsfarther inside the black hole horizon where there is largertruncation error.4 FIG. 9. Density contours (30 equally spaced in the logarithm) for the Emery step problem. The upper (lower) two plotsshow results for resolution 240 ×
80 (480 × ρ min , ρ max ), are (1 . , . × ), (0 . , . × ),(0 . , . × ), and (6 . × − , . × ).
2. Boosted NS
As an additional test of our evolution algorithm, weconsidered a single TOV star with a boost of v = 0 . M ⊙ ). We performed a convergence studyat three resolutions, the lowest of which has approxi-mately 50 points covering the diameter of the star. Themedium (high) resolution has two (three) times the num-ber of points in each dimension, respectively. The AMRhierarchy is identical in all cases, with 7 levels of 2:1refinement, and was determined using truncation errorestimates from the low resolution run. Figure 12 showsthat the constraint violations show the expected second-order convergence to zero. We also compared the perfor-mance of different reconstruction methods (though usingthe HLL flux method throughout). In Fig. 13, we showthe maximum density of the NS as a function of time for various reconstruction methods. Though the driftsand oscillations in density converge away for all meth-ods, we find that WENO5 gives the least density driftcompared to MC and PPM at a given resolution. Thedrift in maximum density with PPM has to do with theway this particular implementation enforces monotonic-ity at extrema, which results in a loss of accuracy (see forexample [102]). Modifying the way the algorithm handlessmooth extrema can reduce this effect. We implementedone such modification (Eqs. 20-23 from [102]), the resultsof which are labeled ‘PPM alt.’ in Fig. 13.
3. Boosted NS flux correction test
As a demonstration of the flux correction algorithm(outlined in Sec. II D) to enforce conservation acrossAMR boundaries, we perform an additional boosted NS5
FIG. 10. Density contours (30 equally spaced in the logarithm) for the 2D Riemann problem with, from left to right, topto bottom: HLL and MC, HLL and WENO5, Marquina and MC, and HLL and MC with mesh refinement. The respectiveminimum and maximum densities, ( ρ min , ρ max ), are (1 . × − , . . × − , . . × − , . . × − , . ×
400 was used. For the final simulation, a refinement region (red box) wasplaced in the middle with equivalent resolution, while the remaining grid has half the resolution ( i.e. this simulation has lowerresolution overall). A CFL factor of 0.5 was used throughout. −3 time (M) || D − D e x a c t || / || D e x a c t || MC Low res.MC Med. res. × (2) MC High res. × (4) WENO High res. × (4) FIG. 11. The L norm of the difference between the numer-ical and exact value of the fluid quantity D (divided by thenorm of the exact solution) for Bondi accretion with MC andWENO5. The low resolution was run with a grid spacing of h = 0 . M BH while the medium and high has twice and 4times the resolution respectively. The results are scaled forsecond-order convergence. test. We use the same conditions as the low resolutionsimulation outlined above in Sec. III C 2 but with a differ-ent AMR hierarchy. In particular, we keep the hierarchyfixed so that the boosted NS will move to areas of succes-sively lower refinement. In Fig. 14, one sees that withoutflux corrections there is a ≈ .
1% loss in fluid rest-massas the NS moves off the highest refinement level, anda ≈ .
8% loss as the NS moves off the next to highestrefinement level. This change in the conserved fluid rest-mass comes from the fact that there is a slight mismatchin fluxes at the mesh refinement boundaries due to trun-cation error. With the flux correction routine activated,this error is eliminated, and the only change in the totalrest-mass is due to the density floor criterion ( i.e., thenumerical atmosphere). As an indication of how the useof flux corrections affects energy and momentum we canalso compare the integrated matter energy and momen-tum as seen by a set of Eulerian observers. The matterenergy density is given by T ab n a n b and the momentumdensity is given by p i = − T ai n a where n a is the timelikeunit normal to the constant t slices. These quantities in-volve combinations of the conserved fluid variables andthe metric and are subject to truncation error, especiallysince in this simulation the NS is allowed to move tolower resolution. In addition, these quantities can varywith time due to gauge effects (though in this case, thevariation due to gauge effects is subdominant to the vari-ation mentioned below), so we use them as an indicationof the effect of flux corrections by comparing them for thesimulations with and without flux corrections to a simi- lar simulation where the AMR tracks the NS and there istherefore essentially no flux across AMR boundaries. Atthe end of the simulation the run without flux correctionshas 0 .
91% less energy compared to a similar simulationwhere the AMR tracks the NS, while the run with fluxcorrections has only 0 .
13% less energy. The run withoutflux corrections has 4 .
0% less momentum (in the boost di-rection) while with flux corrections the comparative lossis 2 . IV. CONCLUSIONS
Numerous scenarios that fall within the purviewof general relativistic hydrodynamics are still mostlyunexplored—especially CO mergers involving neutronstars. There is a rich parameter space, of which largeareas remain uncharted due to uncertainty or potentialvariability in BH and NS masses, BH spin and align-ment, the NS EOS, and other aspects. Beyond the purehydrodynamics problem, the roles of magnetic fields andneutrino physics are just beginning to be explored byvarious groups, and we expect to add support for suchphysics to our code in the future. The potential appli-cations of robust and flexible numerical algorithms forevolving hydrodynamics together with the Einstein fieldequations are manifold. With this in mind, we have im-plemented methods for conservatively evolving arbitraryEOSs, in particular for converting from conserved toprimitive variables without knowledge of derivatives; and7 −5 time (M ) || C a || M Low res.Med. res. × (2) High res. × (3) FIG. 12. The L -norm of the constraint violation ( C a := H a − (cid:3) x a ) in the equatorial plane for a boosted NS simula-tion with v = 0 . M , theADM mass of the NS in its rest frame, and the norm of theconstraints is multiplied by M so as to make it dimensionless. we have implemented numerous reconstruction and fluxcalculation methods that can be used interchangeably tomeet problem specific requirements. Though accuratetreatment of shocks may not be crucial for BH-NS merg-ers (where shocks are not expected to be dynamicallyimportant), the same is not true of NS-NS binaries, espe-cially eccentric ones where the stars may come into con-tact during nonmerger close encounters [103]. We havealso taken care to implement a flux correction algorithmthat preserves the conservative nature of hydrodynamicaladvection across AMR boundaries. Though strict conser-vation is not, strictly speaking, essential (since any non-conservation would be at the level of truncation error), itis an especially appealing property when studying, for ex-ample, CO mergers as potential SGRB progenitors. Af-ter merger, material that did not fall into the black hole— typically on the order of a few percent of the originalNS mass — will fill a large volume making up an ac-cretion disk and potentially unbound material. Thoughaccurately tracking this material is not important for thegravitational dynamics, it is critical for characterizing po-tential EM counterparts to the merger. ACKNOWLEDGMENTS
We thank Adam Burrows, Matt Choptuik, LuisLehner, Scott Noble, Inaki Olabarietta, and Jim Stonefor useful conversations. This research was supported by −4 time (M ) ρ m a x / ρ m a x ( t = ) − WENO−5 Low res.WENO−5 Med. res.MC Low res.PPM Low res.PPM alt. Low res.
FIG. 13. The relative variation of the maximum centraldensity from its initial value ( ρ max /ρ max ( t = 0) −
1) duringa boosted NS simulation with v = 0 . −3 time (M ) ¯ M / ¯ M ( t = ) –1 No flux correc.Flux correc.
FIG. 14. The relative variation of the fluid rest-mass( ¯ M/ ¯ M ( t = 0) −
1) for a boosted NS test with (dotted, redline) and without (solid black line) flux corrections at meshrefinement boundaries. For this test, the mesh hierarchy isfixed so that at t ≈ M the NS has moved from beingcontained entirely on the finest resolution to being containedentirely on the second finest resolution. At t ≈ M the NShas moved from the second to the third finest resolution. Woodhen cluster at Princeton University.
Appendix A: Spectral decomposition of the fluxJacobian
Our conservative formulation of the hydrodynamicalequations (3,4) can be written in vector notation as ∂ t q + ∂ i ( F i ) = S where q is a five dimensional vector of theconserved (in the absence of sources S ) fluid variables q = ( D, S t , S x , S y , S z ) T and the flux F i = ( Dv i , ( S t −√− gP ) v i , S j v i + δ ij √− gP ) T , where the index j in theflux is shorthand for the 3 components ( x, y, z ). Someflux calculation methods such as the Roe solver [77] andthe Marquina flux [78] require the spectral decompositionof the Jacobian ∂ F i ∂ q which we give here. (See [104] for the spectral decomposition for a similar formulation withslightly different conserved variables.) The eigenvaluesare λ ± = αq ( a ± b ) − β i (A1)and λ = αU i − β i (A2)(with multiplicity 3), where a = (1 − c s ) U i , b = c s p (1 − U )[ γ ii (1 − U c s ) − aU i ], q = (1 − U c s ) − , c s is the sound speed, and α, β i , γ ij are metric componentsas in (10). Here and throughout we use i ∈ { x, y, z } torefer to the direction of the flux in the Jacobian withwhich we are concerned, ∂ F i ∂ q . In the following equationswe use the index j as a shorthand for the three spatialcomponents of the eigenvectors (that is, the componentsassociated with S x , S y , and S z ). The indices l and m are fixed by i and the indices n and p are fixed by j asindicated below. The index k is the only index that issummed over. A set of linearly independent right eigen-vectors is given by r ± = (cid:16) , hW [ U k β k − { α ( γ ii − U i U i ) + Aβ i } /B ] , hW ( U j − δ ij A/B ) (cid:17) T , (A3)where A = [ U i c s (1 − U ) ∓ b ] q and B = γ ii − U i ( a ± b ) q , r = (cid:16) κ/ ( HW ( κ/ρ − c s )) , U k β k − α, U j (cid:17) T , (A4) where κ = ∂P∂ǫ , and, r = (cid:16) W U l , hW ( U k β k − α ) U l + hβ l , h ( γ jl +2 W U j U l ) (cid:17) T (A5)where for r , l = y, z, x for i = x, y, z respectively. Theexpression for r can be obtained simply by replacing l with m , where m = z, x, y for i = x, y, z respectively,in the above expression for r . H and W are as definedfollowing (12).We also give the corresponding left eigenvectors.Component-wise, for l ± = ( l D ± , l t ± , l j ± ), l D ± = ∓ f hW V ∓ ξl t ± = ∓ f h ( K − {− γU i + V ∓ ( W ξ − Γ lm ) } + KW V ∓ ξ i /αl j ± = ∓ f h ( γ ln γ mp − γ lp γ mn ) { − KA i − (2 K − V ∓ U i } + (2 K − V ∓ ξW U j i − β j l t ± (A6)where l = y, z, x and m = z, x, y for i = x, y, z respec-tively, and n = y, z, x and p = z, x, y for j = x, y, z re-spectively and Γ lm = γ ll γ mm − γ lm γ lm , ξ = Γ lm − γU i U i , K = (1 − c s ρ/κ ) − , Λ ± = ( a ± b ) q , V ± = ( U i − Λ ± ) / ( γ ii − U i Λ ± ), A i = ( γ ii − U i U i ) / ( γ ii − U i Λ ∓ ), and f − = 2 hW bqξ ( K − γ ii − U i U i ) × [( γ ii − U i Λ + )( γ ii − U i Λ − )] − . Furthermore, l = Wc s ρ ( κ − c s ρ ) (cid:16) h, W/α, W ( U j − β j /α ) (cid:17) , (A7)9and the components of l and l are l D = 0 l t = G lm ( αhξ ) − l j = h δ ji U i G lm + δ jl { γ mm (1 − U i U i ) + γ im U m U i }− δ jm { γ lm (1 − U i U i ) + γ im U l U i } i ( hξ ) − − β j l t (A8)where G lm = ( γ mm U l − γ lm U m ) and for l , l = y, z, x and m = z, x, y for i = x, y, z respectively. The expression for l can be obtained from the above expression for l simplyby interchanging l and m . [1] A. Abramovici et al. Ligo: The laser interferometergravitational wave observatory. Science , 256:325, 1992.[2] H. L¨uck et al.
The geo-600 project.
Class. and Quant.Grav. , 15:1471, 1997.[3] M. Ando et al.
Stable operation of a 300-m laser interfer-ometer with sufficient sensitivity to detect gravitational-wave events within our galaxy.
Phys. Rev. Lett. , 86:3950,2001. astro-ph/0105473.[4] B. Caron et al.
The virgo interferometer.
Class. andQuant. Grav. , 14:1461, 1997.[5] Ramesh Narayan, Bohdan Paczynski, and Tsvi Piran.Gamma-ray bursts as the death throes of massive binarystars.
Astrophys. J. , 395:L83, 1992. astro-ph/9204001.[6] N. Gehrels et al.
A short γ -ray burst apparently asso-ciated with an elliptical galaxy at redshift z = 0.225. Nature (London) , 437:851–854, October 2005.[7] W. H. Lee, E. Ramirez-Ruiz, and J. Granot. A CompactBinary Merger Model for the Short, Hard GRB 050509b.
Astrophys. J. Lett. , 630:L165–L168, September 2005.[8] B. D. Metzger and E. Berger. What is the Most Promis-ing Electromagnetic Counterpart of a Neutron Star Bi-nary Merger?
ArXiv e-prints , August 2011.[9] M. Branchesi, on behalf of the LIGO Scientific Collabo-ration, the Virgo Collaboration, A. Klotz, and M. Laas-Bourez. Searching for electromagnetic counterparts ofgravitational wave transients.
ArXiv e-prints , October2011.[10] The LIGO Scientific Collaboration, Virgo Collabora-tion: J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott,M. Abernathy, T. Accadia, F. Acernese, C. Adams,R. Adhikari, and et al. Implementation and testingof the first prompt search for electromagnetic counter-parts to gravitational wave transients.
ArXiv e-prints ,September 2011.[11] N.K. Glendenning.
Compact Stars: : Nuclear Physics,Particle Physics, and General Relativity . Springer,Berlin, Germany, 2000.[12] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata.Effects of hyperons in binary neutron star mergers.
ArXiv e-prints , October 2011.[13] F. Pannarale, L. Rezzolla, F. Ohme, and J. S. Read.Will black hole-neutron star binary inspirals tell usabout the neutron star equation of state?
ArXiv e-prints , March 2011.[14] J. S. Read, C. Markakis, M. Shibata, K. Ury¯u, J. D. E.Creighton, and J. L. Friedman. Measuring the neutronstar equation of state with gravitational wave observa-tions.
Phys. Rev. D , 79(12):124033, 2009.[15] R. Oechslin and H.-T. Janka. Gravitational Waves fromRelativistic Neutron-Star Mergers with MicrophysicalEquations of State.
Phys. Rev. Lett. , 99(12):121102, 2007.[16] L. K. Tsui, P. T. Leung, and J. Wu. Determinationof the internal structure of neutron stars from gravita-tional wave spectra.
Phys. Rev. D , 74(12):124025, 2006.[17] K. D. Kokkotas and N. Stergioulas. Gravitational Wavesfrom Compact Sources. In A. M. Mour˜ao, M. Pimenta,R. Potting, & P. M. S´a , editor,
New Worlds in As-troparticle Physics: Proceedings of the Fifth Interna-tional Workshop , page 25, 2006.[18] M. Shibata. Constraining Nuclear Equations of StateUsing Gravitational Waves from Hypermassive NeutronStars.
Phys. Rev. Lett. , 94(20):201101, 2005.[19] M. Shibata, K. Taniguchi, and K. Ury¯u. Merger of bi-nary neutron stars with realistic equations of state in fullgeneral relativity.
Phys. Rev. D , 71(8):084021, 2005.[20] J. Faber, P. Grandcl´ement, and F. Rasio. Relativisticcalculations of coalescing binary neutron stars.
Pra-mana , 63:837, 2004.[21] M. Bejger, D. Gondek-Rosi´nska, E. Gourgoulhon,P. Haensel, K. Taniguchi, and J. L. Zdunik. Impactof the nuclear equation of state on the last orbits of bi-nary neutron stars.
Astron. and Astrophys. , 431:297,2005.[22] D. Lai and A. G. Wiseman. Innermost stable circularorbit of inspiraling neutron-star binaries: Tidal effects,post-Newtonian effects, and the neutron-star equationof state.
Phys. Rev. D , 54:3958, 1996.[23] X. Zhuge, J. M. Centrella, and S. L. W. McMillan. Grav-itational radiation from coalescing binary neutron stars.
Phys. Rev. D , 50:6247, 1994.[24] P. J. Montero, O. Zanotti, J. A. Font, and L. Rez-zolla. Dynamics of magnetized relativistic tori oscillat-ing around black holes.
Mon. Not. Roy. Astron. Soc. ,378:1101–1110, July 2007.[25] O. Korobkin, E. B. Abdikamalov, E. Schnetter, N. Ster-gioulas, and B. Zink. Stability of general-relativisticaccretion disks.
Phys. Rev. D , 83(4):043007, February2011.[26] B. D. Farris, Y. T. Liu, and S. L. Shapiro. Binary blackhole mergers in gaseous disks: Simulations in generalrelativity.
Phys. Rev. D , 84(2):024024, July 2011.[27] T. Bode, T. Bogdanovi´c, R. Haas, J. Healy, P. Laguna,and D. Shoemaker. Mergers of Supermassive BlackHoles in Astrophysical Environments.
Astrophys. J. ,744:45, January 2012.[28] V. Paschalidis, Y. T. Liu, Z. Etienne, and S. L.Shapiro. Merger of binary white dwarf–neutron stars:Simulations in full general relativity.
Phys. Rev. D ,84(10):104032, November 2011.[29] Jos A. Font. Numerical hydrodynamics and magneto-hydrodynamics in general relativity.
Living Reviews in Relativity , 11(7), 2008.[30] B. Zink, E. Schnetter, and M. Tiglio. Multipatchmethods in general relativistic astrophysics: Hydrody-namical flows on fixed backgrounds.
Phys. Rev. D ,77(10):103015, May 2008.[31] B. Zink. HORIZON: Accelerated General Relativis-tic Magnetohydrodynamics.
ArXiv e-prints , February2011.[32] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C.Stephens. Relativistic magnetohydrodynamics in dy-namical spacetimes: Numerical methods and tests.
Phys. Rev. D , 72(2):024028, 2005.[33] Tetsuro Yamamoto, Masaru Shibata, and KeisukeTaniguchi. Simulating coalescing compact binaries bya new code (sacra).
Phys. Rev. D , 78:064054, Sep 2008.[34] B. Giacomazzo and L. Rezzolla. WhiskyMHD: a newnumerical code for general relativistic magnetohydro-dynamics.
Class. and Quant. Grav. , 24:235, 2007.[35] M. Thierfelder, S. Bernuzzi, and B. Br¨ugmann. Numer-ical relativity simulations of binary neutron stars.
Phys.Rev. D , 84(4):044012, August 2011.[36] Tanja Bode, Roland Haas, Tamara Bogdanovic, PabloLaguna, and Deirdre Shoemaker. Relativistic Mergersof Supermassive Black Holes and their ElectromagneticSignatures.
Astrophys.J. , 715:1117–1131, 2010.[37] M. Shibata and T. Nakamura. Evolution of three-dimensional gravitational waves: Harmonic slicing case.
Phys. Rev. D , 52:5428, 1995.[38] T. W. Baumgarte and S. L. Shapiro. Numerical in-tegration of Einstein’s field equations.
Phys. Rev. D ,59(2):024007, 1999.[39] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer,M. A. Scheel, and S. A. Teukolsky. Evolving black hole-neutron star binaries in general relativity using pseu-dospectral and finite difference methods.
Phys. Rev. D ,78(10):104015, 2008.[40] M. Anderson, E. W. Hirschmann, S. L. Liebling, andD. Neilsen. Relativistic MHD with adaptive mesh re-finement.
Class. and Quant. Grav. , 23:6503, 2006.[41] D. Garfinkle. Harmonic coordinate method for simulat-ing generic singularities.
Phys. Rev. D , 65(4):044029,2002.[42] F. Pretorius. Numerical relativity using a generalizedharmonic decomposition.
Class. and Quant. Grav. ,22:425, 2005.[43] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, andO. Rinne. A new generalized harmonic evolution sys-tem.
Class. and Quant. Grav. , 23:447, 2006.[44] S. Chawla, M. Anderson, M. Besselman, L. Lehner,S. L. Liebling, P. M. Motl, and D. Neilsen. Mergers ofMagnetized Neutron Stars with Spinning Black Holes:Disruption, Accretion, and Fallback.
Phys. Rev. Lett. ,105(11):111101–+, September 2010.[45] Z. B. Etienne, Y. T. Liu, and S. L. Shapiro. Relativis-tic magnetohydrodynamics in dynamical spacetimes: Anew adaptive mesh refinement implementation.
Phys.Rev. D , 82(8):084031, October 2010.[46] Z. B. Etienne, V. Paschalidis, Y. T. Liu, and S. L.Shapiro. Relativistic MHD in dynamical spacetimes:Improved EM gauge condition for AMR grids.
ArXive-prints , October 2011.[47] Luis Lehner, Carlos Palenzuela, Steven L. Liebling,Christopher Thompson, and Chad Hanna. Intense Elec-tromagnetic Outbursts from Collapsing Hypermassive Neutron Stars. 2011.[48] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, andS. A. Teukolsky. Equation of state effects in black hole-neutron star mergers.
Classical and Quantum Gravity ,27(11):114106–+, June 2010.[49] K. Kiuchi, Y. Sekiguchi, M. Shibata, and K. Taniguchi.Exploring Binary-Neutron-Star-Merger Scenario ofShort-Gamma-Ray Bursts by Gravitational-Wave Ob-servation.
Phys. Rev. Lett. , 104(14):141101, April 2010.[50] K. Kyutoku, H. Okawa, M. Shibata, and K. Taniguchi.Gravitational waves from spinning black hole-neutronstar binaries: dependence on black hole spins andon neutron star equations of state.
Phys. Rev. D ,84(6):064018, September 2011.[51] V. Paschalidis, Z. Etienne, Y. T. Liu, and S. L. Shapiro.Head-on collisions of binary white dwarf-neutron stars:Simulations in full general relativity.
Phys. Rev. D ,83(6):064002, March 2011.[52] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata.Gravitational Waves and Neutrino Emission from theMerger of Binary Neutron Stars.
Phys. Rev. Lett. ,107(5):051102, July 2011.[53] M. Shibata, K. Kiuchi, Y. Sekiguchi, and Y. Suwa.Truncated Moment Formalism for Radiation Hydrody-namics in Numerical Relativity.
Prog. Theor. Phys. ,125:1255–1287, June 2011.[54] B. M¨uller, H.-T. Janka, and H. Dimmelmeier. A NewMulti-dimensional General Relativistic Neutrino Hydro-dynamic Code for Core-collapse Supernovae. I. Methodand Code Tests in Spherical Symmetry.
Astrophys. J.Supp. , 189:104–133, July 2010.[55] P. Cerd´a-Dur´an, J. A. Font, L. Ant´on, and E. M¨uller.A new general relativistic magnetohydrodynamics codefor dynamical spacetimes.
Astron. and Astrophys. ,492:937–953, December 2008.[56] W. H. Lee and W. Klu´zniak. Newtonian Hydrodynamicsof the Coalescence of Black Holes with Neutron Stars. I.Tidally Locked Binaries with a Stiff Equation of State.
Astrophys. J. , 526:178–199, November 1999.[57] S. Rosswog and D. Price. MAGMA: a three-dimensional, Lagrangian magnetohydrodynamics codefor merger applications.
Mon. Not. Roy. Astron. Soc. ,379:915–931, August 2007.[58] J. A. Faber, P. Grandcl´ement, and F. A. Rasio. Mergersof irrotational neutron star binaries in conformally flatgravity.
Phys. Rev. D , 69(12):124036, June 2004.[59] R. Oechslin, H.-T. Janka, and A. Marek. Relativisticneutron star merger simulations with non-zero temper-ature equations of state. I. Variation of binary param-eters and equation of state.
Astron. and Astrophys. ,467:395–409, May 2007.[60] P. Laguna, W. A. Miller, and W. H. Zurek. Smoothedparticle hydrodynamics near a black hole.
Astrophys.J. , 404:678–685, February 1993.[61] S. Rosswog. Conservative, special-relativistic smoothedparticle hydrodynamics.
Journal of ComputationalPhysics , 229:8591–8612, November 2010.[62] M.J. Berger and P. Colella. Local adaptive mesh refine-ment for shock hydrodynamics.
Journal of Computa-tional Physics , 82(1):64, 1989.[63] J. M. Stone, T. A. Gardiner, P. Teuben, J. F. Hawley,and J. B. Simon. Athena: A New Code for AstrophysicalMHD.
Astrophys. J. Supp. , 178:137, 2008. [64] A. S. Almgren, V. E. Beckner, J. B. Bell, M. S. Day,L. H. Howell, C. C. Joggerst, M. J. Lijewski, A. Nonaka,M. Singer, and M. Zingale. CASTRO: A New Com-pressible Astrophysical Solver. I. Hydrodynamics andSelf-gravity. Astrophys. J. , 715:1221–1238, June 2010.[65] D. C. Collins, H. Xu, M. L. Norman, H. Li, and S. Li.Cosmological Adaptive Mesh Refinement Magnetohy-drodynamics with Enzo.
Astrophys. J. Supp. , 186:308–333, February 2010.[66] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zin-gale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Tru-ran, and H. Tufo. FLASH: An Adaptive Mesh Hydro-dynamics Code for Modeling Astrophysical Thermonu-clear Flashes.
Astrophys. J. Supp. , 131:273–334, Novem-ber 2000.[67] K. Beckwith and J. M. Stone. A Second-order GodunovMethod for Multi-dimensional Relativistic Magnetohy-drodynamics.
Astrophys. J. Supp. , 193:6, March 2011.[68] Branson C. Stephens, William E. East, and Frans Pre-torius. Eccentric Black Hole-Neutron Star Mergers.
As-trophys. J. Lett. , 737(1):L5, 2011.[69] W. E. East, F. Pretorius, and B. C. Stephens. Eccentricblack hole-neutron star mergers: effects of black holespin and equation of state.
ArXiv e-prints , November2011.[70] Helmut Friedrich. On the hyperbolicity of ein-stein’s and other gauge field equations.
Communi-cations in Mathematical Physics , 100:525–543, 1985.10.1007/BF01217728.[71] Carsten Gundlach, Jose M. Martin-Garcia, Gioel Cal-abrese, and Ian Hinder. Constraint damping in the Z4formulation and harmonic gauge.
Class. Quant. Grav. ,22:3767–3774, 2005.[72] F. Pretorius. Evolution of Binary Black-Hole Space-times.
Phys. Rev. Lett. , 95(12):121101, 2005.[73] Jos´e A. Font, Mark Miller, Wai-Mo Suen, and MalcolmTobias. Three-dimensional numerical general relativistichydrodynamics: Formulations, methods, and code tests.
Phys. Rev. D , 61:044011, Jan 2000.[74] C. F. Gammie, J. C. McKinney, and G. T´oth. HARM:A Numerical Scheme for General Relativistic Magneto-hydrodynamics.
Astrophys. J. , 589:444–457, May 2003.[75] Miguel Alcubierre, Steven Brandt, Bernd Bruegmann,Daniel Holz, Edward Seidel, et al. Symmetry withoutsymmetry: Numerical simulation of axisymmetric sys-tems using Cartesian grids.
Int.J.Mod.Phys. , D10:273–290, 2001.[76] A. Harten, P.D. Lax, and B.J. van Leer.
SIAM Rev. ,25:35, 1983.[77] F. Eulderink and G. Mellema. General relativistic hy-drodynamics with a ROE solver.
Astron. and Astrophys.Supp. , 110:587, 1995.[78] R. Donat and A. Marquina. Capturing shock reflections:An improved flux formula.
J. Comput. Phys. , 125:42–58, 1996.[79] Chi-Wang Shu and Stanley Osher. Efficient imple-mentation of essentially non-oscillatory shock-capturingschemes, ii.
Journal of Computational Physics , 83(1):32– 78, 1989.[80] E.F. Toro.
Riemann Solvers and Numerical Methods forFluid Dynamics . Springer, Berlin, Germany, 1997.[81] P. Colella and P. R. Woodward. The PiecewiseParabolic Method (PPM) for Gas-Dynamical Simula-tions.
Journal of Computational Physics , 54:174, 1984. [82] Jos´e Ma Mart´ı and Ewald M¨uller. Extension of thepiecewise parabolic method to one-dimensional rel-ativistic hydrodynamics.
Journal of ComputationalPhysics , 123(1):1 – 14, 1996.[83] Guang-Shan Jiang and Chi-Wang Shu. Efficient imple-mentation of weighted eno schemes.
Journal of Compu-tational Physics , 126(1):202 – 228, 1996.[84] A. Tchekhovskoy, J. C. McKinney, and R. Narayan.WHAM: a WENO-based general relativistic numericalscheme - I. Hydrodynamics.
Mon. Not. Roy. Astron.Soc. , 379:469, 2007.[85] Scott C. Noble and Matthew W. Choptuik. Type iicritical phenomena of neutron star collapse.
Phys. Rev.D , 78:064059, Sep 2008.[86] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi.Relativistic equation of state of nuclear matter for su-pernova and neutron star.
Nuclear Physics A , 637:435,1998.[87] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi.Relativistic Equation of State of Nuclear Matter forSupernova Explosion.
Progress of Theoretical Physics ,100:1013, 1998.[88] R.P. Brent.
Algorithms for minimization without deriva-tives . Prentice-Hall, Englewood Cliffs, N.J., 1972.[89] Marsha J Berger and Joseph Oliger. Adaptive meshrefinement for hyperbolic partial differential equations.
Journal of Computational Physics , 53(3):484, 1984.[90] F. Pretorius and M. W. Choptuik. Adaptive mesh re-finement for coupled elliptic-hyperbolic systems.
Jour-nal of Computational Physics , 218:246–274, October2006.[91] Frans Pretorius. Simulation of binary blackhole spacetimes with a harmonic evolution scheme.
Class.Quant.Grav. , 23:S529–S552, 2006.[92] Dieter R. Brill. On the positive definite mass ofthe bondi-weber-wheeler time-symmetric gravitationalwaves.
Annals of Physics , 7(4):466 – 483, 1959.[93] Evgeny Sorkin. On critical collapse of gravitationalwaves.
Classical and Quantum Gravity , 28(2):025011,2011.[94] Gregory B. Cook and Mark A. Scheel. Well-behavedharmonic time slices of a charged, rotating, boostedblack hole.
Phys. Rev. D , 56:4775–4781, Oct 1997.[95] W. Zhang and A. I. MacFadyen. RAM: A RelativisticAdaptive Mesh Refinement Hydrodynamics Code.
As-trophys. J. Supp. , 164:255, 2006.[96] B. Giacomazzo and L. Rezzolla. The exact solutionof the Riemann problem in relativistic magnetohydro-dynamics.
Journal of Fluid Mechanics , 562:223–259,September 2006.[97] B. Giacomazzo, L. Rezzolla, and L. Baiotti. Can mag-netic fields be detected during the inspiral of binaryneutron stars?
Mon. Not. Roy. Astron. Soc. , 399:L164,2009.[98] Luciano Rezzolla and Olindo Zanotti. New relativis-tic effects in the dynamics of nonlinear hydrodynamicalwaves.
Phys. Rev. Lett. , 89:114501, Aug 2002.[99] Ashley F. Emery. An evaluation of several differencingmethods for inviscid fluid flow problems.
Journal ofComputational Physics , 2(3):306 – 331, 1968.[100] Arturo Lucas-Serrano, Jose A. Font, Jose M. Ibanez,and Jose M. Marti. Assessment of a high-resolutioncentral scheme for the solution of the relativistic hydro-dynamics equations.
Astron. Astrophys. , 428:703–715, Black holes,white dwarfs, and neutron stars : the physics of compactobjects . Wiley, New York, 1983.[102] Phillip Colella and Michael D. Sekora. A limiter for ppmthat preserves accuracy at smooth extrema.
Journal ofComputational Physics , 227(15):7069 – 7076, 2008.[103] R. Gold, S. Bernuzzi, M. Thierfelder, B. Bruegmann,and F. Pretorius. Eccentric binary neutron star mergers.
ArXiv e-prints , September 2011.[104] J. M. Ibanez, M. A. Aloy, J. A. Font, J. M. Marti, J. A.Miralles, and J. A. Pons. Riemann Solvers in Gen-eral Relativistic Hydrodynamics.