Numerical relativity simulations in the era of the Einstein Telescope
aa r X i v : . [ g r- q c ] A ug General Relativity and Gravitation manuscript No. (will be inserted by the editor)
Numerical relativity simulations in the era of theEinstein Telescope
Mark Hannam · Ian Hawke
Received: date / Accepted: date
Abstract
Numerical-relativity (NR) simulations of compact binaries are ex-pected to be an invaluable tool in gravitational-wave (GW) astronomy. Thesensitivity of future detectors such as the Einstein Telescope (ET) will placemuch higher demands on NR simulations than first- and second-generationground-based detectors. We discuss the issues facing compact-object simula-tions over the next decade, with an emphasis on estimating where the accuracyand parameter space coverage will be sufficient for ET and where significantwork is needed.
Keywords
Numerical relativity · Black Holes · Neutron Stars · EinsteinTelescope
PACS · Experiments aimed at the first direct detection of gravitational waves are nowunderway with a network of laser-interferometric gravitational-wave (GW)detectors [1,2,3]. At the completion of the current science run the LIGO andVirgo detectors will be upgraded to second generation detectors, with an orderof magnitude improvement in sensitivity, and Advanced LIGO is due to goonline in 2014. It is hoped that the second-generation detectors will yieldenough observations to begin in earnest the field of GW astronomy; see [4] foran overview of the astrophysics and cosmology potential of GW observations.
Mark HannamPhysics Department, University College Cork, Cork, IrelandIan HawkeSchool of Mathematics, University of Southampton, SO17 1BJ, UK
With the goal of achieving a further factor of ten improvement in sensi-tivity, combined with an extension of the detector bandwidth to the range1Hz to 10kHz, a design study for a third generation detector has been sup-ported within the European FP7 framework. If these design goals are met, theEinstein Telescope (ET) would be sensitive to a volume of the universe onemillion times larger than current ground-based detectors, and allow a precisionprobe of gravitational-wave sources throughout the universe to cosmologicaldistances.Extracting the most science from ET observations will require accuratetheoretical knowledge of GW sources — far more accurate than that requiredfor first- and second-generation detectors. One important aspect of sourcemodelling is the numerical solution of Einstein’s equations for the inspiral andmerger of compact bodies (black holes and neutron stars) and the numericalsimulation of the physics of single compact bodies (neutron star oscillations,instabilities and collapse, and core-collapse supernovae amongst others).NR simulations have advanced significantly in the last few years, partic-ularly in describing the last orbits and merger of systems of two black holes(see [5,6,7] for the original breakthrough results, and [8] for an overview ofthe status of current simulations focused on GW detection). The mode-lingof neutron-star binaries has also made significant progress from the first sim-ulations of the Tokyo group ([9]) to the current status as recently reviewedby [10].In this paper we will discuss some of the issues that numerical-relativitysimulations face in producing the results needed for ET data analysis andastrophysics applications.For black-hole binary simulations, the situation is easy to state: since thetotal mass of the system provides an overall scale to the numerical results,the increased bandwidth of ET does not change the physical nature of thesimulations that need to be produced; the black-hole binary parameter spaceremains the same. However, it does change the accuracy requirement of thesimulations, and of the model of the GW parameter-space that can be pro-duced from them.For matter simulations, the situation is quite different. The scales are setby the physics included in the model. The size of the parameter space isdominated by the unknown physics which it is hoped that GW observationswill constrain. In order to make simulations practical, only the “essential”physics will be included in the model. The determination of “essential” isbalanced between current knowledge and expectations of the key physics, andwhat can practically be simulated and detected. Thus the model, resultingparameter space and accuracy requirements will change with the increasedbandwidth of ET.In Section 2 the status of vacuum simulations of binary black holes isdiscussed, including a discussion of the waveform accuracy of current and near-future simulations, the sampling of the parameter space, the construction oftemplate banks and possible future improvements in numerical techniques.Section 3 discusses the situation with matter simulations, assuming that the physical model is “added on” to the spacetime simulation and GWs extractedusing the same techniques. Section 4 summarises the likely status of NR overthe next decade.
Fully general-relativistic simulations of the dynamics of black-hole binary sys-tems during their last orbits and merger have been possible since 2005 [5,6,7]. In these simulations initial data are prescribed for a slice of spacetimethat contains two black holes, usually in orbit, and in many cases with initialparameters chosen, or tuned, such that the black holes are following non-eccentric quasi-circular orbits [11,12,13]. The initial data must also satisfyfour constraint equations in order to be part of a valid solution of Einstein’sequations. The data are advanced to future time slices via a set of evolutionequations; the specific equations depend on the choice of decomposition of thefour-dimensional Einstein equations to 3+1 (space+time) dimensions, and thecoordinates and time-slicing at successive evolution steps depend on the choiceof gauge conditions. We are in principle free to choose the gauge conditions aswe wish (as a result of the coordinate invariance of Einstein’s equations), butin practice are limited to choices that lead to sufficiently stable and accuratenumerical simulations. The issues involved in finding a numerically well-posedand stable formulation of Einstein’s equations, a convenient geometrical rep-resentation of a black-hole spacetime, construction of black-hole-binary initialdata, and numerically and geometrically well-behaved gauge conditions, arediscussed in the textbook [14], and in the review articles [15,16,17].The parameter space of inspiraling black-hole binaries consists of ten pa-rameters: the total mass and mass ratio of the system (or alternatively theindividual masses of the two black holes), the spin vector of each black hole,the eccentricity of the system, and the initial separation of the binary (or al-ternatively the initial phase of the GW signal). Since the total mass providesan overall scale for the solution, this can be removed from the parameter spaceof necessary numerical waveforms. If the numerical late-inspiral-plus-mergerwaveform is combined with an analytic inspiral waveform (for example, cal-culated from the post-Newtonian approximation), the initial separation of thebinary can be made arbitrarily large, and yet one more parameter can be re-moved from the parameter space. Assuming that the process of constructinga “hybrid” waveform is sufficiently robust (and this remains a topic of currentresearch), we are left with a total of eight parameters.To date late-inspiral-merger-ringdown simulations of black-hole binarieshave been performed for systems with mass ratios up to 1:10 (although mostdo not extend beyond 1:4), and systems with a variety of spin magnitudesand orientations, and for several choices of initial eccentricity. A periodicallyupdated summary of simulations that include at least ten GW cycles beforemerger is given in the online version of the review [8]. Note that although alarge number of simulations have been performed, these only account for an extremely sparse sampling of the full parameter space, and do not yet includeany simulations of mass ratios beyond q = m /m = 10, or of spins above a/m ∼ .
92. Efforts to use numerical-relativity simulations to make study theastrophysics of electromagnetic counterparts have also begun [18,19].The issues for black-hole-binary simulations, as they relate to GW astron-omy, are: (1) determining and achieving sufficient accuracy to extract the max-imum physical information from GW observations, which includes producingsimulations that include enough cycles that subsequent analytic-numerical hy-brid waveforms meet these accuracy requirements, (2) producing simulationsof a sufficiently dense sampling of the parameter space, (3) developing methodsand codes that are efficient enough to achieve these goals.In distinguishing the requirements on numerical waveforms between ETand other detectors, for example Advanced LIGO, the key difference is thatsince ET is expected to be roughly an order of magnitude more sensitive thanAdvanced LIGO, the waveform accuracy will also need to be higher to allowthe most science to be extracted from observations. The question for numericalrelativity then is what those accuracy requirements are, and whether they canbe achieved before ET is completed.2.1 Waveform accuracyIn the last few years questions of waveform accuracy requirements have begunto be addressed [20], and it seems that current methods and simulations aresufficiently accurate for the current generation of ground-based detectors, al-though a systematic study has been performed only for the case of equal-massnonspinning binaries [21]. In that work it was shown that, within certain rea-sonable caveats, it will not be possible to use GW observations to distinguishbetween the most accurate numerical waveforms used in the study unless thesignal-to-noise ratio (SNR) is above about 25.SNRs of that magnitude are expected to be rare from the Initial and En-hanced LIGO and Virgo detectors. However, the Einstein Telescope is expectedto be an order of magnitude more sensitive, and SNRs in the range 100-1000may be typical. We first consider whether numerical simulations will be ableto match the accuracy requirements of ET by the time of its construction.Waveform accuracy requirements are discussed in some detail in [20]. Onesimple criterion for waveforms to be sufficiently accurate for all parameterestimation purposes is that the average error in the amplitude and phase(appropriately weighted by the power spectrum of the detector noise [20])satisfy ¯ δφ + ¯ δχ < ρ , (1)where ρ is the SNR, and ¯ δχ is the fractional error in the amplitude and ¯ δφ is the error in the phase. If we assume that most black-hole-binary signals ofinterest for ET will have an SNR of less than 1000, then we require that theaverage phase and fractional amplitude errors be of the order of 10 − . For the waveforms described in [22] (as an example of a state-of-the-art simulation),the maximum phase and amplitude errors are roughly an order of magnitudelarger than the average errors; we expect that the detector-noise-weightedaverage errors will behave similarly.As such, the most accurate numerical simulations to date [23] have totalphase and amplitude errors of 0.02 radian and 0.3% respectively (see the tablein [21] for a convenient comparison of errors between codes), and so we expectthat the average errors are already within the accuracy requirements of ET.The simulation that deals most convincingly with the wave extraction [22]gives average phase and amplitude errors also within the approximate 10 − requirement provided here.So it appears that state-of-the-art numerical simulations are already at,or close to, the accuracy requirements for ET science. Assuming that Moore’sLaw holds, then the increase in computing power by the time ET is operationalwill be 2 / . ≈ any values of the black-holeparameters. As such, the natural approach is to construct an analytic modelof the waveforms, and to use the numerically generated waveforms as inputto calibrate the free parameters of the model. Such models also aim to extendthe numerical waveforms to an arbitrarily large number of GW cycles beforemerger (or, put another way, extend them to arbitrarily low GW frequency).To date there have been two approaches to this problem. One is to startwith an effective-one-body (EOB) model of inspiral-merger-ringdown wave-forms (which in some forms pre-date the existence of full numerical wave-forms [30,31,32,33,34]), and to both calibrate unknown EOB coefficients to numerical results, and to introduce extra “flexibility parameters” that allowyet greater fidelity between the EOB model and the full GR results. EOB-basedmodels for nonspinning waveforms have been pursued by several groups [35,36,37,38,39,40,41,42], and have reached an impressive level of faithfulness tothe full GR waveforms [38,41,42].The second approach has been to construct a purely phenomenologicalansatz motivated by post-Newtonian (PN) predictions for the inspiral, ap-proximate models of the ringdown, and the observed behaviour of the signalin the intermediate plunge and merger phase, as seen for example in [43]. Thismodel, in which each ingredient is motivated by known physics, is extremelyflexible, and has been extended from original work on nonspinning binaries [44,45,46] to binaries with nonprecessing spins [47] and has been shown to alsobe applicable to neutron-star binaries [48]. This procedure is limited only byavailable numerical simulations, although it remains to be seen how easily a“simple” ansatz can be developed for general precessing-spin binaries. Thenonspinning waveform model is constructed from ten mass-ratio-dependentparameters, which can in turn be modelled by second-order polynomials, there-fore requiring at least three simulations to define the model. The EOB models,by contrast, can be calibrated by only one simulation. It remains to be seenhow well this economy of calibration extends to the spinning-binary cases.Both the EOB and pure phenomenological approaches require that thewaveform model for the GW cycles before the beginning of the NR wave-form (i.e., the PN or EOB model) is accurate. This requires comparison ofNR waveforms with their PN counterparts in the regions where they overlap.NR-PN comparisons have been performed for equal-mass nonspinning bina-ries [49,50,51,52], equal-mass binaries with non-precessing spins [53], eccentricbinaries [54], and one configuration of an unequal-mass binary with precessingspins [13]. These results suggest that PN methods perform well up to the pointwhere NR simulations begin. Detailed studies of the accuracy of NR-PN hy-brid waveforms based on these simulations remain to be made, but first studiesof the accuracy of equal-mass hybrids at least suggest that they will fulfil the detection requirements of current ground-based detectors for masses as low as10 M ⊙ [55]. Complementary results for studies of pure PN inspiral waveformssuggest that they in turn will be sufficient for masses up to M ⊙ [56]. How-ever, the wider sensitivity band of ET may require that the hybrid waveformsinclude more accurate PN ingredients, or longer numerical simulations. Thisis a topic that deserves further study in the coming years.2.3 Sampling of the parameter spaceWhichever phenomenological approach turns out to be most practical, it isclear that simulations covering a fine sampling of the black-hole-binary pa-rameter space will be necessary. An exhaustive sampling of the parameterspace is usually presented as a monumental computational challenge. Let us consider the scope of this challenge, and the degree to which it may be metby the time ET might be commissioned.There are eight parameters in the black-hole-binary parameter space. A re-cent phenomenological family of inspiral-merger-ringdown waveforms [47] wasable to faithfully model binaries with non-precessing spins using an analyticansatz in which the coefficients were expressed as third-order polynomials inthe physical parameters of the model (in this case the mass ratio and a param-eter that can be approximately considered as the total spin of the binary). Forsimplicity, let us take the relatively optimistic view that third-order polyno-mials in each of the parameters will be sufficient to model the full parameterspace, and that we therefore require at least four data points in each directionto robustly constrain our model. This means that we will require 4 ≈ , total spin of the binary — with theconsequence that if the individual spins aren’t necessary to model the wave-form, then they will also be difficult to distinguish. This reduces the amountof information that can be obtained from a detection based on this waveformmodel, but also reduces the computational cost of producing the model.As a second example, the works of [57,58] exploit the symmetries of black-hole-binary systems to argue that certain quantities in black-hole mergers (themass, spin and recoil velocity of the final black hole) can be described by amodel that requires only 26 simulations. Their model deals with only quasi-circular inspiral (so the parameter space is reduced to seven, not eight, param-eters), and uses a second-order-polynomial ansatz, but nonetheless requires farfewer simulations than a naive estimation of 3 ≈ ℓ = 2 , m = 2 mode) ofthe waveform, and to date it is only this mode that has been included inanalytic/phenomenological models of black-hole-binary waveforms. Also, twoimportant areas of the parameter space have not yet been accessed by currentsimulations.One is black holes with high spin. All published long simulations to datehave used black-hole initial data that are conformally flat, an assumption thatsimplifies the construction of the data, but limits the spin that can be mod-elled. The behaviour of black-hole-binary systems with near-extremal spins(which may be the most common astrophysically [59,60,61]), have there-fore not yet been studied in detail. Proposals have been made to constructnon-conformally-flat initial data for both excision-based [62,63] and moving-puncture [64] codes, although only the excision data is sufficiently developedfor use in full binary simulations. With workable proposals in place for high-spin data, it is likely that high-spin binaries will be simulated in detail in the next few years, and this region of the parameter space will be well modelledby the time of ET.The situation is different for large mass ratios. In current simulations, thenumerical resolution requirements of the simulation are defined by the smallestblack hole, while the overall scale is set by the total mass. This means thatthe higher the mass ratio ( q = m /m ≥ q = 6 [39], and a short simulation of ≈ . q = 10.Since the simulation cost (both in simulation time and memory usage) scales at best linearly with the mass ratio, and more realistically quadratically, longaccurate simulations of mass ratios q ∼
10 are at the very limit of currentmethods and resources, and q ∼
100 are out of the question. An improvementby a factor of 100 in computer power by the time of ET suggests that massratios q ∼
30 may be possible by that time, but q ∼
100 will require newtechniques; some possibilities are described in Section 2.4.It may ultimately not be necessary to simulate systems with q >
10 in fullgeneral relativity to accurately model the full parameter space. For example,the work in [47] includes simulations up to only q = 3, but models all mass ra-tios by explicitly incorporating an analytically calculated extreme-mass-ratiolimit. The effective-one-body (EOB) procedure also includes by constructionthe extreme-mass-ratio limit [30]. However, in order to verify that the highmass ratios are indeed being faithfully modelled, it will be necessary to per-form at least a few simulations at very high mass ratios. Even if only a fewsuch simulations are required to verify the fidelity of a given model (or set ofmodels), new techniques will be necessary to produce those simulations.2.4 Numerical techniques and room for improvementFor codes that employ the moving-puncture method, numerical methods havenot changed significantly since the first simulations in 2005, but the accuracyand length of the waveforms has improved significantly. The first waveformscovered only a few GW cycles before merger, while the most recent simulationscover up to 30 cycles before merger. The first long simulations that allowedcomparison with post-Newtonian (PN) results had an amplitude uncertaintyof around 10% and an accumulated phase error of around 1 radian beforemerger [49]. The most recent simulations can achieve an amplitude error of0.3% and an accumulated phase error through inspiral, merger and ringdownof 0.02 rad [28].Typically mesh-refinement is used to resolve both the small-scale featuresnear the black holes, the medium-scale physics of the GW signal, and the largescales necessary to push the outer boundary of the computational domain farfrom the source. However, improvements in general code efficiency and theintroduction of higher-order spatial finite-differencing (first sixth-order in [66,50], and later eighth-order in [67,13]) have allowed impressive improvements in both code speed, memory efficiency, and numerical accuracy. More recentresults by two groups [68,29] have demonstrated the efficacy of multipatchmethods, which allow for a far more efficient gridding of the wave extractionregion, which may lead to yet greater accuracy (of the wave extraction) andcomputational efficiency in future simulations. The results of Reisswig et. al. have further achieved Cauchy-Characteristic Extraction (CCE), in which thewaves extracted at some finite spatial distance from the source are then evolvedto null infinity, which is where GW signals can be unambiguously defined [22,69]. These results indeed suggest that this method produces results that aregauge invariant.Simulations with variants of the generalized-harmonic formulation havebeen performed by Pretorius [70,5,71] and by the Caltech-Cornell group [72,52]. The latter’s SpEC code makes use of a multi-domain spectral method,which allows a high degree of accuracy and, like the multi-patch methodsdescribed above, allow high accuracy with minimal computational resourcesin the wave extraction region.These codes represent the current state of the art. However, far greateraccuracy and efficiency will be needed to meet the science goals of ET. Isthere room for improvement?The answer is certainly Yes. A number of options are yet to be fully ex-plored. One that has already received some attention is the use of mixedimplicit-explicit (IMEX) time integration to allow far larger evolution timesteps. A first exploration of this approach has already been made in the con-text of scalar fields on a curved background [73]. Another approach with arelated potential benefit is the use of symplectic integrators [74,75,76,77]. Inthis approach, conserved quantities in the physical system are maintained byconstruction, and this in general allows for far larger time steps. The canonicalexample is the evolution of binaries in Newtonian gravity, where dramatic im-provements in accuracy are possible. General relativistic systems are of coursefar more complex, and there is the added difficulty that energy can be dissi-pated as gravitational radiation, but these issues do not preclude the use ofsymplectic methods, and some research on this topic is already underway.A further improvement in both accuracy and memory efficiency could beachieved by the use of hyperboloidal slices of the spacetime. Here the slices arespacelike near the source, but asymptotically null, meaning that they reach nullinfinity on a finite numerical domain. Since the GW phase is constant on a nullslice, only a finite number of GW cycles need be resolved between the sourceand the outer boundary. This approach would therefore save computational re-sources in the wave zone, while at the same time allowing clean GW extractionat null infinity. Although hyperboloidal initial data for such systems has al-ready been proposed [78,79], stable simulations are more challenging. Progresshas however been made in spherical symmetry and axisymmetry [80], and hy-perboloidal simulations of black-hole binaries may be possible in the next fewyears.In the case of moving-puncture simulations, improvements may also bepossible by better exploiting the natural geometry of black-hole slices in this approach. The method was developed for data that describe the black holes astopological wormholes, but it was later found that the slices quickly asymptoteto cylinders that end at the black-hole throat (also called “trumpets”) [81,82,83]. It may be possible to modify the method to exploit more explicitly thistrumpet form of the data. In particular, if the original concept of the punctureapproach can be resurrected — where the singular nature of the data near theblack holes is described analytically and only small deviations from this form(due to the presence of the second black hole) are described numerically —then it may be possible to achieve dramatic improvements in code accuracyand efficiency.These are a number of avenues for improvement that are already beingexplored, and have the potential by themselves to make possible large numbersof long, accurate simulations possible — but of course we expect that by 2020there will be many other ideas, and it is quite conceivable that by then thestandard methods of the field will be radically different from those today. In contrast to simulations of binary black holes, many calculations of gravita-tional waves have not used full general relativity (for example most calcula-tions of core-collapse as reviewed by e.g. [84] and some binary simulations suchas [85]). This is unsurprising; for scenarios such as core-collapse supernovaethe effect of including the key physics in the model is more important than aperfect model of gravity. However, with the recent successes in evolving thespacetime as detailed in section 2 the use of full general relativity in calcu-lations of gravitational waves is likely to become standard on the timescalesimportant for ET.In addition we note that as the matter terms do not affect the principle partof the spacetime evolution equations there is only one potential complicationfor the spacetime evolution by adding matter. In vacuum evolutions, awayfrom singularities, all quantities remain smooth. As the matter quantities maybecome discontinuous, a simulation that includes such shocks will have aneffect on the spacetime. It is currently unclear what the impact is in GWsimulations, as the only work to look at this ([86]) is in 1 + 1 dimensions.As no effect has been seen in any current simulation we will assume that nofundamental problem will affect matter simulations relevant for ET.From these basic points we will assume that all current GW simulationsincluding matter will be able to benefit from the theoretical and numericaladvances for spacetime and gauge evolution and wave extraction reviewed orproposed in section 2. Therefore the bottleneck in computing accurate, realisticGW signals will be in the physical model of the matter that can practicallybe evolved, and the accuracy of the numerical methods used to evolve it.The parameter space for matter simulations is considerably larger thanfor binary black holes. Firstly the underlying physics of the problem sets themass scale which cannot be removed. Secondly there are a wide variety of GW sources; as well as binary systems there are discs and single star collapse.But most importantly there is the uncertainty in the physics that the modelmust include to be realistic (e.g., neutrino transport plays a crucial role in corecollapse but is much less important in a binary inspiral), and in the details ofthat physics (such as the precise form of the equation of state (EOS) neededto close the perfect fluid equations).The approach so far has been to focus on limited areas in the parameterspace to determine which elements of the physical model affect the qualitativeform of the GWs, and if the GW signal is relatively robust to determine whatqualitative effects minor changes in the parameters make. The key additionalparameters in inspirals are currently believed to be the bulk EOS and themagnetic field; in core collapse neutrinos and composition are more importantthan magnetic fields. Additional effects such as superfluidity, an elastic crust,diffusion, dissipation, or other non-ideal effects may play a role. At presenteven the simplest implementation of these effects is ongoing research.The construction of initial data for binary systems is similar to that forbinary black holes and, as reviewed by [15,10], similar issues of the accuracy ofthe conformal flatness approximation and how to produce quasi-equilibriuminitial data arise. For single star collapse either axi-stationary solutions areperturbed or a pre-calculated model is used, with similar concerns about thephysical accuracy of the initial conditions. In most simulations small parameterstudies show that the qualitative behaviour of the gravitational waves is robust;one fundamental problem will be discussed later.The key question for GWs and in particular for ET is whether numericalrelativity can compute sufficiently accurate signals. As this crucially dependson the physics included in the model to be simulated we will discuss thedifferent input physics separately.3.1 HydrodynamicsDue to the scale and temperature of most interesting GW sources the bulkof the matter can be reasonably modelled as a fluid, and in highly dynamicalsituations such as a binary merger the bulk ideal fluid effects, modelled by therelativistic Euler equations, will dominate. These equations are typically writ-ten in conservation law form [87] as they admit discontinuous (weak) solutionsfor which the conservation law form gives a unique entropy satisfying solution.The microphysics under consideration is determined by the EOS which closesthe system of equations and is typically given as a relation between the pres-sure and the independent thermodynamic and composition variables.The importance of the EOS is obvious, as has been illustrated by a numberof simulations (as reviewed by [10], but see in particular [85,88,89]). Largeand small scale effects both in terms of the stiffness and the effects of heat insimplified form have been covered. However current uncertainty in the EOS islarge, and at present the approaches followed are to restrict to the “best” EOScurrently available or to use a simple parametrized EOS ([90]) to best fit the likely candidates. The effectiveness of the parametrized EOS combined withnumerical simulations has been illustrated by [91], although it should be notedthat the accuracy of the simulations discussed there will not be sufficient forET.Numerical solutions of the relativistic Euler equations use techniques fromthe standard Newtonian Euler equations. In particular it is currently typicalto use High Resolution Shock Capturing (HRSC) methods (see e.g. [92,93]and references therein for reviews). These ensure conservation is preserved(necessary to ensure the correct weak solution is obtained if a shock is present)and ensure that the solution converges near discontinuities. These techniquesare nonlinear in the sense that the information used at any given point dependson the data in the neighbourhood of that point; this is required to bypass theconditions of Godunov’s theorem ([94]) and construct numerical methods thatgive correct solutions that converge faster than first order.The techniques currently used in simulations are extremely similar, evenwhen the spacetime is evolved using spectral methods ([95,96]), multiple patches([97]) or a mesh refined cartesian grid (e.g. [98,99,100]). In nearly every case afinite volume approach combined with a reconstruction-evolution method andapproximate Riemann solvers is used, computing a conservative update overthe neutron star and the surrounding artificial atmosphere. These methodsconverge at best at second order and are computationally expensive com-pared to the methods required for the spacetime. A notable exception is thesmoothed-particle hydrodynamics simulations of e.g. [85,101,102], but as yetthese have not used full general relativity (but note the recent results of [103]).These methods have been used to study purely hydrodynamical scenarios.The numerical accuracy of the resulting GW signals is not always easy toassess, with some work at lower resolutions not claiming accuracy of betterthan 1 orbit ([91]), whilst others at the highest currently available resolutionsclaiming accuracy in key quantities of the order of a percent ([104]). Even insimulations with careful control of the accuracy in the signal the amplitudeaccuracy is nearly an order of magnitude worse than early black hole binarysimulations (compare figure 5 of [104] with figure 3 of [105]). It is also clearthat numerical accuracy is significantly degraded when shocks appear such asin the complex post-merger stage of a binary inspiral. Thus whilst it seemsthat current techniques and computing power are sufficient to make qualitativestatements about gravitational waves and the general effects of changes in theparameter space, detailed quantitative predictions at the level required byET are unlikely to result simply by throwing computing power at currenttechniques, except in certain simple cases.Improvements in numerical techniques to significantly improve the accu-racy of the results, particularly in the complex post-merger regime, are there-fore needed to get the most out of numerical relativity for ET. Extensions ofstandard finite volume techniques to higher order (e.g. [106]) are one obviousavenue to explore. The implementation of more complex but flexible methodssuch as Discontinuous Galerkin finite elements (as implemented in relativityby [107,108]) is another possibility. Finite difference methods have received some attention (e.g. [109,110]) and may be the simplest approach. At presentit seems unlikely that the highly accurate spectral methods will provide asolution due to the problems these methods encounter with discontinuities(although [111] presents one approach in a particular approximation), but hy-brid methods (e.g. [112,113]) may allow for a combination of the best featuresof spectral or simple high order schemes with HRSC methods.Finally it should be noted that the impact of the artificial atmosphereused outside the compact GW source on the accuracy of the results remainsunclear, and methods to avoid the use of an atmosphere (e.g. [114,115]) areat present not sufficiently advantageous to use in relevant simulations. Thepresence of a free surface (in the perfect fluid model) may be an impedimentto the use of more accurate numerical methods, and remains an area wheremore understanding at both numerical and theoretical level is needed. Purelyon numerical grounds the use of finite elements is attractive here, as the gridcan be adapted to the surface of the objects; all that is required is a consistentboundary condition to impose.3.2 Magnetic fieldsSo far the ideal MHD limit has been the focus of a range of studies (for arecent review see [10]), many simply “adding on” the magnetic field to studythe qualitative differences (see e.g. [116,117]), but some studying additionalinstabilities that may arise (e.g. [102,118]). With the increase in the parameterspace due to the strength and topology of the magnetic field the choice of initialconditions becomes ever more important, and this is where one fundamentaldifficulty remains. It is expected that the initial object will be approximatelyaxi-stationary and the precise topology and mixture of toroidal and poloidalcomponents that would make up this initial field has only been studied nu-merically under certain assumptions, see [119,120,121]. Until a self-consistentquasi-stationary solution can be constructed the standard method is to add apurely poloidal magnetic field to the interior; it has been argued ([116]) thatthis is sufficient for the inspiral phase of a magnetized binary.In addition to the conservation constraints, numerical methods should iden-tically preserve the ∇ · B = 0 constraint. There are many possible approaches,with most codes implementing a variant of constrained transport (e.g. [122])and some using hyperbolic divergence cleaning (e.g. [110]). Constrained trans-port methods should ensure the constraint is satisfied to machine precisionwhilst divergence cleaning propagates any errors away; however divergencecleaning is simpler to extend to higher order methods and the presence ofthe neutron star surface and artificial atmosphere may affect the accuracy ofconstrained transport. In all cases the additional constraints to be satisfiedand fields to be evolved lead to lower numerical accuracy; for example [116]shows the numerical error increasing by factors ∼ The transfer of heat and energy through radiation transport changes the localstructure of the fluid, and is known to be crucial in the dynamics of core-collapse supernovae and may play a role in the post-merger dynamics of abinary merger. However the simulation of the full Boltzmann equation in 3 +3+1 dimensional general relativity would seem to be impractical on timescalesrelevant for ET. An estimate ([123]) suggests that even approximate transportmethods require peta scale computing power. A suggestion for modelling thefull problem is given by [124], but on the timescale for ET it seems likely thatapproximation methods such as [125] will be the only practical solution.
The consideration of non-ideal effects may become important post-mergerwhen the matter surrounding the remnant can be both very hot and yet nottoo dense, leading to a plasma where diffusion and resistivity are important.Viscous effects have barely been touched on – for an isolated example see [126].Resistivity effects have been considered and modelled in the magnetosphere,but modelling for neutron stars is only just starting (see [107,108]).The Newtonian limit gives mixed hyperbolic-parabolic systems of equationswith their associated causality problems. The relativistic equivalents have stiffsource terms leading to severe timestep constraints when using standard nu-merical methods. Various numerical schemes have been proposed to bypassthis issue, including the use of IMEX time integrators ([107]). It seems plausi-ble that similar numerical techniques will be adopted as those for large massratio binary black holes.
It is believed that most neutron stars have a region containing superfluidneutron pairs and that magnetized neutron stars may be superconducting insome region. The crust of a mature neutron star is not a fluid but a lattice,best modelled as a relativistic elastic solid inter-penetrated by an additionalfluid. In addition, the possible existence of exotic phases of matter in the innercore leads to another region where multiple particles may be flowing freely and independently. All of these effects involve the local interaction of multiple fluidsor solids (for a review see [127]), with the appearance of additional dynamicsand instabilities.At present no numerical evolutions of any of these effects have been per-formed. Frameworks for evolving multiple fluids have been worked out ([127])but have not been tailored for numerical work that would include shocks. Itseems unlikely that multifluid effects will be visible in numerical simulationsin the inspiral phase, but the additional propagation speeds and instabilitiescould be visible near or post merger. Over the last four decades, numerical relativists have been primarily concernedwith the technical challenges of solving Einstein’s equations numerically: for-mulating Einstein’s equations to be numerically stable, developing appropriategauge conditions, exploring the numerical techniques necessary to accuratelyevolve relativistic spacetimes, and to handle relativistic fluids and magneticfields. Many of these issues have been resolved to a sufficient degree that wecan now extract useful physics from the simulations. In addition, sensitive GWdetectors are now in operation and the analysis of their data requires numer-ical results. The field of numerical relativity is thus shifting from a focus ontheoretical and numerical-analysis issues to questions of astrophysics and GWdata analysis.Over the next decade, this shift will continue, and numerical relativists willinteract more closely with astronomers, astrophysicists, and GW data analysts.Closer collaboration between the NR and DA communities has already begunthrough the Numerical INJection Analysis (NINJA) project, where numericalblack-hole-binary waveforms were injected into simulated Initial LIGO andVirgo detector noise, and then the data were analyzed by a battery of searchand parameter estimation methods, to test their efficacy in dealing with “real”GW signals. The successful first NINJA project [128] is now being followedwith projects to more systematically test search pipelines. Among these arethe suggestion to use actual detector data, which contain the full array of realdetector artifacts and can allow the most conclusive tests of search pipelines,and to extend the original simulated-detector-noise NINJA strategy to mattersources (NS-NS, NS-BH and core collapse).With the advent of Advanced LIGO/Virgo and the space-based detectorLISA [129], the interaction between the NR and DA communities will surelybecome closer, and the possibility of multi-messenger GW astronomy will in-clude astronomers and astrophysicists in these efforts, too. Once GW detec-tions become routine, and numerical simulations computationally cheaper andfaster, it may be typical for a candidate observation to motivate follow-up nu-merical simulations, which in turn lead to more sophisticated data analysisparameter estimation exercises, and comparison with electromagnetic obser-vations. The time scale of this back-and-forth interaction may be as short as weeks or even days. For some regions of the (at least black-hole-binary) pa-rameter space, it is in fact likely that GW astronomers will have ready accessto numerical codes that allow them to “dial up” any waveform they requirefor a particular data-analysis exercise, and the focus of numerical-relativity re-search will be only on those cases (for example, high mass ratios, binaries withcomplex spin interactions, and matter systems) that still present problems.For matter simulations considerable work has been done to explore theparameter space by including as much of the relevant physics as is currentlypractical. However, in order to take full advantage of the additional accuracygiven by ET these simulations will have to improve substantially in certainareas. Firstly the numerical methods currently employed will probably notdeliver the improvements in accuracy required based on the (conservative)increase in computing power projected here. Secondly the community will haveto collaborate closely with the data analysis community in order to determinewhich aspects of the physical model will actually lead to interesting observableeffects, and hence where in the parameter space the numerical simulationsshould be focused. As the discussion in section 3 makes clear it would bepossible to spend all the time and increased computing power in applying andsimulating better physical models and exploring the parameter space. It seemslikely that to take the greatest advantage of ET a narrower focus is needed toproduce the most accurate simulations including the most relevant physics forthe best candidate sources. Acknowledgements
MH was supported by SFI grant 07/RFP/PHYF148, and thanks theUniversity of the Balearic Islands for hospitality while some of this work was completed,and Sascha Husa for numerous discussions. IH was supported by the STFC rolling grantPP/E001025/1 and Nuffield Foundation grant NAL/32622.
References
1. B. Abbott, et al., (2007). arXiv:0711.30412. F. Acernese, et al., Class. Quantum Grav. , S635 (2006)3. S. Hild (for the LIGO Scientific Collaboration), Class. Quantum Grav. , S643 (2006)4. B.S. Sathyaprakash, B.F. Schutz, Living Rev. Rel. , 2 (2009)5. F. Pretorius, Phys. Rev. Lett. , 121101 (2005). DOI 10.1103/PhysRevLett.95.1211016. M. Campanelli, C.O. Lousto, P. Marronetti, Y. Zlochower, Phys. Rev. Lett. , 111101(2006)7. J.G. Baker, J. Centrella, D.I. Choi, M. Koppitz, J. van Meter, Phys. Rev. Lett. ,111102 (2006)8. M. Hannam, Class. Quant. Grav. , 114001 (2009). DOI 10.1088/0264-9381/26/11/1140019. M. Shibata, K. Uryu, Phys. Rev. D61 , 064001 (2000). DOI 10.1103/PhysRevD.61.06400110. J. Faber, Classical and Quantum Gravity (11), 114004 (18pp) (2009). URL http://stacks.iop.org/0264-9381/26/114004
11. H.P. Pfeiffer, et al., Class. Quant. Grav. , S59 (2007). DOI 10.1088/0264-9381/24/12/S0612. S. Husa, M. Hannam, J.A. Gonzalez, U. Sperhake, B. Br¨ugmann, Phys. Rev. D77 ,044037 (2008). DOI 10.1103/PhysRevD.77.04403713. M. Campanelli, C.O. Lousto, H. Nakano, Y. Zlochower, Phys. Rev.
D79 , 084010 (2009)714. M. Alcubierre,
Introduction to 3+1 Numerical Relativity (Oxford University Press,USA, 2008)15. G.B. Cook, Living Rev. Rel. , 5 (2000)16. F. Pretorius, in Physics of relativistic objects in compact binaries: from birth to co-alescence , vol. 359, ed. by M. Colpi, P. Casella, V. Gorini, U. Moschella, A. Possenti(Springer-Verlag, 2009), vol. 359, pp. 269–305. arXiv:0710.133817. S. Husa, Eur. Phys. J. ST , 183 (2007)18. C. Palenzuela, M. Anderson, L. Lehner, S.L. Liebling, D. Neilsen, (2009).arXiv:0905.112119. J.R. van Meter, et al., (2009). arXiv:0908.002320. L. Lindblom, B.J. Owen, D.A. Brown, Phys. Rev.
D78 , 124020 (2008). DOI 10.1103/PhysRevD.78.12402021. M. Hannam, et al., Phys. Rev.
D79 , 084025 (2009)22. C. Reisswig, N.T. Bishop, D. Pollney, B. Szilagyi, (2009). arXiv:0907.263723. M.A. Scheel, et al., Phys. Rev.
D79 , 024003 (2009). DOI 10.1103/PhysRevD.79.02400324. A.M. Sintes, A. Vecchio, (1999). gr-qc/000505825. C. Van Den Broeck, A.S. Sengupta, Class. Quant. Grav. , 1089 (2007). DOI 10.1088/0264-9381/24/5/00526. S. Babak, M. Hannam, S. Husa, B. Schutz, (2008). arXiv:0806.159127. J.I. Thorpe, et al., (2008). arXiv:0811.083328. M. Boyle, et al., Phys. Rev. D78 , 104020 (2008). DOI 10.1103/PhysRevD.78.10402029. D. Pollney, Talk given at NRDA2009, AEI-Potsdam, July 6-9, 200930. A. Buonanno, T. Damour, Phys. Rev. D , 084006 (1999)31. A. Buonanno, T. Damour, Phys. Rev. D , 064015 (2000)32. T. Damour, P. Jaranowski, G. Sch¨afer, Phys. Rev. D , 084011 (2000)33. T. Damour, Phys. Rev. D , 124013 (2001)34. A. Buonanno, Y. Chen, T. Damour, Phys. Rev. D74 , 104005 (2006)35. A. Buonanno, et al., Phys. Rev.
D76 , 104049 (2007)36. T. Damour, A. Nagar, Phys. Rev.
D77 , 024043 (2008). DOI 10.1103/PhysRevD.77.02404337. T. Damour, A. Nagar, E.N. Dorband, D. Pollney, L. Rezzolla, Phys. Rev.
D77 , 084017(2008). DOI 10.1103/PhysRevD.77.08401738. T. Damour, A. Nagar, M. Hannam, S. Husa, B. Br¨ugmann, Phys. Rev.
D78 , 044039(2008). DOI 10.1103/PhysRevD.78.04403939. J.G. Baker, et al., Phys. Rev.
D78 , 044046 (2008). DOI 10.1103/PhysRevD.78.04404640. A.H. Mroue, L.E. Kidder, S.A. Teukolsky, Phys. Rev.
D78 , 044004 (2008). DOI10.1103/PhysRevD.78.04400441. T. Damour, A. Nagar, (2009). arXiv:0902.013642. A. Buonanno, et al., (2009). arXiv:0902.079043. A. Buonanno, G.B. Cook, F. Pretorius, Phys. Rev.
D75 , 124018 (2007). DOI 10.1103/PhysRevD.75.12401844. P. Ajith, et al., Class. Quant. Grav. , S689 (2007). DOI 10.1088/0264-9381/24/19/S3145. P. Ajith, et al., Phys. Rev. D77 , 104017 (2008). DOI 10.1103/PhysRevD.77.10401746. P. Ajith, Class. Quant. Grav. , 114033 (2008). DOI 10.1088/0264-9381/25/11/11403347. P. Ajith, et al., (2009). In preparation48. M. Shibata, K. Kyutoku, T. Yamamoto, K. Taniguchi, Phys. Rev. D79 , 044030 (2009).DOI 10.1103/PhysRevD.79.04403049. J.G. Baker, J.R. van Meter, S.T. McWilliams, J. Centrella, B.J. Kelly, Phys. Rev. Lett. , 181101 (2007). DOI 10.1103/PhysRevLett.99.18110150. M. Hannam, S. Husa, U. Sperhake, B. Br¨ugmann, J.A. Gonzalez, Phys. Rev. D77 ,044020 (2008). DOI 10.1103/PhysRevD.77.04402051. A. Gopakumar, M. Hannam, S. Husa, B. Bruegmann, Phys. Rev.
D78 , 064026 (2008).DOI 10.1103/PhysRevD.78.06402652. M. Boyle, et al., Phys. Rev.
D76 , 124038 (2007). DOI 10.1103/PhysRevD.76.12403853. M. Hannam, S. Husa, B. Bruegmann, A. Gopakumar, Phys. Rev.
D78 , 104007 (2008).DOI 10.1103/PhysRevD.78.104007854. I. Hinder, F. Herrmann, P. Laguna, D. Shoemaker, (2008)55. M. Hannam, et al., (2009). In preparation56. A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, B.S. Sathyaprakash, (2009)57. L. Boyle, M. Kesden, S. Nissanke, Phys. Rev. Lett. , 151101 (2008). DOI 10.1103/PhysRevLett.100.15110158. L. Boyle, M. Kesden, Phys. Rev.
D78 , 024017 (2008). DOI 10.1103/PhysRevD.78.02401759. M. Volonteri, P. Madau, E. Quataert, M.J. Rees, Astrophys. J. , 69 (2005). DOI10.1086/42685860. C.F. Gammie, S.L. Shapiro, J.C. McKinney, Astrophys. J. , 312 (2004)61. S.L. Shapiro, Astrophys. J. , 59 (2005)62. G. Lovelace, R. Owen, H.P. Pfeiffer, T. Chu, Phys. Rev.
D78 , 084017 (2008). DOI10.1103/PhysRevD.78.08401763. G. Lovelace, Class. Quant. Grav. , 114002 (2009). DOI 10.1088/0264-9381/26/11/11400264. M. Hannam, S. Husa, B. Br¨ugmann, J.A. Gonzalez, U. Sperhake, Class. QuantumGrav. , S15 (2007)65. J.A. Gonzalez, U. Sperhake, B. Brugmann, Phys. Rev. D79 , 124006 (2009). DOI10.1103/PhysRevD.79.12400666. S. Husa, J.A. Gonzalez, M. Hannam, B. Br¨ugmann, U. Sperhake, Class. Quant. Grav. , 105006 (2008). DOI 10.1088/0264-9381/25/10/10500667. C.O. Lousto, Y. Zlochower, Phys. Rev. D77 , 024034 (2008). DOI 10.1103/PhysRevD.77.02403468. E. Pazos, M. Tiglio, M.D. Duez, L.E. Kidder, S.A. Teukolsky, (2009). arXiv:0904.049369. J. Winicour, Living Reviews in Relativity (3) (2009). URL
70. F. Pretorius, Class. Quantum Grav. , 425 (2005)71. F. Pretorius, Class. Quantum Grav. , S529 (2006)72. L. Lindblom, M.A. Scheel, L.E. Kidder, R. Owen, O. Rinne, Class. Quantum Grav. , S447 (2006)73. S.R. Lau, H.P. Pfeiffer, J.S. Hesthaven, (2008). arXiv:0808.259774. J.D. Brown, Phys. Rev. D73 , 024001 (2006). DOI 10.1103/PhysRevD.73.02400175. R. Richter, C. Lubich, Class. Quant. Grav. , 225018 (2008). DOI 10.1088/0264-9381/25/22/22501876. R. Richter, Class. Quant. Grav. , 145017 (2009). DOI 10.1088/0264-9381/26/14/14501777. J. Frauendiener, (2008). arXiv:0805.446578. F. Ohme, M. Hannam, S. Husa, N.O. Murchadha, (2009). arXiv:0905.045079. L.T. Buchman, H.P. Pfeiffer, J.M. Bardeen, (2009). arXiv:0907.316380. V. Moncrief, O. Rinne, Class. Quant. Grav. , 125010 (2009). DOI 10.1088/0264-9381/26/12/12501081. M. Hannam, S. Husa, D. Pollney, B. Bruegmann, N. O’Murchadha, Phys. Rev. Lett. , 241102 (2007). DOI 10.1103/PhysRevLett.99.24110282. M. Hannam, S. Husa, F. Ohme, B. Bruegmann, N. O’Murchadha, Phys. Rev. D78 ,064020 (2008). DOI 10.1103/PhysRevD.78.06402083. M. Hannam, S. Husa, N.O. Murchadha, (2009). arXiv:0908.106384. C.D. Ott, Class. Quant. Grav. , 063001 (2009). DOI 10.1088/0264-9381/26/6/06300185. R. Oechslin, H.T. Janka, Phys. Rev. Lett. , 121102 (2007). DOI 10.1103/PhysRevLett.99.12110286. A.P. Barnes, P.G. Lefloch, B.G. Schmidt, J.M. Stewart, Class. Quant. Grav. , 5043(2004). DOI 10.1088/0264-9381/21/22/00387. F. Banyuls, J.A. Font, J.M. Ib´a˜nez, J.M. Mart´ı, J.A. Miralles, Astrophys. J. , 221(1997)88. K. Kiuchi, Y. Sekiguchi, M. Shibata, K. Taniguchi, (2009). arXiv:0904.455189. L. Baiotti, B. Giacomazzo, L. Rezzolla, Phys. Rev. D78 , 084033 (2008). DOI 10.1103/PhysRevD.78.08403390. J.S. Read, B.D. Lackey, B.J. Owen, J.L. Friedman, Phys. Rev.
D79 , 124032 (2009).DOI 10.1103/PhysRevD.79.124032991. J.S. Read, et al., Phys. Rev.
D79 , 124033 (2009). DOI 10.1103/PhysRevD.79.12403392. J.A. Font, Living Reviews in Relativity (7) (2008). URL
93. J.M. Mart´ı, E. M¨uller, Living Reviews in Relativity (7) (2003). URL
94. S.K. Godunov, Mat. Sb. , 271 (1959)95. H. Dimmelmeier, J. Novak, J.A. Font, J.M. Ibanez, E. Muller, Phys. Rev. D71 , 064023(2005). DOI 10.1103/PhysRevD.71.06402396. M.D. Duez, et al., Phys. Rev.
D78 , 104015 (2008). DOI 10.1103/PhysRevD.78.10401597. B. Zink, E. Schnetter, M. Tiglio, Phys. Rev.
D77 , 103015 (2008). DOI 10.1103/PhysRevD.77.10301598. L. Baiotti, I. Hawke, L. Rezzolla, E. Schnetter, Phys. Rev. Lett. , 131101 (2005).DOI 10.1103/PhysRevLett.94.13110199. T. Yamamoto, M. Shibata, K. Taniguchi, Phys. Rev. D78 , 064054 (2008). DOI10.1103/PhysRevD.78.064054100. M. Anderson, et al., Phys. Rev. Lett. , 191101 (2008). DOI 10.1103/PhysRevLett.100.191101101. R. Oechslin, S. Rosswog, F. Thielemann, Phys.Rev. D , 103005 (2002)102. D.J. Price, S. Rosswog, Science (5774), 719 (2006). DOI 10.1126/science.1125201103. S. Rosswog, (2009). arXiv:0907.4890104. L. Baiotti, B. Giacomazzo, L. Rezzolla, Class. Quant. Grav. , 114005 (2009). DOI10.1088/0264-9381/26/11/114005105. J.G. Baker, J. Centrella, D.I. Choi, M. Koppitz, J. van Meter, Phys. Rev. D73 , 104002(2006). DOI 10.1103/PhysRevD.73.104002106. A. Tchekhovskoy, J.C. McKinney, R. Narayan, Mon. Not. Roy. Astron. Soc. , 469(2007). DOI 10.1111/j.1365-2966.2007.11876.x107. C. Palenzuela, L. Lehner, O. Reula, L. Rezzolla, (2008). arXiv:0810.1838108. M. Dumbser, O. Zanotti, (2009). arXiv:0903.4832109. D. Neilsen, E.W. Hirschmann, R.S. Millward, Class. Quant. Grav. , S505 (2006).DOI 10.1088/0264-9381/23/16/S12110. M. Anderson, E. Hirschmann, S.L. Liebling, D. Neilsen, Class. Quant. Grav. , 6503(2006). DOI 10.1088/0264-9381/23/22/025111. S. Bonazzola, L. Villain, M. Bejger, Class. Quant. Grav. , S221 (2007). DOI 10.1088/0264-9381/24/12/S15112. B. Costa, W.S. Don, Journal of Computational and Applied Mathematics (2), 209(2007). DOI DOI:10.1016/j.cam.2006.01.039113. B. Costa, W.S. Don, Journal of Computational Physics (2), 970 (2007). DOIDOI:10.1016/j.jcp.2006.11.002114. W. Kastaun, Phys. Rev. D74 , 124024 (2006). DOI 10.1103/PhysRevD.74.124024115. M.D. Duez, P. Marronetti, S.L. Shapiro, T.W. Baumgarte, Phys. Rev.
D67 , 024004(2003). DOI 10.1103/PhysRevD.67.024004116. Y.T. Liu, S.L. Shapiro, Z.B. Etienne, K. Taniguchi, Phys. Rev.
D78 , 024012 (2008).DOI 10.1103/PhysRevD.78.024012117. B. Giacomazzo, L. Rezzolla, L. Baiotti, (2009). arXiv:0901.2722118. K. Kiuchi, M. Shibata, S. Yoshida, Phys. Rev.
D78 , 024029 (2008). DOI 10.1103/PhysRevD.78.024029119. J. Braithwaite, A. Nordlund, A&A (3), 1077 (2006). DOI 10.1051/0004-6361:20041980. URL http://dx.doi.org/10.1051/0004-6361:20041980 , 1947 (2008)121. J. Braithwaite, MNRAS , 763 (2009)122. G. Toth, J. Comput. Phys. , 605 (2000)123. E. Schnetter, C.D. Ott, G. Allen, P. Diener, T. Goodale, T. Radke, E. Seidel, J. Shalf,in
Petascale Computing: Algorithms and Applications , ed. by D.A. Bader (Chapman& Hall/CRC Computational Science Series, 2007), chap. 24124. B. Zink, (2008). arXiv:0810.5349125. B.D. Farris, T.K. Li, Y.T. Liu, S.L. Shapiro, Phys. Rev.
D78 , 024023 (2008). DOI10.1103/PhysRevD.78.024023126. M.D. Duez, Y.T. Liu, S.L. Shapiro, B.C. Stephens, Phys. Rev.
D69 , 104030 (2004).DOI 10.1103/PhysRevD.69.1040300127. N. Andersson, G.L. Comer, Living Rev. Rel. , 1 (2005)128. B. Aylott, et al., Class. Quant. Grav. , 165008 (2009). DOI 10.1088/0264-9381/26/16/165008129. K. Danzmann, A. R¨udiger, Classical Quantum Gravity , S1 (2003). URL, S1 (2003). URL