Worldtube excision method for intermediate-mass-ratio inspirals: scalar-field model in 3+1 dimensions

Binary black hole simulations become increasingly more computationally expensive with smaller mass ratios, partly because of the longer evolution time, and partly because the lengthscale disparity dictates smaller time steps. The program initiated by Dhesi et al. [Phys. Rev D 104 , 124002 (2021)] explores a method for alleviating the scale disparity in simulations with mass ratios in the intermediate astrophysical range ( 10 − 4 ≲ 𝑞 ≲ 10 − 2 ), where purely perturbative methods may not be adequate. A region (“worldtube”) much larger than the small black hole is excised from the numerical domain, and replaced with an analytical model approximating a tidally deformed black hole. Here we apply this idea to a toy model of a scalar charge in a fixed circular geodesic orbit around a Schwarzschild black hole, solving for the massless Klein-Gordon field. This is a first implementation of the worldtube excision method in full 3+1 dimensions. We demonstrate the accuracy and efficiency of the method, and discuss the steps towards applying it for evolving orbits and, ultimately, in the binary black-hole scenario. Our implementation is publicly accessible in the SpECTRE numerical relativity code.

It is likely that upcoming observing runs by ground-based detectors will continue to record binaries with small massratios.Future ground-based detectors like the Einstein Telescope [10] and Cosmic Explorer [11], featuring an improved low-frequency sensitivity, will be able to detect the capture of stellar-mass black holes (BHs) by intermediate-mass BHs, with mass-ratios down to  ∼ 10 −3 [12].Moreover, spaceborne detectors, like the LISA observatory [13,14], will be sensitive to binaries with mass ratios in the entire range from  ∼ 1 to extreme mass-ratio inspirals with  ∼ 10 −5 [12,[15][16][17].
In anticipation of this remarkable expansion in observational reach, it is important to develop accurate theoretical waveform templates that reliably cover the entire relevant range of mass ratios.Standard Numerical Relativity (NR) methods [18] work well for mass ratios in the range 0.1 ≲  ≤ 1 (see e. g. [19]).However, simulations become progressively less tractable at smaller , and few numerical simulations have been performed at  < 0.1 so far.The root cause is a problematic scaling of the required simulation time with .Fundamentally, one expects the required simulation time to grow in proportion to  −2 , where one factor of  −1 is associated with the number of in-band orbital cycles, and the second factor  −1 comes from the Courant-Friedrich-Lewy (CFL) stability limit on the time step of the numerical simulation, arising from the requirement to resolve the smaller black hole.The state of the art in small- NR is represented by the recent simulations performed at RIT of the last 13 orbital cycles prior to merger of a black-hole binary system with  = 1/128 [20,21].Head-on simulations, where the needed evolution time is orders of magnitudes shorter than for inspirals, are possible at even smaller mass-ratios [22,23].While these simulations represent an important proof of concept, their computational cost is extremely high, and it is presently impossible to explore the full parameter space including spin and eccentricity.
Binaries with extreme mass-ratios, say  ≲ 10 −4 , corresponding to a compact object orbiting a massive black hole in a galactic nucleus, can be modeled with a perturbative expansion in .This "gravitational self-force" (GSF) approach [24,25] incorporates order-by-order in  the small deviations of the motion of the small body away from the geodesic motion that applies for test-bodies.The GSF approach is the only method for modeling extreme-mass-ratio inspirals, and development is ongoing towards waveform models suitable for signal identification and interpretation with LISA [26][27][28][29][30][31].With NR being well-suited to comparable masses and the GSF approach to extreme mass-ratios, the question arises of how to model the intermediate mass-ratio regime.For simple binary systems (of nonspinning black holes in quasi-circular or eccentric inspirals) NR simulations suggest [32,33] that GSF calculations may be sufficiently accurate even at massratios reaching the NR regime."Post-adiabatic" GSF waveforms [31] for non-spinning, quasi-circular binaries have shown those predictions were somewhat over-optimistic in the  > 0.1 range [34], but they have borne out the prediction for smaller mass ratios ≲ 0.1.However, it remains unclear whether the two methods, separately applied, can achieve reliable waveform models of intermediate-mass-ratio inspirals over the full astrophysically relevant parameter space.
In this paper we continue the work of [35] to develop a new approach to the simulation of intermediate-mass-ratio systems, combining NR techniques with black hole perturbation theory.The general idea is to excise a large region around the smaller black hole.Inside this region-a "worldtube" in spacetime-an approximate analytical solution is prescribed for the spacetime metric, arising from the perturbation theory of compact objects in a tidal environment.An NR simulation is set up for the binary, in which the worldtube's interior is excised from the numerical domain, and replaced with the analytical solution.At each time step of the numerical evolution, the numerical solution (outside the worldtube) and analytical solution inside are matched across the worldtube's boundary, in a process that fixes a priori unknown tidal coefficients in the analytical solution, gauge degrees of freedom, and also provides boundary conditions to the NR evolution.The intended effect of this construction is to partially alleviate the scale disparity that thwarts the efficiency of the numerical evolution at small : The smallest length scales on the numerical domain is now that of the worldtube-radius , rather than the scale  2 of the smaller body.As a result, the CFL limit is expected to increase by a factor / 2 ≫ 1, with a comparable gain in computational efficiency.
In Ref. [35], as also in the present work, we consider a linear scalar-field toy model where the small black hole is replaced with a pointlike scalar charge moving on a circular geodesic around a Schwarzschild black hole.Instead of tackling the full Einstein's equations, one solves the less complicated massless linear Klein-Gordon equation for a scalar field.Our previous work [35] decomposed the scalar field into spherical harmonics and solved the resulting 1+1-dimensional (1+1D) partial differential equation for each mode separately.Such a modal decomposition will not be possible in the fully nonlinear BBH case.As a step towards the BBH case, in this paper, we derive and implement a generalized matching scheme in full 3+1D.Our implementation is publicly accessible as part of the SpECTRE platform [36], a new generalrelativistic code developed by the SXS collaboration, which employs a nodal discontinuous Galerkin method with taskbased parallelism.The input file for the simulations presented in this paper is given as supplemental material.(An evolution of the scalar field equation in 3+1D with a point source was performed in [37] using a different method.) The paper is organized as follows.In Section II we describe our scalar-field model, and formulate it as an initial-boundary evolution problem suitable for implementation on SpECTRE.Section III describes the construction of the approximate analytical solution inside the worldtube.In Section IV, we show how the unknown parameters of this local solution can be continuously determined from the evolution data on the worldtube boundary, using a set of ordinary differential equations (in time) derived from the Klein-Gordon equations.The fully specified solution inside the worldtube is then used to formulate boundary conditions for the evolution system.We present the results of our simulations in Section V, and demonstrate a good agreement with both analytical solutions in limiting cases, and numerical results from other simulations.We explore the convergence of our numerical solutions with worldtube size, and show that its rate matches our theoretical expectations.Finally, in Section VI, we summarize our findings and discuss the next steps in our program.We use geometrized units throughout the text with  =  = 1.

II. NUMERICAL FIELD EVOLUTION OUTSIDE THE WORLDTUBE
We place a pointlike particle with scalar charge  on a fixed, geodesic circular orbit around a Schwarzschild black hole of mass .The evolution of the scalar field Ψ is governed by the massless Klein-Gordon equation, Here   is the inverse Schwarzschild metric, and ∇  is the covariant derivative compatible with it.   () is the particle's geodesic worldline parameterized in terms of proper time .In Kerr-Schild coordinates   = (,   ), parametrized by coordinate time , the worldline with orbital radius   and angular velocity  = (/ 3  ) 1/2 is given by where we have fixed the orbital plane and phase without loss of generality.We excise the interior of a sphere with constant Kerr-Schild ), centered on the particle's position, from the numerical domain.We refer to this excision region as the worldtube and elaborate in Sec.IV how boundary conditions are provided to the evolution domain.
Outside the worldtube, the numerical evolution of the scalar-field variable Ψ N ('N ' for 'numerical', to contrast with the analytical solution inside the worldtube, to be introduced below) is governed by the source-free Klein-Gordon equation on the fixed background spacetime: The background spacetime is given in the usual 3+1 split, where  is the lapse,   is the shift and    is the spatial metric on  = const.hypersurfaces.The background spacetime of our simulations is a single Schwarzschild black hole in Kerr-Schild coordinates.The Klein-Gordon equation is transformed into the standard first-order form by introducing the auxiliary variables [38] This introduces two constraint fields [39], which must vanish for any solution to the original, secondorder evolution equation.Following [40], we write the firstorder evolution equations for the vacuum Klein-Gordon equation (3) as The lapse , shift   , spatial metric    , inverse spatial metric    , trace of the extrinsic curvature  :=       and trace of the spatial Christoffel symbol Γ  :=    Γ    appearing in Eqs.(8) depend only on the background Schwarzschild spacetime.Explicitly, they read as follows in Kerr-Schild coordinates: where  = √︃        is the areal radius from the central black hole.The variables  1 and  2 appearing in Eqs.(8) are constraint damping parameters.Compared to the first-order reduction presented in [39], the additional term  1  2     in Eq. (8b) ensures that the system is symmetric hyberbolic for any values of  1 and  2 [40].For  2 , we found that a central Gaussian profile  2 =  − (  ) 2 +  with  = 10,  = 10 −1 / and  = 10 −4 results in a long-term stable evolution for all tested systems.We choose  1 = 0 throughout.
The evolution equations (8) are in the general symmetric hyperbolic form with   := (Ψ, Π, Φ  ) representing the set of first-order variables, enumerated by the indices  and .For the imposition of boundary conditions at a boundary with normal co-vector n , we solve the (left) eigenvalue problem for the eigenvalues  ( â) and eigenvectors  â  , enumerated by the index â.The  ( â) are known as the characteristic speeds, and the parentheses indicate that there is no implicit sum convention on the right hand side of Eq. ( 11).The co-vector n is normalized with respect to the three metric, i.e.    n n  = 1, and we define n =    n  .The characteristic fields  â are obtained by projecting the evolved variables   onto the set of eigenvectors  â  : For the evolution system (8), the characteristic fields are  â = ( 1 ,  2  ,  + ,  − ) with Here,    =    − n n denotes the projection operator orthogonal to n , so that  2  carries only two degrees of freedom.The corresponding characteristic speeds are   1 = − n   (1 +  1 ) ,   2 = − n   and   ± = − n   ± .We note that the fields  ± reduce to the known physical retarded/advanced derivatives   Ψ ±   Ψ in flat space with  2 = 0.The other characteristic fields result from the reduction of the PDE system to first order.Boundary conditions must be specified at the external boundaries of the domain for each characteristic field, if and only if it is flowing into the domain, specifically those with negative characteristic speeds.There are three external boundaries in our domain: one excision sphere within the central black hole, one excision sphere around the scalar charge (the surface of the worldtube), and the outer boundary.
At the black hole excision sphere, all characteristic fields are flowing out of the computational domain into the excised domain, so no boundary conditions need to be applied.For the outer boundary and at the worldtube boundary, the fields  1 ,  2   may require boundary conditions, while  − always requires ones and  + never requires ones.
Boundary conditions for the physical characteristic field  − at the outer boundary are derived from the second-order Bayliss-Turkel radiation condition [41].These boundary conditions are applied with the method of Bjorhus [42].At the worldtube boundary, the local solution inside the worldtube is used to provide boundary conditions for  − , as explained in detail in Section IV.
Boundary conditions for  1 and  2  can be derived by requiring that there are no constraint violations flowing into the domain [43], as described in Appendix A. These constraintpreserving boundary conditions are applied with the method of Bjorhus [42] at the worldtube boundary and at the outer boundary; see Eq. (A11).
The evolution equations ( 8) are solved with SpECTRE [36], which employs a nodal discontinuous Galerkin (DG) scheme in 3+1 dimensions.The domain is built up of several hundred DG elements, each endowed with a tensor product of Legendre polynomials using Gauss-Lobatto quadrature.The elements are deformed from unit cubes to fit the domain structure using a series of smooth maps as illustrated in Fig. 1.Discontinuous Galerkin methods require a choice of numerical flux that dictates how fields are evolved on element boundaries where they are multiply defined [44].Here we employ an upwind flux.
SpECTRE uses dual coordinate frames [45] to solve the evolution equations.The components of the tensors in the evolution Eqs.(8) are constructed in Kerr-Schild coordinates   .We refer to these as the inertial frame because the coordinates are not rotating with respect to the asymptotic frame at spatial infinity.The evolution equations for the inertial components are solved as functions of co-rotating coordinates (t,  ı ) = (t, x, ȳ, z) given by the transformation Tensor components in this frame we denote with a bar, as in  ᾱ β .For more demanding situations (e.g.binary black hole simulations), the transformation   →  ı can take a much more complicated form [46,47].The grid points of the DG domain, as well as the particle position  ı  = (  , 0, 0) are constant in space in these coordinates, which we will refer to as grid coordinates.The internal worldtube solution is evolved in the grid frame directly, which considerably simplifies the formulation of the matching scheme in Sec.IV.A Dormand-Prince time stepper is used to advance the solution of the numerical fields with a global time step.We apply a weak exponential filter to the evolution fields after every time step to ensure stability of the evolution.
The code is parallelized using the heterogeneous task-based parallelism framework Charm++ [48].The inclusion of the worldtube does not adversely impact the parallel efficiency, as its computational cost is negligible compared to even a single DG element evaluation, and no additional communication between cores is introduced.

III. APPROXIMATE SOLUTION INSIDE THE WORLDTUBE
Inside the worldtube, the scalar field is given by an analytical expansion in powers of coordinate distance from the particle's worldline   .We use Ψ A to denote this analytical solution, and we use a formal parameter  = 1 to count powers of the separation between the worldline and the field point.
As in Ref. [35], we split the field Ψ A into a puncture field Ψ P and a regular field Ψ R : Ψ P is an approximate particular solution to the inhomogeneous equation (1), and it will be fully determined in ad- vance; Ψ R is an approximate smooth solution to the homogeneous equation, and it will be determined dynamically through matching Ψ A to Ψ N at the worldtube boundary.
We express both Ψ P and Ψ R in terms of the coordinate distance Δ  :=   − x  , where x  is a reference point on   .For a given field point   at coordinate time , we let x  :=    () be the point on   at the same value of , such that Δ = 0. Tensors evaluated at x  are written with a tilde, as in g .To facilitate matching Ψ A to Ψ N , we ultimately express both Ψ P and Ψ R in the co-rotating grid coordinates (,  ı ) introduced in Eq. ( 14), but most of this section applies in both inertial and co-rotating coordinates.
Unlike in Ref. [35], for Ψ P we use an approximation to the Detweiler-Whiting singular field [49]; this choice ensures that we can calculate the scalar self-force directly from the regular field Ψ R .Covariant expansions of the Detweiler-Whiting singular field are readily available to high order in ; see [50][51][52], for example, with [50] deriving the scalar singular field to the highest order in the literature, O ( 4 ).These covariant expressions contain several ingredients.First among them is Synge's world function (, x) [53], which is equal to half the squared geodesic distance between  and x.Its gradient, σ := ∇ (, x), is a directed measure of distance from x to .The projection of σ tangent to the worldline has magnitude and the projection normal to the worldline has magnitude where ũ is the particle's four-velocity at time .In terms of these quantities, the covariant expansion of Ψ P through order  2 is given by [50,52] Here are contractions of the Weyl tensor    and its derivative evaluated at the reference point x on the particle's worldline.
We now express the covariant expansion (18) in terms of Kerr-Schild coordinates.To achieve this we follow the method in [52], which begins from an expansion of (, x) Differentiating this with respect to x  , we obtain We then use the identity 2 = σ σ to recursively determine the coefficients Ã , B  , C  and so on.This yields, for example, Ã =1 4 g(,) .We now contract σ with the four-velocity, metric and Weyl tensor to get the coordinate expressions for , , C , and C|  as per their definitions ( 16), ( 17), (19), and (20).Our final expression for Ψ P is obtained by substituting all of these results into Eq.( 18) and re-expanding in powers of Δ  .
We write the result in the style of [54]: Here  1 = √︁ ( g + ũ ũ )Δ  Δ  is the leading coordinate approximation to , and P  (Δ  ) is a polynomial in Δ  of homogeneous order .
The form ( 23) is valid in any coordinate system.In the corotating grid coordinates, the reference point on   is x ᾱ =  ᾱ  () = (,   , 0, 0), and the coordinate separation where ũ = (1 − 3/  ) −1/2 .The polynomials P 3 (Δ  ), P 6 (Δ  ), and P 9 (Δ  ) are too long to be included here.Instead we have made them available online as MATHEMATICA code 1 .We now turn to the regular field Ψ R .Because it approximates a smooth homogeneous solution, we can write it as a Taylor series around x  .In the grid coordinates, such an expansion reads , and so on.The coefficients Ψ R ı1 ...ı  in this series contain the full freedom in the approximate solution Ψ A .However, not all of these coefficients are independent; the field equation imposes relationships between them.As shown in Ref. [55], once the field equation is enforced, only the trace-free piece of each Ψ R ı1 ...ı  is left undetermined.An th-order approximate solution Ψ R contains  =0 (2 + 1) = ( + 1) 2 of these undetermined functions.All other functions of  in Ψ R are related to these by ordinary differential equations (ODEs) that result from the field equations.In the next section we show how all the functions Ψ R ı1 ...ı  () can be determined through the combination of (i) matching Ψ R to Ψ N and (ii) solving the ODEs in  that follow from the field equation.

IV. MATCHING METHOD
The idea behind the matching method is straightforward.We numerically solve the scalar wave equation on a Schwarzschild background, excising the worldtube containing the scalar charge from the numerical domain.Inside the worldtube, the solution is given by the analytical approximation Ψ A = Ψ P + Ψ R described above.Outside the worldtube we have the numerical field Ψ N .We demand where Γ = henceforth represents an equality that holds on the (2+1D) worldtube's boundary Γ.We will show that this matching condition, together with the scalar wave equation, fully determine the regular field Ψ R inside the worldtube.This solution, in turn, provides boundary conditions for the evolution of the numerical field, specifically for  − .We formulate the matching scheme in the co-moving grid coordinates  ı introduced in Eq. ( 14).The Euclidean distance to the particle is defined as  := √︁  ı ȷ Δ ı Δ ȷ = √︁    Δ  Δ  .The boundary of the worldtube is located at  = , with normal vector  ı := Δ ı /.We note that  ı is normalized with respect to  ı ȷ , whereas n in Sec.II is normalized with respect to the 3-metric    .
We now introduce the details of our matching scheme for order  = 0, 1 and 2, which are the expansion orders implemented numerically in this work.The matching scheme for an expansion of arbitrary order  is given in Appendix B. We start by re-writing the Taylor expansion in Eq. ( 25) in terms of the quantities  and  ı , and we introduce an analogous expansion for the time derivative of the regular field: where we now drop the order-counting parameter  = 1.The set of coefficients Ψ R 0 (), Ψ R ı (), Ψ R ı ȷ () have one, three and six independent components, respectively, for a total of ten.We will show that all of these can be uniquely determined at each time step from (i) the numerical field Ψ N (,  ı ) at the worldtube boundary, and (ii) the Klein-Gordon equation (3).

A. Worldtube boundary data
At each time step   we enforce the continuity condition both for the field itself and its time derivative.In the following section we will omit explicit expressions which enforce continuity between the time derivative of the regular field Ψ R (,  ı ) and the numerical field   Ψ N (,  ı ) because they are completely analogous to the expressions for the fields themselves.We will utilize symmetric trace-free (STF) tensors, indicated with angular brackets, e.g. ⟨ 1 ...  ⟩ .Note that ; more details about STF tensors are given in Appendix B. Transforming Eq. (27a) to a STF basis using (B7) yields, at order with  ⟨ı⟩ =  ı and  ⟨ı  ȷ⟩ =  ı  ȷ − 1 3  ı ȷ .Equation ( 29) will be used on the left-hand side in Eq. ( 28).
The right-hand side of Eq. ( 28) is obtained by evaluating the puncture field of Eq. ( 23) and its time derivative at the coordinates of the DG collocation points on the worldtube surface, and subtracting them pointwise from the corresponding values of Ψ N (  ,  ı ) and   Ψ N (  ,  ı ).This expression is then projected numerically onto the set of spherical harmonics defined on the worldtube Γ with constant radius , where are the spherical harmonic coefficients of the numerical, regular field Ψ N (  ,  ı ) − Ψ P (  ,  ı ) and Ω is the area element of the flat-space unit 2-sphere.In practice we use real-valued spherical harmonics and evaluate the integral with the Gauss-Lobatto quadrature used by the DG method.Both the spherical harmonics   and the STF normal vector  ⟨ k1 • • •  k ⟩ provide an orthogonal basis for functions on a sphere.They can be transformed into each other using Eqs.(B9), We have thus expressed both sides of the continuity condition (28) in a basis of STF normal vectors, using Eqs.( 29) and (32).Orthogonality of the STF basis allows us to match order by order in the STF expansion: We emphasise that Eqs.(33) contain two distinct sets of coefficients: The Ψ N, R on the left-hand-sides are expansion coefficients on the surface Γ, whereas the Ψ R on the right-handside are the Taylor expansion coefficients of the solution in the interior, Eq. (27a).The continuity conditions for a field expanded to arbitrary order are given in Eq. (B14).For expansion orders  = 0 or  = 1 the second term in Eq. (33a) falls away.The regular field inside the worldtube is then fully determined by the continuity condition and can directly be used to provide boundary conditions for the future evolution.In one dimension, this is equivalent to a linear polynomial in an interval being fully determined by its two endpoints.For  = 2, Eqs.(33) provide only 9 equations for the 10 coefficients of Ψ R (  ,  ı ) because the monopole of the regular numerical field Ψ N, R ⟨0⟩ in Eq. (33a) contributes to both the zeroth-order coefficient Ψ R 0 and the trace of the secondorder coefficient,  ı ȷ Ψ R ı ȷ .More generally, for arbitrary order, the STF expansion on the worldtube, Eq. (B14), provides only the trace-free components of Ψ N, R , so that boundarymatching determines only the trace-free parts of the expansion Ψ R (  ,  ı ) but not its traces.Therefore, for expansions of order  ≥ 2, additional equations are needed to fully determine the regular field inside the worldtube.These are provided by a series expansion of the Klein-Gordon equation, as we describe below in Sec.IV B.
The coefficients of the regular field's time derivative are determined completely analogously, with the continuity condition Ψ N (  ,   ) is evaluated using its evolution equation (8a) and then transformed into the co-moving grid frame by adding the advective term      Ψ N , where    is the instantaneous local grid velocity.The matching conditions for the time derivative of the regular field Ψ  (  ) are then just the time derivative of the matching conditions for Ψ  (  ), Eqs.(33).

B. Klein-Gordon equation
We rewrite the Klein-Gordon equation (3) in grid coordinates where Γ ρ  μ ν Γ ρ μ ν .The metric quantities   and Γ  are expanded in the grid coordinates  ı to the same order  as the regular field at each time step   .For  = 2 these expansions read The expansion coefficients are given by , and similarly for Γ μ 0 , Γ μ ı , and Γ μ ı ȷ .Due to the spherical symmetry of the Schwarzschild spacetime, and our circular-orbit setup, these expansion coefficients are in fact independent of   .
We now expand the Klein-Gordon equation in powers of  by inserting the expansions for Ψ R ,   and Γ  from Eqs. (27a), (36) and (37), respectively, into Eq.(35).The  ( 0 ) piece of the equation reads This ODE provides an additional, independent relation between the expansion coefficients Ψ R 0 , Ψ R ı and Ψ R ı ȷ , which enables us to determine the remaining, trace degree of freedom of the regular field at  = 2. Specifically, combining Eq. ( 38) with the continuity conditions (33), and using We reduce this ODE to first order and use a Dormand-Prince time stepper to advance the zeroth-order coefficient Ψ R 0 and its time derivative to the next time step  +1 , taking the same global time step as the DG evolution.Together with the continuity conditions (33) at time step  +1 , this completely determines all components of the second-order expansion of Ψ R (,  ı ) in Eq. (27a) at  +1 .The coefficients of the numerical, regular field Ψ N, R ⟨ k0 ••• k ⟩ are updated each sub-step.As initial conditions of the ODE (39) we take Ψ R 0 ( 0 ) = Ψ R 0 ( 0 ) = 0.In Appendix B we formulate the generalization of this method to an arbitrary order , and in particular we derive the generalized form of the ODE on Γ.

C. Boundary conditions for Ψ N
Once the expansion of the regular field has been fully determined, it can be used to provide boundary conditions to the DG elements neighboring the worldtube.DG methods commonly formulate boundary conditions between elements using the numerical flux, and these conditions are applied to each of the characteristic fields defined in Eqs.(13).We use the internal solution Ψ A of the worldtube to provide boundary conditions for the characteristic field  − as if the interior of the worldtube were simply another DG element.From the definition of  − in Eq. (13c) and the definitions in Eq. ( 5), we obtain the boundary condition The analytical solution Ψ A (  ) was defined in Eq. ( 15) as the sum of the regular field Ψ R and the puncture field Ψ P , both of which are now fully determined.The time and spatial derivative are simply obtained from The fields Ψ R (  ) and   Ψ R (  ) are given by Eq. ( 27) and its time derivative.The derivative normal to the worldtube boundary is similarly obtained by taking the appropriate spatial derivative of Ψ R in Eq. (27a) analytically.The expression for the puncture field Ψ P (  ) is given in Eq. ( 25), and its time and normal derivative are computed analytically.We evaluate all of these expressions at the grid coordinates  ı of all DG grid points that lie on element faces abutting the worldtube to formulate pointwise boundary conditions.The value of  − (  ) at the boundary is used to apply a correction to the time derivative of the evolution equations using the upwind flux [44].
We initially tried to provide boundary conditions in the above fashion for all characteristic fields entering the numerical domain, including  1 and  2  .However, we found that this caused substantial constraint violations entering the numerical domain at the worldtube boundary.Instead, we use constraint-preserving boundary conditions for  1 and  2  as described in Appendix A.

D. Roll-on function
The initial conditions we use for the simulations are Ψ =   Ψ = 0 for both the DG fields outside the worldtube and the regular field inside it.The puncture field Ψ P added to the regular field in Eq. ( 40) initially creates a discontinuity at the worldtube boundary, due to the unphysical instantaneous appearance of the scalar charge source  = 0. DG methods are very inefficient at resolving discontinuities within elements, due to the Gibbs phenomenon.
To alleviate this, we multiply the puncture field Ψ P with a roll-on function () that smoothly grows from 0 to 1 (up to double precision) between  = 0 and  =  end .We found that this effectively stretches out the initial discontinuity and causes the fields to settle more smoothly to their final values.
The roll-on function ensures a smooth settling of the solution corresponding to the scalar charge slowly being turned on over its first 4 orbits.We found  end = 300 to be a good choice for the simulations with orbital radius   = 5.

E. Error estimates
To estimate the errors that our matching method incurs, we apply the same analysis as we did for the 1 + 1D case in [35].The estimates follow from a Kirchhoff representation of the scalar field.We first consider the field in the numerical domain, outside the tube Γ.Call this region .Inside , our field Ψ N satisfies the same homogeneous field equation as the exact solution Ψ,   ∇  ∇  Ψ N = 0, but it inherits errors that propagate out from Γ.We introduce a retarded Green's function  (,  ′ ) satisfying where  and  ′ denote any two points, primes denote quantities at  ′ , □ :=   ∇  ∇  , and  4 (,  ′ ) :=  4 (   −   ′ ) √ − .If we now take any point  ∈ , then the equations (41) and □Ψ N = 0 imply the identity (42) Integrating this equation over all  ′ ∈  and then integrating by parts, we obtain the Kirchhoff representation Here Σ  ′ is the outward-directed surface element on .For us the relevant portion of  is the tube boundary Γ, where Σ  ′ = O ( 2 ) Ω.As in Eq. ( 31), here Ω is the area element of the unit 2-sphere.
In the integral over Γ, we may replace Ψ N with Ψ A .Our truncated expansion of Ψ A introduces an inherent O ( +1 ) error in Ψ N ( ′ ) and O (  ) error in ∇  ′ Ψ N ( ′ ) on the worldtube.Equation (43) implies that the O (  ) error in ∇  ′ Ψ N ( ′ ) dominates.Accounting for the O ( 2 ) surface element, we see that this creates an O ( +2 ) error in Ψ N ().
An important takeaway from this analysis is that the error in the numerical domain is suppressed by the small spatial size of Γ.As a consequence, the error converges two orders faster than the analogous error in the 1 + 1D problem in [35].
However, we note that this analysis applies only at a fixed location  outside the worldtube.At a point on the worldtube boundary Γ, the errors in Ψ N are inherently O ( +1 ), and the errors in ∇  Ψ N are inherently O (  ).There is no suppression due to the small spatial size of the worldtube in this case.The same is true of the errors at a point outside the worldtube if we consider a point  that is at a fixed multiple of  away from the worldline rather than at a fixed physical location.
We also note that in applications, we require outputs other than Ψ N : the regular field on the particle's worldline and the self-force, for example.The omitted terms in our expansion (25) scale with a power of distance from the worldline, which might make us expect that we incur no error in Ψ R (  ) and   Ψ R (  ) (and therefore in the self-force).However, we can see this is incorrect by referring again to a Kirchhoff representation of the field.Our method enforces the field equation (1) on Ψ A up to an error ∼  −1 (two derivatives of the truncation error in Ψ A ).If we momentarily ignore that error term in the field equation, and if we consider  to be the interior of Γ and repeat the steps that led to Eq. ( 43), then we obtain the Kirchhoff representation If we now take  to be a point   on the worldline and consider the integral over Γ, then we have  (,  ′ ) ∼ 1/ and ∇  ′  (,  ′ ) ∼ 1/ 2 .We can combine this with Σ  ′ ∼  2 and with the errors O ( +1 ) in Ψ A ( ′ ) and O (  ) in ∇  ′ Ψ A ( ′ ) to deduce that the error in If we take a derivative of Eq. ( 44), we find that the error in   Ψ A (  ) is O (  ).These error estimates apply immediately to Ψ R (  ) as well.
It is also straightforward to see that these estimates are not altered by the O ( −1 ) error in the field equation, which we neglected in deriving Eq. (44).That error contributes an error ∼ ∫  −1  (  ,  ′ ) ′ ∼  +1 to Ψ R (  ), consistent with the error from the boundary integral.
In summary, we expect that for an th-order analytical approximation, our method introduces the following errors: Error in where  is a point outside Γ and   is a point on the particle's worldline.Our numerical results in the next section will bear out these predictions.The error in   Ψ R (  ), and hence in the self-force, will be particularly relevant when we allow the system to evolve (as opposed to keeping the particle on a fixed geodesic orbit).We defer discussion of this to the Conclusion.Finally, before proceeding, we note that our error estimate for   Ψ R (  ) might be too pessimistic in some instances.Specifically, time-antisymmetric components, linked to the dissipative pieces of the self-force, might converge more rapidly with .This is because these components arise from the radiative piece of the field, equal to half the retarded solution minus half the advanced solution [56].For these pieces of the field, we can replace the Green's function  in Eq. ( 44) with its radiative piece,  Rad = 1 2 ( Ret −  Adv ). Rad is smooth when its two arguments coincide (because singularities cancel between  Ret and  Adv ), meaning it does not introduce the negative powers of  that  Ret introduces in Eq. ( 44).We therefore might expect that errors scale with a higher power of  in the dissipative components of   Ψ R (  ).That could be extremely beneficial in practice because dissipative effects dominate over conservative ones on the long timescale of an inspiral [57], and dissipative effects must therefore be computed with higher accuracy.However, our numerical experiments in the next section do not entirely bear out this expectation of more rapid convergence, and we leave further investigation of it to future work.

V. RESULTS
We use a central black hole of mass  for the simulations.The excision sphere inside the black hole has radius 1.9.The outer boundary is placed at 400.We use a CFL safety factor of 0.4.The scalar charge is placed on a circular orbit with radius   = 5 with angular velocity The expansion terms of the puncture field converge more quickly with larger orbital radii of the scalar charge.The truncation error of the puncture field and hence of the worldtube solution is therefore particularly large at the relatively small x-axis z-axis orbital radius of   = 5 used in our simulations.This ensures that the scheme is tested in an extreme region, comparable to binary black holes close to merger.Because the error due to the worldtube is comparatively large for a small   , it can be resolved with a lower resolution in the numerical domain, lowering the computational cost of the simulations.
We have implemented the worldtube scheme with the local solution expanded to orders  = 0, 1 and 2. The radius of the worldtube was varied between 0.2 and 1.6.The simulations were run until the field had settled to its steady state solution over the entire domain, which took between 3000 and 7000, depending on the magnitude of the settled error.Figure 2 shows a cut through the equatorial plane of the computational domain.
Figure 3 plots the steady state solution along two lines cut through the domain at late, constant Kerr-Schild times : one along the co-moving -axis connecting the central black hole center and the scalar charge and one along the -axis normal to the charge's orbital plane.The undulations of the scalar field on the -axis correspond to the arms of the spiral in Fig. 2. The missing part of the '-axis' line corresponds to the worldtube; the field increases strongly near the worldtube, because of the scalar charge contained at the center of the worldtube.
We verify the validity and convergence of our simulations in three different ways: First, in Sec.V A, we compare the value of the regular field Ψ R and its spatial derivative  ı Ψ R at the position of the charge to published numerical results obtained using frequency-domain self-force methods.Second, in Sec.V B we compare with the known axially symmetric analytical solution along the -axis, given below in Eq. ( 54).Finally, we perform an internal convergence test along the comoving -axis in Sec.V C. For each simulation in the following sections, the resolution of the DG domain was increased until it no longer affected the steady-state solution.This guaranteed that the error measured was due to the worldtube, not the numerical evolution.A. Regular field Ψ R at the charge's position The value of the regular field for a scalar charge in a circular geodesic orbit in Schwarzschild spacetime has been calculated in self-force literature.We compare the regular field of our simulations with the results of [58], who quote the value Ψ R ref (0) = -0.01023418/for a circular orbit with radius 5.
For expansion orders  = 0 and  = 1, the regular field at the charge position is given directly by the monopole of the numerical field Ψ N, R ⟨0⟩ (  ) in Eq. (33a) (the second term on the right-hand side is absent for  = 0, 1).For  = 2, it is determined by solving the ODE (39) inside the worldtube.
The relative error is the final value of the regular field at the scalar charge, is shown in the top panel of Fig. (4), with each marker representing a simulation.The dashed lines show fits of the data, for each expansion order , to a relation of the form  ∝   , where  (recall) is the worldtube radius.The bottom panel displays the local convergence order, defined through where   are the worldtube radii in our sample, and   are the corresponding errors.We find that the error always decreases with smaller worldtube radius or higher order  of the local solution as expected.Equation ( 46) indicated a convergence order inside the worldtube of  =  + 1 at sufficiently small worldtube radii.For  = 1 and 2 we find that this prediction is confirmed quite well, with global convergence orders measured as ∼ 2.07 and ∼ 3.08, respectively, and local convergence order uniformly close to this value.At  = 0 we measure a global convergence order of ∼ 1.72 and a local convergence order that appears to decrease with the worldtube radius.This suggests that for  = 0 the scheme is not fully in the convergent regime for the values of  we consider; rather, there are still significant contributions from higherorder terms.At smaller worldtube radii , these higher-orderin- contributions become less significant and the local convergence rate approaches the expected value of 1.
We also compare the gradient of Ψ R at the position of the particle, which enters the expression for the self-force acting on the particle due to back-reaction from the scalar field.The value of the radial derivative is given in [58] as 0004149937/ 2 using Schwarzschild coordinates (  ,   ,   ,   ).We are not aware of any works which report the angular derivative of the scalar field at   = 5.Instead, it was computed for us to be    Ψ R ref (0) = −0.01009125769/ 2 using the frequency domain code of [59].The coordinate transformation from Kerr-Schild time  to Schwarzschild time   is given by which makes the conversion between    Ψ R and the Kerr-Schild radial derivative where in the second line we have transformed into the comoving coordinate frame given in Eqs.(14).The reference values of the regular field's gradient at the particle's position are then given by which we use to compare to our simulation values.Figure 5 compares the radial and azimuthal derivative of Ψ R obtained by our worldtube evolutions against these reference values.Shown are the differences from the reference values for orders  = 1 and 2, with power-law fits ∝   .At zeroth order, the regular field is constant across the worldtube so the derivatives can not be computed.The lower panel of Fig. 5 plots the local convergence order  loc defined in (48).
We argued in Eq. ( 47) that the error of the regular field's derivatives at the particle position should scale with the worldtube radius as ∝   .For  = 1, this behavior is confirmed by the local convergence  loc of our simulations with the exception of worldtube radius  = 1.6, which is anomalously lower than expected and skews the global convergence order.51) and (52).Each cross represents the final value of a simulation.The dashed lines are best fits to the power-law relation ∝   .Bottom panel: The local convergence order between the simulations of adjacent worldtube radii as defined in Eq. (48).The simulations at order  = 1 show a convergence order consistent with the predicted rate  = 1.The  = 2 simulations show a higher convergence rate likely due to dominant higher order terms.Some anomalies are not visible.
For  = 2, the radial derivative   Ψ R (linked to the conservative, time-symmetric piece of the self-force) shows a local convergence that approaches the expected order of  2 at smaller worldtube radii.This suggests that the error is just entering the regime where the O (  ) contribution becomes dominant.For the angular derivative   Ψ R (linked to the dissipative, time-antisymmetric piece of the self-force), the local convergence order is larger than 3 for all simulations.This could indicate that the error is still dominated by higherorder contributions at the sampled worldtube radii.Alternatively, it could indicate that dissipative quantities converge more rapidly with  than conservative one, as suggested in Sec.IV E; however, the results for  = 1 do not support that proposal, showing the same convergence rate for   Ψ R as for   Ψ R .We stress that in any case, the convergence is at least as rapid as predicted in Eq. (47).

B. Solution along the 𝑧-axis
The spherical symmetry of the Schwarzschild background allows for the Klein-Gordon Eq. ( 3) to be decomposed into separately evolving spherical harmonic modes Ψ  (, ), where the spherical harmonic decomposition is centered on the black hole (different to the spherical harmonics introduced . The relative error of the scalar field Ψ along the -axis compared to the analytical solution Ψ  given by Eq. ( 54).We show two simulations with worldtube radii 1.6 and 0.6 for each order, 0, 1 and 2. The error decreases with higher order or smaller worldtube radius, as expected.A small region is cut out around  =   = 5, where Eq. ( 54) converges too slowly to be calculated to sufficient accuracy in practice.
in Eq. ( 30), which are centered on the worldtube).On the polar axis ( =  = 0) all modes vanish except the axially symmetric ones, i.e. those with  = 0.These modes are also static and admit simple analytical solutions [35].Along the polar axis these solutions read where   and   are Legendre functions of the first and second kind, respectively,   is the normal vector pointing in the direction of the  coordinate axis, and Θ is the Heaviside function.The full solution along the -axis is then given by where   is normal vector pointing in the coordinate z direction.The expansion (54) converges exponentially in  everywhere except in the neighborhood of  =   , where the convergence is too slow to yield good results in practice.We therefore ignore this region and cut it out of plots when comparing with the analytical solution Ψ  .Figure 6 shows the relative error |Ψ N −Ψ  |/Ψ  between our numerical worldtube solutions Ψ N and Eq. ( 54), computed at late evolution time, after Ψ N has settled into its steady state.The error is fairly constant along the axis.It is immediately clear that smaller  and higher  lead to improved agreement.
To investigate convergence with worldtube radius, we define the  1 -norm Worldtube Radius [M]  54), integrated between  = 10 and  = 100 using the  1 -norm of Eq. (55).Each cross represents the final value of a simulation.The straight lines are a best fit to the power-law relation  ∝   .Bottom panel: The local convergence order between the simulations of adjacent worldtube radii as defined in Eq. ( 48).The simulations with  = 1 and  = 2 show a constant convergence order consistent with the predicted rate  =  + 2. The  = 0 simulations show a higher convergence rate at larger worldtube radii but approach the expected value at smaller radii.
which we use to integrate the error in Fig. 6 between  = 10 and  = 100 for each simulation.Using this norm, the top panel of Fig. 7 plots the relative differences between the analytical solution Eq. ( 54) and numerical solutions Ψ N using various  and  and evaluated at late time in steady state.Each symbol represents the integrated, relative error of a simulation's final value.Also plotted is a best fit of the error convergence ∝   , where  is the worldtube radius and  is the global convergence order.The lower panel of Fig. 7 shows the local convergence order  loc as defined in Eq. (48).
As explained in Section IV E, we expect the convergence order  =  + 2 in the volume outside the worldtube.At order  = 2, the global convergence order is best fit to  = 4.07 which matches the predicted error.Order 1 has a fitted global convergence order of 3.14 and a local convergence order close to this value across the worldtube radii sampled.For the zeroth-order expansion, a global value of 2.33 is calculated, but the local order consistently decreases with smaller worldtube radii, which suggests that the error might still get contributions from higher-order terms at the larger worldtube radii, similar to the zeroth-order expansion in Fig. 4.

C. Solution along the 𝑥-axis
The tests of our method so far compared to previously known data, either at the position of the charge or on the axis.We now evaluate the convergence with worldtube radius in the volume, at locations where no analytic solution is available.To this end, we evaluate our numerical solutions  N along the co-rotating x-axis, which passes through the center of the Schwarzschild black hole and the point charge.The settled field along this axis is shown as the blue curve in Figure 3.The simulation with  = 2 and  = 0.4 is used as a reference solution, denoted Ψ ref , since it has the lowest error inside the worldtube and along the -axis, as demonstrated above.
Figure 8 shows the relative difference with respect to the reference solution between  = 1.9 and  = 100 for two sample simulations at each order.The field along the -axis has more features as it lies in the orbital plane of the charge.The error along the -axis is therefore not quite as smooth as along the -axis; for instance at  ≈ 80 some features are apparent, which coincide to a wave-crest in Ψ N (see Fig. 3).The error decreases with both a higher expansion order  and a decreasing worldtube radius  as is expected.
To quantify the convergence with respect to  we compute the norm Eq. ( 55) integrated along the co-moving x-axis, ∥Ψ N () − Ψ ref () ∥.This difference, normalized, is plotted in Fig. 9, where each marker represents an individual simulation.The straight lines show power law fits ∝   .In the bottom panel we show the local convergence order  loc as defined in Eq. ( 48).The convergence rates in worldtube radius  are close to the expectation from Eq. ( 45),  =  + 2, with global convergence order  equal to 2.25, 3.18, 4.29 for orders 0, 1 and 2, respectively.For  = 0, the local convergence order  loc is steadily decreasing with the worldtube radius  towards the expected value of  = 2, which suggests it is just entering the convergent regime here.The local convergence rate of the  = 2 simulations appears to jump slightly at the smallest worldtube radius sampled, which we attribute to numerical error.48).At first and second order, the simulations reproduce the expected convergence order of  =  + 2. At zeroth order, the local convergence approaches the expected order for smaller worldtube radii, likely indicating that the error still has contributions from higher-order terms in this regime.

VI. CONCLUSIONS
In this paper we continue the work of [35] and explore a novel approach to simulating high mass-ratio binary black holes.A large region (worldtube) is excised from the numerical domain around the smaller black hole, to alleviate the limiting CFL condition due to small grid spacing in this region.The solution inside the worldtube is represented by a perturbative approximation that is determined by the numerical solution on the boundary and in turn provides boundary conditions to the numerical evolution.
We test this method using the toy problem of a scalar charge in circular orbit around a central black hole.The simulations are carried out in 3+1D using SpECTRE, the new discontinuous Galerkin code developed by the SXS collaboration.In order to develop algorithms that generalize to the full GR problem, we do not decompose our solution into spherical harmonics as is the usual approach.We split the solution near the scalar charge into a puncture field, which is fully determined as a local expansion in Sec.III, and a regular field, which is a smooth Taylor series with undetermined coefficients.The expansion coefficients in the regular field are determined by (i) the numerical solution on the worldtube boundary and (ii) the scalar wave equation as described in Sec.IV.Our puncture is constructed from the Detweiler-Whiting singular field,  allowing us to calculate the scalar self-force from our regular field.
We implement the described matching scheme for orders  = 0, 1 and 2 and perform numerical simulations for a circular orbit of radius   = 5 for a variety of worldtube radii.In Sec.IV E we make a theoretical argument for how the error introduced by the excision should converge with the excision radius in-and outside the worldtube.We confirm these results in Sec.V and show that the scheme solves the scalar wave equation with high accuracy even at relatively large worldtube radii.We further validate our method by comparing against known values of the Detweiler-Whiting regular field and its first derivatives on the particle's worldline.
The ultimate goal of the worldtube method is to speed up BBH simulations at large mass-ratios by alleviating the CFL condition.Figure 10 considers the time steps sizes Δ taken by our primary simulations.Plotted is Δ/ vs the worldtube radius .Note that the resolutions of our simulations were adjusted such that for each simulation, numerical truncation error is subdominant compared to the worldtube error, resulting in differences in Δ for simulations with different expansion orders at the same worldtube radius.Nevertheless, it is apparent that for fixed order, the time step is roughly proportional to the worldtube radius.Therefore, the promise that larger worldtube radii allow larger time steps ∝  indeed holds.Ideally, the worldtube error should be comparable to or somewhat smaller than the NR error.Our results show that this can be achieved by either decreasing the worldtube radius or by increasing the expansion order.The former, of course, would lead to a smaller grid-spacing and a more significant CFL condition, whereas the latter has no noticeable performance cost.We have discussed the next steps towards tackling BBH simulations using the worldtube method in [35].
Before tackling the full BBH problem, our next step will be to include the back-reaction of the scalar field onto the charged particle [60].In the present work we have computed the first derivatives of the Detweiler-Whiting regular field, from which we can construct the scalar self-force, but we have so far ignored the effect of that force.Once it is accounted for, the equations of motion for the scalar charge  of bare mass  0 are given by [61]   ∇  (  ) =   Ψ R , (56) Allowing the particle's trajectory to evolve dynamically in this way will be an important step toward the full gravity problem.SpECTRE uses a series of control systems which adjust certain time-dependent parameters of smooth coordinate maps to deform the grid [45].While they usually ensure that excision spheres stay inside a black hole's apparent horizon, we will use these control systems to enable the worldtube to track the inspiraling scalar charge.The control systems are needed for the binary black hole case and such a scalar evolution will ensure they work as expected.
Symmetric, trace-free (STF) tensors are written with angular brackets around the indices.The combination of  STF normal vectors is defined as [78]  ⟨⟩ =  The  ⟨⟩ provide an orthogonal basis for functions on a sphere, and each  ⟨ ⟩ is an eigenfunction of the Laplacian ∇ 2 :=       , satisfying ∇ 2  ⟨⟩ = − ℓ (ℓ+1)  2  ⟨⟩ .For a fixed , the STF tensors  ⟨ ⟩ span the same functions as the set of spherical harmonics   (, ) of rank .This can be seen by expressing the normal vector as   = (sin  cos , sin  sin , cos ), leading to [77]   = Y * ⟨⟩   ⟨⟩ , (B9a) where ⟨⟩  *  Ω, (B9c) We start by expanding the regular scalar field Ψ R (,   ) and its time derivative in a power series around the charge's position to arbitrary order  as shown in Eq. (B10).We will show that all free components of the expansion can be uniquely determined at each time step from (i) numerical data from the worldtube boundary and (ii) the Klein-Gordon equation (3).
The continuity condition at the worldtube boundary is given by Eq. (28).Both sides of this equation can be expressed in a basis of STF normal vectors.The regular field Ψ R (  ,  ı ) on the left is transformed using Eq.(B7) to give

(B12)
As in Sec.IV, the right-hand side of Eq. ( 28) is calculated by projecting the numerical, regular field onto spherical harmonics up to order  to obtain the coefficients  N, R  .This expansion is then transformed to STF normal vectors using Eq.(B9a),

Figure 1 .
Figure 1.Illustration of the computational domain: Shown is the equatorial plane, with height-deformation proportional to the value of the scalar field.The grid lines correspond to the DG-element boundaries of the 3-D numerical evolution.The central blue/green peak represents the region inside the worldtube, where the approximate solution is dominated by the singularity of the scalar field at the point-charge.Left of the peak an excision region is cut out within the horizon of the central black hole.A zoomed-out view of the entire domain is shown in Figure 2.

Figure 2 .
Figure 2. The equatorial plane of the domain, depicting the steadystate solution of the scalar field Ψ N .The scalar charge creates an outward propagating spiral as it orbits the central black hole.Figure 1 shows a tilted perspective zoomed into the center of the same plane with the spiral arms visible in the background.

Figure 3 .
Figure 3.The steady-state solution of the scalar field Ψ N along the co-moving -axis and -axis of the domain.

Figure 4 .
Figure 4. Top panel: The relative error of the regular field at the position of the charge compared to the value computed in [58].Each cross represents the settled error at the final simulation time.The dashed lines are best fits for the relation  ∝   .Bottom panel: The local convergence order between simulations of neighboring worldtube radii.

Figure 5 .
Figure 5. Top panel: The absolute difference between the radial and angular derivatives of the regular field at the position of the particle and the reference values of Eqs.(51) and(52).Each cross represents the final value of a simulation.The dashed lines are best fits to the power-law relation ∝   .Bottom panel: The local convergence order between the simulations of adjacent worldtube radii as defined in Eq.(48).The simulations at order  = 1 show a convergence order consistent with the predicted rate  = 1.The  = 2 simulations show a higher convergence rate likely due to dominant higher order terms.Some anomalies are not visible.

Figure 7 .
Figure 7. Top panel: The relative difference between the numerical, retarded field Ψ N along the -axis and the analytical solution Ψ  given in Eq.(54), integrated between  = 10 and  = 100 using the  1 -norm of Eq. (55).Each cross represents the final value of a simulation.The straight lines are a best fit to the power-law relation  ∝   .Bottom panel: The local convergence order between the simulations of adjacent worldtube radii as defined in Eq. (48).The simulations with  = 1 and  = 2 show a constant convergence order consistent with the predicted rate  =  + 2. The  = 0 simulations show a higher convergence rate at larger worldtube radii but approach the expected value at smaller radii.

Figure 8 .
Figure 8.The relative error along the -axis between  = 10 and  = 100 for two sample simulations of each order.The worldtube is centered at  = 5 and cut out from the plots.

Figure 9 .
Figure 9. Top panel: The relative error  integrated along the comoving x-axis compared to a reference solution Ψ ref with  = 2 and  = 0.4.Each cross represents a simulation, and the straight lines are a best fit of the relation  ∝   .Bottom panel: The local convergence order as defined in Eq. (48).At first and second order, the simulations reproduce the expected convergence order of  =  + 2. At zeroth order, the local convergence approaches the expected order for smaller worldtube radii, likely indicating that the error still has contributions from higher-order terms in this regime.

Figure 10 .
Figure 10.Time steps Δ of the worldtube simulations presented here.