Numerical Relativity in Spherical Polar Coordinates: Evolution Calculations with the BSSN Formulation
Thomas W. Baumgarte, Pedro J. Montero, Isabel Cordero-Carrión, Ewald Müller
NNumerical Relativity in Spherical Polar Coordinates:Evolution Calculations with the BSSN Formulation
Thomas W. Baumgarte,
1, 2
Pedro J. Montero, Isabel Cordero-Carri´on, and Ewald M¨uller Max-Planck-Institute f¨ur Astrophysik, Karl-Schwarzschild-Str. 1, D-85748, Garching bei M¨unchen, Germany Bowdoin College, Brunswick, ME 04011, USA
In the absence of symmetry assumptions most numerical relativity simulations adopt Cartesiancoordinates. While Cartesian coordinates have some desirable properties, spherical polar coordi-nates appear better suited for certain applications, including gravitational collapse and supernovasimulations. Development of numerical relativity codes in spherical polar coordinates has been ham-pered by the need to handle the coordinate singularities at the origin and on the axis, for exampleby careful regularization of the appropriate variables. Assuming spherical symmetry and adopting acovariant version of the BSSN equations, Montero and Cordero-Carri´on recently demonstrated thatsuch a regularization is not necessary when a partially implicit Runge-Kutta (PIRK) method is usedfor the time evolution of the gravitational fields. Here we report on an implementation of the BSSNequations in spherical polar coordinates without any symmetry assumptions. Using a PIRK methodwe obtain stable simulations in three spatial dimensions without the need to regularize the origin orthe axis. We perform and discuss a number of tests to assess the stability, accuracy and convergenceof the code, namely weak gravitational waves, “hydro-without-hydro” evolutions of spherical androtating relativistic stars in equilibrium, and single black holes.
PACS numbers: 04.25.dg, 04.70.Bw, 97.60.Jd, 97.60.Lf
I. INTRODUCTION
The first announcements of successful binary black holesimulations [1–3] marked an important break-throughin numerical relativity and triggered a burst of ac-tivity in the field. While most current simulationsadopt some variation of the BSSN formulation [4–6] to-gether with what have become “standard coordinates”(namely 1+log slicing [7] and the “Gamma-driver” con-dition [8]), different implementations differ in many de-tails. Most current, three-dimensional numerical relativ-ity codes share one feature, though, namely Cartesian co-ordinates. While Cartesian coordinates have many desir-able properties, there are applications, for example grav-itational collapse and supernova calculations, for whichspherical polar coordinates would be better suited. Implementing a numerical relativity code in sphericalpolar coordinates poses several challenges. The first chal-lenge lies in the equations themselves. The original ver-sion of the BSSN formulation, for example, explicitly as-sumes Cartesian coordinates (by assuming that the deter-minant of the conformally related metric be one). Thisissue has been resolved by Brown [10], who introduceda covariant formulation of the BSSN equations that iswell-suited for curvilinear coordinate systems (compare[11]).Another challenge is introduced by the coordinate sin- Cartesian or spherical polar coordinates are not the only twopossibilities, of course. In particular, multi-patch applications,combining some of the advantages of both, may present an at-tractive alternative at least for some applications (see, e.g., [9]for an implementation in numerical relativity). gularities at the origin and the axis, which introduce sin-gular terms into the equations. Although the regularityof the metric ensures that, analytically, these terms can-cel exactly, this is not necessarily the case in numericalapplications, and special care has to be taken in order toavoid numerical instabilities.Several methods have been proposed to enforce reg-ularity in curvilinear coordinates. One possible ap-proach is to rely on a specific gauge, e.g. polar-arealgauge [12, 13], together with a suitable choice of thedynamical variables. Numerous different such methodshave been implemented in spherical symmetry (see, e.g.,[14] for an overview); examples in axisymmetry include[12, 15, 16]. This approach has some clear limitations. Itis not obvious how to generalize these methods to relaxthe assumption of axisymmetry; moreover the restrictionof the gauge freedom prevents adoption of the “standardgauge” that proved to be successful in evolutions withthe BSSN formulation.An alternative method is to apply a regularization pro-cedure, by which both the appropriate parity regular-ity conditions and local flatness are enforced in order toachieve the desired regularity of the evolution equations(see [17–23] for examples). Typically, these schemes in-volve the introduction of auxiliary variables as well asfinding evolution equations for these variables. The re-sulting schemes are quite cumbersome, which may ex-plain why, to the best of our knowledge, no such schemehas been implemented without any symmetry assump-tions.In yet an alternative approach, Cordero-Carri´on et.al. [24] recently adopted a partially implicit Runge-Kutta(PIRK) method (see also [25]) to evolve the hyperbolic,wave-like equations in the Fully Constrained formula-tion of the Einstein equations (see [11]). Essentially, a r X i v : . [ g r- q c ] N ov PIRK methods evolve regular terms in the evolutionequations explicitly, and then use these updated valuesto evolve singular terms implicitly (see [26] and SectionIII A below for details). Following this success, Montero& Cordero-Carri´on [27], assuming spherical symmetry,applied a second-order PIRK method to the full set ofthe BSSN Einstein equations in curvilinear coordinates,and produced the first successful numerical simulationsof vacuum and non-vacuum spacetimes using the covari-ant BSSN formulation in spherical coordinates withoutthe need for a regularization algorithm at the origin (orwithout performing a spherical reduction of the equa-tions, compare [28, 29]).In this paper we present a new numerical code thatsolves the BSSN equations in three-dimensional sphericalpolar coordinates without any symmetry assumptions.The code uses a second-order PIRK method to integratethe evolution equations in time. This approach has theadditional advantage that it imposes no restriction onthe gauge choice. We consider a number of test casesto demonstrate that it is possible to obtain stable androbust evolutions of axisymmetric and non-axisymmetricspacetimes without any special treatment at the origin orthe axis.The paper is organized as follows. In Section II wepresent the basic equations; we will review the covari-ant formulation of the BSSN equations, and will thenspecialize to spherical polar coordinates. In Section IIIwe will briefly review PIRK methods and will then de-scribe other specifics of our numerical implementation.In Section IV we present numerical examples, namelyweak gravitational waves, “hydro-without-hydro” simu-lations of static and rotating relativistic stars, and singleblack holes. Finally we summarize and discuss our find-ings in Section V. We also include two appendices; in Ap-pendix A we describe an analytical form of the flat metricin spherical polar coordinates that provides a useful testof the numerical implementation of curvature quantities,while in Appendix B we list the specific source terms forour PIRK method applied to the BSSN equations.Throughout this paper we use geometrized units inwhich c = G = 1. Indices a, b, . . . denote spacetime in-dices, while i, j, . . . represent spatial indices. II. BASIC EQUATIONSA. The BSSN equations in covariant form
We adopt Brown’s covariant form [10] of the BSSNformulation [4–6]. In particular, we write the conformallyrelated spatial metric ¯ γ ij as¯ γ ij = e − φ γ ij , (1)where γ ij is the physical spatial metric, and e φ a con-formal factor. In the original BSSN formulation the de-terminant ¯ γ of the conformally related metric is fixed to unity, which completely determines the conformal factor.This approach is suitable when Cartesian coordinates areused, but not in more general coordinate systems. Wewill pose a different condition on ¯ γ below, but note al-ready that e φ = (¯ γ/γ ) / . (2)The advantage of this approach is that all quantities inthis formalism may be treated as tensors of weight zero(see also [11]). We also denote¯ A ij = e − φ (cid:18) K ij − γ ij K (cid:19) (3)as the conformally rescaled extrinsic curvature. Slightlydeparting from Brown’s approach we assume this quan-tity to be trace-free, while Brown allows ¯ A ij to have anon-zero trace. In the above expression K ij is the phys-ical extrinsic curvature and K = γ ij K ij its trace.Introducing a background connection ˆΓ ijk (compare[11]) we now define∆Γ ijk = ¯Γ ijk − ˆΓ ijk (4)which, unlike the two connections themselves, transformas a tensor field. We also define the trace of these vari-ables as ∆Γ i ≡ ¯ γ jk ∆Γ ijk . (5)It is not necessary for the background connection to beassociated with any metric. In Section II B below we willspecialize to applications in spherical polar coordinatesand hence will assume that the ˆΓ ijk are associated withthe flat metric in spherical polar coordinates. This as-sumption affects the equations in the remainder of thisSection in only one way, namely, we will assume that theRiemann tensor associated with the connection ˆΓ ijk van-ishes, as is appropriate when the background metric isflat.Finally, we define the connection vector ¯Λ i as a newset of independent variables that are equal to the ∆Γ i when the constraint C i ≡ ¯Λ i − ∆Γ i = 0 (6)holds. The vector ¯Λ i plays the role of the “conformal con-nection functions” ¯Γ i in the original BSSN formulation,but, unlike the ¯Γ i , the ¯Λ i transform as a rank-1 tensorof weight zero (compare exercise 11.3 in [14]). In thefollowing we will evolve the variables ¯Λ i as independentvariables, satisfying their own evolution equation.In order to determine the conformal factor e φ we spec-ify the time evolution of the determinant of the conformalmetric. In this paper we adopt Brown’s “Lagrangian”choice ∂ t ¯ γ = 0 . (7)Defining ∂ ⊥ ≡ ∂ t − L β , (8)where L β denotes the Lie derivative along the shift vector β i , we then obtain the following set of evolution equations ∂ ⊥ ¯ γ ij = −
23 ¯ γ ij ¯ D k β k − α ¯ A ij (9a) ∂ ⊥ ¯ A ij = −
23 ¯ A ij ¯ D k β k − α ¯ A ik ¯ A kj + α ¯ A ij K + e − φ (cid:104) − α ¯ D i ¯ D j φ + 4 α ¯ D i φ ¯ D j φ + 4 ¯ D ( i α ¯ D j ) φ − ¯ D i ¯ D j α + α ( ¯ R ij − πS ij ) (cid:105) TF (9b) ∂ ⊥ φ = 16 ¯ D k β i − αK (9c) ∂ ⊥ K = α K + α ¯ A ij ¯ A ij − e − φ ( ¯ D α + 2 ¯ D i α ¯ D i φ )+4 πα ( ρ + S ) (9d) ∂ ⊥ ¯Λ i = ¯ γ jk ˆ D j ˆ D k β i + 23 ∆Γ i ¯ D j β j + 13 ¯ D i ¯ D j β j − A jk ( δ ij ∂ k α − αδ ij ∂ k φ − α ∆Γ ijk ) − α ¯ γ ij ∂ j K − πα ¯ γ ij S j . (9e)(compare equations (21) in [10]). In the above equations, α is the lapse function, ˆ D i denotes a covariant derivativethat is built from the background connection ˆΓ ijk (andhence, in our implementation, associated with the flatmetric ˆ γ ij in spherical polar coordinates) and the super-script TF denotes the trace-free part. The matter sources ρ , S i , S ij and S denote the density, momentum density,stress, and the trace of the stress as observed by a normalobserver, respectively, and are defined by ρ ≡ n a n b T ab , (10a) S i ≡ − γ ia n b T ab , (10b) S ij ≡ γ ia γ jb T ab , (10c) S ≡ γ ij S ij . (10d)Here n a = ( − α, , ,
0) (11)is the normal one-form on a spatial slice, and T ab is thestress-energy tensor.We compute the Ricci tensor ¯ R ij associated with ¯ γ ij from¯ 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)In all of the above expressions we have omitted termsthat include the Riemann tensor ˆ R ijkl associated with the connection ˆΓ ijk , since these terms vanish for our caseof a flat background.The Hamiltonian constraint takes the form H ≡ K − ¯ A ij ¯ A ij + e − φ ( ¯ R − D i φ ¯ D i φ − D φ ) − πρ = 0 , (13)while the momentum constraints can be written as M i ≡ e − φ (cid:16) √ ¯ γ ˆ D j ( √ ¯ γ ¯ A ij ) + 6 ¯ A ij ∂ j φ −
23 ¯ γ ij ∂ j K + ¯ A jk ∆Γ ijk (cid:17) − πS i = 0 . (14)(see equations (16) and (17) in [10]).We note that when ¯ γ = 1 and ˆΓ ijk = 0, which is suit-able for Cartesian coordinates, the above equations re-duce to the traditional BSSN equations. In the following,however, we will evaluate these equations in spherical po-lar coordinates.Before the above equations can be integrated, we haveto specify coordinate conditions for the lapse α and theshift β i . Unless noted otherwise we will adopt a “non-advective” version of what has become the “standardgauge” in numerical relativity. Specifically, we use the“1+log” condition for the lapse [7] in the form ∂ t α = − αK, (15)and the “Gamma-driver” condition for the shift [8] in theform ∂ t β i = B i (16a) ∂ t B i = 34 ∂ t ¯Λ i . (16b)(compare [27]). These (or similar) conditions play akey role in the “moving-puncture” approach to handlingblack hole singularities in numerical simulations. B. Implementation in spherical polar coordinates
We now focus on spherical polar coordinates, and willassume that the ˆΓ ijk are associated with the flat metricin spherical polar coordinates r , θ , and φ ,ˆ γ ij = η ij = r
00 0 r sin θ . (17)Accordingly, the only non-vanishing components of thebackground connection areˆΓ rθθ = − r ˆΓ rφφ = − r sin θ ˆΓ θφφ = − sin θ cos θ ˆΓ θrθ = r − ˆΓ φrφ = r − ˆΓ φφθ = cot θ. (18)When implementing the above equations in sphericalpolar coordinates, care has to be taken that coordinatesingularities do not spoil the numerical simulation. Thesesingularities appear both at the origin, where r = 0, andon the axis where sin θ = 0. Even for a simple scalarwave, appearances of inverse factors of r and sin θ in theLaplace operator can pose a challenge for a numericalimplementation. In Section III below we discuss a PIRKmethod (see also [24, 27]) that handles these singularitiesvery effectively.An additional challenge in general relativity is thatthese inverse factors of r and sin θ appear through thedynamical variables themselves. Components of the spa-tial metric, for example, scale with powers of r and sin θ ,the inverse metric then scales with inverse powers of thesequantities, and numerical error affecting these terms mayeasily spoil the numerical evolution. It is therefore im-portant to treat these appearances of r and sin θ analyt-ically. We therefore factor out suitable powers of r andsin θ from components of all tensorial objects. We start by writing the conformally related metric ¯ γ ij as the sum of the flat background metric ˆ γ ij and a cor-rection (cid:15) ij (which is not assumed to be small),¯ γ ij = ˆ γ ij + (cid:15) ij . (19)The flat metric ˆ γ ij is given by eq. (17), and we write thecorrection (cid:15) ij in the form (cid:15) ij = h rr rh rθ r sin θh rφ rh rθ r h θθ r sin θh θφ r sin θh rφ r sin θh θφ r sin θh φφ . (20)We similarly rescale the extrinsic curvature ¯ A ij as¯ A ij = a rr ra rθ r sin θa rφ ra rθ r a θθ r sin θa θφ r sin θa rφ r sin θa θφ r sin θa φφ , (21)and the connection vector ¯Λ i as¯Λ i = λ r λ θ /rλ φ / ( r sin θ ) . (22)We treat the shift β i and B i similarly, and finally rewritethe evolution equations (9) for the coefficients h ij , a ij and λ i etc.We can compute the connection coefficients (4) from∆Γ ijk = 12 ¯ γ il (cid:16) ˆ D j ¯ γ kl + ˆ D k ¯ γ jl − ˆ D l ¯ γ jk (cid:17) . (23)Since ˆ D i ˆ γ jk = 0 we can compute the derivatives of thespatial metric ˆ D i ¯ γ jk = ˆ D i (cid:15) jk (24) In an alternative approach, one could represent the metric in anorthonormal frame, so that the correct powers of r and sin θ areabsorbed in the unit vectors. in terms of the coefficients h ij . Direct calculation usingthe flat connection (18) yieldsˆ D r ¯ γ rr = h rr,r ˆ D r ¯ γ rθ = rh rθ,r ˆ D r ¯ γ rφ = r sin θh rφ,r ˆ D r ¯ γ θθ = r h θθ,r ˆ D r ¯ γ θφ = r sin θh θφ,r ˆ D r ¯ γ φφ = r sin θh φφ,r ˆ D θ ¯ γ rr = h rr,θ − h rθ ˆ D θ ¯ γ rθ = r ( h rθ,θ + h rr − h θθ )ˆ D θ ¯ γ rφ = r sin θ ( h rφ,θ − h θφ )ˆ D θ ¯ γ θθ = r ( h θθ,θ + 2 h rθ )ˆ D θ ¯ γ θφ = r sin θ ( h θφ,θ + h rφ )ˆ D θ ¯ γ φφ = r sin θh φφ,θ ˆ D φ ¯ γ rr = h rr,φ − θh rφ ˆ D φ ¯ γ rθ = r ( h rθ,φ − cos θh rφ − sin θh θφ )ˆ D φ ¯ γ rφ = r sin θ ( h rφ,φ + sin θh rr + cos θh rθ − sin θh φφ )ˆ D φ ¯ γ θθ = r ( h θθ,φ − θh θφ )ˆ D φ ¯ γ θφ = r sin θ ( h θφ,φ + sin θh rθ + cos θh θθ − cos θh φφ )ˆ D φ ¯ γ φφ = r sin θ ( h φφ,φ + 2 sin θh rφ + 2 cos θh θφ ) (25)The (flat) covariant derivative of the connection vector¯Λ i can similarly be expressed in terms of the λ i asˆ D r ¯Λ r = ∂ r λ r ˆ D θ ¯Λ r = ∂ θ λ r − λ θ ˆ D φ ¯Λ r = ∂ φ λ r − sin θλ φ ˆ D r ¯Λ θ = 1 r ∂ r λ θ ˆ D θ ¯Λ θ = 1 r (cid:0) ∂ θ λ θ + λ r (cid:1) ˆ D φ ¯Λ θ = 1 r (cid:0) ∂ φ λ θ − cos θλ φ (cid:1) ˆ D r ¯Λ φ = 1 r sin θ ∂ r λ φ ˆ D θ ¯Λ φ = 1 r sin θ ∂ θ λ φ ˆ D φ ¯Λ φ = 1 r sin θ (cid:0) ∂ φ λ φ + sin θλ r + cos θλ θ (cid:1) (26)Using the above expressions, we can compute the Riccitensor (12) as follows. In the first term on the right-handside of (12) we write the second covariant derivative of¯ γ ij as a sum of first partial derivatives of the quantitiesˆ D i ¯ γ ij and (flat) connection terms multiplying the ˆ D i ¯ γ ij ,ˆ D k ˆ D l ¯ γ ij = ∂ k ( ˆ D l ¯ γ ij ) (27) − ( ˆ D m ¯ γ ij )ˆΓ mlk − ( ˆ D l ¯ γ mj )ˆΓ mik − ( ˆ D l ¯ γ im )ˆΓ mjk . We then insert the expressions (25) into the first term onthe right-hand side and evaluate all derivatives explicitly,so that these terms can be written in terms of secondpartial derivatives of the coefficients h ij . Once this stephas been completed, we add those remaining terms forwhich the flat background connection (18) is nonzero.The resulting equations are rather cumbersome, andit is easy to introduce typos in the numerical code. Thenumerical examples of Section IV are excellent tests ofthe code. In Appendix A we describe another analyticaltest that we have found very useful to check our imple-mentation of curvature quantities.As a final comment we note that the condition (7) de-termines the time evolution of the determinant ¯ γ of theconformally related metric, but not its initial value. Thelatter can be chosen freely in this scheme, in particu-lar it does not need to be chosen equal to that of thebackground metric ˆ γ (unlike in the original BSSN formu-lation). For some of our numerical simulations, however,in particular for the rotating star simulations of SectionIV B 2, we found that rescaling the conformally relatedmetric so that its determinant becomes ˆ γ = r sin θ im-proved the stability of the simulation, so that it requireda smaller coefficient η in the Kreiss-Oliger dissipationterm (34) below. III. NUMERICAL IMPLEMENTATIONA. PIRK methods
The origin of the numerical instabilities in curvilinearcoordinate systems are related to the presence of stiffsource terms in the equations, e.g. factors of 1 /r or1 / sin ( θ ) that become arbitrary large close to the ori-gin or the axis. In the following we will refer to theseterms as “singular terms”. PIRK methods evolve allother, i.e. regular, terms in the evolution equations ex-plicitly, and then use these updated values to evolve thesingular terms implicitly. This strategy implies that thecomputational costs of PIRK methods are comparableto those of explicit methods. The resulting numericalscheme does not need any analytical or numerical inver-sion, but is able to provide stable evolutions due to itspartially implicit component. We refer to [26] for a de-tailed derivation of PIRK methods (up to third order),and limit our discussion here to a simple description ofthe second-order PIRK method that is implemented inour code.Consider a system of partial differential equations (cid:26) u t = L ( u, v ) ,v t = L ( u ) + L ( u, v ) , (28)where L , L and L are general non-linear differentialoperators. We will denote the corresponding discretizedoperators by L , L and L , respectively. We will furtherassume that L and L contain only regular terms, andhence will update these terms explicitly, whereas the L operator contains the singular terms and will therefore be treated partially implicitly. Note that L is assumedto depend on u only. In the case of the BSSN equationsthis holds for almost all variables; the one exception canbe treated as discussed in the paragraph below equation(B6) in Appendix B, where we provide the exact form ofthe source terms.In our second-order PIRK scheme we update the vari-ables u and v from an old timestep n to a new timestep n + 1 in two stages. In each of these two stages, we firstevolve the variable u explicitly, and then evolve the vari-able v taking into account the updated values of u for theevaluation of the singular L operator. For the system ofequations (28), the first stage u (1) = u n + ∆ t L ( u n , v n ) ,v (1) = v n + ∆ t (cid:20) L ( u n ) + 12 L ( u (1) ) + L ( u n , v n ) (cid:21) , (29)is followed by the second stage u n +1 = 12 (cid:104) u n + u (1) + ∆ t L ( u (1) , v (1) ) (cid:105) ,v n +1 = v n + ∆ t (cid:2) L ( u n ) + L ( u n +1 )+ L ( u n , v n ) + L ( u (1) , v (1) ) (cid:105) . (30)In the first stage, u is evolved explicitly; the updatedvalue u (1) is used in the evaluation of the L operator forthe computation of v (1) . In the second stage, u is againevolved explicitly, and the updated value u n +1 is used inthe evaluation of the L operator for the computation ofthe updated values v n +1 .Our PIRK method is stable as long as the timestep islimited by a Courant condition; see eq. (33) below.We include all singular terms appearing in the sourcesof the equations in the L operator. Firstly, the confor-mal metric components, h ij , the conformal factor, φ , thelapse function, α , and the shift vector, β i , are evolvedexplicitly (as u is evolved in the previous PIRK scheme);secondly, the traceless part of the extrinsic curvature, a ij ,and the trace of the extrinsic curvature, K , are evolvedpartially implicitly, using updated values of α , β i , φ and h ij ; then, the λ i are evolved partially implicitly, usingthe updated values of α , β i , φ , h ij , a ij and K . Finally, B i is evolved partially implicitly, using the updated val-ues of the previous quantities. Lie derivative terms andmatter source terms are always included in the explicitlytreated parts. In Appendix B, we give the exact form ofthe source terms included in each operator. B. Numerical grid
We adopt a centered, fourth-order finite differencingrepresentation of most spatial derivatives. For each gridpoint, the finite-differencing stencil therefore involves thetwo nearest neighbors in each direction (see Fig. 1).An exception from our fourth-order differencing are ad-vective derivatives along the shift, for which we use a
FIG. 1: A schematic representation of our cell-centered gridstructure in spherical polar coordinates, for one fixed valueof φ . Grid points, marked by the crosses, are placed at thecenter of grid cells, so that no grid point ends up at the center( r = 0) or on the axes (sin θ = 0 or sin θ = π ). Our interiorgrid, bordered by solid lines in the figure, covers the region0 ≤ r ≤ r max and 0 ≤ θ ≤ π (as well as 0 ≤ φ ≤ π ). Assuggested by the two highlighted stencils, our fourth-orderdifferencing scheme requires two levels of ghost zones outsideof the interior grid, indicated by the dotted lines. second-order (one-sided) upwind scheme. Because of thesecond-order time evolution, and the second-order advec-tive terms, our scheme is overall second-order accurate,even though for some cases with vanishing shift we havefound that the error appears to be dominated by thefourth-order terms.We adopt a cell-centered grid, as shown schematicallyin Fig. 1. Specifically, we divide the physical domaincovered by our grid, 0 < r < r max , 0 < θ < π and 0 <φ < π into N r × N θ × N φ cells with uniform coordinatesize∆ r = r max /N r , ∆ θ = π/N θ , ∆ φ = 2 π/N φ . (31)Because of our fourth-order finite differencing scheme weneed to pad the interior grid with two layers of ghostzones. Except at the outer boundary, each ghost zonecorresponds to some other zone in the interior of the grid(with some other value of θ and φ ), so that these ghostszones can be filled by copying the corresponding valuesfrom interior grid points.As a concrete example, consider a grid point with an-gular coordinates θ and φ , say, in the innermost radialzone (highlighted by a (blue) filled circle in Fig. 1). Eval-uating the partial derivative with respect to r at this center axisvectors r - + θ + - φ - -tensors rr + + rθ - - rφ + - θθ + + θφ - + φφ + +TABLE I: Parity conditions for components of vectors andtensors as implemented in our coordinate-based code. Com-ponents of vectors and tensors have to be multiplied with thecorresponding sign when they are copied into ghost zones atthe center or the axis. point requires two grid points that, formally, have neg-ative radii. We can fill these two required ghost pointsby finding the corresponding points in the interior of thegrid, which have angular coordinates π − θ and φ + π .Similarly, evaluating a derivative with respect to θ fora point with angular coordinates ( θ, φ ) next to the axis(see the (red) filled square in Fig. 1) requires ghost pointsthat can be filled by finding the corresponding grid pointswith azimuthal angle φ + π in the interior of the grid.For scalar functions the corresponding function valuescan be copied immediately, but for components of vec-tors or tensors, expressed in spherical polar coordinates, apossible relative sign has to be taken into account. Essen-tially, this occurs because, in spherical polar coordinates,the unit vectors may point into the opposite physical di-rection when we identify a ghost zone with an interiorpoint, i.e. when we go from ( θ, φ ) to ( π − θ, φ + π ) or( θ, φ + π ). We list these relative sign changes, as imple-mented in our coordinate-based code, in Table I.We also require two sets of two ghost zones for φ , whichcan be filled directly using periodicity.At the outer boundary we also require two ghost zones,as suggested by the (red) squared stencil in Fig. 1. Weimpose a Sommerfeld boundary condition, which is anapproximate implementation of an outgoing wave bound-ary condition, to fill these ghost zones. In our coordinate-based code we implement this condition by tracing anoutgoing radial characteristic from each of the outerboundary grid points back to the previous time level.We then interpolate the corresponding function to theintersection of that characteristic and the previous timelevel, and copy that interpolated value, multiplied by asuitable fall-off in r , into the boundary grid point. Weassume a fall-off with r − for all metric variables (i.e. h ij , a ij , φ and K ) as well as the lapse α , but a r − fall-off forthe shift β i as well as λ i .The PIRK method of Section III A is stable as longas the time step ∆ t is limited by a Courant-Friedrichs-Lewy condition. In order to evaluate this condition wefirst find the smallest coordinate distance ∆ min betweenany two grid-points in our cell-centered, spherical polargrid. This minimum distance is approximately∆ min = min(∆ r, (∆ r/ θ, (∆ r/
2) sin(∆ θ/ φ ) . (32)We then set ∆ t = C ∆ min , (33)where we have chosen a Courant factor C = 0 . f KO = η t (cid:18) (∆ r ) ∂ f∂r + (∆ θ ) ∂ f∂θ + (∆ φ ) ∂ f∂φ (cid:19) (34)to the right hand side of the evolution equation for eachdynamical variable f . Here η is a dimensionless coeffi-cient which we have chosen between 0 (for some of ourshort-time evolutions) and 0.001 for the rotating neutronstar simulation in Sect. IV B 2. IV. NUMERICAL EXAMPLESA. Weak gravitational waves
As a first test of our codes we consider small-amplitudegravitational waves on a flat Minkowski background. Fol-lowing Teukolsky [31] we construct an analytical, linearsolution for quadrupolar ( (cid:96) = 2) waves from a function F ( r, t ) = A ( r ± t ) e − ( r ± t ) /λ , (35)where the constant A is related to the amplitude of thewave and λ to its wavelength (see also Section 9.1 in[14]). We set λ = 1, by which all length scales becomedimensionless. We will consider axisymmetric ( m = 0)and non-axisymmetric ( m = 2) modes separately.
1. Axisymmetric waves
We first consider axisymmetric m = 0 waves. Sincethese solutions are independent of the coordinate φ , wemay choose N φ as small as possible (which is N φ = 2in our code) without loss of accuracy. We also choose asmall amplitude of A = 10 − , so that deviations from theanalytic solution, which is accurate only to linear orderin A , are dominated by our finite-difference error, andnot by terms that are higher-order in A .In the following we show results for a numerical gridwith (40 N, N,
2) grid points, where N = 1, N = 2 or FIG. 2: Snapshots of the metric coefficient h rr for an axisym-metric m = 0 small-amplitude gravitational wave at differentinstances of time. For this simulation we used a grid of size(160 , ,
2) and imposed the outer boundary at r max = 8 . r in the (arbitrary) direction θ = 1 .
61 and φ = 4 .
71. Differences between the numericalresults (marked by crosses) and the analytical solution (solidlines) are smaller than the width of the lines in this graph. N = 4, and imposing the outer boundary at r = 8 . β i = 0 instead of theGamma-driver condition (16).In Fig. 2 we show snapshots of the metric function h rr at different instances of time for our highest-resolutionsimulation with N = 4. For each time, we include thenumerical results as crosses, as well as the analytical so-lution as a solid line. The differences between the nu-merical results and analytical solution are well below theresolution limit of this graph, so that the two cannot bedistinguished in this Figure.In Fig. 3 we show a convergence test for these waves.Specifically, we compute the L -norm of the differencebetween the analytical solution h rr and the analyticalsolution, || ∆ h rr || = 1 V (cid:18)(cid:90) ( h num rr − h ana rr ) dV (cid:19) / , (36)where V is the coordinate volume of the numerical grid.In Fig. 3 we show these norms as a function of time for N = 1, N = 2 and N = 4. The norms are rescaledwith a factor N ; the convergence of the resulting errorcurves indicates that, at these early times, the error isdominated by the fourth-order differencing of the spatialderivatives. In spherical polar coordinates, the Courant FIG. 3: The norm of the error in the quantity h rr as a func-tion of time for a small-amplitude, axisymmetric gravitationalwave. We show results for simulations with a grid of size(40 N, N, N = 1, N = 2 and N = 4, with the outerboundary imposed at r = 8 .
0. At these early times, the errorappears to be dominated by the fourth-order differencing ofthe spatial derivatives. condition (33) limits the time step to such small valuesthat the second-order errors associated with our PIRKmethod are smaller than the fourth-order error of ourspatial derivatives (for vanishing shift).
2. Nonaxisymmetric waves
Non-axisymmetric gravitational waves represent arare example of an analytical, time-dependent, three-dimensional, albeit weak-field solution to the Einsteinequations. Clearly, this solution represents a stringenttest for our code. In Fig. 4 we show results for an m = 2 wave, againfor an amplitude A = 10 − . As in Fig. 2, we graphsolutions for h rr as functions of r at different instances oftime. Again, our numerical solution (marked by crosses)can hardly be distinguished from the analytical solution(shown as solid lines). For a code in Cartesian coordinates, even a spherically symmetricsolution represents a stringent test, because the symmetry is notreflected by the numerical grid. In our code, however, numericalexpressions simplify for spherical or axisymmetric solutions, sothat they do not test every aspect of the code.
FIG. 4: Snapshots of the metric coefficient h rr for a non-axisymmetric m = 2 small-amplitude gravitational wave atdifferent instances of time. For this simulation we used agrid of size (40 , ,
64) and imposed the outer boundary at r max = 4 .
0; we show data as a function of r in the direction θ = 1 .
62 and φ = 3 .
19. Numerical results are marked by thecrosses, while the analytical solution is shown as the solid line.
B. Hydro-without-hydro
As a test of strong-field, but regular solutions we con-sider spacetimes containing relativistic stars. In general,this requires evolving the stellar matter self-consistentlywith the gravitational fields, for example by solving theequations of relativistic hydrodynamics. Since this is be-yond the scope of this paper, we here adopt the “hydro-without-hydro” approach suggested by [32]. In this ap-proach, which can also be described as an “inverse-Cowling approximation”, we leave the matter sourcesfixed, and evolve only the gravitational fields. In thisway, it is possible to assess the stability of a spacetimeevolution code, and its capability of accurately evolvingstrong but regular gravitational fields in spacetimes withstatic matter, without having to worry about the hydro-dynamical evolution. These simulations serve as both atestbed and a preliminary step towards fully relativistichydrodynamical simulations of stars. In this Section weconsider static and uniformaly rotating stars separately.
1. Spherical neutron stars
We first consider non-rotating relativistic stars, de-scribed by the Tolman-Oppenheimer-Volkoff (TOV) so-lution [33, 34]. We focus on a polytropic TOV star with
FIG. 5: Snapshots of the conformal exponent φ and the lapse α for a rapidly rotating star (see text for details). We showboth functions both at the initial time, and at a late time t = 318 M . We also show both functions along rays in twodifferent directions, one very close to the equator, the otherpointing close to the pole. Both profiles remain very similarto their initial data throughout the evolution. polytropic index Γ = 2, and with a gravitational massof about 85% of the maximum-allowed mass. For thismodel, the central density is about 40% of that of themaximum mass model. We evolved this star with the1+log slicing condition for the lapse (15), but kept theshift fixed to zero. Because the spacetime is sphericallysymmetric, we could choose both N θ and N φ as small aspossible ( N θ = N φ = 2) without loss of accuracy.Even for very modest grid resolutions in the radial di-rection (e.g. N r = 40, with the outer boundary imposedat four times the stellar radius), we found that the grav-itational fields settle down into an equilibrium that issimilar to the initial data. After this initial transition,which is caused by the finite-difference error, the stel-lar surface as well as the outer boundaries (see [32]), thesolution remains stable.
2. Rotating neutron stars
The evolution of the spacetime of a rapidly rotating rel-ativistic star is a more demanding test than the previousone, as it breaks spherical symmetry and instead involvesaxisymmetric non-vacuum initial data in the strong grav-ity regime. The initial data used for this test are the nu-merical solution of a stationary and axisymmetric equi-librium model of a rapidly and uniformly rotating rel- ativistic star [35], which is computed using the Lorenecode [36].We consider a uniformly rotating star with the sameΓ = 2 polytropic equation of state as the non-rotatingmodel of Sect. IV B 1. Our particular model has the samecentral rest-mass density as that non-rotating model, butrotates at 92% of the allowed mass-shedding limit (for astar of that central density); expressed in terms of thegravitational mass M , the corresponding spin period isapproximately 157 M . The ratio of the polar to equa-torial coordinate radii for this model is 0 .
7. For thissimulation we adopted both the 1+log condition for thelapse (15) and the Gamma-driver condition for the shift(16).For this test we adopted a grid of size (48 , , . M , which equals fourtimes the equatorial radius. In Fig. 5 we show the initialand late-time profiles of the conformal exponent φ andthe lapse α , both in a direction close to the equator andclose to the axis. Evidently, both functions remain veryclose to their initial values throughout the evolution, asthey should. C. Schwarzschild
In this Section we present results for two different sim-ulations involving Schwarzschild black holes. In SectionIV C 1 we evolve a Schwarzschild black hole in a “trum-pet” geometry [37–39], which, in the limit of infinite res-olution, is a time-independent solution to the Einsteinequations given our slicing conditions (15). In SectionIV C 2 we adopt wormhole initial data and follow the co-ordinate transition to a trumpet geometry.
1. Trumpet initial data
Maximally sliced trumpet data [38] represent a time-independent slicing of the Schwarzschild spacetime thatsatisfies our slicing condition (15). The solution canbe expressed analytically in isotropic coordinates, albeitonly in parametrized form [40]. In this Section we adoptthese trumpet data as initial data, so that, in the con-tinuum limit, the solution should remain independent oftime.For trumpet data the conformal factor diverges at r = 0. While, on our cell-centered grid, functions arenever evaluated directly at the origin, derivatives in theneighborhood of the singularity at the origin are clearlyaffected by the singular behavior of the conformal fac-tor. However, the great virtue of the “moving-puncture”gauge conditions (15) and (16) is that these errors onlyaffect the neighborhood of the puncture, and do not spoilthe evolution globally [2, 3, 37, 41]. In the following wewill demonstrate these properties in our code using spher-ical polar coordinates.0 FIG. 6: The difference between the maximum of the radialcomponent of the shift vector, β r max , and its initial value β r max | t=0 , as a function of time. For these simulations weused a grid of size (160 N, ,
2) for N = 1, 2, 4 and 8, andimposed the outer boundary at r max = 16 . M . We rescaleall differences with N , so that the convergence of these linesdemonstrates second-order convergence. For the simulations presented in this section weadopted a numerical grid of size size (160 N, ,
2) for N = 1, 2, 4 and 8, with the outer boundary imposedat r max = 16 . M .In Fig. 6 we show results for the maximum of the radialcomponent β r of the shift vector as a function of time.Specifically, we show the difference between these maxi-mum values β r max and their initial values β r max | t =0 . Sinceour trumpet data represent a time-independent solutionto the Einstein equations and our slicing and gauge con-ditions (15) and (16), these differences should convergeto zero as the grid resolution is increased. In Fig. 6 wemultiply the differences with N ; the convergence of theresulting lines therefore demonstrates second-order con-vergence of the simulation. Apparently the error in thesesimulations is dominated by the second-order advectiveterms.We also note that the outer boundary introduces errorterms that depend on both the grid resolution and thelocation of the outer boundary. Since the latter does notdecrease when we increase the grid resolution, the codeconverges more slowly in regions that have come intocausal contact with the outer boundary. We thereforeinclude in Fig. 6 only sufficiently early times, before thelocation of the shift’s maximum is affected by the outerboundary. FIG. 7: Initial data and final profiles of the conformal ex-ponent φ , the lapse function α , and the shift β ˆ r , showingthe coordinate transition from wormhole initial data to time-independent trumpet data. The (blue) long-dashed lines rep-resent the initial data at t = 0, the (red) dashed lines show ournumerical results at time t = 79 M , and the (black) solid linesshow the analytical trumpet solution [40]. The initial dataappear double-valued because we graph this functions as afunction of the areal radius R (see text for details). For thesesimulations we adopted a grid size (10240 , ,
2) with the outerboundary imposed at r = 256 M . (In these graphs we did notinclude the innermost two grid points, which are affected bythe singular behavior of the puncture.)
2. Wormhole initial data
We now turn to evolutions of wormhole initial data,representing a horizontal slice through the Penrose dia-gram of a Schwarzschild black hole. For these data, theconformal factor is given by ψ = 1 + M r , (37)the conformally related metric is flat, ¯ γ ij = η ij , and theextrinsic curvature vanishes, ¯ A ij = 0 = K . Instead ofchoosing the Killing lapse and Killing shift, which wouldleave these data time-independent, we choose, at t = 0,a “pre-collapsed” lapse [8] α = ψ − (38)and a vanishing shift, β i = 0. We then evolve thelapse and the shift with the 1+log condition (15) andthe Gamma-driver condition (16).Since these initial data do not represent a time-independent solution to the Einstein equations together1 FIG. 8: The maximum of the radial shift β r as a functionof time. We show results for different grid sizes (1280 N, , N = 1, 2, 4 and 8, with the outer boundary imposed at256 M . After a brief transition from the initial data β r = 0,the shift settles down into a new equilibrium. For relativelycoarse grid resolutions the shift experiences a slow drift, butthis drift disappears as the grid resolution is increased. with our gauge conditions, we observe a non-trivial timeevolution that represents a coordinate evolution. Forthe “non-advective” 1+log condition (15), this coordi-nate transition results in the maximally sliced trumpetgeometry of Section IV C 2 [39]. In Fig. 7 we show thiscoordinate transition for the conformal exponent φ , thelapse α and the shift β r .We note that some care has to be taken when the nu-merical and analytical results are compared. The an-alytical solution of [40] assumes ¯ γ ij = η ij . We alsochoose ¯ γ ij = η ij in our initial data, but this relationis not necessarily maintained during the time evolution,so that the numerical and analytical solutions may berepresented in different spatial coordinate systems (buton the same spatial slice). In order to compare the twosolutions we therefore graph all quantities as a functionof the gauge-invariant areal radius R . Since for wormholedata each value of R > M corresponds to two values ofthe isotropic radius r , the initial data in Fig. 7 appeardouble-valued. For these comparisons with the analyt-ical solution we also graph the orthonormal componentof the shift β ˆ r rather than the coordinate component β r itself. Fig. 7 clearly shows the coordinate transition fromwormhole initial data to the trumpet equilibrium solu-tion.In Fig. 8 we show the maximum of the radial shift β r as a function of time. After a brief period of a coordinate FIG. 9: Profiles of violations of the Hamiltonian constraint(13) at time t = 79 M . As in Fig. 8 we show results forgrid sizes (1280 N, ,
2) for N = 1, 2, 4 and 8, with the outerboundary imposed at 256 M . All results are rescaled with N ;the convergence of the resulting lines demonstrates second-order convergence of our code. transition the shift settles down into a new equilibrium.We show results for grid sizes (1280 N, ,
2) for N = 1, 2,4 and 8, with the outer boundary imposed at 256 M . Thegraph shows that differences between the different resultsdecrease rapidly as the grid resolution is increased. Forour coarser grid resolutions the shift still experiences aslow drift after the initial transition, but this drift de-creases as the grid resolution is increased.Finally, in Fig. 9, we show profiles of the violations ofthe Hamiltonian constraint (13) at time t = 79 M . In thisgraph all results are rescaled with N ; the convergence ofthe resulting lines demonstrates that the numerical errorin these simulations is again dominated by the second-order implementations of the advective shift terms, andpossibly the time evolution. V. DISCUSSION
In this paper we demonstrate that a PIRK method canbe used to solve the Einstein equations in spherical po-lar coordinates without any need for any regularizationat the origin or on the axis. Specifically, we integratea covariant version of the BSSN equations in three spa-tial dimensions without any symmetry assumptions. Tothe best of our knowledge, these calculations representthe first successful three-dimensional numerical relativitysimulations using spherical polar coordinates. We con-2sider several test cases to assess the stability, accuracyand convergence of the code, namely weak-field “Teukol-sky” gravitational waves, “hydro-without-hydro” simula-tions of static and rotating relativistic stars, and singleblack holes.Spherical polar coordinates have several advantagesand disadvantages over Cartesian coordinates. At leastin single-grid calculations, spherical polar coordinates al-low for a more effective allocation of the numerical gridpoints for applications that involve one center of mass,for example gravitational collapse of single stars or su-pernovae. This is true even for uniform grids, which weadopt in this paper, but curvilinear coordinate systemsalso facilitate the use of non-uniform grids (e.g. a log-arithmic radial coordinate) to achieve a high resolutionnear the origin while keeping the outer boundary suffi-ciently far.Spherical polar coordinates have another strong advan-tage over Cartesian coordinates. In simulations of super-novae or gravitational collapse, for example, the shape ofthe stellar objects is not well represented by Cartesiangrids. This mismatch between the symmetry of the ob-ject and the grid creates direction-dependent numericalerrors, which are observed to trigger m = 4 modes thatgrow in time. Since spherical polar coordinates mimicthe symmetry of collapsing stars more accurately, we ex-pect that this problem can at least be reduced with thesecoordinates.However, spherical polar coordinates also have disad-vantages. One of these disadvantages is of practical na-ture: the equations in spherical polar coordinates includemany more terms than those in Cartesian coordinates,which makes the numerical implementation more cum-bersome and error prone. Spherical polar coordinatesalso introduce coordinate singularities that traditionallyhave created many numerical problems; but these prob-lems can be avoided when using a PIRK method.Perhaps the most severe disadvantage of spherical po-lar coordinates is caused by the Courant limitation on thetime step. As shown in eq. (33), the close proximity ofgrid points close to the origin limits the size of the timesteps ∆ t to increasingly small values as the resolutionis increased. In three-dimensional simulations, ∆ t de-creases approximately with the product ∆ t ∝ ∆ r ∆ θ ∆ φ .This is a severe disadvantage compared to Cartesian co-ordinates where typically ∆ t ∝ ∆ x i . However, this prob-lem is not unique to numerical relativity, and instead iswell-known from dynamical simulations in spherical po-lar coordinates in any field. Accordingly, several differ-ent approaches to either solving or reducing this problemhave been suggested.One possible approach is to reduce the grid resolutionin the angular directions, N θ and N φ , close to the ori-gin. However, for many applications the angular depen-dence of the solution may be independent of the radius,so that this approach might severely limit the accuracyof the results. It may also be possible to replace thePIRK method in a sphere around the origin with a com- pletely implicit scheme, so that the time step there is nolonger limited by the Courant condition (33). Similarimplicit/explicit (IMEX) “split-by-region” schemes havebeen suggested, for example, in [42] in the context ofspectral schemes. Finally, the “Yin-Yang” method sug-gested in [43] mitigates the restrictions imposed by theCourant condition (33) as follows. Note that the small-est physical distance between grid points, which in turnlimits the time step ∆ t , occurs next to the axis. In theYin-Yang method, the unit sphere is therefore coveredby two different grids that are rotated by an angle of 90degrees with respect to each other. Each one covers onlya region around its equator, thereby avoiding the mostsevere limitation on the time step next to the axis, butcombined both grids cover the entire unit sphere.Despite the small time step, however, we have beenable to complete all simulations presented in this papereven with a serial code – in fact, some of our simulationswere performed on a laptop computer. Acknowledgments
TWB and ICC gratefully acknowledge support fromthe Alexander-von-Humboldt Foundation, TWB wouldalso like to thank the Max-Planck-Institut f¨ur Astro-physik for its hospitality. This work was supportedin part by the Deutsche Forschungsgemeinschaft (DFG)through its Transregional Center SFB/TR7 “Gravita-tional Wave Astronomy”, and by NSF grant PHY-1063240 to Bowdoin College.
Appendix A: A numerical test for curvaturequantities in spherical polar coordinates
In spherical polar coordinates, in particular in the ab-sence of any symmetry assumptions, the numerical im-plementation of curvature quantities involves a signifi-cant number of terms that can easily introduce mistakes(see Section II B). One way of testing this part of the nu-merical code is to compare with known analytical solu-tions, for example for the Schwarzschild metric. However,most analytical solutions feature symmetries (e.g. spher-ical symmetry for Schwarzschild) that simplify the prob-lem in the spherical polar coordinates of our code. Asa consequence, many terms vanish identically for thesesolutions, so that not all terms in the code are tested. Inthis Appendix we describe a simple test that is also ana-lytical, but is neither spherically nor axially symmetric,and hence a very stringent test.Starting with the flat metric in Cartesian coordinateswe introduce a coordinate transformation of each coordi-nate x i that only depends on that coordinate itself; the3resulting metric then takes the form¯ γ ij = f ( x ) 0 00 g ( y ) 00 0 h ( z ) , (A1)where f ( x ), g ( y ) and h ( z ) are arbitrary functions. Trans-forming this metric into spherical polar coordinates leadsto a metric for which, in general, all coefficients are non-zero and depend on the coordinates in potentially com-plicated ways.In Cartesian coordinates, the only non-vanishingChristoffel symbols are∆Γ xxx = ¯Γ xxx = f (cid:48) ( x )2 f , (A2)∆Γ yyy = ¯Γ yyy = g (cid:48) ( y )2 g , (A3)∆Γ zzz = ¯Γ zzz = h (cid:48) ( z )2 h (A4)where the prime denotes a derivative with respect to theargument. Given that the ∆Γ ijk transform like tensors,we can obtain the corresponding coefficients in sphericalpolar coordinates with a simple coordinate transforma-tion. For sufficiently general functions f ( x ), g ( y ) and h ( z ), all 18 components of ∆Γ iij in spherical polar coor-dinates will be non-zero. This yields analytical expres-sions for the connection coefficients (23) that can thenbe compared with numerical results.Similarly, the connection functions are given by3 . in Λ i = (cid:18) f (cid:48) ( x )2 f , g (cid:48) ( y )2 g , h (cid:48) ( z )2 h (cid:19) (A5)in Cartesian coordinates, and can be transformed intospherical polar coordinates with a simple coordinatetransformation.Finally, all components of the Ricci tensor in spheri-cal polar coordinates should converge to zero, since themetric (A1) is still flat. In Fig. 10 we show numericalexamples for f ( x ) = 1 + 0 . x ,g ( y ) = 1 + 0 . y , (A6) h ( z ) = 1 + 0 . z . All components of ¯ R ij are non-zero, but converge to zeroas the grid resolution is increased. In the graph we rescaleall results with N , so that the convergence of the result-ing quantities indicates fourth-order convergence of ourimplementation of the Ricci tensor, as expected. Appendix B: Detailed source terms included in thePIRK operators for the evolution equations
We evolve the evolution eqs. (9), (15)-(16) using asecond-order PIRK method. In this Appendix we pro-vide details on how we split the right-hand sides of these
FIG. 10: Values of the norms of the components of the Riccitensor, || R ij || , for the flat metric (A1) with functions (A6),evaluated using grid sizes (16 N, N, N ) for N = 1, 2, 4 and8. All values are rescaled with N , so that the convergence ofthese results indicates fourth-order convergence of our imple-mentation. equations into the explicit and partially implicit opera-tors.We start each time step by evolving the conformal met-ric components, h ij , the conformal factor φ , the lapsefunction, α , and the shift vector, β i , explicitly, i.e., allthe source terms of the evolution equations of these vari-ables are included in the L operator of the second-orderPIRK method.We then evolve the traceless part of the extrinsic cur-vature, a ij , and the trace of the extrinsic curvature, K ,partially implicitly. More specifically, the corresponding L and L operators associated with the evolution equa-tions for a ij and K in terms of the original BSSN variable¯ A ij , related to a ij through eq. (21), are L
2( ¯ A ij ) = e − φ (cid:104) − α ¯ D i ¯ D j φ + 4 α ¯ D i φ ¯ D j φ + 4 ¯ D ( i α ¯ D j ) φ − ¯ D i ¯ D j α + α ¯ R ij (cid:105) TF , (B1) L
3( ¯ A ij ) = −
23 ¯ A ij ¯ D k β k − α ¯ A ik ¯ A kj + α ¯ A ij K, (B2) L K ) = − e − φ ( ¯ D α + 2 ¯ D i α ¯ D i φ ) + α ¯ A ij ¯ A ij , (B3) L K ) = α K . (B4)The λ i are evolved partially implicitly, using the updatedvalues of α , β i , φ , h ij , a ij and K . In terms of the originalBSSN variable ¯Λ i , related to λ i through eq. (22), the4operators are L i ) = ¯ γ jk ˆ D j ˆ D k β i + 13 ¯ D i ¯ D j β j − α ¯ γ ij ∂ j K − A jk ( δ ij ∂ k α − αδ ij ∂ k φ − α ∆Γ ijk ) , (B5) L i ) = 23 ∆Γ i ¯ D j β j . (B6)We note that the evaluation of the Ricci tensor ¯ R ij ineq. (B1) requires updated values of Λ i before they becomeavailable. It is possible to either replace these updatedvalues with old values, or to update the Λ i provisionallyin a purely explicit step, use these values in eq. (B1), but then overwrite these values after the Λ i are updatedpartially implicitly. We have used the latter approach inthe simulations presented in this paper.Finally, the B i are evolved partially implicitly, usingthe updated values of the previous quantities, L B i ) = 34 ∂ t ¯Λ i , (B7) L B i ) = 0 . (B8)Matter source terms and Lie derivative terms are alwaysincluded in the explicitly treated parts. [1] F. Pretorius, Phys. Rev. Lett. , 121101/1 (2005).[2] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo-chower, Phys. Rev. Lett. , 111101/1 (2006).[3] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, andJ. van Meter, Phys. Rev. Lett. , 111102/1 (2006a).[4] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor.Phys. Suppl. , 1 (1987).[5] M. Shibata and T. Nakamura, Phys. Rev. D , 5428(1995).[6] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D ,024007/1 (1998).[7] C. Bona, J. Mass´o, E. Seidel, and J. Stela, Phys. Rev.Lett. , 600 (1995).[8] M. Alcubierre, B. Br¨ugmann, P. Diener, M. Koppitz,D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D , 084023 (2003).[9] E. Schnetter, P. Diener, E. N. Dorband, and M. Tiglio,Class. Quantum Grav. , S553 (2006).[10] J. D. Brown, Phys. Rev. D , 104029/1 (2009).[11] S. Bonazzola, E. Gourgoulhon, P. Grandcl´ement, andJ. Novak, Phys. Rev. D , 104007/1 (2004).[12] J. M. Bardeen and T. Piran, Phys. Reports , 205(1983).[13] M. W. Choptuik, Phys. Rev. D , 3124 (1991).[14] T. W. Baumgarte and S. L. Shapiro, Numerical relativity:Solving Einstein’s equations on the computer (CambridgeUniversity Press, Cambridge, 2010).[15] C. R. Evans, in
Dynamical Spacetimes and NumericalRelativity , edited by J. M. Centrella (Cambridge Univer-sity Press, Cambridge, 1986), pp. 3–39.[16] S. L. Shapiro and S. A. Teukolsky, Phys. Rev. D , 2739(1992).[17] M. W. Choptuik, E. W. Hirschmann, S. L. Liebling, andF. ˜Pretorius, Class. Quantum Grav. , 1857 (2003).[18] O. Rinne and S. J.M., Class. Quantum Grav. , 1143(2005).[19] M. Alcubierre and J. Gonz´alez, Comp. Phys. Comm. , 76 (2005).[20] M. Ruiz, R. Takahashi, M. Alcubierre, and D. N´u˜nez,Gen. Rel. Grav. , 1705 (2008a).[21] O. Rinne, Class. Quantum Grav. , 048003 (2009). [22] M. Alcubierre and M. Mendez, Gen. Rel. Grav. , 2769(2011).[23] E. Sorkin, Phys. Rev. D , 084062 (2010).[24] I. Cordero-Carri´on, P. Cerd´a-Dur´an, and J. Ib´a˜nez, Phys.Rev. D , 044023 (2012).[25] U. M. Ascher, S. J. Ruuth, and B. T. R. Wetton, SIAMJ. Numer. Anal. , 797 (1995).[26] I. Cordero-Carri´on and P. Cerd´a-Dur´an, arXiv preprint1211.5930 (2012).[27] P. J. Montero and I. Cordero-Carri´on, Phys. Rev. D ,124037/1 (2012).[28] D. Garfinkle, C. Gundlach, and D. Hilditch, Class. Quan-tum Grav. , 075007 (2008).[29] S. Bernuzzi and D. Hilditch, Phys. Rev. D , 084003(2010).[30] H.-O. Kreiss and J. Oliger, Global atmospheric researchprogramme publications series (1973).[31] S. A. Teukolsky, Phys. Rev. D , 745 (1982).[32] T. W. Baumgarte, S. A. Hughes, and S. L. Shapiro, Phys.Rev. D , 087501/1 (1999).[33] R. C. Tolman, Phys. Rev. , 364 (1939).[34] J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. ,374 (1939).[35] S. Bonazzola, E. Gourgoulhon, M. Salgado, and J. A.Marck, Astron. Astrophys. , 421 (1993).[36] URL .[37] M. Hannam, S. Husa, D. Pollney, B. Bruegmann, andN. O’Murchadha, Phys. Rev. Lett. , 241102/1 (2007).[38] M. Hannam, S. Husa, N. ´O. Murchadha, B. Br¨ugmann,J. A. Gonz´alez, and U. Sperhake, J. Phys. Conf. Series , 012047/1 (2007).[39] M. Hannam, S. Husa, F. Ohme, B. Br¨ugmann, and N. ´O.Murchadha, Phys. Rev. D , 064020/1 (2008).[40] T. W. Baumgarte and S. G. Naculich, Phys. Rev. D ,067502/1 (2007).[41] J. D. Brown, Phys. Rev. D , 044018/1 (2008).[42] A. Kanevsky, M. H. Carpenter, D. Gottlieb, and J. S.Hesthaven, J. Comp. Phys. , 1753 (2007).[43] A. Wongwathanarat, N. J. Hammer, and E. M¨uller, As-tron. Astrophys.514