An Implementation of DF-GHG with Application to Spherical Black Hole Excision
Maitraya K Bhattacharyya, David Hilditch, K Rajesh Nayak, Sarah Renkhoff, Hannes R. Rüter, Bernd Bruegmann
AAn Implementation of DF-GHG with Application to Spherical Black Hole Excision
Maitraya K Bhattacharyya , , David Hilditch , K Rajesh Nayak , ,Sarah Renkhoff , Hannes R R¨uter , and Bernd Br¨ugmann Indian Institute of Science Education and Research Kolkata, Mohanpur 741246, India Center of Excellence in Space Sciences India, Mohanpur 741246, India Centro de Astrof´ısica e Gravita¸c˜ao - CENTRA,Departamento de F´ısica, Instituto Superior T´ecnico - IST,Universidade de Lisboa - UL, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal Max Planck Institute for Gravitational Physics (Albert Einstein Institute), 14476 Potsdam-Golm, Germany Theoretical Physics Institute, University of Jena, 07743 Jena, Germany (Dated: January 29, 2021)We present an implementation of the dual foliation generalized harmonic gauge (DF-GHG) for-mulation within the pseudospectral code bamps . The formalism promises to give greater freedomin the choice of coordinates that can be used in numerical relativity. As a specific application wefocus here on the treatment of black holes in spherical symmetry. Existing approaches to black holeexcision in numerical relativity are susceptible to failure if the boundary fails to remain outflow.We present a method, called DF-excision, to avoid this failure. Our approach relies on carefullychoosing coordinates in which the coordinate lightspeeds are under strict control. These coordinatesare then combined with the DF-GHG formulation. After performing a set of validation tests in asimple setting, we study the accretion of large pulses of scalar field matter on to a spherical blackhole. We compare the results of DF-excision with a naive setup. DF-excision proves reliable evenwhen the previous approach fails.
I. INTRODUCTION
Free-evolution formulations of GR for numerical rela-tivity (NR) are built with a number of requirements inmind. Foremost in this list is that the specific PDE prob-lem to be solved must be well-posed. The easiest way toguarantee well-posedness of the initial value problem isto try and render the equations hyperbolic so that text-book theorems may be applied. This, in turn, requires achoice of gauge. Considering the popular harmonic gaugechoice (cid:3) X α = 0 we see already that that such a choicerequires a choice of coordinates. But, in case we alreadyhave a choice of coordinates in mind that do not satisfythis condition, the latter may be problematic. It turnsout that what is really required for hyperbolicity is a sen-sible choice of tensor basis. If, with this tensor basis fixedwe change coordinates it turns out that in many cases theequations remain hyperbolic. This strategy is regularlyused within the SpEC numerical relativity code [1, 2] totreat compact binary systems with coordinates that areapproximately corotating with the system, but alwayswith a single foliation of spacetime by a time coordi-nate T . To overcome this restriction one may turn to thedual-foliation (DF) formalism which, as first presentedin [3], allows us to employ a tensor basis associated withcoordinates X α = ( T, X i ) whilst actually working in co-ordinates x α = ( t, x i ). The DF formalism has been usedin a number of places in the literature [4–11] for mathe-matical analysis and is under active investigation for thetreatment of future null infinity.In this paper, we present the first implementation ofthe dual-foliation generalized harmonic gauge (DF-GHG)formulation of GR, which was made in our pseudospec-tral code bamps [12, 13]. In performing the implementa- tion we have made a number of validation tests, a few ofwhich are presented below. But to try and demonstratethe potential of the formalism, we concentrate primarilyon the specific use case of black hole excision. The numer-ical binary black hole breakthrough [14–16] rests, looselyspeaking, on the backs of two different approaches fortreating the strong-field region, black hole excision andthe moving-puncture method. Each has strengths andweaknesses. Excision, as suggested by Unruh to Thorn-burg [17] and developed by many authors, see [18–26] fora selection, relies on the idea that nothing can escapefrom the black hole region, so it should be possible tosimply remove that region from the computational do-main without affecting the domain of outer communica-tion whatsoever. This has the advantage that the mostviolent spacetime region is not treated, and the remain-ing solution may be reasonably expected to be smooth.There are, however, two important requirements to over-come. First, the intuitive idea that nothing can escape needs to be encoded in a formal sense within the equa-tions. This is not trivial because if the excision boundaryflaps around wildly and fails to remain an outflow bound-ary, then we need to give boundary conditions. Even inthe Minkowski spacetime it is possible to introduce anexcision boundary that satisfies the first condition, bysimply taking a sphere and expanding it radially at thespeed of light. Secondly therefore, we must guaranteethat the physical domain is not discarded at the speedof light. For this we need to insure that a small partof the black hole region stays within the computationaldomain. Assuming that the apparent horizon remainsinside spatial slices of the event horizon, this could bedone by making sure that the apparent horizon remainson the grid. A final, less fundamental, but nevertheless a r X i v : . [ g r- q c ] J a n desirable property is to control the coordinate positionof the apparent horizon within the domain.Within the SpEC code these necessary conditions forexcision are enforced by choosing spatial coordinates x i with a control system [27] that monitors the position ofthe apparent horizon and drives the coordinates in a de-sirable direction. This approach is very effective in prac-tice, but as far as we are aware is not guaranteed neverto fail, even in spherical symmetry. It also has the tech-nical disadvantage that because the notion of apparenthorizon is quasilocal, it is not obvious that textbook well-posedness results can be applied directly. In this paperwe use the DF-GHG formulation with coordinates care-fully chosen for excision. Although we work in the spher-ical setting we believe that it may eventually be possibleto use the key ingredients of our method in a more gen-eral context, subsuming our coordinate choice within thecontrol system setup. Unsurprisingly the core point ofour coordinate choice is the use of an area-locking ra-dial coordinate which, combined with insights from thedynamical horizons framework [28, 29] guarantees thefirst two properties mentioned above. In the near-futurethe development presented here also has the importantuse that it will allow us to generalize our earlier per-turbative work [30] on the spherical scalar field withina fixed Schwarzschild background to treat perturbationsrobustly in the fully nonlinear setting, which the naiveapproach of [31] was incapable of. These results will bereported upon elsewhere.We begin in Section II with an overview of the DF-GHG formulation. In Section III we then describe, atthe continuum level, each of the coordinate choices thatwe test in our implementation. In Section IV we give abrief overview overview of the bamps code, before pre-senting our results in Section V. Finally we conclude inSection VI. Geometric units are used throughout. II. THE DUAL FOLIATION FORMULATION
In this section, we provide a brief summary of thedual foliation (DF) formalism. Readers interested ina more detailed approach to the topic may look at [3]and [5]. The principal idea behind the DF approach isto consider two coordinate systems defined in the sameregion of spacetime x µ = ( t, x i ) and X µ = ( T, X i ),hereby referred to as the lower case and upper case co-ordinates respectively. They are shown in Fig. 1. Asa matter of convention, Greek indices go over spaceand time, Latin indices a, b, c, d, e stand for abstract in-dices, whereas i, j, k, l, m, p represent spatial componentsin the x µ basis, and when underlined stand for spatialcomponents in the X µ basis.The two time coordinates t and T define two foliationsof spacetime, the lower case and upper case foliation re-spectively. In practical applications of the DF formalismwe aim to exploit good properties of each coordinate sys-tem. As mentioned in the introduction, in our specific FIG. 1. The DF approach: A spacetime with two dif-ferent slicings and two coordinate systems, the upper casecoordinates (
T, X i ) and the lower case coordinates ( t, x i ). N a and n a denote the timelike unit normal vectors for thetwo slices, the inner product of which is the Lorentz fac-tor W = − ( N a n a ). V a = W (N) ⊥ ba n b and v a = W ⊥ ba N b denote the two boost vectors. setting, this will mean choosing the upper case coordi-nates (and their associated tensor basis) to be generalizedharmonic (cid:3) X α = H α , which is then used to guaranteesymmetric hyperbolicity of the field equations we solve.In later sections we will see that we can then choose thelower case coordinates x µ in a variety of ways, includingchoices that are useful for black hole excision.In the lower case foliation, we can define the lapse,normal vector, time vector, projection operator and shiftvector as α = ( −∇ a t ∇ a t ) − , n a = − α ∇ a t,t a ∇ a t ≡ , ⊥ ab = δ ab + n a n b ,β a = ⊥ ba t b , β i = − αn a ∇ a x i . (1)Similar quantities may be defined in the upper case foli-ation A = ( −∇ a T ∇ a T ) − , N a = − A ∇ a T,T a ∇ a T ≡ , (N) ⊥ ab = δ ab + N a N b ,B a = (N) ⊥ ba T b , B i = − AN a ∇ a X i . (2)The projection operator ⊥ ab with two indices downstairsis the natural induced metric γ ab on the lower case fo-liation and a similar result follows for the upper casefoliation where the naturally induced metric is denotedby (N) γ ab . The covariant derivative associated with γ ab is denoted by D and the corresponding connection is de-noted by Γ. For the upper case spatial metric (N) γ ab ,the associated covariant derivative is (N) D and the corre-sponding connection is (N) Γ.The relationship between the upper case and the lowercase unit normal vector is given by N a = W ( n a + v a ) , n a = W ( N a + V a ) , (3)where W is called the Lorentz factor and is defined as W = − ( N a n a ) = 1 √ − v i v i = 1 (cid:112) − V i V i , (4)and V a and v a are the upper case and lower case boostvectors defined as v a = 1 W ⊥ ba N b , V a = 1 W (N) ⊥ ba n b . (5)The Jacobian matrix, defined as J µµ ≡ ∂X µ /∂x µ can bedecomposed in the 3 + 1 form n α J αα N α = − W, n α J iα ≡ π i ,J αi N α = W v i , J ii ≡ φ ii . (6)In matrix form, the Jacobian can be represented as J = (cid:18) A − W ( α − β i v i ) απ i + β i φ ii − A − W v i φ ii (cid:19) , (7)with the inverse being J − = (cid:18) α − W ( A − B i V i ) A Π i + B i Φ ii − α − W V i Φ ii (cid:19) . (8)Note that the quantities π i and Π i can be written interms of the lapse, shift and boost vectors π i = W V i − W A − B i , Π i = W v i − W α − β i . (9)Another important result we will need is that for a firstorder evolution system in upper case coordinates of theform ∂ T u = ( A A p + B p ) ∂ p u + A S , (10)where u is the state vector, A p are the principal matrices,and S contains the source terms, can be rewritten interms of the lower case coordinates as( + A V ) ∂ t u = αW − (cid:16) A p ( ϕ − ) pp − ( + A V )Π p (cid:17) ∂ p u + αW − S , (11)where ϕ ii = (N) γ iµ J µi is called the projected Jacobianand A V ≡ A i V i . A. DF in the generalized harmonic formulation
In this subsection, we look at the generalized harmonicformalism employed using the dual foliation approach.Our discussion will closely follow [5] but with the addition of the γ parameter, which because of the subtle asymp-totics on hyperboloidal slices was earlier hard-coded tovanish. We start essentially with the first order GHGequations of [32] which are in turn based on earlier workof Garfinkle [33]. With the γ parameter turned back on,they read ∂ T g µν = (1 + γ ) B i ∂ i g µν + AS ( g ) µν ,∂ T Φ iµν = B j ∂ j Φ iµν − A∂ i Π µν + γ A∂ i g µν + AS (Φ) iµν ,∂ T Π µν = γ γ B i ∂ i g µν + B i ∂ i Π µν − A (N) γ ij ∂ i Φ jµν + AS (Π) µν , (12)where the source terms are given by S ( g ) µν = − Π µν − γ A − B i Φ iµν ,S (Φ) iµν = − γ Φ iµν + 12 N α N β Φ iαβ Π µν + (N) γ jk N α Φ ijα Φ kµν ,S (Π) µν = 2 g αβ (cid:16) (N) γ ij Φ iαµ Φ jβν − Π αµ Π βν − g δγ Γ µαδ Γ νβγ (cid:17) − (cid:18) ∇ ( µ H ν ) + γ Γ αµν C α − γ g µν Γ α C α (cid:19) − N α N β Π αβ Π µν − N α (N) γ ij Π αi Φ jµν + γ (cid:104) δ α ( µ N ν ) − g µν N α (cid:105) C α − γ γ A − B i Φ iµν , (13)where we haveΓ αµν ≡ (N) γ i ( µ | Φ i | ν ) α − (N) γ iα Φ iµν + N ( µ Π ν ) α − N α Π µν , (14)and the equations are subject to both the reduction con-straints C iµν = ∂ i g µν − Φ iµν , (15)and the GHG constraints C µ = g αβ Γ µαβ + H µ = 0 . (16)The functions H µ are called the gauge source functionswhich are functions of the coordinates and the metric.Now, considering these evolution equations to be in thestandard form of Eqn. (10), we can easily obtain the A p matrices which is a slight modification from that givenin [5] A p = γ A − B p γ δ pi − δ pi γ γ A − B p − (N) γ pj . (17)To avoid repeating the calculation in [5], we observethat by inclusion of the γ terms within the coordinatechange (11) is straightforwardly done by modifying theform given in [5] with the Sherman-Morrison formula.Doing so we arrive at the lower case time evolution equa-tions ∂ t g µν = (cid:18) β p − αv p + γ αB p W ( A + B V ) ( ϕ − ) pp (cid:19) ∂ p g µν + αW − s ( g ) µν ,d t Φ iµν = (cid:16) β p δ ji − αv p δ ji + αW v i ( g − ) pj (cid:17) d p Φ jµν + αW − g pi (cid:16) γ ∂ p g µν − ∂ p Π µν (cid:17) + αW − s (Φ) iµν ,∂ t Π µν = β p ∂ p Π µν − αW ( g − ) pi d p Φ iµν + αW − s (Π) µν − γ (cid:18) αv p − γ αB p W ( A + B V ) ( ϕ − ) pp (cid:19) ∂ p g µν . (18)We use a shorthand notation which abbreviates the con-traction with the projected Jacobian d µ Φ iµν = ϕ ii ∂ µ Φ iµν , (19)and write B V = B p V p . The boost metric is g ij = γ ij + W v i v j . (20)In terms of the upper case sources, the lower case sourcescan be written as s ( g ) µν = S ( g ) µν γ A − B V ,s (Φ) iµν = S (Φ) iµν + W V i (cid:16) V j S (Φ) jµν − γ S ( g ) µν + S (Π) µν (cid:17) ,s (Π) µν = γ S ( g ) µν γ A − B V + W (cid:16) V j S (Φ) jµν − γ S ( g ) µν + S (Π) µν (cid:17) . (21)We take s i to be an arbitrary spatial vector of unit mag-nitude with respect to it ( g − ) ij s i s j = 1 and define aprojection operator orthogonal to s i by q ⊥ ij = γ ij − ( g − ) ik s k s j . (22)The characteristic variables of the system are given by u ˆ0 µν = g µν ,u ˆ Biµν = q ⊥ ji Φ jµν + W q ⊥ ji v j (Π µν − γ g µν ) ,u ˆ ± µν = Π µν ∓ W (cid:112) v s ) ( g − ) ij s j Φ iµν − γ g µν , (23)with the corresponding characteristic speeds β s − α v s + γ B p ( ϕ − ) pp s p A + γ B V ,β s − α v s ,β s ± α (cid:112) v s ) . (24) B. DF scalar field
For completeness, in this section, we compute the fieldequations for the scalar field project employing the DFformalism. This calculation is directly analogous to thatin the last subsection. We start with the first order formof the scalar field equations written in the upper casecoordinates ∂ T Φ = B i ∂ i Φ + AS (Φ) ,∂ T χ i = B j ∂ j χ i + A∂ i Π + γA∂ i Φ + AS ( χ ) i ,∂ T Π = B i ∂ i Π + A (N) γ ij ∂ j χ i + AS (Π) . (25)where the source terms are given by S (Φ) = Π ,S ( χ ) i = A − χ j ∂ i B j + A − Π ∂ i A − γχ i ,S (Π) = K Π + A − χ i (N) γ ij ∂ j A − (N) γ ij (3) Γ kij χ k . (26)The first order system is of the form given in Eqn. (10),therefore we can construct the principal matrices as A p = γδ pi δ pi (N) γ pj . (27)As mentioned in [5], when the Lorentz factor W isbounded, it is possible to invert the coefficient (cid:0) + A V (cid:1) ,which gives (cid:0) + A V (cid:1) − = − γW V i (N) g ji W V i − γ ( W − W V j W , (28)where (N) g ji = (N) γ ji + W V j V i . Using this informationin Eqn. (11), we arrive at the evolution equation for thescalar field variables in the lower case coordinates ∂ t Φ = ( β p − αv p ) ∂ p Φ + αW − s (Φ) ,d t χ i = (cid:16) β p δ ji − αv p δ ji + αW v i ( g − ) pj (cid:17) d p χ j + αW − g pi ( γ∂ p Φ + ∂ p Π) + αW − s ( χ ) i ,∂ t Π = β p ∂ p Π + γαv p ∂ p Φ + αW ( g − ) pi d p χ i + αW − s (Π) . (29)Here again we use a shorthand notation which abbrevi-ates contraction with the projected Jacobian d µ χ i = ϕ ii ∂ µ χ i . (30)In terms of the upper case sources the lower case sourcesbecome s (Φ) = S (Φ) ,s ( χ ) i = S ( χ ) i + W V i (cid:16) V j S ( χ ) j − γS (Φ) − S (Π) (cid:17) ,s (Π) = − γS (Φ) − W (cid:16) V j S ( χ ) j − γS (Φ) − S (Π) (cid:17) . (31)The characteristic variables of the system are given by u ˆ0 = Φ ,u ˆ Bj = q ⊥ ij χ i − W q ⊥ ij v i Π ,u ˆ ± = − Π ∓ W (cid:112) v s ) ( g − ) ij s j χ i − γ Φ , (32)with the corresponding characteristic speeds β s − α v s , β s − α v s , β s ± α (cid:112) v s ) . (33)Here s again denotes an arbitrary unit vector which isnormalized against the boost metric and is spatial withrespect to n a . III. DF JACOBIANS
In this paper, we are going to present the first numer-ical tests with DF-GHG with an aim to not only changethe spatial coordinates [27] but also the foliation. DF-GHG is implemented in 3d but for now we focus on spher-ical tests. To demonstrate that everything in the code iscorrect, we implement a list of Jacobians, some analyticand some that require the evolution of additional fields.As a sanity check, the simplest Jacobian that we imple-ment is the identity Jacobian t = T, x i = X i , (34)which of course gives the correct result that we wouldexpect in a run without DF when the same gamma pa-rameters are chosen for the job. Although these Jaco-bians are primarily built for use in spherically symmetricspacetimes, the implementation itself is made in our fully3d code. This has the twin advantages that, using theCartoon method [25, 34] for symmetry reduction, we candevelop and turn around simulations very quickly, butsimultaneously end up with code that can be used in amore general context. In the following subsections weconsider: Analytic Jacobians:
In these tests the two sets of co-ordinates are related by given closed form expres-sions. The new aspect is that in the past all sim-ulations were performed under the simplifying as-sumption T = t . Vanishing shift Jacobian:
Close to the threshold ofblack hole formation in vacuum there are indica-tions [31] that popular choices of generalized har-monic coordinates form coordinate singularities. Itis known that asymptotically flat spacetimes canalways be foliated using a vanishing shift, whichthis Jacobian choice enforces.
Areal radius Jacobian:
In spherical symmetry thereis a close relationship between the geometric radialcoordinate and the null expansion. In this Jacobian we exploit this relationship to build coordinates inwhich (as long as we excise close enough to the ap-parent horizon) the excision boundary remains out-flow for sure, and for which the apparent horizon isguaranteed to stay in the computational domain.
DF-excision Jacobian:
This Jacobian is an adjust-ment to the previous setup in which we use a so-lution to the eikonal equation to get tight controlalso over the incoming coordinate light speeds.Alternative choices will be presented in future work.
A. Analytic Jacobians
First we consider the analytic Jacobian described bythe following relations t = T, x i = f ( t, r ) X i . (35)Here we choose f ( t, r ) such that at t = 0 and for largeradius, the upper case and the lower case coordinatesmatch with each other. We do not yet have provisions forthe applying outer boundary conditions in the DF case,so at large radius we require the coordinates to changeback to GHG where the usual GHG boundary conditionsin bamps can be applied. A choice for f which satisfiesthese conditions is given by f ( t, r ) = 1 + t A e − ( r − r ) e − ( t − t ) , (36)where the Gaussian is centered such that its values ap-proximately reach machine precision or less near theouter boundary. Likewise, we consider another Jacobianwhich is described by the following relations t = f ( t, r ) T, x i = X i . (37) B. Vanishing shift Jacobian
The vanishing shift Jacobian which keeps the lowercase shift zero at all times. Such a choice of coordinatesmay be useful when performing simulations of gravita-tional collapse. First we choose t = T, (38)which makes some other quantities trivial, that is W = 1 , α = A, V i = 0 , v i = 0 . (39)With these choices, we can write down the first ofEqn. (9) as π i = − A − B i . (40)This also simplifies the evolution equation for φ ii inEqn. (7) which can be obtained using Cartan’s magicformula [3] ∂ t φ ii = − D i B i + L β φ ii . (41)The first term in the right hand side of the above equa-tion can be considered a source term, because by addi-tion of the reduction constraints, given in Eqn. (15), allfirst derivatives of metric components can be replacedby evolved variables, whereas the second term should beideally zero since we want the lower case lapse to be zero.However, since we do not yet have outer boundary condi-tions in the lower case coordinates, we will employ a tran-sition function approach. In this approach, we choose β i = Ω( r ) B i , (42)where Ω is zero at small radii and transitions to one atlarge radii. This allows us to apply the standard GHGboundary conditions for the outer boundary.The first source term in Eqn. (41) can be written as ∂ i B i = J ki ∂ k B i . (43)The upper case spatial derivatives of the upper case shiftcan then be written down in terms of the lapse, shift,extrinsic curvature and the Christoffel symbols, the ex-pressions for which are given below [35] ∂ m B l = Γ lm + B l Γ m + AK lm − B n Γ lmn − B l B n Γ mn , (44)where further we use the expressions for the extrinsiccurvature K ij = − A Γ ij . (45)The Lie derivative term in Eqn. (41) is only non-vanishing for the subpatch where the transition happensand for all outer subpatches. It can be written down as L β φ ii = Ω B k ∂ k φ ii + φ ij ∂ i β j , (46)where ∂ i β j = ( ∂ r Ω)Θ i B j + Ω( ∂ i B j ) . (47)where Θ i are functions of the angular coordinates whichare same for both the upper case and lower case coordi-nates and are related to the Cartesian coordinates by therelation x j = r Θ j , X j = R Θ j , Θ i = δ ii Θ i . (48)Putting all of this together, we can construct the requiredJacobian from which the inverse Jacobian can be com-puted numerically J = (cid:32) − B i + Ω B j δ jj φ ij φ ii (cid:33) . (49) C. Areal radius Jacobian
In this setup, we choose the Jacobian such that thelower case radial coordinate is the areal radius. Withthis choice, we can show that the position of the appar-ent horizon is located at the zero crossing of the out-going radial coordinate lightspeed. Now, since the po-sition of the apparent horizon in spherical symmetry inthese coordinates can only increase as the simulation pro-gresses [28, 29], if the apparent horizon appears on thegrid at the beginning of the simulation, it must do so atlater times. Consequently, as a result of the weak cos-mic censorship conjecture, the event horizon stays on thenumerical domain at all times. This ensures a successful‘excision’ strategy.Consider the upper case foliation whose spatial lineelement is given by ds = L dR + (N) γ T R d Ω . (50)Here L is called the length scalar, which is to the 2 + 1split what the lapse is to the 3 + 1 case. The relation-ship between the upper case and the lower case radialcoordinate is given by r = η ( r, (N) γ T ) R, (51)where η is a function chosen such that for small r thelower case radial coordinate becomes the areal radius co-ordinate whereas for large values of r , it becomes thestandard radial coordinate as given by GHG. This en-sures that normal GHG outer boundary conditions canbe applied for the system. A possible choice of η is of theform η = (cid:112) (N) γ T χ ( r ) + 1(1 − χ ( r )) , (52)where χ ( r ) is any suitable transition function varyingfrom zero to one with increasing r . In principle, a hy-perbolic tangent function would serve the purpose butfor reasons of rapid convergence, we choose a low orderpolynomial function which transitions at the penultimatesubpatch. The functional form of χ ( r ) which transitionsfrom one to zero between r = r and r = r can be givenby χ ( r ) = , r < r , − a ( r − r ) − a ( r − r ) , r ≤ r ≤ r , , r > r . (53)where a = − / ( r − r ). The derivatives of η are givenby ˜ ∂ r η = (cid:112) (N) γ T χ (cid:48) ( r ) − χ (cid:48) ( r ) , ˜ ∂ ( (N) γ T ) η = χ ( r )2 √ (N) γ T . (54)Here the tilde on the partial derivatives means that thederivative must be taken keeping the other argument con-stant. We shall now construct the various components ofthe inverse Jacobian by noting that in this case α = A, V i = 0 , W = 1 . (55)This information can be used to construct the ( J − ) i components of the inverse Jacobian. The spatial compo-nents of the inverse Jacobian as given in Eqn. (8) follow-ing the relation given belowΦ ji = ∂ i x j = Θ j ∂ i r + r∂ i Θ j , (56)using the relationship between r and Θ j as specified inEqn. (48). The two terms of the above equation can beevaluated using the fact that ∂ i r = Θ i η + ( r ˜ ∂ ( (N) γ T ) η∂ i (N) γ T ) /ηκ , (57)where κ ( r, (N) γ T ) ≡ − r ˜ ∂ r ηη , (58)where we use Eqn. (51) and the fact that ∂ i R = Θ i . Forthe second term, we have ∂ i X j = δ ji = ( ∂ i R )Θ j + R∂ i Θ j , (59)which gives ∂ i Θ j = ηr (cid:16) δ ji − Θ i Θ j (cid:17) . (60)We will now calculate the ( J − ) i component which canthen be used to construct the time-space part of the in-verse Jacobian in Eqn. (8)( J − ) i = ∂ T x i = Θ i ∂ T r, Π j = Θ j ∂ T r − B i Φ ji A . (61)Now, the upper case time derivative of r can be computedfrom Eqn. (51) in a straightforward manner ∂ T r = ( r/η ) ˜ ∂ ( (N) γ T ) η∂ T (N) γ T κ . (62)For the sake of completeness, we also provide the up-per case time and spatial derivatives of (N) γ T which areneeded to construct the above quantities. An expressionfor (N) γ T can be written down in terms of the lapse, thedeterminant of the Cartesian form of the metric and thelength scalar as (cid:112) − g sph = AL √ q, (63)where g sph is the determinant of the metric in sphericalcoordinates and q is the determinant of the two metricgiven by (cid:18) R (N) γ T R sin θ (N) γ T (cid:19) . (64) Using the fact that (cid:112) − g sph = R sin θ √− g cart , (65)we obtain an expression for (N) γ T where the reference tothe Cartesian form of the metric is suppressed for thesake of brevity (N) γ T = √− gAL . (66)From here, the upper case time and spatial derivativesof (N) γ T can be obtained in a straightforward manner byusing the derivatives of √− g , A and L . Using standardresults from the literature [35], we can compute time andspatial derivatives of the square root of the determinantof the metric ∂ T √− g = √− g Γ µ µ , ∂ i √− g = √− g Γ µiµ , (67)and also for the lapse ∂ T A = A (Γ − B m Γ m ) , ∂ i A = A (Γ i − B m Γ im ) . (68)The derivatives of the upper case length scalar L can becomputed from its definition L − = (N) γ ij Θ i Θ j . (69)To see that the apparent horizon is located at the zerocrossing of the outgoing radial coordinate lightspeed inarea locking coordinates, we consider the expression forthe expansion which can be written as [35] H = 1 L (cid:18) R + 1 (N) γ T ∂ R (N) γ T (cid:19) − (N) K θθ , (70)where (N) K ij is the extrinsic curvature in the upper casefoliation. A similar expression of course holds in an ar-bitrary foliation. We have H ∝ (cid:0) ∂ T + C R + ∂ R (cid:1) R (N) γ T , (71)where C R + is called the outgoing radial coordinate light-speed and is defined as C + = − B R + A/L where hereand in the following we suppress the label R . Now intro-ducing area locking coordinates (˚ T , ˚ R = R √ (N) γ T ), wecan write H ∝ (cid:16) ∂ ˚ T + c ˚ R + ∂ ˚ R (cid:17) ˚ R , ∝ c ˚ R + ˚ R. (72)From the above expression, we see that in the case of theapparent horizon, where the expansion is zero, c ˚ R + = 0as ˚ R is greater than zero. D. Dual frame excision Jacobian
As an addition to the previous strategy, which ensuresthe correct sign of the outgoing radial coordinate light-speed c + at the inner boundary of the simulation pro-vided that we excise close enough to the apparent hori-zon, we would like to exactly control the incoming radialcoordinate lightspeed c − , at least near the black hole.If c − can be set to − T, X i ) are the generalized harmonic coordinateswhereas the lower case coordinates ( t, x i ) are defined bythe Jacobian to be described shortly. The angular coor-dinates are kept to be the same in both cases. The re-lationship between the upper case and lower case radialcoordinate is kept same as the previous strategy. Fur-thermore, a new coordinate ˚ v is introduced r = η ( r, (N) γ T ) R, ˚ v = ˚ T + r, (73)where ˚ T = √ (N) γ T T . The only condition that we imposeon ˚ v is that it be a null coordinate, that is, it satisfies theeikonal equation g ab ∇ a ˚ v ∇ b ˚ v = 0 , (74)The eikonal equation inside the transition regionwhere ˚ T = t can be expanded using the expression forthe coordinate lightspeeds along the radial direction c ± = − β r ± α/l, (75)keeping in mind that l − = γ rr :1 α (1 + c + )(1 + c − ) = 0 . (76)From the above expression, it can be clearly seen thatwhen c + is not equal to − c − takes the value of − t = η ( r, (N) γ T ) T, = ˚ T χ ( r ) + 1(1 − χ ( r )) T, (77)Using Eqn. (73), we arrive at the final relation betweenthe upper case and the lower case time coordinate t = ˚ vχ − rχ + (1 − χ ) T. (78)It is clear from the above expression that we have toevolve derivatives of ˚ v and T to obtain the different com-ponents of the inverse Jacobian. Instead of the compo-nents of the Jacobian, we can choose to evolve an equiv-alent set of quantities which are known as the ‘opticalJacobians’ V − i ≡ − ∂ i ˚ v, E − ≡ N µ ∂ µ ˚ v. (79) In terms of these new variables, it is straightforward toshow that the eikonal equation in Eqn. (74) can be rewrit-ten as E − = (N) γ ij V − i V − j . (80)It is clear from the above expressions that the evolu-tion equation for E − can be completely dropped in fa-vor of V − i . However, we choose to keep them since V − j and E − satisfy the eikonal equation, which can be usedto construct a constraint monitor.We ask the reader to refer to [5] for a complete deriva-tion for the equations of motion for the optical Jacobiansand only provide a brief summary of the final equationshere. The evolution equations in the upper case coordi-nates can be written as a set of advection equations suchthat this subsystem is minimally coupled to the first or-der GHG system ∂ T V − i = ( B j − AS j − ) ∂ j V − i + AS ( V − ) i ,∂ T ln E − = ( B j − AS j − ) ∂ j ln E − + AS ( E − ) , (81)where S j − = E − − V j − and the source terms are given by S ( V − ) i = A − V − j ∂ i B j + S j (3) Γ kij V − k − A − E − ∂ i A,S ( E − ) = K S − S − − L S − ln A. (82)Since our objective is to evolve the optical Jacobians inthe lower case ‘Cartesian’ coordinates, we must transformthe evolution equations in Eqn. (81) using Eqns. (10)and (11). To do this, we compute ( ϕ − ) ii which can bewritten down in terms of Φ ii , Π i and V i , (cid:0) ϕ − (cid:1) ii = Φ ii + Π i V i . (83)The complete principal matrix A p associated with theGHG variables and our new variables can be expressedas A p = (cid:18) Λ
00 Λ (cid:19) , (84)where Λ = γ A − B p γ δ pi − δ pi γ γ A − B p − (N) γ pj , (85)and Λ = (cid:18) − S p − − S p − (cid:19) . (86)Using the above expressions and Eqn. (11), the evolutionequations in the lower case coordinates can be written as ∂ t V − i = ( β j − αs j − ) ∂ j V − i + αW − s ( V − ) i ,∂ t ln E − = ( β j − αs j − ) ∂ j ln E − + αW − s ( E − ) , (87) r . . . . . L ( r ) Data with scalar field Pure Schwarzschild case .
98 1 .
99 2 .
00 2 .
01 2 .
02 2 .
03 2 .
04 2 . r t i m e AHloc output EHloc output . . . . . M A D M FIG. 2. Left: An example of constrained solved initial data with the blue line representing data corresponding to a non-zeroscalar field and the green line representing the pure Schwarzschild case. The inset plot shows the values of M ADM which aregenerated during the iterative solve. Right: A comparison between the output of the event horizon locator and apparent horizonlocator for a simulation with a lapse perturbation in the Schwarzschild spacetime. The deviation in the two outputs at latertimes demonstrates the event horizon locator trying to ‘find’ the horizon. Here this effect is exaggerated because we chose apoor initial guess for the position of the event horizon on purpose. where s i − = (cid:0) ϕ − (cid:1) ii S i − W (cid:16) W v j ( ϕ − ) ji S i − (cid:17) + v i , (88)and the lower case source terms are related to the uppercase sources as s ( V − ) i = (cid:16) W v j (cid:0) ϕ − (cid:1) jj S j − (cid:17) − S ( V − ) i ,s ( E − ) = (cid:16) W v j (cid:0) ϕ − (cid:1) jj S j − (cid:17) − S ( E − ) . (89)It is straightforward to obtain the equations of motionfor the other two variables ∂ t R = J , ∂ t ˚ v = αW E − − b i W V i − b i ϕ ji V − j , (90)where b i = − αv i + β i . Finally, we can use the aboveinformation to construct the different components of theoptical Jacobian. The Φ ji component are evaluated inthe same way as the previous case. Now( J − ) = ∂ T t, = (cid:16) AE − − B j V − j (cid:17) χ + ˚ v∂ r χ∂ T r − χ∂ T r − r∂ r χ∂ T r − T ∂ r χ∂ T r + (1 − χ ) , ( J − ) i = ∂ T x i = Θ i ∂ T r, ( J − ) i = ∂ i t = − V − i χ + ˚ v∂ i χ − ( ∂ i r ) χ − r∂ i χ − T ∂ i χ, (91)where ∂ i χ = ( ∂ r χ )Θ i Φ ii . (92) The ∂ T r term in these equations can be simplified usingEqn. (62). To initialize the evolved quantities at the be-ginning, we propose a choice which leads to the Jacobianbeing identity initially. A choice of the evolved quantitiesis V − i = − ∂ i r, ˚ v = r, T = 0 . (93)The choice of E − is not independent but follows from theeikonal equation.Another important point to note is that we employ thecartoon method to compute the y and z derivatives usingKilling vectors. The formula for doing this is providedbelow ∂ y V − i = h ( x ) (cid:16) δ xi δ jy − δ yi δ jx (cid:17) V − j ,∂ z V − i = − h ( x ) (cid:16) δ xi δ jz − δ zi δ jx (cid:17) V − j . (94)Note that h ( x ) = 1 for the on-axis case and h ( x ) = 1 /x otherwise.Lastly, we briefly describe the constraint preserv-ing outer boundary conditions, the constraint beingEqn. (80). At the outer boundary, we choose ˙ V − i tobe equal to zero. This requires a choice of ˙ E − which isgiven by ˙ E − = − V − i V − j E − (N) γ ik (N) γ jl ∂ T (N) γ kl . (95) IV. CODE SETUP
In this section we describe our numerical setup, initialdata and post-processing tools.0
A. Code overview
The bamps code [12, 31, 36, 37] is built for large scale,parallel numerical evolutions of hyperbolic systems. Sev-eral different approximation schemes are implemented,including DG schemes [36], but here we use exclusivelya multidomain pseudospectral method to solve our firstorder symmetric hyperbolic PDEs described in the previ-ous sections. Each individual numerical domain is calleda subpatch. Within each subpatch spatial derivativesare approximated using Chebyshev polynomials imple-mented, as usual, by matrix multiplication. Data arecommunicated between patches using a penalty methodwhich is applied to the incoming characteristic variablesat each subpatch boundary. Our domain always has asmooth timelike outer boundary at a fixed radial coordi-nate r . Because of this we need to apply boundary condi-tions. These need to be constraint preserving, to controlundesirable gauge effects, and to control the physical be-havior at the boundary [38]. For now, to avoid introduc-ing too many new complications into the code at once,we choose Jacobians that transition to the identity in aneighborhood of the outer boundary. This allows us torecycle our boundary conditions for the GHG formula-tion, essentially those of Rinne [39], directly. For evolu-tion in time we use a fourth order Runge-Kutta method.Because we will be treating spherical spacetimes we usethe Cartoon method [25, 34] to suppress two spatial di-mensions. With this reduction our tests are very fast,the longest taking just a few minutes on a large desktopmachine. We have tested the implementation by evolv-ing our spherical data with the full 3d setup and obtainperfectly consistent results, and so do not discuss theseslower computations further. For a deeper technical de-scription of the code we direct the reader to [31]. B. Initial data
Our system involves a scalar field minimally coupled tothe metric. To evolve such a system, we must first solvefor constraint preserving initial data which can then beevolved using a combination of DF-GHG and DF scalarfield projects. We shall provide the necessary coupledordinary differential equations for the sake of complete-ness. Consider the coordinates ( r, θ, φ ) in which the lineelement of the spatial metric can be written as ds = l ( r ) dr + r d Ω , (96)where l is the lower case length scalar. The form of theextrinsic curvature follows in a straightforward manner K ij = K rr ( r ) 0 00 r K T ( r ) 00 0 r sin θK T ( r ) . (97)We can now use this to obtain the Hamiltonian and themomentum constraints which are given below respec- tively4 K T K − K T + 2 (cid:0) rl (cid:48) + l − l (cid:1) r l = 8 π (cid:18) Φ (cid:48) l + Π (cid:19) , rK (cid:48) T + 3 K T − K ) r = 8 π ΠΦ (cid:48) , (98)where K is the trace of the extrinsic curvature, Φ is thescalar field and Π is related to the time derivative of thescalar field as Π = − α (cid:0) ∂ t Φ − β i ∂ i Φ (cid:1) . (99)From the Hamiltonian and momentum constraints we ar-rive at the ODEs that we can solve using a Runge-Kuttamethod dldr = l r (cid:0) − r K T l K + 3 r K T l + 4 πr l Π − l + 4 πr Φ (cid:48) + 1 (cid:1) ,dK T dr = − K T + 4 πr ΠΦ (cid:48) + Kr . (100)This ODE is solved using an iterative method, withthe trace of the extrinsic curvature taken to be that inSchwarzschild with the M ADM mass taken to be one. Thevalues of l ( r ) obtained in the first iteration is then usedto construct the new ADM mass, defined by M ADM = 12 r (cid:0) l ( r ) − (cid:1) , r → ∞ (101)This is continued until the difference between the lastand the second last evaluation of the ADM mass meetsa tolerance level. The spatial metric quantities can thenbe reconstructed using Eqn. (96), while the coordinatelightspeed C + constructed from C + = r − M ADM r + 2 M ADM , (102)can be used to reconstruct the lapse and the shift α = l + C + l , β r = 1 − C + . (103)An example of the initial data solver in action is shownin the left plot of Fig. 2. C. Apparent and event horizon finders
We require diagnostic tools for post-processing to en-sure that the excised region of spacetime remains insidethe black hole event horizon at all times during the nu-merical evolution. For this purpose we use two tools, theapparent horizon, defined locally on a given hypersurfaceand the event horizon which is a global property of thespacetime.1The apparent horizon is defined as the outermostmarginally outer trapped surface on a given spatial hy-persurface, that is, it is defined by the vanishing of theexpansion parameter of the outgoing null geodesics. Inspherical symmetry, the condition for the apparent hori-zon is given by [35] H = 1 l (cid:18) r + 1 γ T ∂ r γ T (cid:19) − K θθ = 0 , (104)where the metric is represented in the lower case basis.This is implemented in the AHloc feature of bamps . Analternative and simpler way to find the apparent horizonin the area locking coordinates is that the zero crossingof the lower case outgoing coordinate lightspeed c + cor-responds to the position of the apparent horizon.We will now describe the implementation of a newevent horizon finder EHloc for bamps . The event horizonin general is a 2 + 1 null surface which is the boundary ofthe black hole region from which no future pointing nullgeodesics can escape to null infinity J + [40]. Hence, oneway to obtain approximations of event horizons in nu-merical spacetimes is to integrate null geodesics forwardin time all over the numerical domain. The disadvantageof this method is that simulations can only be performedfor a finite time and hence it is not straightforward to findescaping null geodesics [41]. A more efficient algorithmis to integrate outgoing null geodesics or null surfacesbackwards in time, since then the event horizon acts asan attractor of null geodesics [41, 42].The geodesic method for integrating backwards in timeis considered to be the most accurate method and prob-lems mentioned in the literature like tangential driftingare not seen in practice [43]. Hence, this is the methodwe have used for EHloc .Unlike
AHloc , it is essential for
EHloc to be run in post-processing when the black hole is no longer ringing butis rather Schwarzschild. In such cases, we can start fromthe last numerical slice integrate the geodesic equationbackwards d x α dλ + Γ αβγ dx β dλ dx γ dλ = 0 , (105)where λ is the affine parameter and x α is the 4-position ofthe geodesic. The initial conditions for the geodesic areso chosen that it is outgoing. In spherical symmetry, thegeodesic equation can be represented by a set of coupledordinary differential equations [44] d Π r dt = − α ,r + ( α ,r Π r − αK rr Π r Π r )Π r + β r,r Π r − αγ rr,r Π r Π r ,drdt = α Π r − β r . (106)Here α , β i are the lapse and shift respectively, K ij is theextrinsic curvature and γ ij is the inverse of the spatialmetric. All quantities mentioned here are represented in the lower case basis. The intermediate variable Π r isrelated to the momentum p r = dx r /dλ asΠ r ≡ p r (cid:112) γ ij p i p j . (107)A Runge-Kutta integrator is used for performing the timestepping while the data is loaded and then interpolatedusing Chebyshev functions. A brief description of thegrid and interpolation setup is given below.The data is written on every patch at the Gauss-Lobatto points x α = − cos (cid:18) πβN − (cid:19) , (108)where N is the number of points on each grid and β =0 , . . . , N −
1. Chebyshev polynomials are used to performthe spectral interpolation T n ( x ) = cos( n cos − x ) . (109)These polynomials are defined in the interval [ a, b ] by achange of variable: y ≡ x − ( b + a ) ( b − a ) . (110)The coefficients of interpolation a , . . . , a n − are foundout by solving: T ( x ) . . . T N − ( x )... . . . T ( x n − ) T N − ( x N − ) a ... a N − = u ... u N − , (111)where u , . . . , u N − are the given values at the N points.The spatial derivatives are computed by a matrix multi-plication [13]: ( ∂ x u ) α = N − (cid:88) k =0 D αk u k , (112)where D αβ is the Gauss-Lobatto derivative matrix givenby D αβ = − N − +16 , α = β = 0 , q α ( − α + β q β ( x α − x β ) , α (cid:54) = β, − x β − x β ) , α = β = 1 , . . . , N − , N − +16 , α = β = N − , (113)where q α = 2 at the boundary points and q α = 1 else-where.In practice, we do not compute the diagonal terms ofthe derivative matrix but use the identity which gives thederivative matrix better stability against rounding errors D αα = − N − (cid:88) k =0 ,k (cid:54) = α D αk . (114)2The time interpolation on the data is performed usinga linear interpolation algorithm and a judicious choiceof the number of points needs to be taken into account.The number of data steps loaded into memory also affectsperformance. However, both of these problems are hard-ware specific and hence we do not go into details here.As sanity checks, we test the event horizon finder witha multi-patch simulation of the Schwarzschild spacetimeand another with a gauge perturbation. These resultshave also been compared with the output of the appar-ent horizon finder and seen to be in good agreement. Thegauge perturbation case with both EHloc and
AHloc out-puts are shown in the right plot of Fig. 2.
V. NUMERICAL RESULTSA. Tests with analytic Jacobians
We begin our numerical experiments by first testingout our implementation of the analytic Jacobians. As ini-tial data, we choose the metric components to be those ofthe Schwarzschild spacetime in Kerr-Schild coordinates.We perform simulations for both the analytic Jacobianskeeping the function f to be f = 1 + 0 . t e − ( r − e − ( t − . (115)The numerical domain for these simulations are from r ∈ [1 . , .
8] and they are performed at 4 different resolu-tions starting from 20 patches, 11 points and increasingthe number of points by 10 in each case. We finally plotthe harmonic constraints C α = H α + g βγ Γ αβγ , (116)with radius at four different resolutions and find thatthe constraints converge with increasing resolution. Aplot of such a convergence test is provided in Fig. 3. Wealso perform tests of the time derivatives of the harmonicconstraints given by F α (cid:39) ∂ N H α + g βγ ∂ N Γ αβγ − Γ βγα ∂ N g βγ , (117)where (cid:39) denotes equality up to the combinations of thereduction constraints and ∂ N ≡ N α ∂ α (see [13, 32] fordetails). We find that they also converge with increasingresolution. B. Tests with the vanishing shift Jacobian
Another numerical experiment we perform is to keepthe lower case shift to be zero by a suitable choice of Ja-cobian. This experiment is performed on the Minkowskispacetime by adding a Gaussian gauge wave of the form α = 1 + G (cid:48) e − w (cid:48) ( r − r ) , (118) where the parameters of the perturbation are given by G (cid:48) = 0 . , r = 0 , w (cid:48) = 1 . (119)As has been demonstrated in the calculations of sec-tion III, the upper case and the lower case coordinatesmatch in the outermost subpatch of the simulation whilethe penultimate subpatch serves as the transition region.This can be seen clearly in the plots in Fig. 4 where theupper case shift is represented by the blue curve and thelower case shift is represented by the green curve. Thelower case shift is successfully kept to zero in patchesinside the transition zone, while outside the transitionzone, it is seen to agree with the upper case shift. Aconvergence test is also performed considering the reduc-tion constraints, the harmonic constraints and the timederivatives of the harmonic constraints and we see con-vergence with increase in resolution, as is expected. C. Tests with the areal radius Jacobian
We now perform simulations of a massless scalar fieldminimally coupled to general relativity in spherical sym-metry. As a first try, we evolve the spacetime in gener-alized harmonic coordinates [32] using the old excisionsetup, that is, there are no boundary conditions placedat the inner boundary which is expected to be outflowat the beginning of the simulation. We also monitor thesigns of the two radial coordinate lightspeeds as given inEqn. (75) at the inner boundary of the simulation at alltimes. We set the scalar field data initially to be of thefollowing profileΦ =
Cr e − ( r − r ) /σ , Π = − Crσ ( r − r ) e − ( r − r ) /σ . (120)The specific parameters which are chosen for the run are C = 0 . , r = 11 . , σ = 1 . (121)Using these parameters, we then solve for the spatial met-ric, extrinsic curvature, lapse and shift in the initial datausing the method prescribed in section IV B.A plot of the outgoing radial coordinate lightspeed C + ,as can be seen in the bottom row, left of Fig. 5 shows thatit assumes a positive sign at the inner boundary for sometime during the simulation. This indicates that the ex-cision strategy has failed as the inner boundary has notremained an outflow boundary during those times. Thisexperiment clearly demonstrates that for certain configu-rations of the matter content, the existing excision strat-egy is unsuccessful.We now perform the same experiment, but this time weswitch on DF and the areal radius Jacobian as describedin section III. As described before, the use of the arealradius ensures that the position of the apparent horizoncan only monotonically increase with time, which ensures3 r − − − − − H a r m o n i cc o n s t r a i n t ( C T )
11 21 31 41 r . . . . . . . L a p s e lowercase lapse uppercase lapse
90 100 1100 . . . . FIG. 3. Left: A convergence test performed with the time component of the harmonic constraints when the analytic Jacobiangiven in Eqn. (37) is considered. The numbers in the legend correspond to the number of points per patch considered. Right:A comparison between the lower case and the upper case lapse at time t (cid:39) r − . − . . . . s h i f t r − . − . . . . uppercase shift lowercase shift r − . − . − . − . . . . . FIG. 4. A snapshot of the upper case shift and the lower case shift at three different times given by t (cid:39) .
7, 13 . . M ). In the first figure, we see as expected irrespective of the upper case shift, the lower case shift is vanishingly smallinside the transition region. In the second figure, we demonstrate the effect of the transition region on the lower case shift. Inthe third figure, we see the form of the lower case shift mostly outside the transition region where it is expected to agree withits upper case counterpart. that if it is initially located inside the numerical domain,it shall do so at all times. Like before, the lightspeedsat the inner boundary are monitored at all times. It canbe clearly seen, from the green line in the bottom left ofFig. 5 that the value of C + at the inner boundary remainsnegative at all times during the simulation thereby indi-cating that the new excision strategy is successful. Al-though not seen in the plot, the C − lightspeed, while itsvalue fluctuates, remains negative at all times for boththe DF and the non-DF case. This fact can be seen fromEqn. (75) which shows that if C + remains negative at alltimes, so must the value of C − . We also perform simula-tions with different quantities of scalar field content and find many other cases where the new method proves tobe successful where the old one does not.With the parameters which have been provided above,we perform simulations at three different resolutions,having 29, 37 and 45 points per patch and plot the re-duction constraints which are defined by (15) and can becomputed in the lower case coordinates from C iαβ = (cid:0) ϕ − (cid:1) ki ∂ k g αβ + V i Π αβ − Φ iαβ , (122)as a function of space for a given value of time. We seeconvergence with increase in resolution, as can be seen inthe top left plot of Fig. 5.4 r − − − − R e du c t i o n c o n s t r a i n t ( C X TT )
29 37 45 time . . . . . . . . E H / A H p o s i t i o n AH position EH position time − . − . − . − . − . . . c + without dfex with dfex r t AH position EH position − . . . . . . . FIG. 5. Top row, left: A convergence test performed with the
XT T component of the reduction constraints at t = 30 M performed at three different resolutions, the number of points per patch being mentioned in the legend. These simulations areperformed with the GHG, DF and scalar field projects. Top row, right: The position of the apparent horizon and the eventhorizon as a function of time in one of these simulations with DF, GHG and scalar field. The position of the apparent horizon,also where c + crosses zero is seen to increase monotonically as a function of time. Bottom row, left: A plot showing the valueof the outgoing radial coordinate lightspeed c + at the inner boundary of two simulations, one performed with DF-excisionand another performed without it. The simulation performed with DF-excision switched on has a negative value of c + at theinner boundary throughout the simulation while the non-DF simulation fails in maintaining that. Bottom row, right: The timeevolution of the scalar field along with the position of the event horizon and the apparent horizon. We also employ our event horizon and apparent hori-zon finders to track the location of the horizons. A su-perposition of the output of the two finders is provided atthe top right of Fig. 5. As expected the apparent horizongrows monotonically with time from 2 M to (cid:39) . M .The event horizon also shows a monotonic behavior inthese coordinates.Another experiment we perform involves evolving theSchwarzschild spacetime with a lapse perturbation α = α (cid:48) + H (cid:48) e − w (cid:48) ( r − r ) ,∂ i α = ∂ i α (cid:48) − H (cid:48) w (cid:48) ( r − r ) e − w (cid:48) ( r − r ) x i r , (123)where α (cid:48) is the natural lapse associated with the Schwarzschild metric in Kerr-Schild coordinates. Thespecific parameters which are chosen for this experimentare H (cid:48) = 1 , r = 10 , w (cid:48) = 1 . (124)The lapse perturbation is shown in the inset plot of Fig. 6.We observe the zero crossing of the outgoing radial co-ordinate lightspeed C + as this corresponds to the posi-tion of the apparent horizon. At the beginning of thesimulation, this stays at 2 M and continues to remain sothroughout the entire duration of the simulation. Thisis shown by the blue and green lines in the left plot ofFig. 6 which correspond to the lightspeed at the begin-ning and at the end of the simulation. This is indeed5the desired behavior since the position of the apparenthorizon should not change as there is no physical pertur-bation. D. Tests with the DF-excision Jacobian
Finally, we perform numerical experiments with thedual foliation eikonal Jacobian. The numerical setupagain consists of a massless scalar field minimally cou-pled to general relativity. Control of the C + radial coor-dinate lightspeed is borrowed from the treatment in theprevious section. In this section, our goal is also to con-trol the ingoing radial coordinate lightspeed C − to beidentically − C = 0 . , r = 15 , σ = 1 . (125)Our principal objective in this experiment is to grow theapparent horizon as much as possible by letting accret-ing scalar field fall into the black hole horizon. For thisspecific choice of parameters, we see that the apparenthorizon position grows from (cid:39) M to 4 . M thereby reg-istering (cid:39) C − is seen to vary freely while our methodensures that the lower case c − is strictly kept to be equalto − VI. CONCLUSIONS
In this paper, we have presented the first implementa-tion of the DF-GHG formulation, together with the DF-scalar field. The implementation was made within the bamps code. We performed a battery of tests involvingseveral Jacobians, but with an emphasis on black holeexcision. Although the tests performed are in sphericalsymmetry as proof of concept and also for reasons of effi-ciency, the whole implementation itself was made in the full 3 + 1 setting. In addition to this, we introduced ourevent horizon finding code
EHloc .To test the newly written DF-GHG project, we haveperformed elementary tests with two analytic Jacobians,in one of which we consider two different foliations forthe upper case and lower case coordinates. After this,we tested the vanishing shift Jacobian, another importantcase in which the lower case shift is kept zero at all times,which we expect to be helpful while considering cases ofgravitational collapse and black hole formation.Finally, we considered the two most important Jaco-bians for our black hole excision work, the areal radiusJacobian and the DF-excision Jacobian. In the areallocking case, we saw that when the lower case coordi-nates are made to include the areal radius, the apparenthorizon is located at the zero-crossing of the outgoingradial coordinate lightspeed c + . By basic results for dy-namical horizons, these coordinates also have the specialproperty that the position of the apparent horizon cannotdecrease. Thus if the apparent horizon is initially locatedon the numerical grid, it stays so throughout the simula-tion. Assuming the weak cosmic censorship conjecture,we then ensure that the event horizon, being located out-side the apparent horizon, also remains on the numericalgrid. In the case of the DF-excision Jacobian, we care-fully control the ingoing radial coordinate lightspeed c − to be equal to − c − to be aconstant everywhere on the numerical grid, barring thetransition region and outside, ensures that the redshiftor blueshift that potentially arises out of using ‘artificial’coordinates is avoided. We performed a series of tests onthe Minkowski spacetime with lapse perturbations or per-turbed Schwarzschild spacetimes by employing a combi-nation of the DF, DF-GHG and DF-scalarfield projects.These tests were validated by performing several con-vergence tests, which demonstrate clean spectral conver-gence.The excision setup presented here is of course highly specialized when compared with the full control systemapproach used in the SpEC code. That said it pro-vides precisely the functionality needed for our near-termwork, and has the advantage that we use only pointwise,rather than quasilocal manipulation of our variables inconstructing the Jacobians. Therefore, at the continuumlevel, basic theorems can be trivially applied to our for-mulation. To avoid coupling through derivatives betweenthe Jacobian and evolution equations, which would re-quire a more careful mathematical analysis, it was cru-cial that we could replace first derivatives of the metricusing the reduction constraints. This works because theexpansion contains at most one derivative of the met-ric. Since this fact remains true even in the absence ofspherical symmetry we hope, eventually, to generalize theDF-excision strategy to the full 3 + 1 setting. For now itis unclear whether or not this will pan out, since there isa qualitative difference between 1d and 3d excision thatcannot be overlooked. But if successful the generalization6 . . . . . . r . . . . . c + t = 0 M t = 15 M r − − − − − − − R e du c t i o n c o n s t r a i n t ( C X TT )
29 33 37 41 45 r . . . L a p s e FIG. 6. Left: A demonstration of the fact that the position of the apparent horizon does not change when hit by a lapseperturbation. In the figure, the outgoing radial coordinate lightspeed c + is shown at two different times and it is seen that thezero crossing remains at 2 M throughout the simulation. The inset plot shows the lapse perturbation in the initial data. Right:A convergence plot of the XT T component of the reduction constraints performed with five different resolutions of the samelapse perturbation simulation. The legend shows the number of points per patch. time . . . . . A H p o s i t i o n r − . − . − . − . − . − . − . − . C − / c − uppercase C − lowercase c − r − − − H a r m o n i cc o n s t r a i n t ( C T )
11 21 31 41
80 85 90 − . − . FIG. 7. Left: A demonstration showing the growth of the apparent horizon by (cid:39) C − /c − in both the upper case and thelower case coordinates with the lower case result shown to be − would provide an improved moving excision strategy forbinary black holes within a pseudospectral code.In order to perform our numerical tests, appropriateboundary conditions at the outer boundary must be pro-vided. At present, we do not yet have outer boundaryconditions in the code for the DF projects. To over-come this problem at the outer boundary, we ensure thatthe Jacobian transitions into the identity Jacobian atthe outermost subpatch. This is achieved by using a low order polynomial transition function, which worksremarkably well in practice when the transition is placedat the two ends of the penultimate subpatch, and we ex-pect that various alternative configurations would behavesimilarly. A desirable alternative would be to implementouter boundary conditions in the code that take care ofthe full DF infrastructure, including the management oftwo time coordinates. Work on this will be reported on inthe near future. An immediate goal is to use the methods7developed here to study systematically, in the sphericalcontext, the transition from the linear regime we studiedin [30] to the case with arbitrary non-linear perturba-tions. ACKNOWLEDGMENTS
We are grateful to Thanasis Giannakopoulos and Is-abel Su´arez Fern´andez and for helpful discussions andfeedback on the manuscript. MKB and KRN acknowl-edges support from the Ministry of Human ResourceDevelopment (MHRD), India, IISER Kolkata and theCenter of Excellence in Space Sciences (CESSI), In- dia, the Newton-Bhaba partnership between LIGO In-dia and the University of Southampton, the NavajbaiRatan Tata Trust grant and the Visitors’ Programme atthe Inter-University Centre for Astronomy and Astro-physics (IUCAA), Pune. CESSI, a multi-institutionalCenter of Excellence established at IISER Kolkata isfunded by the MHRD under the Frontier Areas ofScience and Technology (FAST) scheme. DH grate-fully acknowledges support offered by IUCAA, Pune,where part of this work was completed. The workwas partially supported by the FCT (Portugal) IF Pro-gram IF/00577/2015, Project No. UIDB/00099/2020 andPTDC/MAT-APL/30043/2017. [1] Mark A. Scheel, Harald P. Pfeiffer, Lee Lindblom,Lawrence E. Kidder, Oliver Rinne, and Saul A. Teukol-sky, “Solving Einstein’s equations with dual coordinateframes,” Phys. Rev. D , 104006 (2006), gr-qc/0607056.[2] SpEC - Spectral Einstein Code, .[3] David Hilditch, “Dual Foliation Formulations of Gen-eral Relativity,” arXiv e-prints , arXiv:1509.02071 (2015),arXiv:1509.02071 [gr-qc].[4] David Hilditch and Milton Ruiz, “The initial boundaryvalue problem for free-evolution formulations of Gen-eral Relativity,” Class. Quant. Grav. , 015006 (2018),arXiv:1609.06925 [gr-qc].[5] David Hilditch, Enno Harms, Marcus Bugner, HannesR¨uter, and Bernd Br¨ugmann, “The evolution of hyper-boloidal data with the dual foliation formalism: Mathe-matical analysis and wave equation tests,” Class. Quant.Grav. , 055003 (2018), arXiv:1609.08949 [gr-qc].[6] Andreas Schoepe, David Hilditch, and Marcus Bugner,“Revisiting Hyperbolicity of Relativistic Fluids,” Phys.Rev. D97 , 123009 (2018), arXiv:1712.09837 [gr-qc].[7] David Hilditch and Andreas Schoepe, “Hyperbolicity ofdivergence cleaning and vector potential formulations ofgeneral relativistic magnetohydrodynamics,” Phys. Rev.D , 104034 (2019), arXiv:1812.03485 [gr-qc].[8] Edgar Gasperin and David Hilditch, “The Weak NullCondition in Free-evolution Schemes for Numerical Rel-ativity: Dual Foliation GHG with Constraint Damping,”Class. Quant. Grav. , 195016 (2019), arXiv:1812.06550[gr-qc].[9] Miguel Duarte and David Hilditch, “Conformally flatslices of asymptotically flat spacetimes,” Class. Quant.Grav. , 145018 (2020), arXiv:1909.06135 [gr-qc].[10] Edgar Gasperin, Shalabh Gautam, David Hilditch, andAlex Va˜n´o Vi˜nuales, “The Hyperboloidal Numerical Evo-lution of a Good-Bad-Ugly Wave Equation,” Class.Quant. Grav. , 035006 (2020), arXiv:1909.11749 [gr-qc].[11] Shalabh Gautam, Alex Va˜n´o-Vi˜nuales, David Hilditch,and Sukanta Bose, “Summation by Parts and TruncationError Matching on Hyperboloidal Slices,” arXiv e-prints, arXiv:2101.05038 (2021), arXiv:2101.05038 [gr-qc].[12] Bernd Br¨ugmann, “A pseudospectral matrix method fortime-dependent tensor fields on a spherical shell,” J. Comput. Phys. , 216–240 (2013), arXiv:1104.3408[physics.comp-ph].[13] David Hilditch, Andreas Weyhausen, and BerndBr¨ugmann, “Pseudospectral method for gravitationalwave collapse,” Phys. Rev. D93 , 063006 (2016),arXiv:1504.04732 [gr-qc].[14] Frans Pretorius, “Evolution of binary black hole space-times,” Phys. Rev. Lett. , 121101 (2005), gr-qc/0507014.[15] John G. Baker, Joan Centrella, Dae-Il Choi, MichaelKoppitz, and James van Meter, “Gravitational waveextraction from an inspiraling configuration of mergingblack holes,” Phys. Rev. Lett. , 111102 (2006), gr-qc/0511103.[16] Manuela Campanelli, Carlos O. Lousto, Pedro Mar-ronetti, and Yosef Zlochower, “Accurate evolutions oforbiting black-hole binaries without excision,” Phys. Rev.Lett. , 111101 (2006), gr-qc/0511048.[17] J. Thornburg, “Coordinates and boundary conditionsfor the general relativistic initial data problem,” Class.Quantum Grav. , 1119–1131 (1987).[18] Miguel Alcubierre and Bernd Br¨ugmann, “Simple exci-sion of a black hole in 3+1 numerical relativity,” Phys.Rev. D , 104006 (2001), gr-qc/0008067.[19] Edward Seidel and Wai-Mo Suen, “Towards a singularity-proof scheme in numerical relativity,” Phys. Rev. Lett. , 1845–1848 (1992), gr-qc/9210016.[20] P. Anninos, G. Daues, J. Mass´o, E. Seidel, and W.-M.Suen, “Horizon boundary condition for black hole space-times,” Phys. Rev. D , 5562–5578 (1995).[21] G. B. Cook, M. F. Huq, S. A. Klasky, M. A. Scheel, A. M.Abrahams, A. Anderson, P. Anninos, T. W. Baumgarte,N. T. Bishop, S. R. Brandt, J. C. Browne, K. Camarda,M. W. Choptuik, R. R. Correll, C. R. Evans, L. S. Finn,G. C. Fox, R. G’omez, T. Haupt, L. E. Kidder, P. La-guna, W. Landry, L. Lehner, J. Lenaghan, R. L. Marsa,J. Mass’o, R. A. Matzner, S. Mitra, P. Papadopoulos,M. Parashar, L. Rezzolla, M. E. Rupright, F. Saied,P. E. Saylor, E. Seidel, S. L. Shapiro, D. Shoemaker,L. Smarr, W.-M. Suen, B. Szil’agyi, S. A. Teukolsky,M. H. P. M. van Putten, P. Walker, J. Winicour, andJ. W. York, Jr., “Boosted three-dimensional black holeevolutions with singularity excision,” Phys. Rev. Lett. , 2512–2516 (1998). [22] Jonathan Thornburg, “A 3+1 computational schemefor dynamic spherically symmetric black hole space-times – II: Time evolution,” (1999), gr-qc/9906022, gr-qc/9906022.[23] D. Shoemaker, K. Smith, U. Sperhake, P. Laguna,E. Schnetter, and D. Fiske, “Moving black holes via sin-gularity excision,” Class. Quantum Grav. , 3729–3743(2003), gr-qc/0301111.[24] Gioel Calabrese, Luis Lehner, Oscar Reula, Olivier Sar-bach, and Manuel Tiglio, “Summation by parts and dis-sipation for domains with excised regions,” Class. Quan-tum Grav. , 5735–5758 (2004), gr-qc/0308007.[25] Frans Pretorius, “Numerical relativity using a general-ized harmonic decomposition,” Class. Quant. Grav. ,425–451 (2005), arXiv:gr-qc/0407110.[26] U. Sperhake, B. Kelly, P. Laguna, K. L. Smith, andE. Schnetter, “Black-hole head-on collisions and gravi-tational waves with fixed mesh-refinement and dynamicsingularity excision,” Phys. Rev. D , 124042 (2005),gr-qc/0503071.[27] Daniel A. Hemberger, Mark A. Scheel, Lawrence E. Kid-der, Bela Szilagyi, Geoffrey Lovelace, et al. , “Dynami-cal Excision Boundaries in Spectral Evolutions of BinaryBlack Hole Spacetimes,” Class.Quant.Grav. , 115001(2013), arXiv:1211.6079 [gr-qc].[28] Sean Hayward, “Spin-coefficient form of the new laws ofblack hole dynamics,” Class. Quantum Grav. , 3025(1994), gr-qc/9406033.[29] Abhay Ashtekar and Badri Krishnan, “Dynamical hori-zons and their properties,” Phys. Rev. D , 104030(2003), gr-qc/0308033.[30] Maitraya K. Bhattacharyya, David Hilditch, K. Ra-jesh Nayak, Hannes R. R¨uter, and Bernd Br¨ugmann,“Analytical and numerical treatment of perturbed blackholes in horizon-penetrating coordinates,” Phys. Rev. D , 024039 (2020), arXiv:2004.02558 [gr-qc].[31] David Hilditch, Andreas Weyhausen, and BerndBr¨ugmann, “Evolutions of centered Brill waves witha pseudospectral method,” Phys. Rev. D96 , 104051(2017), arXiv:1706.01829 [gr-qc].[32] Lee Lindblom, Mark A. Scheel, Lawrence E. Kidder,Robert Owen, and Oliver Rinne, “A new generalized har-monic evolution system,” Class. Quant. Grav. , S447–S462 (2006), gr-qc/0512093. [33] D. Garfinkle, “Harmonic coordinate method for simu-lating generic singularities,” Phys. Rev. D , 044029(2002), gr-qc/0110013.[34] M. Alcubierre, S. R. Brandt, B. Br¨ugmann, D. Holz,E. Seidel, R. Takahashi, and J. Thornburg, “Symmetrywithout symmetry: Numerical simulation of axisymmet-ric systems using Cartesian grids,” Int. J. Mod. Phys. D , 273–289 (2001), gr-qc/9908012.[35] Miguel Alcubierre, Introduction to 3+1 Numerical Rela-tivity (Oxford University Press, Oxford, 2008).[36] Marcus Bugner, Tim Dietrich, Sebastiano Bernuzzi, An-dreas Weyhausen, and Bernd Br¨ugmann, “Solving 3Drelativistic hydrodynamical problems with WENO dis-continuous Galerkin methods,” Phys. Rev.
D94 , 084004(2016), arXiv:1508.07147 [gr-qc].[37] Hannes R. R¨uter, David Hilditch, Marcus Bugner, andBernd Br¨ugmann, “Hyperbolic Relaxation Method forElliptic Equations,” Phys. Rev.
D98 , 084044 (2018),arXiv:1708.07358 [gr-qc].[38] Olivier Sarbach and Manuel Tiglio, “Continuum anddiscrete initial-boundary value problems and einstein’sfield equations,” Living Reviews in Relativity (2012),arXiv:1203.6443 [gr-qc].[39] Oliver Rinne, “Stable radiation-controlling boundaryconditions for the generalized harmonic Einstein equa-tions,” Class. Quant. Grav. , 6275–6300 (2006),arXiv:gr-qc/0606053.[40] Stephen W. Hawking and George F. R. Ellis, The largescale structure of spacetime (Cambridge University Press,Cambridge, England, 1973).[41] Jonathan Thornburg, “Event and apparent horizon find-ers for 3 + 1 numerical relativity,” Living Rev. Relativity(2006), [Online article], gr-qc/0512169.[42] P. Diener, “A new general purpose event horizon finderfor 3D numerical spacetimes,” Class. Quantum Grav. ,4901–4917 (2003), gr-qc/0305039.[43] Michael I. Cohen, Harald P. Pfeiffer, and Mark A. Scheel,“Revisiting Event Horizon Finders,” Class.Quant.Grav. , 035005 (2009), arXiv:0809.2628 [gr-qc].[44] Andy Bohn, Lawrence E. Kidder, and Saul A. Teukol-sky, “Parallel adaptive event horizon finder for nu-merical relativity,” Phys. Rev. D94