{'fraction_non_alphanumeric': 0.035983315694453094, 'fraction_numerical': 0.01070783788831476, 'mean_word_length': 4.345535219079312, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Excising a boosted rotating black hole with overlapping grids

We use the overlapping grids method to construct a fourth order accurate discretization of a first order reduction of the Klein-Gordon scalar field equation on a boosted spinning black hole blackground in axisymmetry. This method allows us to use a spherical outer boundary and excise the singularity from the domain with a spheroidal inner boundary which is moving with respect to the main grid. We discuss the use of higher order accurate energy conserving schemes to handle the axis of symmetry and compare it with a simpler technique based on regularity conditions. We also compare the single grid long term stability property of this formulation of the wave equation with that of a different first order reduction.


I. INTRODUCTION
Black hole excision has become an important technique in numerical relativity. First proposed by Unruh [1], excision consists in placing an outflow inner boundary which eliminates the black hole singularity from the domain. Sometimes combined with singularity-avoiding slicings, it is used in all of the current long-term black hole evolutions. (For review articles see Refs. [2,3,4,5], some more recent work includes [6,7,8].) Excision has proved most successful with black holes at a fixed coordinate location. For example, in studies of orbiting compact objects (black holes and/or neutron stars), long runs have been achieved in co-rotating frames, where dynamic gauge conditions attempt to keep the black holes at fixed locations. For more general orbits, and for simulations over many orbital time scales, it is anticipated that the ability to move the black holes on the computational grid through a stable excision algorithm will be crucial.
Excision methods for moving black holes typically require the extrapolation of data onto the trailing edge of the black hole [9,10,11,12,13]. We recently reported on a new excision method using simultaneous, multiple coordinate patches, and implemented these numerically using overlapping grids [14]. In this method each boundary, whether an outer boundary or an excision boundary, moving or static with respect to a main coordinate system, is fixed to at least one coordinate system. The different coordinate patches overlap just like the charts of an atlas. Information is communicated from one grid to the other using only interpolation, without any decomposition into ingoing and outgoing variables. This excision method has a number of advantages, including: (1) the possibility of choosing coordinate patches which conform to each individual boundary, giving smooth numerical boundaries and thus simplifying the boundary treatment; (2) the ability to adapt the coordinates near the black hole to the horizon geometry, allowing the excision surface to be placed relatively far from the singularity, and thereby "excising the excisable"; (3) the use of simple data structures and the implementation of standard methods for parallelization and adaptive mesh refinement due to the fact that individual grids are logically Cartesian.
In our previous work we demonstrated this excision method by evolving a Klein-Gordon massless scalar field on a boosted Schwarzschild background. Numerical tests showed that the scheme was stable and second order convergent for very large boost parameters (v = 0.98c). Long-term convergence test of spherical waves in Minkowski space also showed that the interfaces between the overlapping grids did not introduce growing numerical errors. Moreover, our tests indicated that the overlapping grids technique is robust, in the sense that stability was not dependent on some of the fine details of the implementation. For example, stability was independent of the interpolation method, the number of grids used, the physical overlap size of the grids (which was kept fixed while testing convergence), and their relative resolutions.
Overlapping grids have also been used recently in evolutions of the full Einstein equations. Thornburg [15], for example, used six spherical coordinate patches to excise a Kerr black hole (a = 0.6). The patches are designed such that coordinate lines in one direction overlay each other, so that the interpolation required is effectively one dimensional. Using the BSSN formulation, Thornburg was able to evolve the Einstein equations for 1500 M , and it appeared that the instabilities were related to the outer boundary treatment. He also used overlapping patches in an apparent horizon solver, which is now available in Cactus [16].
More recently, Anderson and Matzner [17] have performed simulations of black holes with overlapping grids. They used two stereographic patches to cover the spherical excision region, and solved the standardġ-K ADM equations in a constrained evolution. They were able to evolve boosted black holes (v = 0.5c) across the computational domain, and runs of single, stationary black holes ran for more than 700 M .
Finally, in a related approach, Reula, Tiglio and Lehner [18] use multiple coordinate patches designed such that all boundary points on neighboring patches are aligned. Rather than overlapping, the grids just touch. The communication between the touching grids is based on the Carpenter-Gottlieb-Abarbanel's method [19], which requires the computation of the characteristic variables. They successfully evolved a Klein-Gordon scalar field on a Kerr background spacetime and are currently working towards fully relativistic evolutions.
In this paper we expand upon our previous work, and choose to study improvements in the following directions: (1) we generalize the background spacetime geometry to allow for spinning black holes, and investigate the effect of high spin (a = 0.99) on the system's stability; (2) we upgrade our derivative and interpolation operators to give global fourth order accuracy; (3) we introduce a spherical patch for the outer boundary, removing all boundary corners.
While in [14] we put strong emphasis on the use of energy conserving schemes based on difference operators satisfying the summation by parts rule, here we allow ourselves to explore discretizations for which a stability proof based on the discrete energy method is not immediately available. At times, such discretizations can be much simpler, particularly in the higher order accurate case.
The paper is organized as follows: In Sec. II we outline the main ingredients of the overlapping grids method. Sec. III is a generalization of Sec. III.B of [14] to the spinning case with the inclusion of a potential in the wave equation. The various coordinate systems used and the regularization of the equations on the axis are discussed in Sec. IV. Before describing in Sec. VI the details of the discretization, in Sec. V we produce group velocity diagrams which illustrate how information propagates through the domain. A high resolution convergence test using the forcing solution method is described in Sec. VII. This paper contains three appendices.

II. OVERLAPPING GRIDS
When solving hyperbolic initial-boundary value problems numerically, it is often difficult, if not impossible, to accurately represent the entire domain, boundary included, with a single grid. For example, at least two coordinate patches are required to cover the entire surface of a 3-sphere without coordinate singularities. When solving for the spacetime representing a binary black hole collision, one may want to use a spherical domain with the outer boundary placed sufficiently far away, in which the black hole singularities are excised. To represent the interior of a sphere with two (or one after the merger) smaller spheres removed, more than one coordinate system is needed, specially if one wants the coordinates to be adapted to the boundaries of the problem.
The overlapping grids methods provides a simple and flexible solution to such problems. Our scheme is based on that described in Starius [20] and Ref. [21]. To demonstrate how the algorithm works we consider the first order FIG. 1: This illustrative example consists of two overlapping grids. The one on the left is the main grid. The other one is adapted to the physical boundary of the problem, which may be moving with respect to the main grid. The gridpoints marked with a solid dot are updated via interpolation. The overlapping grids algorithm is described in the body of the paper.
linear hyperbolic system where u is a vector valued function representing tensor field components, in a domain Ω(t) with smooth boundary ∂Ω(t). The problem includes appropriate initial and boundary data. In general, system (1) is specified in a coordinate system {t, x} which is not and cannot be adapted to the boundary of the problem. We will refer to this boundary as the physical boundary, and we allow for it to be moving with respect to the main coordinate system. We introduce a coordinate system {t ′ , x ′ } adapted to ∂Ω(t), in the sense that the boundary surface can be represented by x ′i = const., for some i. The equations in this coordinate system will take a different form where the components of u ′ are related to the components of u via tensor transformation laws. The initial and boundary data are also transformed. To maintain simultaneity of patches we restrict the coordinate transformation by demanding that t ′ = t. The problem is discretized by introducing a grid for each coordinate system. Any grid point of the main grid which lies beyond a x ′i = const. line (or surface) are dropped, as illustrated in Fig. 1. The two grids overlap and the physical overlap size is kept constant when performing grid refinement tests.
The numerical computation of the right hand side of (1) at a particular gridpoint x ij requires information from a number of gridpoints in each coordinate direction. With a fourth order accurate centered difference operator, for example, one needs two gridpoints in each direction. Gridpoints of the main grid, at which the right hand side cannot be evaluated, are updated via interpolation from the other grid. All components of u ′ are interpolated onto the main grid and the tensor law transformation is used to evaluate the components of u. A similar procedure is done at the non-physical boundary of the second grid. Points that require interpolation are marked with a solid dot in the figure. The numerical treatment of the physical boundary is as that of a single grid boundary.
In our code we use n-th order Lagrange interpolation, which in two dimensions is given by If f is sufficiently smooth, the interpolating function is a n-th order approximation of f (x, y). Note that if the physical boundary moves with respect to the main grid, gridpoints may have to be dropped from or added to the main grid and the set of points which require interpolation needs continuous updating. This is done at the end of each full time step of the time integrator of choice (e.g. fourth order Runge-Kutta). When grid points are added to the main grid, the grid adapted to the boundary is able to provide accurate data for these points.
The overlapping grids method requires artificial dissipation for stability [22]. We use sixth order dissipation, see Eq. (52), which has a seven point stencil in each direction. This means that in our code we actually interpolate three, rather than two gridpoints.
Here we consider the case in which a scalar field propagates on a boosted spinning black hole background. Two boundaries are introduced: an inner and an outer boundary. The first one represents the excision surface and is purely outflow. It requires no boundary data. The second one is introduced for computational reasons. We need to have a bounded spatial domain because of limited computational resources. To handle the two boundaries we introduce two additional coordinate patches. One patch is adapted to the outer boundary and one patch is co-moving with the black hole, and fixed to the inner excision boundary. We choose cylindrical coordinates for the main coordinate patch. The black hole is boosted with velocity β along the axis of symmetry with respect to this coordinate system. Spheroidal coordinates are used on the second patch, such that the location of the event horizon is at a constant coordinate value, and spherical coordinates on the outer patch. We require that all data in all three coordinate systems be simultaneous. By adapting these coordinates to the black hole horizon, we may excise the spherical grid at the event horizon for all values of the boost parameter. This allows for an efficient use of the excision technique, as we can excise the excisable. The layout of the grids is illustrated in Fig. 2. For the fourth order accurate case, which is the minimum order that we require, Eq. (3) on an uniform grid takes the form The boundary treatment and the discretization near the axis of symmetry are described in Sec. VI.

III. A FIRST ORDER REDUCTION OF THE WAVE EQUATION
In this section we write down the wave equation on a general curved background and discuss a particular first order reduction. In Appendix C we consider a different first order reduction, which has different features, and compare the numerical stability of the two formulations. For a definition of strong and symmetric hyperbolicity we refer the reader to [21,23]. The rest of this section is a generalization of Sec. III.B of [14].
The equation of motion for a scalar field propagating on a curved background (M, g) is given by the second order wave equation where ∇ denotes the covariant derivative associated with the Lorentz metric g and V (Φ) is a potential. In terms of the tensor density γ µν = √ −gg µν , where g = det(g µν ), the wave equation can be written as If we introduce the auxiliary variables T = ∂ t Φ and d i = ∂ i Φ, we can rewrite Eq. (6) as a first order system, An attractive feature of this particular first order reduction is that the constraint variables propagate trivially, namely ∂ t C i = 0. This ensures that any solution of (7-9) which satisfies the constraints initially, will satisfy them at later times, even in the presence of (static) boundaries.
The characteristic speeds in an arbitrary direction n, with | n| = 1, are given by the eigenvalues of the matrix These are s ± = (γ tn ± (γ tn ) 2 − γ tt γ nn )/(−γ tt ) = β n ± α √ h nn and s 0 = 0 with multiplicity equal to the spatial dimension of the problem, where α is the lapse function, β i the shift vector, and h ij is the induced 3metric on the t = const. slices in the Arnowitt-Deser-Misner (ADM) decomposition (C1). Hyperbolicity requires that the characteristic speeds be real, namely, (γ tn ) 2 − γ tt γ nn = det(h ij ) h nn ≥ 0 for any n, which will be true as long as the t = const. hypersurfaces are spacelike.
One can verify that where ξ and η are functions of t and x, is the most general symmetric matrix that satisfies HA i = (HA i ) T , i = 1, 2, 3. When positive definite, which will be the case if and only if ∂ t is timelike and ξ > 0, η > 0, it represents the most general symmetrizer of system (7)(8)(9). Using the fact that γ tt < 0 (because the t = const. slices are spacelike) and γ ij = one can show that the symmetrizer is positive definite if and only if the vector field ∂ t is timelike. The symmetrizer can be used to construct an energy and obtain energy estimates. The time derivative of is given by where n i is the outward pointing unit normal to ∂Ω. If V (Φ) is quadratic in Φ, e.g. V = 1 2 m 2 Φ 2 , ∂ t is time-like, and maximally dissipative boundary conditions are used, we have an energy estimate. Note that in this case (12) corresponds to the choice ξ = m 2 /2 and η = 1/2 in (11). If, furthermore, the background admits a time-like vector field k and we use a coordinate system adapted to it, ∂ t = k, the components of γ µν will be time independent and we have, ignoring boundary terms, energy conservation.
The integrand of the surface term can be written as where λ ± = γ n ± γ tn and are the orthonormal characteristic variables of HA n . To simplify the notation we have introduced the quantities γ n = δ µν γ µn γ νn ,γ µn = γ µn /γ n and γ i ⊥ . The latter satisfies δ ij γ i ⊥ γ j ⊥ = 1 and δ ij γ i ⊥ γ jn = 0. To express the primitive variables in terms of the characteristic variables we invert Eqs. (15)- (16), We use Eqs. (15)- (18) in the boundary conditions to prescribe exact data to the incoming characteristic variable w (+λ+;n) . The variable Φ is a zero speed characteristic variable for any direction n and requires no boundary data.
We also assume axisymmetry, which implies that there exists a spacelike Killing field ψ = ψ µ ∂ µ = ∂ φ . By adopting coordinate systems adapted to the Killing field, we have that the metric components are independent of the φ coordinate. Since we are only interested in axisymmetric solutions of the wave equation, i.e., solutions which do not depend on φ, the variable d φ , which represents ∂ φ Φ, can be eliminated from the system.
In the next section we define the various coordinate systems used in our overlapping grids scheme and specialize Eqs. (7)-(9) to these coordinates.

IV. THE KERR METRIC IN KERR-SCHILD COORDINATES
We are interested in the case in which the background is given by a rotating Kerr black hole. We will use Kerr-Schild coordinates [24,25]. After recalling basic properties of these coordinates we explicitly compute the coefficients γ µν needed in the evolution system (7)-(9) and determine the regularized equations on the axis. The Kerr-Schild metric components are given by where and r is determined implicitly by or explicitly by where ζ 2 = x 2 + y 2 + z 2 . The inverse metric can be written as where ℓ µ ≡ η µν ℓ ν = g µν ℓ ν is a null vector.
The quantities M and a are constants, M representing the mass and M a the angular momentum of the black hole as measured from infinity. We restrict ourselves to the case a 2 < M 2 . We recall that the event horizon is located at r = r + = M + (M 2 − a 2 ) 1/2 , the Cauchy horizon at r = r − = M −(M 2 −a 2 ) 1/2 , and the stationary limit surface, the set of points in which the Killing field κ = ∂/∂t becomes null, is given by r = M + (M 2 − a 2 cos 2 θ) 1/2 . Another set of points in which κ is null is given by r = M − (M 2 − a 2 cos 2 θ) 1/2 . The disc x 2 + y 2 ≤ a 2 , z = 0 corresponds to r = 0. The ring x 2 + y 2 = a 2 , z = 0 is a curvature singularity. For later convenience we introduce the quantity ρ 2 BL (r, θ) ≡ r 2 + a 2 cos 2 θ.

A. Boosted cylindrical coordinates
The main coordinate system is obtained by performing a Lorentz boost, followed by a transformation to cylindrical coordinates. Under a Lorentz boost, i.e., in the new coordinatest . At timet the singularity is located atx 2 +ȳ 2 = a 2 ,z = −βt.
Note that in these coordinates the time derivative of √ −g vanishes.
In this coordinate system the symmetrizer associated with this first order reduction does not lead to a conserved energy (∂tγμν = 0), except for β = 0. The region in which the system is symmetric hyperbolic is determined by the set of points in which ∂t is timelike, namely In the unboosted case the system is symmetric hyperbolic outside the stationary limit surface, and only strongly hyperbolic inside (see Fig. 4). For β = 0 the region of lack of symmetric hyperbolicity increases. This issue was explored in greater detail in [14] for the a = 0 case. On the axis of symmetry (ρ = 0) equation (8) needs to be regularized. We want to express it in a form which avoids "0/0". This can be done by taking the limitρ → 0 in the equations. For this purpose it is convenient to introduce the rescaled quantities which have a finite limit forρ → 0 provided that r ≥ r 0 > 0. The right hand side of (8) atρ = 0 becomes In Sec. V we analyze the group velocity of the system and produce plots which give a graphical representation of how information propagates throughout the domain.

B. Co-moving spherical coordinate system
To excise the black hole we introduce a spherical coordinate system {t ′ , r ′ , θ ′ , φ ′ } which is related to the unboosted Cartesian coordinates {t, x, y, z}, the coordinates in which the black hole is at rest, by Note that for β = 0 the time coordinate coincides with the Kerr-Schild time and not that of the Boyer-Lindquist coordinates [26]. In addition, the azimuthal angle φ ′ does not coincide with the Boyer-Lindquist φ BL , but it is related to it via The inverse coordinate transformation of (28) is given by The coordinates are adapted to the event horizon in the sense that its location (r ′ = r + ) is time independent and setting t ′ =t ensures simultaneity with the main coordinate system. The components of the inverse metric can be constructed from where Note that for a = 0 there is a region outside the event horizon in which the system is not symmetric hyperbolic, even when β = 0.
As we did in the cylindrical case on the axis of symmetry (θ ′ = 0 or θ ′ = π) we need to regularize the equations. Introducing the quantities and taking the limit θ ′ → θ 0 , where θ 0 = 0, π, the right hand side of (8) on the axis becomes where the components ofγ µ ′ ν ′ are understood to be evaluated at θ ′ = 0, π. Note that in general the characteristic speeds in the rotating case are higher than those of the non-rotating case. For example, in the radial direction at (r, θ) = (M, π/2) they can be as large as as 4/3 in the extremal unboosted case a = M , β = 0. In the nonrotating case, instead, the same characteristic speeds are bounded by 1 outside the event horizon.

C. Outer spherical coordinate system
In order to have a smooth outer boundary we introduce a spherical coordinate system {t,r,θ,φ}, which is related to the boosted cylindrical coordinate system through the time independent transformation ρ =r sinθ (34) z =r cosθ and its inverser which ensure simultaneity. In this coordinate system we have √ −ḡ s =r 2 sinθ As for the co-moving spherical coordinate system, to have the evolution equation in a form that avoids "0/0" it is convenient to introduce the rescaled quantities We get

D. Transformation of variables
The interpolation procedure which is used to communicate information between the grids requires the coordinate transformations and the tensor transformation laws to relate the components of the evolved fields. The boosted cylindrical and co-moving spherical coordinate systems are related bȳ and the inverse transformation The field Φ transforms like a scalar and the fields T and d i transform like components of 1-forms. In this case we haveT where (T , dρ, dz) and (T ′ , d r ′ , d θ ′ ) are the fields on the cylindrical and spherical grids, respectively. We also need the transformations between the boosted cylindrical and the outer spherical coordinate systems, Eqs. (34)- (35), and and dr = sinθdρ + cosθdz (42) dθ =r cosθdρ − sinθdz (43)

V. GROUP VELOCITY DIAGRAMS
To have a graphical representation of how information propagates through the domain it is useful to construct group velocity or light cones diagrams. These diagrams are obtained by looking at the group velocity of the system at selected points of the domain. We briefly describe how the procedure works. More details can be found in [27,28,29].
Consider a system of quasilinear partial differential equations in two spatial dimensions The characteristic speeds, or phase velocities, in the direction n = (n 1 , n 2 ) are given by the eigenvalues of A i n i . We focus on one of them, which we denote by λ(t, x, u, n). The group velocity is given by its gradient with respect to n, i.e., v g (t, x, u, n) = ∂λ ∂n 1 , ∂λ ∂n 2 In order to produce light cones plots at a certain time t 0 , about a solution u 0 , we select a number of uniformly spaced points of the domain, x i , and plot the parametric function where n(τ ) = (cos τ, sin τ ) with τ ∈ [0, 2π) and α is a constant introduced for aesthetic reasons. It avoids that the cones either overlap or are too small. We produced plots for the linear system (7)-(9) in the cylindrical coordinate system. Fig. 3 illustrates the light cone structure in the non-rotating case, whereas Fig. 4 represents the light cones in the rapidly spinning case, and Fig. 5 shows the rapidly spinning, highly boosted case. It is interesting to see how in Fig. 5 the light cones tilt. While the light cones that are being approached by the black hole are almost unperturbed, those behind are tilted in such a way that the information is forced to follow the black hole. This "dragging" effect is clearly noticeable in our numerical simulations.

VI. DISCRETIZATION
To discretize system (7)-(9) we replace the partial derivative ∂ i with the finite difference operator D i , without expanding derivatives of products. This leads to the semi-discrete system The brown dot is the ring singularity. The outer region between the event horizon and the green line, the stationary limit surface, represents the ergoregion. Between the two green lines the vector field ∂t is spacelike and, therefore, the system is only strongly hyperbolic.
where we have not explicitly written out the gridfunction indeces. As was shown in [14], this type of discretization leads to discrete energy conservation provided that ∂ t is a time-like Killing field [47].
In the interior we use the centered fourth order accurate finite difference operators = (−u i+2,j + 8u i+1,j − 8u i−1,j + u i−2,j )/(12h 1 ) We now discuss the discretization near the axis of symmetry and near the physical boundaries.

A. Axis of symmetry
The difference operators (50) and (51) have a 5 point stencil. Gridpoints that are too close to the axis require special treatment. If a gridpoint is close to but not on the axis of symmetry we use Eqs. (47)-(49) combined with the regularity conditions to evaluate difference operators in the direction normal to the axis. With respect to the coordinate normal to the axis Φ, T , and d A are even and d n is odd, where n is normal to the axis and A is tangent. If a gridpoint lies on the axis we use the regularized equations combined with the regularity conditions. Specifically, in the boosted cylindrical coordinate case we use ∂tT = γtzD (z)T + 2γtρT + D (z) (γtzT ) + 2γρρD (ρ) reg dρ + 2γρzdz + D (z) (γzzdz) + ∂tγttT + ∂tγtzdz + − dV dΦ /(−γtt).
In the co-moving spherical grid on the axis of symmetry (θ = 0 and θ = π) we use the approximation whereas in the outer spherical grid we use The operator D (n) reg represents the centered fourth order accurate operator computed using the regularity conditions. Note that this discretization is not energy conserving. This is discussed further in Appendix B. Clearly, the derivative operator in the direction parallel to the axis needs to be modified near the physical boundary and boundary data needs to be provided at the outer boundary. This is discussed in the next subsection.

B. Boundary conditions
We always place the excision surface r ′ = r ′ min between between the Cauchy and the event horizon, in which case no boundary conditions are required. The fourth order accurate difference operator near the excision surface in the normal direction to it, however, needs to be modified. The modification is given in Appendix A. A similar modification of the difference operator is required near the outer boundary of the outer spherical grid,r =r max . In addition, here we give data to the incoming characteristic variables in the normal direction, including at those points which lie on the axis, where the normal is chosen to be parallel to the axis (see Fig. 4 of [14]). We do not couple the ingoing to the outgoing characteristic variables. Our numerical implementation is based on Olsson's method [30].

C. Artificial dissipation
It is known that overlapping grids require artificial dissipation for stability [22]. Therefore to the right hand side of the discretized system we add a term of the form We modify this operator near the physical boundaries according to where we have only kept the relevant gridpoint index. The modification near i = N is obtained from the above with the replacements i → N − i and D ± → −D ∓ . This modification was derived by combining the extrapolation conditions h 3 D 3 + u −j = 0 and h 3 D 3 − u N +j , j = 1, 2, 3 with (52) whenever gridpoints outside the domain are needed. The result is then multiplied by h to ensure that the modification is a third order correction, which does not lower the global accuracy of the scheme [31]. Near and on the axis of symmetry dissipation (in the normal direction) is computed exploiting the regularity conditions of the fields. We were not able to modify the dissipative operator such that a discrete energy estimate holds and without lowering the global accuracy of the scheme.

D. The discrete energy method
Unlike for the second order accurate case considered in [14], we are not able to show that our discretization satisfies a discrete energy estimate, not even on a single grid and with homogeneous boundary data. We should stress that the inability to obtain a discrete energy estimate is by no means a proof of instability. The discrete energy method is only a sufficient condition for stability.
In Sec. VII numerical experimentation is used to determine the stability of our scheme.

E. Choice of time step and dissipation parameter
We obtain the fully discrete system by integrating the semi-discrete system of ordinary differential equations with fourth order Runge-Kutta. Whenever explicit finite difference schemes are used to approximate hyperbolic problems, the ratio between the time step size k and the mesh size h = min{h i }, the Courant factor, cannot be greater than a certain value [32], which is inversely proportional to the characteristic speeds of the system. In Fig. 6 we numerically estimate allowable values for the Courant factor by examining the 2D wave equation written in first order form, ∂ t u 0 = ∂ i u i , ∂ t u i = ∂ i u 0 . We plot the Courant limits for fourth order Runge-Kutta as a function of the artificial dissipation parameter, assuming fourth order, centered differencing for the spatial derivatives (50)-(51) and sixth order dissipation (52). A similar plot for the second order accurate case can be found in Fig. 7 of [14]. From the plot we see that the value σ ≈ 0.0075, while damping the high frequency modes, does not force us to reduce the Courant factor. This is the value of dissipation parameter that we use in our simulations. It is about three times smaller than the correspondent value for a second order accurate approximation (σ ≈ 0.02).

VII. NUMERICAL EXPERIMENTS
In this section we outline the test carried out to verify that our overlapping grid code is fourth order convergent. We faced some unexpected issues when testing for long term stability. These are discussed in Appendix C.
Let u be the exact solution of the continuum problem at time t and v h the solution of the fully discrete approximation at time step n, such that t = kn, obtained with a mesh size h. The convergence rate is computed as where · h is a discrete L 2 -norm. Neglecting round-off errors one should expect that lim h→0 Q h = 4 for a fourth order accurate scheme. To use this equation we need an exact solution of the continuum problem. We use the same forcing solution technique described in [14]. Let us rewrite the partial differential equation as L(u) = 0 and let w be an arbitrary function. If w is inserted into the equation, in general, it will produce a non vanishing right hand side, Clearly, the modified equationL(u) ≡ L(u) − F = 0 has w as an exact solution and the convergence of the code can be tested using Eq. (53). We chose w(t, r, θ) = sin(t + r) cos(nθ), where {t, r, θ} are the spheroidal coordinates of the rest frame and n is an integer. This is an exact solution of as long as F is given by Both w and F are scalar quantities. The evolution equation (8) is modified according to where g = det(g µν ). On the axis of symmetry we use the limits lim θ→0 sin nθ sin θ = n and lim θ→π sin nθ sin θ = (−1) n+1 n. The result of our convergence test for M = 1, β = −0.85, a = 0.99, and n = 2, confirms that we have fourth order convergence on each grid and is shown in Fig. 7. The simulation was stopped when the inner spherical grid was about to touch the outer spherical grid.
Furthermore, Fig. 8 illustrates that failing to add artificial dissipation to the scheme can lead to numerical instabilities.

VIII. CONCLUSION
The difficulties associated with the construction of a stable scheme for the wave equation on a black hole background are also present in simulations of fully nonlinear general relativistic systems. Therefore, the ability to handle a scalar field propagating around a boosted spinning black hole is one of the first necessary steps towards the construction of a successful, long term stable binary black hole collision numerical code from which gravitational wave templates can be extracted. In this work we showed that the overlapping grids method is compatible with fourth order accuracy, provided that sufficiently high order interpolation and artificial dissipation are used. We also noticed that simple schemes at times work better than more complicated techniques based on conservation of discrete quantities (see, for example, the treatment of the axis in Appendix B).
In Appendix C we observed that different hyperbolic first order reductions, which have different levels of hyperbolicity (symmetric or only strongly hyperbolic) inside (or near) the excision region, have different stability properties. This was unexpected since in the continuum limit the treatment of the excision region, if contained inside the black hole, should not affect the solution outside and therefore should not introduce any instabilities [48]. We noticed that for the T -formulation the behavior of the solution is highly sensitive to the discretization (changing the way the dissipative operator is modified at the inner boundary, for example, can have noticeable effects on the growth rate of the error). It would be interesting to combine the Π-and T -formulations by using the former in the inner spheroidal patch and the latter in the other patches, closely resembling the "interpolating formulation" introduced in [33].
In this work, as in [14], we only considered fully first order systems. Recently, second order in space systems have generated more interest, both at the continuum and discrete level [34,35,36,37,38,39,40,42]. We intend to apply the overlapping grid method to such systems and, most importantly, to the fully nonlinear dynamical case, in which non trivial issues arise. Unlike in the toy model cases that we have investigated so far, in the nonlinear case one has to face the problem of tracking a suitable outflow surface containing the black hole singularity and the dynamical generation of a grid adapted to this surface. We are currently working on a number of related points.
Overlapping grids represent possibly the simplest and most flexible technique capable of accurately handling smooth and time dependent boundaries within finite differencing. It is our hope that it will become a useful tool for the the binary black hole problem.

Acknowledgments
We would like to thank Eric Hirschmann for many interesting discussions. This research was supported in part by NSF and NASA under grants NSF-PHY-0244699, NSF-PHY-0326311 and NASA-NAG5-13430 to Louisiana State University and the Horace Hearne Jr. Labratory for Theoretical Physics. GC was also supported by a Marie Curie Intra-European Fellowship within the 6th European Community Framework Program.

APPENDIX A: MODIFICATION OF THE 4TH ORDER ACCURATE OPERATOR
The fourth order accurate centered difference operator is given by When no boundaries are present (or if the gridfunctions are periodic), we have where (u, v) h = h j u j v j . In [21] it is shown that, by appropriately modifying the difference operator and the scalar product near and at the boundary, one can recover the summation by parts rule [43] (u, We are interested in a globally fourth order accurate scheme. For this purpose one could use a sixth order accurate operator at the interior with a third order accurate modification near the boundaries, or a fourth order accurate operator in the interior with a third order accurate modification near the boundaries. Whereas the first operator satisfies the summation by parts rule with respect to a diagonal scalar product, the second one requires a non-diagonal scalar product Since the global accuracy is the same, in this work we choose the second option as it involves a smaller stencil and therefore requires fewer computational resources. Its modification near the boundary is given by The expression for the σ ij coefficients can be found in [21], where stability proofs for linear hyperbolic problems without corners can be found.

APPENDIX B: HIGHER ORDER ACCURATE DISCRETIZATIONS NEAR THE AXIS OF SYMMETRY
To illustrate the difficulties that arise with higher order discretizations of axisymmetric systems we consider the polar wave equation written in first order form where P | ρ=0 = 0. This system admits the following conserved energy in the sense that its time derivative gives only boundary contributions From the vanishing of P at ρ = 0 we have that Hence we consider the following semidiscrete approximation and discrete energy Let us assume that D is a finite difference operator satisfying the summation by parts rule (A2) with respect to a diagonal scalar product. The time derivative of (B9) gives In order to have an energy estimate consistent with the one of the continuum we need to ensure that the last term vanishes, i.e., the difference operator at the axis has to satisfy 2αDP 0 h = σ 0 D(ρP ) 0 (B10) for some α. As was pointed out in [14], if a first order accurate one sided difference operator is used on the axis, then (B10) is satisfied for α = σ 0 /2 = 1/4.
Using Maple we were able to obtain a (non unique) locally second order accurate modification of the difference operator which leads to discrete energy conservation. However, we had to exploit the regularity conditions at the axis (P is odd and T and ρP are even functions of ρ). The semidiscrete system where D is the standard fourth order accurate centered approximation of ∂ ρ , conserves the following discrete en- Numerical experiments done at high resolutions suggest that the overall order of accuracy of the scheme is three. We now compare this discretization with a fourth order accurate approximation, obtained by simply exploiting the regularity conditions of the fields. We use the approximation (B6-B8) where D is the standard centered fourth order accurate finite difference operator, appropriately modified at the outer boundary and use the fact that T and ρP are symmetric across the axis. We see that the time derivative of the discrete energy (B9) is Clearly, there is no value of α that leads to cancellations of the terms near the axis. Furthermore, if corrections of the form T 2 i h 2 or P 2 i h 2 are added to the energy, products like T i P i+2 or P i T i+2 appear in the estimate. Also, any mixed product P i T j h 2 would give rise to P 2 and T 2 terms. This discretization does not conserve the energy (B9). In particular for α = 1/4 one is left with the axis terms 1 6 (T 0 P 2 + 2T 1 P 1 )h .
Clearly, the fourth order accurate discretization based on regularity conditions, although not energy conserving, is much simpler. The introduction of ghost zones makes it numerical implementation trivial. Fig. 9 illustrates how, for a particular problem and with some particular initial data, the energy of the fourth order accurate operator oscillates about the value of the energy conserving scheme.
Furthermore, no long term growth of the error was observed with the simpler fourth order approximation. In fact, our experiments indicate that, in general, the errors are much smaller than with the energy conserving scheme (B11). This is due to the lower order of local accuracy near the axis of symmetry and is illustrated in Fig. 10.
In our overlapping grid code we adopted the simpler fourth order accurate approximation near the axis.

APPENDIX C: A COMPARISON BETWEEN TWO FIRST ORDER REDUCTIONS
In [14] we performed long term stability tests to check that the interpolation process did not introduce unwanted growth of the error. There we used M = a = β = 0 and had only the spherical inner patch and the cylindrical patch. Here we repeat the same test, but with nontrivial values for the black hole mass and spin. Interestingly, this test revealed that the approximation for the initial-boundary value problem for the reduction (7)- (9), which we will refer to as T -reduction from now on, suffers from exponential growth of the error whenever the domain contains a region in which the system is only strongly hyperbolic. These errors are concentrated near the excision region and, at the resolutions that we typically use, become evident only after about a hundred crossing times.
To exclude the possibility of the overlapping grid method being responsible for this growth, we eliminated all grids, except the inner one. The growth was still there. Furthermore, the growth was present also when the second order accurate, energy conserving method for axisymmetric systems developed in [14] and the modification of the dissipative operator constructed in [44] were used.
We recall that the first order formulation (7)- (9) is symmetric hyperbolic only outside the ergoregion (outside the black hole in the non rotating case). It is important to realize that the use of the discrete energy method to obtain stable discretization relies on the fact that the energy is a positive definite quadratic form of the main variables. Not having symmetric hyperbolicity on the entire domain, we have no guarantee of obtaining a stable discretization.
This prompted us to investigate a different first order reduction, which, unlike the T -reduction, is symmetric hyperbolic everywhere. The energy associated with the symmetrizer, however, is not conserved and therefore we have no preferred way of discretizing the system. After describing this alternative reduction we experimentally compare the stability properties of the two formulations.
As in Sec. III we start with the wave equation around a curved background Using the ADM decomposition of the metric [45] where h = det(h ij ). We introduce the auxiliary variables Π = (β i ∂ i Φ − ∂ t Φ)/α and d i = ∂ i Φ. Defining the wave equation can be written in the form We will refer to this first order reduction of the wave equation as the Π-reduction. In the particular case of the Kerr-Schild metric we have The principal part is given by where β n = β i n i and h jn = h ji n i . The characteristic variables and speeds in the direction n are where A is a vector orthogonal to n. The inverse transformation is The most general symmetrizer is where ξ = ξ(t, x) > 0 and η = η(t, x) > 0. This reduction, unlike the T -reduction, has an everywhere positive definite symmetrizer (it is symmetric hyperbolic) if and only if the t = const. slices of spacetime are space-like. Using the regularity conditions (the fields Φ, Π, d r are even and d θ is odd) and the fact that on the axis ∂ θ α = 0 and αh ij Γ θ ij d θ = −αh θθ ∂ θ d θ we obtain the following regularized equation on the axis We compare single grid, long term stability of the Tand Π-formulations using the same second order accurate finite difference operator modified as in [14], where one sided difference operators were used at the physical boundaries and at the axis. The dissipative operator is modified according to prescription given in [44] at the physical boundaries and regularity is used at the axis. No special grouping of variables is done for the Π-formulation. Long term stability tests were done using the forcing solution of Sec. VII with n = 2. At the outer boundary r = r max we give data to the incoming characteristic variables of A n . (In the T -formulation we actually give data to the incoming characteristic variables of HA n . One can show that both methods control the boundary term arising in the energy estimate [46].) Our tests revealed that the Π-formulation, despite not being energy conserving, outperformes the T -formulation in terms of stability and long term stability. Fig. (11) clearly shows that if the excision boundary is placed well inside the event horizon the T -discretization is unstable, whereas the Π-discretization is stable. Interestingly, in the highly spinning case, the rate of growth of the error in the T -formulation does not seem to increase with resolution. Note that by excising at r min = M = 1, the region of lack of symmetric hyperbolicity is actually larger in the non-spinning case.
In an attempt to understand the different behavior of the T -and Π-reductions of the wave equation, we perform a Laplace-Fourier analysis of an analogous constant coefficient toy model initial-boundary value problem.
Consider the wave equation around Minkowski in Cartesian coordinates where η µ ′ ν ′ = diag{−1, +1, +1, +1}. Under the following change of coordinates where β i are constant, the wave equation takes the form We consider two different types of reductions: the Treduction and the Π-reduction Notice that the constraints C j ≡ d j − ∂ j φ = 0, which are introduced in the reduction process, propagate differently in the two formulations. Whereas in the first case we have ∂ t C j = 0, in the second case we have ∂ t C j = β i ∂ i C j . In both cases the φ variable decouples from the system. We will drop it in the analysis that follows.
The T -reduction has the following symmetrizer and characteristic speeds 0, β n ± 1 .
The Π-reduction has a simpler symmetrizer and the characteristic speeds are given by β n , β n ± 1 .
Most importantly, whereas (C16) is positive definite for any β i , (C14) is positive definite if and only if δ ij β i β j < 1. The last condition is equivalent to the requirement that the vector field ∂ t be time-like. Assume β 1 > 1 and consider the quarter space problem x > 0 with periodic solutions in y and z and no boundary conditions at x = 0. With the Π-reduction energy estimates can be obtained in the outflow case using the energy method. For the T -formulation a Laplace-Fourier analysis yields the following eigenvalue problem (a system of ordinary differential equations in the variable x) If an eigenvalue with ℜ(s) > 0 exists, then the initialboundary value problem is ill-posed in any sense [21]. The eigenvalues are Because of the requirement that the solution belongs to L 2 (0, +∞), we must discard those eigenvalues which have positive real part. Since ℜ(λ ± ) > 0 if ℜ(s) > 0, the problem is not obviously ill-posed.
This analysis reinforces our suspicion that the observed exponential frequency dependent growth is merely due to the discretization. More work is needed to exactly establish the cause of the growth. It is possible that applying the Laplace transform method to the semi-discrete problem may shed some light.  [8,14] × [0, π], 120 × 120 for the outer spherical grid. The lack of stability is evident. On the other hand, as shown in the inset, when artificial dissipation is added to the overlapping grids scheme (σ = 0.02) stability is restored. FIG. 9: We compare the energy conserving discretization (B11) with the simpler fourth order accurate one based on regularity. We plot the discrete energy (B9) with α = 1/4 as a function of time. The grid consist of 200 gridpoints and the domain is 0 ≤ ρ ≤ 20. We use 4th order Runge-Kutta with a small Courant factor, ∆t/∆ρ = 1/10 to be close to the semidiscrete approximation. We use initial data T = sin 6 (ρ) for 0 ≤ ρ ≤ π and zero elsewhere. The energy of the simpler approximation oscillates about the value of the energy of the energy conserving scheme. We use 50 gridpoints and the domain 0 ≤ ρ ≤ 5. The Courant factor is ∆t/∆ρ = 1. We use the forcing solution method with exact solution Φ = ρ sin ρ sin t and forcing term F = − sin t(sin ρ+3ρ cos ρ)/ρ. The errors of the simpler fourth order accurate approximation are clearly much smaller than those of the energy conserving scheme.