SENR/NRPy+: Numerical Relativity in Singular Curvilinear Coordinate Systems
aa r X i v : . [ g r- q c ] M a r SENR/NRPy+: Numerical relativity in singular curvilinear coordinate systems
Ian Ruchlin, Zachariah B. Etienne,
1, 2 and Thomas W. Baumgarte Department of Mathematics, West Virginia University, Morgantown, West Virginia 26506, USA Center for Gravitational Waves and Cosmology, West Virginia University,Chestnut Ridge Research Building, Morgantown, West Virginia 26505, USA Department of Physics and Astronomy, Bowdoin College, Brunswick, Maine 04011, USA (Dated: March 30, 2018)We report on a new open-source, user-friendly numerical relativity code package calledSENR/NRPy+. Our code extends previous implementations of the BSSN reference-metric for-mulation to a much broader class of curvilinear coordinate systems, making it ideally suited tomodeling physical configurations with approximate or exact symmetries. In the context of modelingblack hole dynamics, it is orders of magnitude more efficient than other widely used open-source nu-merical relativity codes. NRPy+ provides a Python-based interface in which equations are writtenin natural tensorial form and output at arbitrary finite difference order as highly efficient C code,putting complex tensorial equations at the scientist’s fingertips without the need for an expensivesoftware license. SENR provides the algorithmic framework that combines the C codes generatedby NRPy+ into a functioning numerical relativity code. We validate against two other established,state-of-the-art codes, and achieve excellent agreement. For the first time—in the context of mov-ing puncture black hole evolutions—we demonstrate nearly exponential convergence of constraintviolation and gravitational waveform errors to zero as the order of spatial finite difference deriva-tives is increased, while fixing the numerical grids at moderate resolution in a singular coordinatesystem. Such behavior outside the horizons is remarkable, as numerical errors do not converge tozero near punctures, and all points along the polar axis are coordinate singularities. The formula-tion addresses such coordinate singularities via cell-centered grids and a simple change of basis thatanalytically regularizes tensor components with respect to the coordinates. Future plans includeextending this formulation to allow dynamical coordinate grids and bispherical-like distribution ofpoints to efficiently capture orbiting compact binary dynamics.
I. INTRODUCTION
The Laser Interferometer Gravitational-Wave Obser-vatory (LIGO) and Virgo Scientific Collaboration’s directdetections of six gravitational wave events—five blackhole binary mergers [1–5] and one neutron star binarymerger [6, 7]—have ushered in the age of gravitationalwave astrophysics and multimessenger astronomy. Asthe signal-to-noise ratios of gravitational wave detectionsgrow with increased interferometer sensitivity, the needto continually improve our theoretical models of thesephenomena is critical, as new physics may otherwise bemissed.Just over a decade ago, breakthroughs in numericalrelativity [8–10] opened the door to simulating the in-spiral, merger, and ringdown phases of black hole bina-ries in vacuum. Black hole simulations—and, indeed, allcompact binary simulations—span many orders of mag-nitude both in length scale and timescale. Making themcomputationally tractable for reliable gravitational wavepredictions requires that the underlying numerical gridstructure be tuned to optimally sample the space. Or-dered meshes that map to Cartesian grids are most prac-tical, as this greatly simplifies algorithms for high-orderapproximations of spatial derivatives in Einstein’s equa-tions.There are currently two basic approaches in numeri-cal relativity to setting up numerical grids for compactbinary simulations. The most popular is to apply adap- tive mesh refinement (AMR), where the grids consist ofnested Cartesian coordinate boxes at different, discretenumerical resolutions. This enables the application ofhighest numerical resolution where it is most needed:the strongly curved spacetime fields inside and aroundcompact objects. The most widely used AMR infras-tructure is provided by the open-source Cactus/Carpetcode [11–14] within the Einstein Toolkit [15, 16] (ETK).One downside to this approach is that the compact ob-jects of interest are typically round, not rectangular,which results in a highly inefficient distribution of gridpoints inside and near compact objects (see, e.g., [17]for a more detailed analysis). In addition, the suddenchange in grid resolution at the refinement boundariescan produce spurious reflections in sharp gauge modesor high-frequency gravitational wave features when us-ing the moving puncture formalism, resulting in poorconvergence in gravitational waveforms extracted fromthese numerical models [18, 19].The current alternative to adaptive mesh refinementin numerical relativity codes, pioneered by the Sim-ulating eXtreme Spacetimes (SXS) Collaboration [20]in their Spectral Einstein Code (SpEC [21], andits in-development successor SpECTRE [22]), is tosmoothly juxtapose a large number of curvilinear three-dimensional grid patches, each with a smooth one-to-onemapping to a Cartesian grid. Evaluating derivatives withthese grids requires computation of Jacobians, which,alongside managing the dynamics of the grid structureitself, remains a significant contributor to their codes’computational costs for compact binary inspirals that areaimed at generating gravitational wave predictions.Both approaches involve control systems that adjustthe grids to track compact objects as they orbit, al-though the SXS control system is far more complex be-cause their formalism [23] additionally requires for stabil-ity that black hole interiors be very carefully excised fromthe computational domain. Negative side effects of thealgorithmic complexity for both methods include a steeplearning curve for new users and difficulty in interpretingnumerical errors.We address these drawbacks by developing a new codethat builds on an innovative rescaling approach for solv-ing Einstein’s equations in spherical coordinates. Theapproach is designed to take full advantage of symme-tries in the underlying configuration, requiring as fewnumerical grid points as possible and unlocking the desk-top as a powerful tool for numerical relativity. Our workbuilds on strategies that, in the context of spherical coor-dinates, rescale tensors component by component so thatthe detrimental effects of coordinate singularities on nu-merics are completely removed [24–29]. Treating all co-ordinate singularities analytically, the equations can thenbe integrated numerically without encountering instabili-ties. We generalize this approach to a much broader classof static orthogonal coordinate systems by absorbing thecoordinate singularities out of the tensor components andinto a noncoordinate basis. The rescaling strategy is im-plemented in the context of the BSSN reference-metricformulation to enable highly efficient puncture black holeevolutions using the moving puncture approach [9, 10, 30]in a broad class of spherical-, cylindrical-, and Cartesian-like coordinate systems without special integration meth-ods or introducing gaps in the numerical grid.We implement this approach within a new, open-sourcecode package called SENR/NRPy+ [17]. At its core,SENR/NRPy+ aims to be as algorithmically simple anduser friendly as possible, all while being highly efficient.In short, SENR/NRPy+ aims to minimize both hu-man and computational expense while maximizing sci-ence outcomes.SENR/NRPy+ is built upon the philosophy that thedistribution of points on the numerical grid should takemaximum advantage of approximate symmetries in thephysical system. Compact object systems of interest ingravitational wave astronomy typically possess a high de-gree of angular symmetry, making spherical- and, moregenerally, cylindrical-like coordinate systems ideal candi-dates for efficient sampling.NRPy+ (“Python-based code generation for numeri-cal relativity and beyond”) is designed to convert theBSSN reference-metric formulation of the Einstein equa-tions, in a broad class of orthogonal coordinate systems,from Einstein-like notation directly into C code. It op-erates without the need for expensive, proprietary com-puter algebra systems like Mathematica or Maple. As itsname suggests, NRPy+ is based entirely in Python anddepends only on the standard Python computer algebra package SymPy [31] for symbolic algebra, which is widelyavailable on supercomputing clusters.SENR (“the Simple, Efficient Numerical Relativitycode”) incorporates C codes generated by NRPy+ toform a complete, OpenMP-parallelized [32] numerical rel-ativity code. Its skeletal structure makes the algorithmicunderpinnings of numerical relativity codes transparentto the user.We verify SENR/NRPy+ by direct comparisons withtwo other numerical relativity codes that are both well es-tablished in the literature. In the context of strongly per-turbed Minkowski spacetime (a version of the robust sta-bility test [33–36]), we achieve roundoff-level agreementwith the Baumgarte et al. [29] code, which evolves theBSSN equations in spherical coordinates at fixed fourth-order finite difference accuracy. We also demonstrate ex-cellent agreement between the results of SENR/NRPy+and the ETK in the context of simulating a single, dy-namical black hole. Then, we perform simulations ofsingle- and double-black-hole spacetimes, demonstratingthat the finite difference truncation error converges tozero with increasing grid resolution at the expected rate.Perhaps most importantly, we show for the first timethat—in the context of moving puncture evolutions—thetruncation error in our finite differencing scheme con-verges to zero nearly exponentially outside punctureblack hole horizons with linear increase in the finite differ-ence order, keeping the numerical grids fixed at moderateresolution.The fact that we observe nearly exponential conver-gence is remarkable for two reasons. First, the numericalgrids chosen for these simulations possess coordinate sin-gularities at all points where r sin( θ ) = 0. Therefore,near-exponential convergence in regions within causalcontact of these singularities demonstrates that our ten-sor rescaling algorithm and cell-centered grids completelyeliminate convergence problems related to these coordi-nate singularities. Second, a puncture black hole exhibitsnonsmooth fields at the site of the puncture, meaning weshould have no a priori expectation of near-exponentialnumerical convergence outside the black hole, either. Weattribute the observed convergence to the fact that thecharacteristics of the physical (nongauge) fields, in thevicinity of the puncture, point towards the puncture. When finite difference truncation error dominates, we expectnumerical error to scale approximately as (cid:12)(cid:12) C n u ( n +1) ( ξ ) (cid:12)(cid:12) (∆ x ) n ,where n = N FD is the finite difference order, ξ is in the neigh-borhood of the point at which we evaluate the derivative of thefunction u , and C n ∼ / n for a centered stencil on a numeri-cal grid with uniform spacing ∆ x . In the case that (cid:12)(cid:12) u ( n +1) ( ξ ) (cid:12)(cid:12) isbounded (e.g., u ( r ) = sin r ), pure exponential convergence of thefinite difference derivative error is observed as n is increased and∆ x is held fixed (again, assuming that truncation error dom-inates). When simulating gravitational fields, (cid:12)(cid:12) u ( n +1) ( ξ ) (cid:12)(cid:12) cangrow as n ! (e.g., u ( r ) = 1 / (1 − r ), r = 1), reducing the rate ofexponential convergence at the finite difference orders we typi-cally choose ( N FD ∈ { , , , , } ). We refer to this behavior as nearly exponential convergence . With sufficient resolution inside of the horizon, the er-rors resulting from finite differencing across the puncturesingularity become trapped near the puncture, and arenot able to escape and contaminate the simulation atlarge [37, 38].Future SENR projects will involve further extendingthe formalism to handle dynamical, bispherical-like co-ordinate systems, so that compact binary dynamics maybe modeled with minimum computational expense. Inthis work we demonstrate near-exponential convergenceof gravitational waves with increased finite differencingorder in the context of head-on collisions of punctureblack holes, and will build on this success to tackle theorbital black hole binary problem as our next step.This paper is organized as follows. In Sec. II, wepresent the reference-metric formulation of the BSSNand gauge evolution equations. In Sec. III, we outlinethe tensor component rescaling procedure that makes itpossible to evolve gravitational fields in singular coor-dinate systems. In Sec. IV, we describe the structureof the SENR/NRPy+ code, including the implementa-tion of grid structures, coordinate system options, diag-nostics, and boundary conditions. In Sec. V, we firstdemonstrate that SENR/NRPy+ agrees with the resultsof other established numerical relativity codes. Then, weshow that numerical errors converge to zero at the ex-pected rates for a nonspinning black hole with varyinggrid resolution or finite difference order, and in differentcoordinate systems and gauges. Finally, in the contextof head-on collisions of two nonspinning puncture blackholes, we demonstrate near-exponential convergence ofthe gravitational waveforms with increasing finite differ-ence order. We conclude and present plans for futurework in Sec. VI.Throughout this paper geometrized units are adopted,in which c = 1 and G = 1. Latin indices ( i, j, k, . . . )denote spatial degrees of freedom and obey the Einsteinsummation convention. II. BSSN AND GAUGE EVOLUTIONEQUATIONS
In this section, we describe our strategy for solv-ing Einstein’s equations, based on the tensor-weight-zero BSSN [39–41] formulation of Brown [26]. Ournumerical implementation extends the rescaling ap-proach, developed by Baumgarte et al. [42] for spher-ical coordinates, to a broader class of singular curvi-linear coordinate systems. In this paper, we focuson spherical-, cylindrical-, and Cartesian-like coordinatesystems, though the method can be easily extended tomany others. For any of our coordinate grids, the refer-ence metric ˆ γ ij represents the flat space metric compo-nents expressed in a coordinate basis. In Sec. III, we useˆ γ ij in the rescaling procedure to define a noncoordinatebasis, in terms of which which all tensor components areexplicitly free of coordinate singularities. We assume that the background is independent of the coordinate time t ,so that ∂ t ˆ γ ij = 0. All hatted quantities are associatedwith the reference metric.The reference metric is used to decompose the confor-mal metric ¯ γ ij into a correction about the flat background¯ γ ij = ˆ γ ij + ε ij , (1)where the components in ε ij are not necessarily smalland contain the metric fields that are evolved on ournumerical grids. The conformal metric is related to thephysical spatial metric γ ij though a conformal rescaling γ ij = e φ ¯ γ ij . (2)Taking the determinant of (2) we observe that the con-formal factor e φ can be expressed as e φ = (cid:18) γ ¯ γ (cid:19) / , (3)where γ ≡ det( γ ij ) and ¯ γ ≡ det(¯ γ ij ) are the metric de-terminants. Quantities associated with the conformalmetric ¯ γ ij are barred. For example, the conformal co-variant derivative operator ¯ D i is defined with respect tothe conformal metric. The inverse conformal metric ¯ γ ij is defined to satisfy ¯ γ ik ¯ γ kj = δ ij , (4)where δ ij is the Kronecker delta tensor.The conformal rescaling (2) is not yet unique. In con-trast to the original BSSN formulation, in which ¯ γ wasset to unity, we adopt Brown’s “Lagrangian” choice [26] ∂ t ¯ γ = 0 , (5)so that ¯ γ remains equal to its initial value. In particular,this implies that both γ and ¯ γ are allowed to transform asdeterminants, i.e., as scalar densities of weight two. Ac-cording to (3) the conformal factor e φ then transformsas a scalar, rather than a scalar density, and all othertensorial objects in our formalism similarly transform astensors of weight zero (see also [24, 26]). We choose¯ γ = ˆ γ ≡ det(ˆ γ ij ) in the initial data for all applicationsin this paper. Since both determinants remain indepen-dent of time, they remain equal to each other throughoutthe evolution.It is well known that Christoffel symbols do not trans-form covariantly between coordinate systems. However,the difference of two sets of Christoffel symbols is tenso-rial. We define the tensor∆ ijk ≡ ¯Γ ijk − ˆΓ ijk , (6)whose indices are raised and lowered with the conformalmetric. It is useful to construct a vector by taking thetrace ∆ i ≡ ¯ γ jk ∆ ijk . (7)In addition, the tensor-weight-zero conformal connectioncoefficient three-vector ¯Λ i is evolved independently, andsatisfies the initial constraint C i ≡ ¯Λ i − ∆ i = 0 . (8)The conformal, trace-free part of the extrinsic curva-ture is denoted¯ A ij = e − φ (cid:18) K ij − γ ij K (cid:19) , (9)where K ij is the physical extrinsic curvature and K = γ ij K ij is the mean curvature.Defining the hypersurface-normal derivative operator ∂ ⊥ ≡ ∂ t − L β , (10)where L β is the Lie derivative along the shift vector β i ,the BSSN evolution system in vacuum is written as [26] ∂ ⊥ ε ij = 23 ¯ γ ij (cid:0) α ¯ A kk − ¯ D k β k (cid:1) + 2 ˆ D ( i β j ) − α ¯ A ij , (11a) ∂ ⊥ ¯ A ij = −
23 ¯ A ij ¯ D k β k − α ¯ A ik ¯ A kj + α ¯ A ij K + e − φ (cid:8) − α ¯ D i ¯ D j φ + 4 α ¯ D i φ ¯ D j φ +4 ¯ D ( i α ¯ D j ) φ − ¯ D i ¯ D j α + α ¯ R ij (cid:9) TF , (11b) ∂ ⊥ W = − W (cid:0) ¯ D k β k − αK (cid:1) , (11c) ∂ ⊥ K = 13 αK + α ¯ A ij ¯ A ij − e − φ (cid:0) ¯ D i ¯ D i α + 2 ¯ D i α ¯ D i φ (cid:1) , (11d) ∂ ⊥ ¯Λ i = ¯ γ jk ˆ D j ˆ D k β i + 23 ∆ i ¯ D j β j + 13 ¯ D i ¯ D j β j − A ij ( ∂ j α − ∂ j φ ) + 2 ¯ A jk ∆ ijk − α ¯ γ ij ∂ j K . (11e)In the above, “TF” denotes the trace-free part of theexpression in brackets, the conformal factor is evolvedas W = e − φ (following, e.g., [28, 43] to ensure smootherspacetime fields near puncture black holes), and the com-ponents of the conformal Ricci tensor are calculated by¯ R ij = −
12 ¯ γ kl ˆ D k ˆ D l ¯ γ ij + ¯ γ k ( i ˆ D j ) ¯Λ k + ∆ k ∆ ( ij ) k + ¯ γ kl (cid:16) mk ( i ∆ j ) ml + ∆ mik ∆ mjl (cid:17) . (12)The trace-free condition ¯ γ ij ¯ A ij = 0, which can be vio-lated by numerical error, is enforced dynamically by theterm proportional to ¯ A ii in Eq. (11a) [26]. In this paper,we restrict ourselves to vacuum spacetimes so that all ofthe matter source terms vanish.The evolution system is completed with specificationof the gauge: the lapse function α and the shift vector β i . Unless otherwise stated, we employ the advective 1+loglapse condition [44] ∂ α = − αK (13)and the advective Gamma-driver shift condition [30] ∂ β i = B i , (14a) ∂ B i = 34 ∂ ¯Λ i − ηB i . (14b)Here, B i is an auxiliary vector, η is a (dimensionful)damping parameter [45], and the (noncovariant) advec-tive derivative operator is defined as (see Case No. 8 inTable I of [46]) ∂ ≡ ∂ t − β j ∂ j . (15)The BSSN equations coupled to these gauge conditionsare known together as the moving puncture approach [9,10, 30]. A total of 24 fields are evolved. III. TENSOR RESCALING
In the previous section, we described the reference-metric formulation of the BSSN evolution equations inthe moving puncture approach (11), (13), and (14).These equations form the foundation for solving Ein-stein’s equations in any coordinate system we like, andwe leverage this coordinate freedom to take maximum ad-vantage of near symmetries in physical systems. While,from an analytic perspective, there is no problem solvingEinstein’s equations in arbitrary coordinate systems, nu-merical solutions diverge if the chosen coordinate systempossesses a coordinate singularity.Well-known examples of coordinate singularities in-clude the points at ρ = 0 in cylindrical coordinates orat r sin( θ ) = 0 in spherical coordinates. Tensor compo-nents that are regular everywhere in a Cartesian basiswill inherit these singularities during the change of basisto those coordinates. The fact that these singularitiesare a consequence of the coordinates themselves, and notof the underlying tensor fields, means that they can inprinciple be scaled out of the tensor components in away that is consistent with the adopted reference-metricformulation of BSSN.The goal of the tensor rescaling is to analytically ab-sorb singular terms into the noncoordinate basis. Tensorcomponents are naturally regular with respect to the co-ordinates when expressed in terms of the noncoordinatebasis.In Sec. III A, we show a simple example of rescaling arank-1 tensor with one vector component and a rank-2tensor with one two-dual-vector component in the case ofordinary spherical coordinates. We generalize the rescal-ing procedure to arbitrary coordinate distributions inSec. III B. A. Spherical rescaling examples
In this section, we refer to objects in Cartesian coor-dinates with indices i and j , and to objects in sphericalcoordinates with indices a and b .Consider a finite vector field V with components V i in Cartesian coordinates y ( i ) = ( x, y, z ) and a Cartesiancoordinate basis ∂∂ y ( i ) V = V i ∂∂ y ( i ) ≡ V x ∂∂ x + V y ∂∂ y + V z ∂∂ z . (16)An index in parentheses labels the individual coordinatefunctions or basis vectors, not vector components. Forsimplicity, and without loss of generality, we considera vector that possesses only one nontrivial component( V y = V z = 0) V = V x ∂∂ x . (17)The component V x is a smooth, finite function of theCartesian coordinate values.In Cartesian coordinates, the natural set of noncoor-dinate basis vectors e ( i ) coincides with the coordinatebasis, and therefore V = V x e ( x ) . (18)Now, we transform the vector to ordinary, uniform spher-ical coordinates r ( a ) = ( r, θ, ϕ ), related to the Cartesiancoordinates by x = r sin( θ ) cos( ϕ ) , (19a) y = r sin( θ ) sin( ϕ ) , (19b) z = r cos( θ ) . (19c)This coordinate relationship characterizes the familiarJacobian matrix with components ∂y i ∂r a = sin( θ ) cos( ϕ ) r cos( θ ) cos( ϕ ) − r sin( θ ) sin( ϕ )sin( θ ) sin( ϕ ) r cos( θ ) sin( ϕ ) r sin( θ ) cos( ϕ )cos( θ ) − r sin( θ ) 0 (20)and the inverse Jacobian matrix ∂r a ∂y i satisfies ∂y i ∂r a ∂r a ∂y j = δ ij , ∂y i ∂r b ∂r a ∂y i = δ ab . (21)An application of the ordinary derivative chain rule yieldsthe transformation formula for the coordinate basis vec-tors ∂∂ y ( i ) = ∂r a ∂y i ∂∂ r ( a ) . (22)The Cartesian components V i are transformed from thespherical coordinate components V a using the inverse Ja-cobian V i = V a ∂y i ∂r a . (23) It is this fact, that a tensor’s components transform ina way that is inverse to the basis, which preserves thegeometric meaning of a tensor field represented in anycoordinate system. The resulting vector components inthe spherical coordinate basis are V = V a ∂∂ r ( a ) = V x ∂r a ∂x ∂∂ r ( a ) = V x sin( θ ) cos( ϕ ) ∂∂ r + V x cos( θ ) cos( ϕ ) r ∂∂ θ − V x sin( ϕ ) r sin( θ ) ∂∂ ϕ . (24)The coordinate singularities that have been introducedby the Jacobian become obvious when the componentsare evaluated at the origin r = 0 and along the polaraxis sin( θ ) = 0.We absorb this undesirable behavior into the basis vec-tors as follows. For the noncoordinate spherical basis e ( a ) , choose orthogonal vectors e ( r ) = ∂∂ r , (25a) e ( θ ) = 1 r ∂∂ θ , (25b) e ( ϕ ) = 1 r sin( θ ) ∂∂ ϕ . (25c)In terms of these, the vector components become V = V x sin( θ ) cos( ϕ ) e ( r ) + V x cos( θ ) cos( ϕ ) e ( θ ) − V x sin( ϕ ) e ( ϕ ) , (26)which are manifestly regular over the entire domain.To demonstrate that this strategy is more generallyapplicable, consider a rank-2 tensor W with only onenontrivial component in a Cartesian coordinate two-dual-vector basis W = W xx d x ⊗ d x , (27)where ⊗ is the tensor product and d is the exterior deriva-tive operator acting on the scalar coordinate functions toproduce one-forms d y ( i ) (i.e. dual vectors). The trans-formation rule dual to (22) is d y ( i ) = ∂y i ∂r a d r ( a ) . (28)Two contractions with the inverse Jacobian matrixtransforms the tensor components to the spherical co-ordinate basis W = W xx sin ( θ ) cos ( ϕ ) d r ⊗ d r + 2 W xx r sin( θ ) cos( θ ) cos ( ϕ ) d r ⊗ d θ − W xx r sin ( θ ) sin( ϕ ) cos( ϕ ) d r ⊗ d ϕ + W xx r cos ( θ ) cos ( ϕ ) d θ ⊗ d θ − W xx r sin( θ ) cos( θ ) sin( ϕ ) cos( ϕ ) d θ ⊗ d ϕ + W xx r sin ( θ ) sin ( ϕ ) d ϕ ⊗ d ϕ . (29)Notice that the components of W vanish when evaluatedat, say, r = 0 and θ = π , which amounts to a coordi-nate singularity that destroys information about the ten-sor’s value in other bases (e.g., Cartesian) at these points.This singular behavior is dual to that seen above, wherethe component values of V in the spherical coordinatebasis became unbounded at certain locations.Again, this is ameliorated by an alternative choice ofbasis. The noncoordinate spherical basis one-forms aredefined as the dual to the basis vectors (25) e ( r ) = d r , (30a) e ( θ ) = r d θ , (30b) e ( ϕ ) = r sin( θ ) d ϕ . (30c)In this noncoordinate basis, the tensor components be-come W = W xx sin ( θ ) cos ( ϕ ) e ( r ) ⊗ e ( r ) + 2 W xx sin( θ ) cos( θ ) cos ( ϕ ) e ( r ) ⊗ e ( θ ) − W xx sin( θ ) sin( ϕ ) cos( ϕ ) e ( r ) ⊗ e ( ϕ ) + W xx cos ( θ ) cos ( ϕ ) e ( θ ) ⊗ e ( θ ) − W xx cos( θ ) sin( ϕ ) cos( ϕ ) e ( θ ) ⊗ e ( ϕ ) + W xx sin ( ϕ ) e ( ϕ ) ⊗ e ( ϕ ) . (31)Now, for any point in the spherical domain, it can beshown that all of the above components vanish simulta-neously if and only if W xx = 0.For both V and W —as well as higher rank tensors ingeneral—these arguments extend to arbitrary Cartesiantensors, allowing all components to be nontrivial. B. The general rescaling procedure
The rescaling examples reviewed in the previous sec-tion can be generalized by treating the noncoordinate ba-sis as a collection of matrix operators (and the dual basisas the associated inverse operators), and using them toproject the singularities out of tensor components repre-sented in a coordinate basis. Since the difference betweenbases is only a coordinate transformation, and the BSSNformulation we adopt is covariant, we are free to applythis strategy to remove coordinate singularities from ten-sorial expressions in the formulation. This section re-views the general procedure.We denote the noncoordinate basis vectors by e i ( j ) ,where the index ( j ) lists the individual basis vectors and i labels the components of a particular vector with respectto the underlying coordinate basis. By definition, the setof basis vectors { e ( i ) } is linearly independent and spansthe tangent space at every point in the spatial hypersur-face. We restrict our consideration to time-independent,orthonormal bases. There exists a dual basis e ( i ) j satisfy-ing e ( i ) k e k ( j ) = e ( k ) j e i ( k ) = δ ij . (32) The components of the flat background reference met-ric, represented in a coordinate basis, are related alge-braically to the basis dual vectors viaˆ γ ij = δ kl e ( k ) i e ( l ) j . (33)In this way, the reference metric is constructed as the“product” of basis vectors, or, equivalently, that the basisconstitutes the “square root” of the reference metric [47].We treat this relationship as the definition of the nonco-ordinate basis components in terms of the known flatspace reference-metric components in the correspondingcoordinate basis. The noncoordinate basis, defined assuch, is an orthonormal basis. The components e ( i ) j aresometimes referred to as the “scale factors” of the refer-ence metric, and they contain the singularities associatedwith the coordinates. The components of other tensorsare made regular with respect to the coordinates by fac-toring out the scale factors in the appropriate way. Forexample (1) ε ij = h kl e ( k ) i e ( l ) j (34)and (8) ¯Λ i = ¯ λ j e i ( j ) . (35)The rescaled components are recovered by the inverserelationships h ij = ε kl e k ( i ) e l ( j ) (36)and ¯ λ i = ¯Λ j e ( i ) j . (37)The rescaled tensor components h ij , ¯ λ i , and so on,are regular with respect to the coordinate singularities.(Note that h ij and ¯ λ i in spherical coordinates corresponddirectly to the functions of the same name in [29].) TheBSSN evolution equations (11) are tensorial, and aretherefore independent of the choice of basis.The rescaled fields are those that are integrated anddifferentiated numerically on the coordinate grid, de-scribed next, in Sec. IV. This rescaling procedure enablesus to achieve stable and convergent solutions in a broadclass of singular coordinate systems, as we demonstratein Sec. V. IV. THE SENR/NRPY+ CODE
Numerical relativity codes built to evolve 3+1 initialvalue formulations of Einstein’s equations generally con-tain thousands of lines of code just to express the neededequations for initial data, time evolution, and diagnos-tics. Early incarnations were largely coded by hand, ex-acerbating the already laborious and time-consuming de-bugging process. Most numerical relativity groups havemigrated to automatic code generation, typically relyingon closed-source, proprietary computer algebra systemslike Maple or Mathematica to directly convert tensorialexpressions typed by hand directly in Einstein-like nota-tion into, e.g., highly optimized C code. Kranc [48] is oneexample of a very nice open-source, Mathematica-basedpackage for converting Einstein’s equations—written inEinstein-like notation—into optimized C code.Proprietary packages like Mathematica or Maple re-quire expensive licenses that some users simply cannot af-ford, creating a barrier to entry for potential developers.Further, most numerical relativity simulations are per-formed with supercomputing systems, on which licensesfor, e.g., Mathematica or Maple are not available—againdue to the high licensing cost. As a workaround, numeri-cal relativists will often generate their code locally, usingMathematica or Maple, and then transfer it to the super-computing system—just another way that these licensescan inconvenience users.NRPy+ (“Python-based code generation for numeri-cal relativity and beyond”) aims to address these issues.It is the first open-source [49, 50], non-Mathematica- orMaple-based code generation package for tensorial ex-pressions written in Einstein notation. NRPy+ is writ-ten entirely in Python and depends only on the stan-dard Python computer algebra package SymPy [31] forsymbolic algebra, which is widely available on supercom-puting clusters.If we wish to solve Einstein’s equations in a new coor-dinate system with NRPy+, we need only define the cor-responding reference metric in terms of its scale factors.Using these as input, NRPy+ generates Einstein’s equa-tions in these coordinates and outputs OpenMP-capableC code that is highly optimizable (SIMD vectorized) bycompilers, resulting in a tremendous performance boostcompared to simple serial implementations. NRPy+ alsoleverages SymPy’s ability to eliminate common subex-pressions from complicated algebraic expressions, mini-mizing the number of floating point operations per ex-pression evaluation.SENR (“the Simple, Efficient Numerical Relativitycode”) is a complete, OpenMP-parallelized [32] numeri-cal relativity code, incorporating the C codes generatedby NRPy+ wherever complicated tensorial expressionsare needed. Its skeletal structure makes the algorithmson which numerical relativity codes are based transpar-ent to the user.The division of labor between SENR and NRPy+ isoutlined in Table I, providing a convenient launchingpoint for later subsections that expand on this structure. Both Python 2.7+ and 3.0+ are supported.
A. Computational grid structures
Each of the 24 evolved fields defined in Sec. II are sam-pled by a discrete computational grid, represented as anumerical array storing the function value at each gridpoint. We define a uniformly sampled unit cube gridwith coordinate labels x ( i ) = (cid:0) x , x , x (cid:1) ≡ ( x1 , x2 , x3 ),where x represents the first spatial degree of freedom inEinstein notation and x1 represents the first coordinateas it appears in SENR/NRPy+, and so on for the othertwo coordinates. These coordinates correspond to therescaled tensor basis, in which coordinate singularitieshave been removed. Thus we perform numerical integra-tions and finite difference operations on this grid usingordinary, uniformly spaced stencils. The uniform coordi-nates are mapped to the nonuniformly distributed Carte-sian coordinates y ( i ) = y ( i ) (cid:0) x ( j ) (cid:1) ≡ ( y1 , y2 , y3 ), chosento exploit the near-symmetries of the physical system ofinterest. Tensors in the y ( i ) coordinates exist in a Carte-sian coordinate basis with trivial symmetry and parityconditions (i.e., no inner boundaries).The user specifies the number of grid points N i dedi-cated to each coordinate direction x i , fixing the uniformgrid cell spacing ∆ x i = 1 N i . (38)The current method requires N i to be even and N i ≥ N FD (see Sec. IV A 2), which determines the number of “ghostzone” points N G = N FD / N FD / β i ∂ i ) are ap-proximated by upwinded finite differences, which employasymmetrical stencils with N FD / N FD / − N T i = N i + 2 N G . (39)The total number of points on the uniform grid is simply N T = Y i N T i . (40)The grid points themselves are located at x ( i ) ( j ) = ∆ x i (cid:18) j − N G + 12 (cid:19) , (41)where the grid index j ∈ { , , . . . , N T i − } . This pro-duces a grid that guarantees functions are never evalu-ated on the coordinate singularities at, e.g., r = 0, θ = 0, TABLE I. The division of labor between the SENR C code and the NRPy+ Python code, with a link to the relevant sectiondetailing each task. To perform a simulation, NRPy+ is first run to automatically generate C files containing necessary initialdata, evolution, and diagnostic equations, coupled to highly optimized finite difference codes as needed for spatial derivatives.SENR contains all of the infrastructure needed to make use of these C codes in the context of a full numerical simulation,complete with highly efficient time evolution algorithms, boundary conditions, diagnostics, and checkpointing capabilities.Task Description SectionCoordinates The curvilinear coordinates are defined in terms of the uniform coordinates within
NRPy+ .Only the scale factors (the square root of each diagonal reference-metric component) must bedefined; all hatted quantities in Eq. (11) are evaluated directly from these scale factors. In addi-tion, mappings from the chosen curvilinear coordinate basis to spherical and Cartesian coordinatebases are defined. The former is necessary for transforming initial data (currently expressed ex-clusively in spherical coordinates) to the desired uniform coordinates. The latter is necessary totransform from the chosen basis to evaluate the ADM integrals (expressed in Cartesian coordi-nates for convenience in interpreting the linear and angular momentum components). IV A 1Initial Data Various initial configurations are included inside
NRPy+ , including multiple black hole Brill-Lindquist [51] initial data, conformally curved UIUC [52] initial data for single Kerr black holes,as well as analytical static trumpet initial data [53]. In addition, there are several choices for theinitial gauge conditions. As described above, all initial data are currently written in a sphericalbasis, and converted to the desired curvilinear coordinate basis by
NRPy+ . SENR reads theinitial data code generated by
NRPy+ to define each of the 24 BSSN fields evolved on our grids(“grid functions”) at the initial time. IV CBoundaryConditions
SENR fills the ghost zone points with data from the grid interior for each of the evolved gridfunctions. Specialized boundary condition routines are written by hand for each type of boundarycondition, whether they be inner (e.g., θ <
NRPy+ constructs finite difference stencils at the user-specified accuracy order on the uniform grid. Upwinded derivatives are enabled by default on allshift-advection terms (as is typical; see e.g., [30, 54, 55]). IV A 2EvolutionEquations
NRPy+ constructs and outputs the evolution equation right-hand-side C codes, expressing allspatial derivatives of the grid functions as finite differences. The C codes include reading allneeded data from memory. IIDiagnostics
NRPy+ outputs the needed C codes for the BSSN constraints, the ADM integrands, and thespherically symmetric horizon finder. Additionally,
SENR includes a code that interpolates gridfunction data onto a Cartesian grid to be evaluated using the large suite of diagnostic utilitieswithin the ETK, such as generic horizon finding and gravitational wave extraction. Diagnosticsare periodically evaluated by
SENR and stored to disk. IV BNumericalGrids
SENR allocates memory for coordinate grids and evolved grid functions given the number ofgrid points and the coordinate definitions. IV ATimeIntegration
SENR computes the largest-allowed CFL time step from the proper distance, determined by thereference metric (output from
NRPy+ ) and the chosen grid. It iterates the grid functions tothe next time step, evaluating the evolution equation C codes. Generally, RK4 is used for timeintegrations. IV D or θ = π in spherical coordinates. Points for which allthree coordinates satisfy 0 < x i < x i < x i >
1. Coordinate options
In the present work, we demonstrate BSSN evolutionin three different classes of coordinate system: Cartesian-like, cylindrical-like, and spherical-like. These are dis-tinguished by the number and character of their innerboundary conditions. Boundary conditions must be ap-plied on all faces of our numerical grid cube, whether the face maps to another face (e.g., in the case of the ϕ = 0 and ϕ = 2 π faces in spherical coordinates) orcorresponds to the outer boundary (e.g., at r = r max inspherical coordinates). The details of how to specify theinner and outer boundary ghost zone points are discussedin Sec. IV E.Finite difference stencils are evaluated on the uniform x ( i ) grid, which possesses a one-to-one mapping to thenonuniformly sampled Cartesian coordinates y ( i ) . Forexample, in spherical-like coordinates, the nonuniformgrid is related to the uniform grid via y1 = x = f ( x1 ) sin ( π x2 ) cos (2 π x3 ) , (42a) y2 = y = f ( x1 ) sin ( π x2 ) sin (2 π x3 ) , (42b) y3 = z = f ( x1 ) cos ( π x2 ) . (42c)The coordinate distribution is generalized by the function f which is required to be invertible, odd with respectto the origin, and at least twice differentiable. Theseproperties of f determine the symmetry conditions forthe inner boundaries of the grid, which are describedin detail for spherical-like coordinates and the generalcase in Sec. IV E. The choice f ( x1 ) = r max x1 reducesthis to ordinary, uniform spherical coordinates extend-ing out to radius r max /M >
0. It is often useful toadopt a logarithmically distributed radial coordinate ofthe form f ( x1 ) = A sinh( x1 /w ), where A, w > f ( x1 ) = a x1 + b x1 + c x1 , where appropriate choicesfor the coefficients a , b , and c lead to increasing the rela-tive coordinate density on a spherical shell, which is idealfor sampling a black hole horizon or neutron star crust.The rescaling procedure developed in this paper alsoallows for the angular coordinates to be redistributedin a similar way, but we restrict our discussion to ra-dial rescalings for the sake of simplicity. To be clear,by “radial,” we refer to { x1 , x2 , x3 } in Cartesian-like,to { x1 , x3 } (cylindrical radius and height) in cylindrical-like, and to { x1 } (radius) in spherical-like coordinates.To construct the noncoordinate basis for use in thetensor component rescaling, we start with the flat metricin the y ( i ) Cartesian coordinate basis. Next, transform itto the coordinate basis of the uniform grid x ( i ) using theJacobian matrix. We identify ˆ γ ij with the flat metric inthe x ( i ) basis. Then, the noncoordinate basis componentsare defined by Eq. (33). When the coordinate system isorthogonal, as is the case for all examples in Sec. V, thenonly three of the nine basis components are nonzero.The coordinate system distributions and the corre-sponding basis components are summarized in Table II.
2. Numerical representation of spatial derivatives
The user-specified n = N FD sets the order of the finitedifference stencil approximation, so that the truncationerror scales as E FD ∼ O (cid:0) ∆ x n (cid:12)(cid:12) ∂ n +1 x u (cid:12)(cid:12)(cid:1) (43)for each of the 24 dynamical fields u = { ε ij , ¯ A ij , W, K, ¯Λ i , α, β i , B i } . We demonstratethat finite-difference truncation errors converge to zeroat the prescribed rates in Sec. V B 1. To calculatethe finite difference stencil coefficients, NRPy+ in-verts the corresponding linear system of Taylor series coefficients at user-specified order, akin to invertingthe Vandermonde matrix for Lagrange polynomialinterpolation [56]. Adopting a simple syntax, NRPy+automatically replaces all spatial derivatives that appearin expressions with the appropriate finite differenceapproximation, at the desired order. Such spatialderivatives appear throughout the right-hand sides ofthe evolution system, Eqs. (11), (13), and (14), anddiagnostics—including the BSSN constraint equationsand ADM integrals (see Sec. IV B).We use Kreiss-Oliger dissipation [57, 58] to diffuse un-resolved, high-frequency modes that can reduce the con-vergence order. A standard high-order derivative opera-tor L KO acting on the grid function as L KO u = − ǫ KO ( − n/ n ∆ t (cid:0) ∆ x i (cid:1) n ∂ ni u (44)is added to the right-hand sides of the evolution equa-tions (11), (13), and (14). Note that L KO is not atensorial derivative in general, and its inclusion in theevolution equations violates spatial covariance. How-ever, the coefficient of artificial dissipation is chosen suchthat the contribution of L KO u vanishes in the contin-uum limit. The dimensionless Kreiss-Oliger parameter ǫ KO = ǫ KO (cid:0) x i (cid:1) is allowed to vary smoothly over space,and typically approaches ǫ KO = 0 .
99 in the weak field re-gion. In particular, we often use a spherically symmetrictransition function ǫ KO ( r ) = ǫ KO0 (cid:20) erf (cid:18) r − r KO w KO (cid:19) + 1 (cid:21) , (45)where ǫ KO0 , r KO , w KO > r is a radial coordinate. We generally set this function to beless than 10 − near the origin, so that its nonsmoothnessat the origin is made irrelevant relative to nonsmoothnesscaused by roundoff error. In particular, we usually set ǫ KO0 = 0 . r KO /M = 2, and w KO /M = 0 . B. Diagnostics
We employ a variety of diagnostics to monitor the ac-curacy of our calculations, as well as to probe the physicalproperties of the simulated spacetimes. Diagnostic rou-tines fall broadly into two categories: diagnostics gener-ated in NRPy+, and diagnostic routines within the ETK.
1. Diagnostics generated by NRPy+
The following describes the constraints, the ADM inte-grals, and a spherically symmetric horizon finder, whichare the diagnostics written in Python in NRPy+. Theycontain spatial derivatives of the evolved fields, whichare approximated by the automatically generated finitedifference stencils. The resulting C code is evaluated bySENR during data output, after the time step iteration.0
TABLE II. Summary of coordinate system choices. Finite difference operations take place on the uniformly sampled unit cubegrid x ( i ) = ( x1 , x2 , x3 ), and are mapped to the nonuniformly sampled Cartesian grid y ( i ) = ( y1 , y2 , y3 ). The nontrivial scalefactors constitute the basis, which is used to rescale tensor quantities. The function f allows the coordinates to be redistributedon the nonuniform grid, and the prime mark indicates differentiation with respect to the function argument. Although notshown here, our method allows for the more general case of a different redistribution function for each independent coordinate.Coordinates Definitions Scale Factors y1 y2 y3 e ( x1 ) x1 e ( x2 ) x2 e ( x3 ) x3 Cartesian-like f ( x1 ) f ( x2 ) f ( x3 ) f ′ ( x1 ) f ′ ( x2 ) f ′ ( x3 )Cylindrical-like f ( x1 ) cos (2 π x2 ) f ( x1 ) sin (2 π x2 ) f ( x3 ) f ′ ( x1 ) f ( x1 ) f ′ ( x3 )Spherical-like f ( x1 ) sin ( π x2 ) cos (2 π x3 ) f ( x1 ) sin ( π x2 ) sin (2 π x3 ) f ( x1 ) cos ( π x2 ) f ′ ( x1 ) f ( x1 ) f ( x1 ) sin ( π x2 ) In terms of the BSSN variables (see Sec II), the Hamil-tonian constraint takes the form [29]
H ≡ K − ¯ A ij ¯ A ij + e − φ (cid:0) ¯ R − D i φ ¯ D i φ − D φ (cid:1) = 0 , (46)where ¯ R = ¯ γ ij ¯ R ij , and the momentum constraint is M i ≡ e − φ (cid:18) ˆ D j ¯ A ij + 2 ¯ A k ( i ∆ j ) jk + 6 ¯ A ij ∂ j φ −
23 ¯ γ ij ∂ j K (cid:19) = 0 . (47)The Hamiltonian, momentum, and conformal connectioncoefficient (8) constraints are monitored throughout thesimulation as a measure of numerical accuracy. In ad-dition, the ADM surface integrals for total mass M ADM ,linear momentum P i ADM , and angular momentum J i ADM also serve as diagnostics. The integrands are evaluatedon a spherical surface on the boundary of the spatial hy-persurface at spatial infinity. Numerically, the integralsare approximated by two-dimensional Riemann sums ona spherical surface which is near the outer boundary, ide-ally in the weak field region. Supposing that the space-time is asymptotically flat, and that the spacetime metric g µν approaches the Minkowski metric η µν at least as fastas g µν − η µν = O (1 /r ) when r → ∞ , then the ADMintegrals take the form [42] M ADM = lim r →∞ π ˛ γ ij ( ∂ i γ kj − ∂ k γ ij ) √ γ d S k , (48a) P i ADM = lim r →∞ π ˛ (cid:0) K ij − γ ij K (cid:1) √ γ d S j , (48b) J i ADM = lim r →∞ [ ijk ]8 π ˛ y j ( K kl − γ kl K ) √ γ d S l . (48c) The term ¯ A ik ∆ jjk in Eq. (47) of this paper is missing from themomentum constraint in Eq. (17) of [26] and Eq. (14) of [29]. Inthe notation of their respective articles, the expression in [26] iscorrected by the substitution g → g/ ˜ g , and in [29] by ¯ γ → ¯ γ/ ˆ γ .Be mindful of a parenthesis size mismatch in Eq. (14) of [29]. The vector components d S i play the role of the outward-oriented surface element induced at spatial infinity, y i are the components of a Cartesian coordinate vector,and [ ijk ] is the totally antisymmetric Levi-Civita sym-bol. Note that the ADM integrals are not covariant aswritten, and they must be evaluated in asymptoticallyCartesian coordinates. This allows us to easily interpretthe directionality of the P i ADM and J i ADM components.In the special case of a spherically symmetric configu-ration, the expansion of outgoing null geodesics Θ takesthe simplified form [59]Θ( r ) = 4¯ γ θθ ∂ r φ + ∂ r ¯ γ θθ e φ ¯ γ θθ √ ¯ γ rr − K θθ ¯ γ θθ . (49)The coordinate radius of the apparent horizon r H is de-fined to satisfy Θ( r H ) = 0. Numerically, SENR evaluatesΘ at every grid point using finite difference stencil C codegenerated by NRPy+. Then, it searches for the pair ofneighbors that straddle Θ = 0, and linearly interpolatesbetween those two points to approximate r H .
2. Diagnostics provided by the Einstein Toolkit
A translation layer for the ETK is implemented inSENR, where the fields ¯ γ ij , e φ , ¯ A ij , and K are inter-polated onto a Cartesian grid and then converted tothe ADM quantities γ ij and K ij in the Cartesian ba-sis. These data are fed into the ETK to unlock a widevariety of diagnostic tools [60–62], including the horizonfinder thorn AHFinderDirect [63] and the Ψ gravita-tional waveform extraction thorn WeylScal4 [64, 65]. Themeasured Ψ contains information relating to the gravi-tational wave strain in the transverse-traceless gauge andthe weak field region via [42]Ψ = ¨ h + − i ¨ h × , (50)where h + and h × are the gravitational wave strain am-plitudes of the “plus” and “cross” polarization states,respectively, and the dots denote time derivatives.1 C. Initial data
NRPy+ implements initial data for zero- ( γ ij = ˆ γ ij ),one-, and two-black-hole spacetimes. Single Kerr blackhole initial data are available in UIUC conformally curvedcoordinates [52], Schwarzschild trumpet coordinates [53],and boosted Schwarzschild black holes in isotropic coor-dinates [66]. Two-black-hole initial data take the formof initial black holes at rest (Brill-Lindquist [51]). Allimplemented initial data solve the Hamiltonian (46) andmomentum (47) constraints exactly. Expressions for allinitial data evolved in this work are presented alongsidetheir results in Sec. V.Typically, these initial data types are most naturallyrepresented in the spherical coordinate basis. The Ja-cobian matrix is used to transform the initial data fromuniform spherical coordinates to the desired uniformlysampled x ( i ) grid. Finally, rank-1 and -2 tensors aretransformed to the noncoordinate basis, according to theprocedure in Sec. III. D. Time integration
The evolution equations (11), (13), and (14) are allfirst-order-in-time partial differential equations that maybe written in the form ∂ t u ( t ) = L ( u ( t ) , t ) , (51)where u = { ε ij , ¯ A ij , W, K, ¯Λ i , α, β i , B i } is a vector com-posed of the 24 evolved fields. As can be seen fromEqs. (11), (13), and (14), the differential operator L de-pends on multiple components of u , as well as their firstand second spatial derivatives. All spatial derivatives ofthe evolved fields in L are calculated using finite differ-ences on the uniformly sampled x ( i ) grid.To advance the grid function in time u ( t ) → u ( t + ∆ t ),we adopt the fourth-order Runge-Kutta (RK4) See Ref. [67]. See also [68] for details on the stability propertiesof the fully explicit Runge-Kutta methods at various order. Thewidely used RK4 method is conditionally stable for this appli-cation, although PIRK4 might allow for similar accuracy withlarger timesteps. method [69] k = L ( u ( t ) , t ) , (52a) u = u ( t ) + ∆ t k , (52b) k = L (cid:18) u , t + ∆ t (cid:19) , (52c) u = u ( t ) + ∆ t k , (52d) k = L (cid:18) u , t + ∆ t (cid:19) , (52e) u = u ( t ) + ∆ tk , (52f) k = L ( u , t + ∆ t ) , (52g) u ( t + ∆ t ) = u ( t ) + ∆ t k + 2 k + 2 k + k ) . (52h)Being fourth order means that the error associatedwith the time stepping, at a fixed time, scales as E RK4 ∼ O (cid:0) ∆ t (cid:1) . Immediately after evaluating the RK4steps given by Eqs. (52b), (52d), and (52f), boundaryconditions are applied to u , u , and u , respectively. Ifthe boundary conditions are time dependent, then theyare applied at t + ∆ t/ t + ∆ t on the third substep. At the end of thefull RK4 iteration, boundary conditions are applied to u at time t + ∆ t . Finally, the algebraic correction¯ γ ij → (cid:18) ˆ γ ¯ γ (cid:19) / ¯ γ ij (53)is applied to the metric components to enforce the La-grangian specification constraint (5), where ¯ γ = ˆ γ is theconformal metric determinant on the initial slice.When applying this standard, explicit RK4 algorithm,the Courant-Friedrichs-Lewy (CFL) condition [69] mustbe satisfied. For a numerical grid with coordinates x ( i ) ,SENR finds the smallest proper distance ∆ s min alongeach of the independent coordinate directions∆ s min = min (cid:0) √ ¯ γ ∆ x1 , √ ¯ γ ∆ x2 , √ ¯ γ ∆ x3 (cid:1) , (54)where ∆ x i is the uniform grid spacing between adjacentpoints in the x i direction, and ¯ γ ii is evaluated at x i . Thenthe time step is ∆ t = C ∆ s min , (55)where the Courant factor is set to C = 0 . t ∝ ∆ x i , for cylindrical co-ordinates ∆ t ∝ (cid:0) ∆ x i (cid:1) , and for spherical coordinates∆ t ∝ (cid:0) ∆ x i (cid:1) . The higher-order dependence in the caseof cylindrical and spherical coordinates is due to the fo-cusing of grid points along the symmetry axis or near the2origin. This CFL restriction can be softened significantlyby clever choice of coordinate redistribution function f (see Sec. IV A 1). E. Boundary conditions
As described in Sec. IV A, SENR/NRPy+ maps a uni-formly sampled unit cube grid to a nonuniformly dis-tributed Cartesian coordinate grid, chosen to efficientlysample the space. Our grids are cell centered, so nopoints exist precisely on any of the faces of the unit cube.However, for the points that are nearest to the faces, butinside the cube, the finite difference stencils for spatialderivatives will reach outside of the cube. To ensure thatthese stencils correspond to valid data, we add a collec-tion of points to a shell region exterior to the cube calledthe “ghost zone.” If N G ghost zone points are neededoutside each boundary, then this would increase the to-tal number of points in the grid from N x1 × N x2 × N x3 to( N x1 + 2 N G ) × ( N x2 + 2 N G ) × ( N x3 + 2 N G ).Prior to evaluation of the right-hand sides of the BSSNequations, these ghost zone points must be filled. Someghost zone points, including those at ϕ < ϕ > π on spherical- or cylindrical-like coordinate grids, map topoints inside the cube; we call these inner boundaries.The remaining ghost zone points map back to points out-side the interior of the cube (e.g., r > r max in spherical-like coordinate grids); we call these outer boundaries.Outer boundary ghost zone points may be filled in ac-cordance with the desired outer boundary condition. Al-though the widely used Sommerfeld outer boundary con-dition is also implemented, we find the simple quadraticextrapolation condition to be quite effective on spherical-like coordinate grids u ( x ) ≈ u ( x − ∆ x ) − u ( x − x ) + u ( x − x ) , (56)where u is any of the evolved fields and x is the coordinate x ( i ) perpendicular to the boundary. As with any approx-imate boundary condition, this condition produces un-wanted ingoing modes that contaminate the interior ofthe simulation. In practice, logarithmically spaced radialcoordinates enable us to push the outer boundary out ofcausal contact with the origin for as long as we care tosimulate.Inner boundary conditions depend on the coordinatesystem, and must account for intrinsic periodic, axial,and radial symmetries. The coordinate redistributionfunction f (see Sec. IV A 1) is required to be odd, whichensures that ghost zone points across inner boundariescoincide exactly with other points on the grid interior,respecting the desired symmetries.In the case of scalar functions, these symmetry condi-tions simply copy the appropriate values of the functionfrom the grid interior to its ghost zone partner. Vec-tors and higher rank tensors, however, are sensitive tochanges of sign in the basis when evaluated across innerboundaries. In that case, an appropriate change of sign must be copied into the ghost zone along with the func-tion value itself. We refer to these changes of sign as parity conditions .In the following, we will show by example the symme-try and parity conditions specific to the spherical coor-dinate topology. Then a generic algorithm for assigningghost zone values in arbitrary coordinates is described.
1. Spherical boundary conditions example
To derive the boundary conditions appropriate for a co-ordinate system, first express those coordinates in termsof a system whose boundaries are well understood. Tothis end, we choose ordinary Cartesian coordinates onthe domain −∞ < x, y, z < ∞ . Since each coordinateis unbounded from both above and below, there are noinner boundaries. In addition, every point on the com-putational grid has a unique ( x, y, z ) label.Now consider ordinary spherical coordinates, which arerelated to the Cartesian coordinates in the usual way[Eq. (19)]. The spherical coordinate domain is boundedby 0 < r < ∞ , 0 < θ < π , and 0 < ϕ < π . In thiscase, there is only one outer boundary, corresponding to r → ∞ ; what remains are the five inner boundaries.To find the symmetry conditions for these inner bound-aries, first recall that a scalar function g has a particu-lar value at some location, regardless of the underlyingcoordinate choice. Next, evaluate the function in thespherical coordinate ghost zone, identify the correspond-ing Cartesian coordinate values, and then find the pointin the spherical grid interior that corresponds to thosesame Cartesian coordinates. This links each point in theghost zone to a unique point in the grid interior.By this procedure, the five symmetry conditions forthe spherical coordinate inner boundaries are found tobe g ( − r, θ, ϕ ) = g ( r, π − θ, π + ϕ ) , (57a) g ( r, − θ, ϕ ) = g ( r, θ, π + ϕ ) , (57b) g ( r, π + θ, ϕ ) = g ( r, π − θ, π + ϕ ) , (57c) g ( r, θ, − ϕ ) = g ( r, θ, π − ϕ ) , (57d) g ( r, θ, π + ϕ ) = g ( r, θ, ϕ ) . (57e)These correspond to radial symmetry about the origin[Eq. (57a)], axial symmetry about the north [Eq. (57b)]and south [Eq. (57c)] poles, and periodic symmetryaround the axis in the negative [Eq. (57d)] and positive[Eq. (57e)] orientations. For any point in the ghost zone,there exists a combination of the inner boundary sym-metry rules (57) that maps that point to either the gridinterior or the outer boundary. Note that these symme-tries refer only to the coordinate distribution, and notthe evolved fields, which are allowed to be completelyasymmetrical.As mentioned above, these symmetry conditions aresufficient for filling the inner boundary ghost zone pointsof a scalar function. For the case of vectors and higher3rank tensors, however, parity conditions must also betaken into account.The needed parity conditions are found by comparingthe basis vectors in the ghost zone to their counterpartsin the grid interior. Again, we express the basis vec-tors in terms of a basis in which the parity conditionsare well understood. Since the Cartesian basis has noinner boundaries, the Cartesian basis has only the triv-ial parity conditions, which is to say that there are nochanges of sign. Start with the noncoordinate sphericalbasis e ( a ) (25) which is expressed in terms of the spheri-cal coordinate basis ∂∂ r ( a ) . Then, transform the sphericalcoordinate basis to the Cartesian basis (22) e ( r ) = ∂x i ∂r ∂∂ x ( i ) , (58a) e ( θ ) = 1 r ∂x i ∂θ ∂∂ x ( i ) , (58b) e ( ϕ ) = 1 r sin( θ ) ∂x i ∂ϕ ∂∂ x ( i ) . (58c)Remembering e ( i ) = ∂∂ x ( i ) in Cartesian coordinates andcontracting with the inverse of the Jacobian matrix (20)results in e ( r ) ( r, θ, ϕ ) = sin( θ ) cos( ϕ ) e ( x ) + sin( θ ) sin( ϕ ) e ( y ) + cos( θ ) e ( z ) , (59a) e ( θ ) ( r, θ, ϕ ) = cos( θ ) cos( ϕ ) e ( x ) + cos( θ ) sin( ϕ ) e ( y ) − sin( θ ) e ( z ) , (59b) e ( ϕ ) ( r, θ, ϕ ) = − sin( ϕ ) e ( x ) + cos( ϕ ) e ( y ) , (59c)where we indicate on the left-hand side the explicit func-tional dependence of the spherical noncoordinate basisvectors on ( r, θ, ϕ ). To find the parity conditions, com-pute the dot product between a basis vector in the ghostzone and its partner in the grid interior, according toEq. (57). If the coordinate system is properly con-structed, then this dot product should always evaluateto ±
1, where the negative case indicates that the basischanges sign in the ghost zone.
TABLE III. Cylindrical- and spherical-like parity conditionsfor rank-1 and rank-2 tensors. “Radial” refers to the parityacross r = 0, “Axial” refers to the parity across ρ = 0 orsin( θ ) = 0, and “Periodic” refers to the parity across ϕ = 0or ϕ = 2 π . Compare with Table I in [29].Coordinates Component Radial Axial PeriodicCylindrical-like ρ − + ϕ − + z + + ρρ + + ρϕ + + ρz − + ϕϕ + + ϕz − + zz + +Spherical-like r − + + θ + − + ϕ − − + rr + + + rθ − − + rϕ + − + θθ + + + θϕ − + + ϕϕ + + + For the current example, we find e ( r ) ( − r, θ, ϕ ) · e ( r ) ( r, π − θ, π + ϕ ) = − , (60a) e ( θ ) ( − r, θ, ϕ ) · e ( θ ) ( r, π − θ, π + ϕ ) = +1 , (60b) e ( ϕ ) ( − r, θ, ϕ ) · e ( ϕ ) ( r, π − θ, π + ϕ ) = − , (60c) e ( r ) ( r, − θ, ϕ ) · e ( r ) ( r, θ, π + ϕ ) = +1 , (60d) e ( θ ) ( r, − θ, ϕ ) · e ( θ ) ( r, θ, π + ϕ ) = − , (60e) e ( ϕ ) ( r, − θ, ϕ ) · e ( ϕ ) ( r, θ, π + ϕ ) = − , (60f) e ( r ) ( r, π + θ, ϕ ) · e ( r ) ( r, π − θ, π + ϕ ) = +1 , (60g) e ( θ ) ( r, π + θ, ϕ ) · e ( θ ) ( r, π − θ, π + ϕ ) = − , (60h) e ( ϕ ) ( r, π + θ, ϕ ) · e ( ϕ ) ( r, π − θ, π + ϕ ) = − , (60i) e ( r ) ( r, θ, − ϕ ) · e ( r ) ( r, θ, π − ϕ ) = +1 , (60j) e ( θ ) ( r, θ, − ϕ ) · e ( θ ) ( r, θ, π − ϕ ) = +1 , (60k) e ( ϕ ) ( r, θ, − ϕ ) · e ( ϕ ) ( r, θ, π − ϕ ) = +1 , (60l) e ( r ) ( r, θ, π + ϕ ) · e ( r ) ( r, θ, ϕ ) = +1 , (60m) e ( θ ) ( r, θ, π + ϕ ) · e ( θ ) ( r, θ, ϕ ) = +1 , (60n) e ( ϕ ) ( r, θ, π + ϕ ) · e ( ϕ ) ( r, θ, ϕ ) = +1 . (60o)The same procedure can be used to determine thesymmetry and parity conditions for cylindrical-like co-ordinates. The parity conditions for cylindrical- andspherical-like coordinates across each inner boundary aresummarized in Table III.Operationally, vector and tensor component values arecopied to the ghost zone using the symmetry condition.4For every basis vector attached to each tensor compo-nent, compensate for any changes in sign in the basis byapplying the parity condition.
2. The general boundary condition procedure
The procedure described in the previous section is gen-eralized to construct an explicit routine for filling in-ner and outer boundary ghost zone points in other co-ordinate systems. Here we review the routine, whichwas developed for SENR/NRPy+ to validate the spe-cial, coordinate-specific inner and outer boundary condi-tion routines developed for cylindrical- and spherical-likegrids.Starting with the mapping between Cartesian coordi-nates and the coordinates of interest [e.g., Eq. (19)], theroutine automatically classifies the ghost zones as inneror outer boundaries, maps each inner ghost zone point toits partner on the grid interior, constructs the noncoordi-nate basis vectors, and determines the parity conditions.The symmetry and parity conditions are compiled intoa list at the start of a simulation, and the routine laterruns through the list to apply the boundary conditionsas needed.The ghost zone points are filled, in order, outward fromthe interior grid. This is of particular importance at theouter boundary, where the one-sided quadratic extrapo-lation stencil equation (56) requires that outer boundaryghost zone points be filled in an outward-going direction.In addition, for cases of high symmetry and large N FD ,it is possible for the finite difference stencil to be widerthan the grid interior, which means that a ghost zonepoint could be mapped to another ghost zone point. As-signing the boundary values in an outward fashion avoidsattempts to copy to the ghost zone from uninitializedmemory. V. RESULTS
In this section we present a number of validation andverification tests performed on the SENR/NRPy+ code.In Sec. V A, we compare its results with two other nu-merical relativity codes, and demonstrate that, e.g., dif-ferences in results between SENR/NRPy+ and the codeof Baumgarte et al. [29] are at the level of roundoff errorin the case of ordinary spherical coordinate evolutions ofa strongly perturbed Minkowski spacetime. In Sec. V B,we demonstrate in the contexts of a single, nonspinningblack hole and a double black hole head-on collision thatthe finite difference truncation error converges to zerowith increasing grid resolution at the expected rate, andthat increasing finite difference order with fixed grid res-olution results in near-exponential convergence of the er-ror.These tests also act to showcase the efficiency ofSENR/NRPy+ against other open-source numerical rela- tivity codes in the context these physical scenarios, as alltests were performed on desktop-scale computers, usingat most only about 1 . A. Code comparisons
Here we directly compare SENR/NRPy+ results withtwo other established BSSN evolution codes. These testsverify that all of the evolution equations and diagnos-tics are implemented correctly, and that the simulationsdo indeed contain black hole horizons with the expectedproperties.
1. Robust stability test: Roundoff-level agreement betweenSENR/NRPy+ and BMCCM
This section compares SENR/NRPy+ to the spherical-polar, spatially fourth-order finite differenced BSSN codeof Baumgarte et al. [29] (hereafter BMCCM). BMCCMincludes many features beyond the scope of this pa-per [70–74]; here we focus on initial data that representa version of the robust stability test bed [33–36], whichinvolves a strong random perturbation about flat space-time. At each grid point, the random number genera-tor drand48 [75] is seeded with a unique (but constant)integer tied to the grid index. Then, each of the gridfunctions are populated, in turn, with Minkowski initialdata ( e φ = 1, ε ij = 0, ¯ A ij = 0, K = 0, α = 1, β i = 0,and B i = 0) plus a random value picked from the uni-form distribution [ − . , . is stable in the presence of coordinate sin-gularities [68], we also implement PIRK2 in SENR todirectly compare results with the established BMCCMcode.We compare between SENR/NRPy+ and BMCCM the L norm of the Hamiltonian violation over the entire grid,distilling an entire grid’s worth of data down to a singlenumber at each time. Double precision floating pointarithmetic maintains approximately 16 digits of signifi-cance in any mathematical operation, limiting the extentto which it can be said that two algorithms are in numer-5 D i g it s o f A g r ee m e n t Light-crossing TimeComparisonCalibration 1Calibration 2
FIG. 1. Digits of agreement between SENR/NRPy+ and BM-CCM evolving random initial data, measuring the L normof the Hamiltonian constraint over the entire grid. The initialdata are a random perturbation about flat spacetime of maxi-mum magnitude 0 .
02. The calibration runs compare BMCCMwith itself randomly perturbed, initially, in the least signifi-cant digit. Both codes terminate due to NaNs at exactly thesame iteration. ical agreement. Over the course of a simulation, errorsat the 16th significant digit will gradually rise—an un-avoidable phenomenon when using finite-precision arith-metic known as roundoff error. If two codes are shown toproduce results that agree to within roundoff error, theyare functionally identical, so that if one code has beenproven robust, then the other code possesses identicalrobustness.The comparison is calibrated by first running BMCCMwith the random initial data. Then, it is run again withidentical initial data, except that the least significantdigit at every point is reset to a random number. Thistiny initial difference grows over time due to roundofferror. Given the strong, discontinuous nature of the ran-dom perturbations away from flat space, infinities andNaNs eventually develop, causing the simulations to ter-minate. For good measure, we perform the calibrationa second time, starting with an identical perturbation ofMinkowski, but re-randomizing the least significant digitfor each function at every grid point. The results divergefrom the base case in a similar fashion, shown as dashedlines in Fig 1.Finally, we run the initial data of the BMCCM basecase through SENR/NRPy+ on an identical grid. Theresulting differences are illustrated as a solid line in Fig 1.SENR/NRPy+ maintains significance at least as wellas the BMCCM calibrations, indicating the two codesagree to within roundoff error. Thus the results of theSENR/NRPy+ and BMCCM codes are numerically in-distinguishable.
2. SENR/NRPy+ and the Einstein Toolkit: Comparison ofnonspinning puncture black hole evolutions
In this section, results from SENR/NRPy+ are com-pared with the ETK in the context of nonspinning blackhole evolutions, tracking the apparent horizon radii. Theinitial data represent a single wormhole slice of theSchwarzschild black hole spacetime in isotropic coordi-nates. These data are conformally flat ( ε ij = 0), max-imally sliced ( K = 0), and exist at a moment of timesymmetry ( ¯ A ij = 0). The conformal factor is the solutionto the flat space Laplace equation when the asymptoticflatness condition lim r →∞ e φ = 1 is imposed at infinity e φ = 1 + M r , (61)where r is the isotropic radial coordinate distance to thepuncture located at the origin. The constant of integra-tion is chosen such that the total ADM energy (48a) ofthe slice is M ADM = M . As initial gauge conditions, weuse the popular “pre-collapsed” lapse [30] α = e − φ (62)and vanishing shift β i = 0 and B i = 0. The shift evo-lution damping parameter η (14b) influences the coordi-nate radius of the black hole apparent horizon. We foundthat our choice of η = 0 . /M results in a horizon thatquickly settles down to a static state.In the ETK simulation, we make use of the open-sourceCactus/Carpet [13, 14] AMR infrastructure to place asingle wormhole black hole at the origin, surrounded by8 levels of mesh refinement. The AMR grids are Carte-sian and adopt a Cartesian basis. Each refinement levelcontains the same number of points (excluding ghost andAMR buffer zones), and the grid spacing doubles eachtime a refinement boundary is crossed (starting from theorigin and moving outward). The grid outer boundary isset to r max /M = 128. There are 32 uniformly spaced gridpoints across each refinement level in each coordinate di-rection, giving 16 points across the horizon when at itssmallest (on the initial slice). On the finest refinementlevel, the grid spacing along each of the coordinate direc-tions is ∆ r min /M = 0 . f ( x1 ) = r max sinh ( x1 /w )sinh (1 /w ) , (63)with all tensorial variables expressed in the spherical ba-sis. The grid parameters are tuned to match both the6outer boundary of the ETK simulation, as well as theminimum grid spacing at the black hole. The ETK AMRgrid has 32 × / √ ≈ N x1 = 148and N x2 = N x3 = 2 (the minimum number of angulargrid points required by our boundary condition mod-ule), let w = 0 . r max /M = 128. This results in ∆ r min /M ≈ . N FD = 6 and RK4 time integra-tion in SENR/NRPy+ to match ETK.Despite tuning the grid resolutions and basic numeri-cal evolution strategies to be consistent across codes, thechosen shift condition in both codes is not covariant [notethe partial derivatives appearing in Eqs. (14) and (15),and see discussion in Brown [26] for how to make thisshift condition covariant]. Thus we should not in gen-eral expect results that agree between the codes, as theydo not adopt the same coordinate system. For spheri-cally symmetric spacetimes, however, the partial deriva-tives ∂ i in Eqs. (14) and (15) can be replaced with thecovariant derivative ˆ D i in both spherical and Cartesiancoordinates, showing that, in this special case, the gaugeconditions are geometrically identical.In Fig. 2, we monitor the coordinate radius of a punc-ture black hole as an indicator of the spacetime fieldand shift dynamics (obviously, the shift condition di-rectly influences the coordinate radius of a punctureblack hole). Notice that the apparent horizon radii mea-sured by SENR/NRPy+ and ETK match extremely wellover time , starting from the expected initial coordinateradius of the wormhole throat r H /M = 0 . r H /M ≈ . r H / M t / M SENR/NRPy+Einstein Toolkit FIG. 2. Single black hole wormhole initial data evolved inthe standard gauge with sixth-order finite differencing: com-parison between SENR/NRPy+ and the ETK. Evolutionof the apparent horizon coordinate radius as measured inSENR/NRPy+ by finding the root of the null expansion (49),and in the ETK using the AHFinderDirect thorn [63]. expected.
B. Convergence tests
Section V B 1 shows that finite difference truncationerrors converge to zero with increasing grid resolutionat the expected rates in SENR/NRPy+. We show inSec. V B 2 that the truncation error converges to zeronearly exponentially with linear increase in the finite dif-ference order, keeping the numerical grids held fixed. Weexplore the convergence behaviors of physical quantitiesextracted from the evolved fields. Then, in Sec. V B 3we evolve the dynamical head-on collision from rest oftwo nonspinning black holes, and confirm that the ring-down of the merged black hole gravitational waveformmatches the analytical prediction both in frequency andamplitude.
1. Convergence of puncture black hole evolutions
The BSSN equations are solved on the uniformly sam-pled x ( i ) grid, which directly maps to the solution on thenonuniformly sampled, Cartesian y ( i ) grid. In particular,we use the linear coordinate redistribution function f ( x1 ) = r max x1 , (64)and similarly for x2 or x3 , depending on the choice ofcoordinates. Numerical errors in solving these equationsstem largely from the finite difference representation ofspatial derivatives (i.e., truncation errors of these deriva-tives are typically dominant). Finite difference operatorseffectively fit a polynomial to a function sampled at a7fixed number of neighboring points, so that the deriva-tive of the polynomial acts as an approximation to theexact derivative. Truncation error—i.e., the error causedby approximating functions with finite polynomials of de-gree D —drops as the sample rate to some power N FD that is related to D . In our finite difference schemes, N FD = D .In this section, we confirm that in fact the error dropsin proportion to our underlying uniform grid spacing E FD ∼ O (cid:0) ∆ x N FD (cid:1) , in the context of a single, nonspin-ning puncture black hole evolved using the BSSN formal-ism, in which N FD ∈ { , , } . We monitor the Hamil-tonian constraint violation. With uniform resolution, itbecomes very expensive to simultaneously resolve the ap-parent horizon and push the outer boundary far from thepuncture. This limits the total run time before error fromthe outer boundary contaminates the horizon.As an alternative approach, we perform convergencetests in which the grid resolution is held fixed and thefinite difference order is allowed to vary. In this way,we increase the degree of the finite difference polynomialthat is fit to each function at each point. Therefore, weexpect the convergence rate to be approximately expo-nential, provided roundoff error is sufficiently small andthe underlying functions are smooth (see, e.g., [77] foradditional discussion of “exponential” convergence).These tests adopt the same isotropic wormhole initialdata and initial gauge conditions as in Sec. V A 2. Weset the outer boundary to r max /M = 10. The fastestwaves on the grid (related to the 1+log lapse condition;see, e.g., [19, 46]) propagate, with speed √ α , inwardfrom the outer boundary and outward from the puncture.These simulations end at t/M ≈ N i / N FD , where N i ∈ { , , , } represents the number of gridpoints (in the nonangular coordinate directions) in eachrun. We let N FD = 2 in the Cartesian case, N FD = 4in the cylindrical case, and N FD = 6 in the sphericalcase. The lack of convergence at r/M & anticonvergent inside the hori-zon ( r/M . . outside the horizon,hitting roundoff error at eighth order. We have observednonconvergent behavior propagating outside the horizononly in cases when the total number of points inside thehorizon are set to be so small that finite difference sten-cils outside the horizon touch grid points immediatelysurrounding the puncture.More alarming than the puncture itself, the sharplapse wave (a “gauge shock” [78, 79]) that propagates -12-10-8-6-4-2 0 0 2 4 6 8 10 l og H a m ilt on i a n C on s t r a i n t ( N i / ) x / M 200254322410-12-10-8-6-4-2 0 0 2 4 6 8 10 l og H a m ilt on i a n C on s t r a i n t ( N i / ) ρ / M 200254322410-12-10-8-6-4-2 0 0 2 4 6 8 10 l og H a m ilt on i a n C on s t r a i n t ( N i / ) r / M 200254322410 FIG. 3. Single black hole wormhole initial data evolved in thestandard gauge: convergence of truncation errors to zero ofHamiltonian constraint violation in linearly distributed Carte-sian ( top panel ), cylindrical ( middle panel ), and sphericalcoordinates ( bottom panel ). The legends indicate the num-ber of grid points N i in each (nonangular) direction. The datashown here are scaled by the factor ( N i / N FD . Data aremeasured along a radial line at t/M ≈
5. Cartesian, cylindri-cal, and spherical coordinate evolutions are performed with N FD = 2, 4, and 6, respectively. We adopt the linear redistri-bution function parameter r max /M = 10 [Eq. (64)]. -12-10-8-6-4-2 010 -3 -2 -1 l og H a m ilt on i a n C on s t r a i n t ρ / MFD 2FD 4FD 6FD 8-12-10-8-6-4-2 010 -3 -2 -1 l og H a m ilt on i a n C on s t r a i n t r / MFD 2FD 4FD 6FD 8 FIG. 4. Single black hole wormhole initial data evolvedin the standard gauge: convergence of truncation errors tozero of Hamiltonian constraint violation in sinh-cylindrical( top panel ) and sinh-spherical ( bottom panel ) coordi-nates. Numerical grids are held fixed at moderate resolu-tion, and only finite difference order is increased. Hamilto-nian constraint violation is measured along a radial line at t/M ≈ η = 2 /M . We adopt redistribution parameters r max /M = 1000 and w = 0 . N x1 = 200 and N x2 = N x3 = 2 points, and thecylindrical-like grid uses N x1 = 200, N x2 = 2, and N x3 = 400. outward from the horizon in moving puncture simula-tions is another source of nonconvergence. Remarkably,exponential-like convergence is restored in the wake ofthis gauge wave pulse; in Fig 4, the pulse has reached r/M ≈
425 at the time of measurement. Restorationof convergence after the gauge pulse is very difficult toachieve on Cartesian AMR grids, as sharp outgoing wavesare partially reflected off of refinement boundaries. (See,e.g., [19] for discussion of how this problem might be mit-igated.)Based on these results, we anticipate much cleanerexponential-like convergence in black hole evolutions for which this gauge pulse does not exist. Next we explorejust such a case: evolutions of trumpet black hole initialdata.
2. Convergence of evolved static trumpet initial data
The trumpet solution represents a time-independentslicing of the Schwarzschild spacetime [53]. In this sec-tion, we adopt the trumpet solution to show that the nu-merical evolution outside the horizon is completely dom-inated by truncation error.The trumpet data are conformally flat ( ε ij = 0), anddescribe a single black hole with mass M . With thechoice f = R = M in [53], the trumpet conformalfactor is e φ = r Mr (65)and the nonvanishing extrinsic curvature terms are K = M ( M + r ) , (66a)¯ A rr = − M M + r ) , (66b)¯ A θθ = ¯ A ϕϕ sin ( θ ) = 2 M r M + r ) . (66c)Using an alternative to the standard 1+log conditiongiven by Eq. (13), the lapse is evolved according to acondition consistent with staticity ∂ α = − α (1 − α ) K . (67)For the shift vector evolution equation, we desire onlythat the right-hand sides vanish analytically (althoughnumerical error is expected to result in specious evolu-tion). To this end, we adopt the nonadvecting Gamma-driver condition ∂ t β i = B i , (68a) ∂ t B i = 34 ∂ t ¯Λ i − ηB i . (68b)The initial lapse and shift take on the forms α = rM + r , (69a) β r = M r ( M + r ) . (69b)We use damping parameter η = 0 . /M , although theresults are not very sensitive to its particular value be-cause the evolution begins and remains in a quasistaticstate.Analytically, the trumpet solution with these gaugeconditions is static, but numerical errors result in un-wanted evolution of the fields away from the initial data.9We perform numerical evolutions of these data on fixednumerical grids, subject to these gauge evolution equa-tions, to confirm that evolution away from the initial datadisappears nearly exponentially with increased finite dif-ference order.We choose a fixed, spherical-like coordinate grid with N x1 = 128 and N x2 = N x3 = 2, where the radial pointsare distributed according to Eq. (63) with w = 0 . r max /M = 1000 (the location of the outer boundary).We perform numerical evolutions on these fixed grids atfinite difference orders N FD ∈ { , , , } .After t/M ≈ W , and the gauge functions α and β r , are compared withtheir initial values. The relative differences are shown inFig. 5 for varying N FD . As in the puncture evolution ofthe previous section, we observe nonconvergent numeri-cal errors inside the horizon ( r/M . . r/M = 10 and100.
3. Gravitational waves from a head-on collision
In Sec. V B 1, we demonstrated that nearly exponen-tial convergence outside a puncture black hole horizon(and inside the region causally influenced by the outerboundary) is restored in the wake of a sharp gauge wave.In [18, 19], it is posited that this sharp wave causes non-convergent errors in moving puncture black hole binarysimulations on AMR grids, due to reflections off refine-ment boundaries. These nonconvergent errors have a di-rect impact on the convergence of the gravitational wavesin these simulations.In this section, we explore the convergence of gravita-tional wave signals from a head-on Brill-Lindquist blackhole collision on spherical-like coordinate grids, keepingthe numerical grids fixed at a moderate resolution andvarying only finite difference derivative order, choosing N FD ∈ { , , , , } .Brill-Lindquist initial data [51] represent nonspinningblack holes starting from rest. The initial data are con-structed from a superposition of isotropic wormhole slicesof the Schwarzschild spacetime. Though the formulationholds for an arbitrary number of black holes, here we useonly two.The wormhole slice of the Schwarzschild spacetimewith isotropic radial coordinate r is conformally flat( ε ij = 0) on the initial Cauchy surface. This hyper-surface is maximal ( K = 0) and exists at a moment oftime symmetry ( ¯ A ij = 0). The conformal factor is the -12-10-8-6-4-2 0 l og | ∆ W / W | -12-10-8-6-4-2 0 l og | ∆ α / α | FD 2FD 4FD 6FD 8-12-10-8-6-4-2 010 -4 -3 -2 -1 l og | ∆ β r / β r | r / M FIG. 5. Trumpet black hole initial data evolved for t/M ≈ W ( top panel ), lapse α ( middle panel ), and radial component of the shift β r ( bot-tom panel ), versus distance from the origin. The trumpetblack hole is centered on the origin. solution to the flat-space Laplace equation, which allowsfor a direct superposition of single wormhole conformalfactors given in Eq. (61) e φ = 1 + M + r + + M − r − , (70)where r ± = p r + b ∓ br cos( θ ) (71)is the isotropic radial coordinate distance from the coor-dinate origin to the puncture with mass parameter M ± ± b/M above (+)or below ( − ) the origin. This configuration is axisym-metric with respect to the polar axis, though, as with allof the other cases, we perform our simulation in full 3+1dimensions.For the runs presented here, we focus on an equal-mass case in which M ± = M/
2, so that the ADM massintegral (48a) gives M ADM = M . We confirm that theADM momentum (48b) and angular momentum (48c)integrals vanish, as expected. We also chose b/M = 0 . β i = 0 and B i = 0 for the initial shift. The shift evolution dampingparameter is set to η = 2 /M , which results in an equi-librium remnant horizon radius of r H /M ≈ .
4. We usea spherical-like grid with redistribution function equa-tion (63), so that x1 corresponds to a radial coordi-nate. We choose r max /M = 1000 and w = 0 .
125 with N x1 = 400, N x2 = 64, and N x3 = 2 (where N x3 is the axisof symmetry for the collision). Evolving wormhole initialdata with the ordinary 1+log lapse condition Eq. (13)results in a sharp, nonconvergent gauge pulse that prop-agates outward from the puncture. We find, however,that excellent convergence is restored in the wake of thepulse as it moves towards the outer boundary.The black holes are both nonspinning and are releasedfrom rest, so they remain on the polar axis during in-fall and collide head-on. They merge to form a single,strongly perturbed black hole, which quickly rings downto a stationary state as the time-changing quadrupole inthe horizon is radiated away in the form of gravitationalwaves. The result is that the waveform after merger be-haves as an exponentially damped harmonic oscillator.In spin-weight-two spheroidal harmonics, the fundamen-tal gravitational wave mode in the ℓ = 2 harmonic hascomplex frequency ωM ≈ . − . i [80]. Withonly a constant amplitude A f and phase offset φ f act-ing as fitting parameters, the magnitude of the expectedringdown signal ℜ (cid:0) Ψ rd4 (cid:1) = A f exp( − . t ) cos(0 . t + φ f ) (72)is plotted in the upper panel of Fig. 6, atop the ampli-tude of the real part of the dominant ℓ = 2 mode ofΨ measured in our simulation. Excellent agreement isobserved between the expected ringdown signal and theresults from our simulation over more than six decades inamplitude. Further, the symmetry of the head-on colli-sion is expected to result in gravitational waves that arein a pure + polarization state, which we confirm by mea-suring the imaginary part of Ψ to be zero to roundofferror [see Eq. (50)].The bottom panel of Fig. 6 demonstrates that differ-ences in waveforms at adjacent finite difference orders(keeping the spherical-like coordinate grids fixed at mod-erate resolution) converge nearly exponentially with in-creased finite difference order. Notice that after the peakgravitational wave signal has passed, | FD 6 − FD 8 | and | FD 8 − FD 10 | are at times influenced by roundoff error,as evidenced by their suddenly stochastic behavior. We confirmed this feature by repeating the simulations with long double (80-bit) floating-point precision.The peak amplitude occurs at retarded time of approx-imately t ret /M ≈
18, when the differences between wave-forms at adjacent finite difference orders are near theirpeaks. At this time, this mode of Ψ gains about an or-der of magnitude more precision with each increment offinite difference order. The particular rate of exponentialconvergence depends on the grid spacing.The wave carries away energy from the black hole, soto what degree is convergence in the waveform reflectedin the Hamiltonian (energy) constraint? Figure 7 plotsthe Hamiltonian constraint at the time in which the peakgravitational wave signal crosses the gravitational waveextraction radius ( t/M ≈ . r ext /M ≈ . r/M & t ret /M ≈ | FD 8 − FD 10 | and | FD 6 − FD 8 | curves are easily distinguishable inFig. 6. We conclude that even if the Hamiltonian con-straint is dominated by roundoff error, the gravitationalwaveforms may still be convergent.Based on these results, we conclude thatSENR/NRPy+ would be an excellent tool for studyingperturbed black holes. VI. CONCLUSION
In this work, we extend the reference-metric formula-tion of the BSSN equations pioneered by [24, 26, 28, 29] tohandle Cartesian-, cylindrical-, and spherical-like numer-ical grids. At the heart of this strategy, a noncoordinatebasis is adopted to remove from tensorial variables all co-ordinate singularities that arise from the choice of certainbases. Treating these singularities analytically, we suc-cessfully evolve black holes using the moving puncturesapproach [9, 10, 30] without resorting to special inte-gration techniques, and without encountering numericalinstabilities.We announce a new numerical relativity code packagecalled SENR/NRPy+, which implements this approach.It is fully open source, open development, and nonpro-prietary. NRPy+, written entirely in Python, convertstensorial expressions and their derivatives in Einsteinnotation to optimized C code, representing derivativeswith suitable finite difference approximations. Our cur-rent implementation supports Cartesian-like, spherical-like and cylindrical-like coordinates, but our methodscan be generalized easily for other orthogonal coordinatesystems. SENR contains all of the basic numerical al-1 -13-11-9-7-5 0 50 100 150 200 l og | ∆ R e ( Ψ ) | t ret / M |FD 2 - FD 4||FD 4 - FD 6||FD 6 - FD 8||FD 8 - FD 10|-12-10-8-6-4 l og | R e ( Ψ ) | Measured WaveformExpected Ringdown
FIG. 6. Analysis of gravitational waves from the head-on colli-sion of two puncture black holes resulting from Brill-Lindquistinitial data evolved with the standard gauge conditions. Thecommon axis shows the retarded time t ret ≡ t − r ext , where thegravitational wave extraction radius is r ext /M ≈ .
8. (
Toppanel ) Dominant ( ℓ = 2 , m = 0) mode of Ψ , |ℜ (Ψ ) | ℓ =2 ,m =0 ,plotted atop (cid:12)(cid:12) Ψ rd4 (cid:12)(cid:12) , the quasinormal mode ringdown ampli-tude and frequency expected from black hole perturbationtheory [Eq. (72)]. Numerical results are plotted at finite dif-ference order N FD = 10. ( Bottom panel ) Convergence ofabsolute differences in |ℜ (Ψ ) | ℓ =2 ,m =0 at adjacent finite dif-ferencing orders, keeping the spherical grid fixed at moderateresolution. gorithms needed for a numerical relativity code, makinguse of the C codes generated by NRPy+ where complextensorial expressions are required. To the best of ourknowledge, this is the first open-source numerical rela-tivity code that lets the user select from a broad rangeof curvilinear coordinate systems.The SENR/NRPy+ implementation of the BSSNequations is validated against two other well estab-lished numerical relativity codes. In the context of aMinkowski spacetime with strong random perturbations,we achieve roundoff-level agreement with the BMCCMcode, which evolves the BSSN equations in spherical co-ordinates at fixed fourth-order finite difference accuracy.We similarly observe excellent agreement between the re-sults of SENR/NRPy+ and the ETK’s McLachlan BSSNthorn [15] in the context of a single puncture black holeevolution. We also show that, for both single and dou-ble black hole spacetimes, the finite difference truncation -16-14-12-10-8-6-4-2 010 -2 -1 l og H a m ilt on i a n C on s t r a i n t r / MFD 2FD 4FD 6FD 8FD 10 FIG. 7. Brill-Lindquist head-on collision of two punctureblack holes: Hamiltonian constraint violation at t/M ≈ . r ext /M ≈ .
8, the radius at which gravitationalwaves are extracted in Fig. 6. Spherical grids are fixed atmoderate resolution; only finite differencing order is varied.The black vertical line denotes the approximate location ofthe peak gravitational wave strain. error converges with increasing grid resolution at the ex-pected rate. In addition, we demonstrate exponentialconvergence of the error with increasing finite differenceorder, keeping the grid resolution constant.A number of physical diagnostic quantities are im-plemented in SENR/NRPy+, including constraint vio-lations and ADM integrals. For additional tools, we de-veloped an ETK compatibility layer within SENR thatinterpolates quantities in the chosen curvilinear coordi-nate basis to the Cartesian basis, and onto a Cartesiangrid. In this way, SENR can take direct advantage ofthe large suite of ETK-based diagnostic utilities, includ-ing, e.g., apparent horizon finders and gravitational wavediagnostics. These diagnostics are applied to head-on col-lisions of two nonspinning black holes to show that theabsolute difference between gravitational waveforms con-verges exponentially at successive finite difference order.We conclude that our extended formalism forBSSN on arbitrary coordinate grids as implementedSENR/NRPy+ provides an outstanding tool for analyz-ing perturbed black holes without approximation. Fur-ther, the spherical-like coordinate systems adopted areideal for gravitational wave extraction and analysis. Wenext plan to add coordinate system dynamics and explorebispherical-like coordinate geometries, so that black holebinaries may be modeled with extreme efficiency. Allsimulations displayed in this work can be performed onaging desktop computers (except for a few of the high res-olution Cartesian runs, which were executed on a desk-top with additional RAM), and given the ability of thesecoordinate systems to exploit near-symmetries near com-pact objects, we anticipate that SENR/NRPy+ may bethe first code to unlock the desktop as a powerful tool for2fully general relativistic gravitational wave astrophysics.
ACKNOWLEDGMENTS
We wish to thank Vassilios Mewes, Manuela Campan-elli, and Yosef Zlochower for many useful discussions. Wethank Isabel Cordero-Carri´on and Pedro Cerd´a-Dur´anfor their correspondence regarding the stability of Runge-Kutta integration. We also thank Saul Teukolsky forconversations concerning the expected rates of conver-gence of our finite difference approximations. Fundingfor computer equipment was provided in part by NSFEPSCoR Grant No. OIA-1458952. Some computations were performed on single nodes of West Virginia Univer-sity’s Spruce Knob high-performance computing cluster,funded in part by NSF EPSCoR Research InfrastructureImprovement Cooperative Agreement No. 1003907, thestate of West Virginia (WVEPSCoR via the Higher Ed-ucation Policy Commission), and West Virginia Univer-sity. Some computational resources were also providedby the Extreme Science and Engineering Discovery En-vironment (XSEDE) Jetstream cloud computing infras-tructure at Indiana University, via the XSEDE Cam-pus Champion for West Virginia University account TG-TRA140024. This work was also supported by NSFGrants No. 1402780 and No. 1707526 to Bowdoin Col-lege. [1] B. P. Abbott, R. Abbott, T. D. Abbott, M. R.Abernathy, F. Acernese, K. Ackley, C. Adams,T. Adams, P. Addesso, R. X. Adhikari, andet al., Physical Review Letters , 061102 (2016),arXiv:1602.03837 [gr-qc].[2] B. P. Abbott, R. Abbott, T. D. Abbott, M. R.Abernathy, F. Acernese, K. Ackley, C. Adams,T. Adams, P. Addesso, R. X. Adhikari, andet al., Physical Review Letters , 241103 (2016),arXiv:1606.04855 [gr-qc].[3] The LIGO Scientific Collaboration, the Virgo Collab-oration, B. P. Abbott, R. Abbott, T. D. Abbott,F. Acernese, K. Ackley, C. Adams, T. Adams, P. Ad-desso, and et al., Phys. Rev. Lett. , 221101 (2017),arXiv:1706.01812 [gr-qc].[4] B. P. Abbott, R. Abbott, T. D. Abbott, F. Ac-ernese, K. Ackley, C. Adams, T. Adams, P. Ad-desso, R. X. Adhikari, V. B. Adya, andet al., Physical Review Letters , 141101 (2017),arXiv:1709.09660 [gr-qc].[5] The LIGO Scientific Collaboration, the Virgo Col-laboration, B. P. Abbott, R. Abbott, T. D. Ab-bott, F. Acernese, K. Ackley, C. Adams, T. Adams,P. Addesso, and et al., ArXiv e-prints (2017),arXiv:1711.05578 [astro-ph.HE].[6] B. P. Abbott, R. Abbott, T. D. Abbott, F. Ac-ernese, K. Ackley, C. Adams, T. Adams, P. Ad-desso, R. X. Adhikari, V. B. Adya, andet al., Physical Review Letters , 161101 (2017),arXiv:1710.05832 [gr-qc].[7] B. P. Abbott, R. Abbott, T. D. Abbott, F. Ac-ernese, K. Ackley, C. Adams, T. Adams, P. Ad-desso, R. X. Adhikari, V. B. Adya, and et al.,Astrophysical Journal Letters , L12 (2017),arXiv:1710.05833 [astro-ph.HE].[8] F. Pretorius, Physical Review Letters , 121101 (2005),gr-qc/0507014.[9] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Physical Review Letters , 111101 (2006),gr-qc/0511048.[10] J. G. Baker, J. Centrella, D.-I. Choi,M. Koppitz, and J. van Meter,Physical Review Letters , 111102 (2006),gr-qc/0511103. [11] T. Goodale, G. Allen, G. Lanfermann, J. Mass´o,T. Radke, E. Seidel, and J. Shalf, in Proceedings of the 5th International Conference on High Performance Computing for Computational Science ,VECPAR’02 (Springer-Verlag, Berlin, Heidelberg, 2003)pp. 197–227.[12] E. Schnetter, S. H. Hawley, and I. Hawke,Classical and Quantum Gravity , 1465 (2004),gr-qc/0310042.[13] “Cactus Computational Toolkit,” .[14] “Carpet: Adaptive mesh refinement for the Cactusframework,” .[15] F. L¨offler, J. Faber, E. Bentivegna, T. Bode, P. Di-ener, R. Haas, I. Hinder, B. C. Mundim, C. D. Ott,E. Schnetter, G. Allen, M. Campanelli, and P. Laguna,Classical and Quantum Gravity , 115001 (2012),arXiv:1111.3344 [gr-qc].[16] “Einstein Toolkit Consortium Homepage,” .[17] “SENR/NRPy+ website,” http://math.wvu.edu/~zetienne/SENR/ .[18] Y. Zlochower, M. Ponce, and C. O.Lousto, Phys. Rev. D , 104056 (2012),arXiv:1208.5494 [gr-qc].[19] Z. B. Etienne, J. G. Baker, V. Paschalidis, B. J. Kelly,and S. L. Shapiro, Phys. Rev. D , 064032 (2014),arXiv:1404.6523 [astro-ph.HE].[20] “SpEC: Spectral Einstein Code website,” https://black-holes.org .[21] B. Szil´agyi, International Journal of Modern Physics D , 1430014 (2014),arXiv:1405.3693 [gr-qc].[22] L. E. Kidder, S. E. Field, F. Foucart, E. Schnet-ter, S. A. Teukolsky, A. Bohn, N. Deppe,P. Diener, F. H´ebert, J. Lippuner, J. Miller,C. D. Ott, M. A. Scheel, and T. Vincent,Journal of Computational Physics , 84 (2017),arXiv:1609.00098 [astro-ph.HE].[23] F. Pretorius, Classical and Quantum Gravity , 425 (2005),gr-qc/0407110.[24] S. Bonazzola, E. Gourgoulhon, P. Grandcl´ement,and J. Novak, Phys. Rev. D , 104007 (2004),gr-qc/0307082.[25] M. Shibata, K. Ury¯u, and J. L. Friedman,Phys. Rev. D , 044044 (2004), gr-qc/0407036. [26] J. D. Brown, Phys. Rev. D , 104029 (2009),arXiv:0902.3652 [gr-qc].[27] E. Gourgoulhon, (Springer-Verlag Berlin Heidelberg, 2012).[28] P. J. Montero and I. Cordero-Carri´on, Phys. Rev. D , 124037 (2012),arXiv:1204.5377 [gr-qc].[29] T. W. Baumgarte, P. J. Montero, I. Cordero-Carri´on,and E. M¨uller, Phys. Rev. D , 044026 (2013).[30] M. Alcubierre, B. Br¨ugmann, P. Diener, M. Kop-pitz, D. Pollney, E. Seidel, and R. Takahashi,Phys. Rev. D , 084023 (2003), gr-qc/0206072.[31] A. Meurer, C. P. Smith, M. Paprocki, O. ˇCert´ık, S. B.Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K.Moore, S. Singh, T. Rathnayake, S. Vig, B. E. Granger,R. P. Muller, F. Bonazzi, H. Gupta, S. Vats, F. Johans-son, F. Pedregosa, M. J. Curry, A. R. Terrel, v. Rouˇcka,A. Saboo, I. Fernando, S. Kulal, R. Cimrman, andA. Scopatz, PeerJ Computer Science , e103 (2017).[32] L. Dagum and R. Menon,IEEE Comput. Sci. Eng. , 46 (1998).[33] B. Szil´agyi, R. G´omez, N. T. Bishop, and J. Winicour,Phys. Rev. D , 104006 (2000), gr-qc/9912030.[34] B. Szil´agyi, B. Schmidt, and J. Winicour,Phys. Rev. D , 064015 (2002), gr-qc/0106026.[35] M. Alcubierre, G. Allen, C. Bona, D. Fiske, T. Goodale,F. S. Guzm´an, I. Hawke, S. H. Hawley, S. Husa,M. Koppitz, C. Lechner, D. Pollney, D. Rideout,M. Salgado, E. Schnetter, E. Seidel, H.-a. Shinkai,D. Shoemaker, B. Szil´agyi, R. Takahashi, and J. Wini-cour, Classical and Quantum Gravity , 589 (2004),gr-qc/0305023.[36] “Apples With Apples: Robust Stability Test,” .[37] J. D. Brown, Phys. Rev. D , 044018 (2008),arXiv:0705.1359 [gr-qc].[38] J. D. Brown, Phys. Rev. D , 084042 (2009),arXiv:0908.3814 [gr-qc].[39] T. Nakamura, K. Oohara, and Y. Kojima,Progress of Theoretical Physics Supplement , 1 (1987).[40] M. Shibata and T. Nakamura,Phys. Rev. D , 5428 (1995).[41] T. W. Baumgarte and S. L. Shapiro,Phys. Rev. D , 024007 (1999), gr-qc/9810065.[42] T. W. Baumgarte and S. L. Shapiro, Numerical Relativ-ity: Solving Einstein’s Equations on the Computer (Cam-bridge University Press, New York, NY, USA, 2010).[43] P. Marronetti, W. Tichy, B. Br¨ugmann, J. Gonz´alez,and U. Sperhake, Phys. Rev. D , 064010 (2008),arXiv:0709.2160 [gr-qc].[44] C. Bona, J. Mass´o, E. Seidel, and J. Stela,Physical Review Letters , 600 (1995), gr-qc/9412071.[45] E. Schnetter, Classical and Quantum Gravity , 167001 (2010),arXiv:1003.0859 [gr-qc].[46] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I.Choi, Phys. Rev. D , 124011 (2006), gr-qc/0605030.[47] S. Carroll, Spacetime and Geometry: An Introduction to General Relativity (Addison Wesley, 2004).[48] S. Husa, I. Hinder, and C. Lechner,Computer Physics Communications , 983 (2006),gr-qc/0404023.[49] “The 2-Clause BSD License,” https://opensource.org/licenses/BSD-2-Clause . [50] “The GNU General Public License v3.0,” .[51] D. R. Brill and R. W. Lindquist,Phys. Rev. , 471 (1963).[52] Y. T. Liu, Z. B. Etienne, and S. L. Shapiro,Phys. Rev. D , 121503 (2009).[53] K. A. Dennison and T. W. Baumgarte,Classical and Quantum Gravity , 117001 (2014),arXiv:1403.5484 [gr-qc].[54] J. Thornburg, Classical and Quantum Gravity , 3665 (2004),gr-qc/0404059.[55] M. Chirvasa and S. Husa, ArXiv e-prints (2008),arXiv:0812.3752 [gr-qc].[56] N. Macon and A. Spitzbart,The American Mathematical Monthly , 95 (1958).[57] H. Kreiss and J. Oliger, Methods for the Approximate So-lution of Time Dependent Problems , Global AtmosphericResearch Programme (GARP): GARP Publication Se-ries, Vol. 10 (GARP Publication, 1973).[58] M. Alcubierre,
Introduction to 3+1 Numerical Relativity ,International series of monographs on physics (OxfordUniv. Press, Oxford, 2008).[59] J. Thornburg, Phys. Rev. D , 104007 (1999),gr-qc/9801087.[60] O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnet-ter, Phys. Rev. D , 024018 (2003), gr-qc/0206008.[61] L. B. Szabados, Living Reviews in Relativity , 4 (2004).[62] E. Schnetter, B. Krishnan, and F. Beyer,Phys. Rev. D , 024028 (2006), gr-qc/0604015.[63] J. Thornburg, Classical and Quantum Gravity , 743 (2004),gr-qc/0306056.[64] J. Baker, M. Campanelli, and C. O. Lousto,Phys. Rev. D , 044001 (2002), gr-qc/0104063.[65] D. Pollney, C. Reisswig, E. Schnetter, N. Dor-band, and P. Diener, Phys. Rev. D83 , 044045 (2011),arXiv:0910.3803 [gr-qc].[66] I. Ruchlin, J. Healy, C. O. Lousto, and Y. Zlochower,Phys. Rev. D , 024033 (2017).[67] I. Cordero-Carri´on, P. Cerd´a-Dur´an, and V. Mewes, pri-vate communication.[68] I. Cordero-Carri´on and P. Cerd´a-Dur´an, ArXiv e-prints(2012), arXiv:1211.5930 [math-ph].[69] W. H. Press, S. A. Teukolsky, W. T. Vetterling, andB. P. Flannery, Numerical Recipes 3rd Edition: The Artof Scientific Computing , 3rd ed. (Cambridge UniversityPress, New York, NY, USA, 2007).[70] T. W. Baumgarte, P. J. Montero, andE. M¨uller, Phys. Rev. D , 064035 (2015),arXiv:1501.05259 [gr-qc].[71] T. W. Baumgarte and P. J. Mon-tero, Phys. Rev. D , 124065 (2015),arXiv:1509.08730 [gr-qc].[72] T. W. Baumgarte and C. Gundlach,Physical Review Letters , 221103 (2016),arXiv:1603.04373 [gr-qc].[73] C. Gundlach and T. W. Baum-garte, Phys. Rev. D , 084012 (2016),arXiv:1608.00491 [gr-qc].[74] A. J. Miller and T. W. Baumgarte,Classical and Quantum Gravity , 035007 (2017),arXiv:1607.03047 [gr-qc].[75] C. S. Roberts, The Bell System Technical Journal , 2053 (1982).[76] I. Cordero-Carri´on and P. Cerd´a-Dur´an, SEMA SIMAISpringer Series, (2014). [77] J. Boyd, Chebyshev and Fourier Spectral Methods: Second Revised Edition ,Dover Books on Mathematics (Dover Publications, 2001).[78] M. Alcubierre, Classical and Quantum Gravity , 607 (2003),gr-qc/0210050. [79] M. Alcubierre, Classical and Quantum Gravity , 4071 (2005),gr-qc/0503030.[80] E. Berti, V. Cardoso, and A. O. Starinets,Classical and Quantum Gravity26