Numerical convergence of model Cauchy-Characteristic Extraction and Matching

Gravitational waves provide a powerful enhancement to our understanding of fundamental physics. To make the most of their detection we need to accurately model the entire process of their emission and propagation toward interferometers. Cauchy-Characteristic Extraction and Matching are meth-ods to compute gravitational waves at null inﬁnity, a mathematical idealization of detector location, from numerical relativity simulations. Both methods can in principle contribute to modeling by providing highly accurate gravitational waveforms. An underappreciated subtlety in realising this potential is posed by the (mere) weak hyperbolicity of the particular PDE systems solved in the characteristic formulation of the Einstein ﬁeld equations. This shortcoming results from the popular choice of Bondi-like coordinates. So motivated, we construct toy models that capture that PDE structure and study Cauchy-Characteristic Extraction and Matching with them. Where possible we provide energy estimates for their solutions and perform careful numerical norm convergence tests to demonstrate the eﬀect of weak hyperbolicity on Cauchy-Characteristic Extraction and Matching. Our ﬁndings strongly indicate that, as currently formulated, Cauchy-Characteristic Matching for the Einstein ﬁeld equations would provide solutions that are, at best, convergent at an order lower than expected for the numerical method, and may be unstable. In contrast, under certain conditions, the Extraction method can provide properly convergent solutions. Establishing however that these conditions hold for the aforementioned characteristic formulations is still an open problem.


I. INTRODUCTION
The use of present and improved future gravitational wave (GW) detectors such as the ground-based advanced LIGO and Virgo [1,2], Kagra [3], the Cosmic Explorer [4], and the Einstein Telescope [5], as well as the space-borne detectors like LISA [6], Taiji [7], and Tian-Qin [8], is already deepening and will continue to deepen our understanding of gravity.To harness their power, gravitational waveforms of high fidelity are necessary for the parameter estimation process [9].These waveform models are used to compare against observational data, and thus to infer the underlying properties of the sources.
In the modeling process, a GW detector is typically assumed to lie infinitely far away from the source in an asymptotically flat spacetime.After emission, GWs propagate toward future null infinity at the speed of light where they can be detected.There are different methods to calculate the detected gravitational radiation; see [10] for a review.A widely used technique is the extrapolation of the Weyl scalar ψ 4 to null infinity, which, in suitable coordinates [11], can be used to predict the effects of GWs on a detector.For this method, the initial boundary value problem (IBVP) that results in the production of GWs is solved, in a bounded numerical domain with outer boundary r out .The value of ψ 4 is then calculated on worldtubes of different radii, and consequently used to fit a polynomial expansion, which is finally taken to calculate ψ 4 at the limit r → ∞.This approach however is susceptible to systematic extrapolation errors, which can be avoided using the Cauchy-characteristic extraction (CCE) method [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29].CCE, as illustrated in Fig. 1, allows one to calculate the gravitational radiation at future null infinity I + directly.In CCE, information from the Cauchy domain is passed through a worldtube T to the characteristic one, and propagated to I + .
Since the evolution of the Cauchy and the characteristic setups are decoupled, this method can still be affected by errors related to artificial boundary conditions imposed on the outer boundary of the Cauchy domain.One way to control these errors is to enlarge the Cauchy domain, such that the artificial boundary conditions on r out do not affect the solution within a smaller part of the domain, specifically up to the worldtube r ext .The characteristic setup is then used to extract the solution from r ext to I + .However, the initial data provided on the initial null hypersurface for the characteristic initial boundary value problem (CIBVP) very often are not compatible with the IBVP solution outside r ext , so there is still a type of error that has not been removed [30].
A way to avoid this type of error as well is to evolve the Cauchy and the characteristic setups simultaneously, and allow information to flow via T both ways.This approach is often called Cauchy-characteristic matching (CCM) and can in principle provide gravitational waveforms of high fidelity [20,31,32].Recently, there has r out r ext

T I +
c h a r a c t e r i s t i c Cauchy i + i 0 FIG. 1.The conformal diagram of part of Minkowski spacetime, where spatial i 0 , future temporal i + and future null I + infinities are shown.In gravitational waveform modeling, the detector is assumed to be at I + .As is common in this context, we abuse the terminology "Cauchy" to denote the IBVP for a domain foliated by spacelike hypersurfaces that do not extend to i 0 , but are truncated to some radius r * , with r * = rout for CCE and r * = rext for CCM.By "characteristic" we refer to the CIBVP for GR.For Cauchy-characteristic extraction (CCE) and matching (CCM) information between the Cauchy and characteristic domains is communicated via the worldtube T .
been an important effort to improve and optimize characteristic codes that can perform CCE and CCM [23,29].However, novel results regarding the hyperbolicity of general relativity (GR) in the characteristic setup [33,34] pose pressing questions as to whether CCE and CCM as currently performed for GR can indeed deliver on their promise to produce highly accurate waveforms.Within the framework of numerical relativity, the CIBVP for the Einstein field equations (EFE) is typically constructed using a Bondi-like coordinate system [20,35,36].In CCE and CCM, the resulting Bondilike CIBVP is combined with the IBVP with a strongly or symmetric hyperbolic formulation of GR, such as the generalized harmonic formulation [37][38][39].Characteristic codes for GR based on a Bondi-like formulation have been widely used in numerical simulations.For instance, the first long-term evolution of a single black hole by the Binary Black Hole Grand Challenge Alliance [40], was achieved in such a setup.Many more codes have solved a Bondi-like CIBVP providing, for instance studies of relativistic stars [41,42], and gravitational collapse [43][44][45][46][47][48][49].This success has led to the belief that these Bondi-like CIBVPs are well posed.
A partial differential equation (PDE) problem is well posed if it has a unique solution that depends continu-ously on the given data in an appropriate norm.Let us for example consider the first order, linear, constant coefficient system with real-valued variables and coefficients where u is the state vector, ∂ p denote spatial derivatives, B p are the principal matrices and the product Bu denotes source terms.This system is weakly hyperbolic (WH) if the principal symbol P s ≡ B p s p has real eigenvalues for all unit spatial covectors s p , and is called strongly hyperbolic (SH) if furthermore P s is uniformly diagonalizable for all s p .In addition, if all the principal matrices B p can be brought simultaneously in a symmetric form, then the system is called symmetric hyperbolic (SYMH).
Let us now consider the initial value problem (IVP) for the above system.Continuous dependence of the solution u at all times t on the initial data f can be understood as being able to provide an estimate of the form with real constants K ≥ 1 and α ∈ R, and || • || denoting a suitable norm.The degree of hyperbolicity of the PDE system determines well-posedness of the problem, with the IVP for SH systems being (strongly) well posed in the L 2 norm of the system.On the contrary, the IVP of a system that is only WH but not SH, is ill posed in the L 2 norm.However, under certain conditions it may be weakly well posed in a different norm that is not equivalent to L 2 .More details on weak well-posedness can be found in [50,51].The essential difference between strong and weak well-posedness is that the first, in contrast to the second, is maintained in the presence of lower order perturbations (such as source terms).
The hyperbolicity analysis of free evolution schemes used in numerical relativity can be performed by bringing them to the form (1), via linearization and the constant coefficient approach.In [33] such an analysis was presented for the Bondi-Sachs free evolution system.Existence and uniqueness of solutions to the Bondi-like CIBVP were shown in [52,53].In [54], a specific subsystem of the full system was found symmetric hyperbolic, and consequently certain norm estimates were provided.However, the hyperbolicity analyses presented in [33,34] which consider the full PDE system (in the linear, constant coefficient approximation), show that it is only weakly hyperbolic.This result is in contrast to that of [54] and implies that a Bondi-like CIBVP is ill posed in L 2 (in the linear, constant coefficient approximation).It is therefore remarkable that several characteristic codes can perform long simulations and produce physically sensible results.This fact, in combination with the existence of symmetric hyperbolic formulations of GR that involve third order metric derivatives and use Bondi-like gauges [55][56][57][58], suggests that there may still be some norm that is not equivalent to L 2 , which can be used to demonstrate well-posedness of the standard Bondi-like CI(B)VP with energy estimates.By standard here we mean that GR is formulated only via the EFE so that only up to second order metric derivatives are involved, rather than including the Bianchi identities as in [55][56][57][58].However, finding such a norm and showing even weak well-posedness for these second order Bondilike CI(B)VPs is still an open question.
Here, we concern ourselves with the following question: what are the best error estimates we can obtain for CCE and CCM with present setups and how can we test them numerically?To address this question we focus on simple toy models that have a structure similar to that of GR in the gauges typically employed.More specifically, we employ a characteristic model that is WH in a similar way to free evolution Bondi-like schemes, whereas the Cauchy model is SYMH like for instance GR in harmonic gauge.For completeness and comparison, we also consider cases with SYMH CIBVPs and WH IBVPs.The WH system is such that its IVP is weakly well posed.The paper has the following structure.In Sec.II we introduce the models we use for our analysis and provide energy estimates for each one of them.We also discuss under which conditions we can obtain an energy estimate for CCE and CCM.In Sec.III we discuss the discretization of the models, and perform numerical experiments to address convergence of solutions to discrete approximations of the previously discussed norms.We summarize our findings in Sec.IV and present our conclusions in Sec.V.

II. THE MODELS
In this section we construct toy models and provide energy estimates for their solutions.For simplicity, we focus on first order, linear, constant coefficient, homogeneous and real-valued PDEs.We study energy estimates both for SYMH and WH models for the IBVP and the CIBVP, as well as for their CCE and CCM.The structure of the WH models are motivated by popular Bondi-like formulations of GR, and follows the model used in [33].Working with SYMH and not just SH systems allows us to mimic the best case scenario in a current GR CCM implementation, matching between a SYMH and a WH system for the IBVP and CIBVP respectively.Furthermore, estimates including boundary terms are simpler to obtain for SYMH systems than for systems that are only SH.The estimates presented in the section form the basis of our numerical tests presented in Sec.III.
A. The IBVP Our Cauchy PDE model system is with ϕ 1 , ψ v1 the right-moving fields, and ψ 1 the leftmoving.This system is SYMH when the boxed term of Eq. (3a) is included and only WH otherwise.To obtain an energy estimate for the solution to the IBVP of the SYMH system (3) we follow the standard procedure.Let us start with the L 2 norm After acting on it with ∂ t , using the right-hand side (rhs) of Eq. ( 3) to replace ∂ t u 1 and collecting total derivatives, we obtain where the total ∂ z terms vanish since the data are assumed to be periodic in the z direction throughout.Evaluating the total ∂ ρ terms and integrating in t ∈ [t i , t f ] we obtain where we have rearranged terms such that the left-hand side (lhs) includes only the solution and the rhs only given data.This shows that the solution is completely controlled by the given data in the norm described by the lhs of Eq. ( 5).
Let us next consider the WH version of (3).First, let us note that we fail to obtain an energy estimate in the L 2 norm for the WH case, due to the boxed term in Eq. (3a).Instead, motivated by the analysis presented in Sec.III.B of [33], we start with the "lopsided" or q norm where q is due to the adoption of the notation of [50] for weak well-posedness and "lopsided" refers to the fact that only the angular derivative of a certain variable is included in the norm, and not of all variables.This norm leads to where now, due to the WH structure of the system the term 2ψ v1 ∂ z ϕ 1 does not form a total derivative.The lopsided modification of the original L 2 norm was performed exactly in a way that allows us to control this term.To see this, let us first integrate in t and drop the total derivative term ∂ z ψ 2 1 to obtain 2 the latter reads as which after using Grönwall's inequality yields where we have rearranged terms such that the lhs includes only the solution, and the rhs only the given data, showing that the solution is completely controlled by the given data in the q norm (6).It is important to remember that this result is subject to the structure of the source terms, and that for generic source terms one cannot provide such an energy estimate.This is the essential difference between strong and weak well-posedness.Since the q norm involves a specific derivative in its integrand, we find it useful to employ a norm for which one can show energy estimates for the SYMH IBVP, which also contains derivatives.We denote by H 1 the following norm: Let us consider the first order system that is formed after we enlarge the SYMH system (3), by including ∂ ρ u 1 and ∂ z u 1 as variables.This system has the same principal structure as the original SYMH (3), and an energy estimate for its IBVP can be obtained starting with the H 1 norm (8).Following the same steps that led to the estimate (5) in the L 2 norm, one can show that where again the lhs involves only the solution, and the rhs the given data.We can use this result to assess the well-posedness of the IBVP of the SYMH system (3), in a norm that includes derivatives.Notice that in comparison to the q norm (7), the H 1 norm (9) contains the same derivatives on all variables and so is symmetric, rather than lopsided.

B. The CIBVP
The characteristic model we consider is where the fields ϕ 2 , ψ v2 are right moving, and propagate with speed unity and ψ 2 is left moving with speed 1/2.The system has two equations with derivatives only intrinsic to the characteristic surfaces N u and one that includes also a derivative transverse to them.The system is SYMH when the boxed term in Eq. (10a) is included, and only WH otherwise.If we consider the coordinate transformation the SYMH version of the system (10) is the characteristic version of the Cauchy system (3), provided that and similarly for the WH one.We consider the CIBVP for the characteristic domain D 2 shown in Fig. 3, that consists of u ∈ [u i , u f ], x ∈ [x min , x max ], and z ∈ [0, 2π) with z a periodic direction.The given data are the left-moving field ψ 2 (u i , x, z) on an initial characteristic N ui , its value ψ 2 (u, x max , z) on the right boundary T xmax and the values of the rightmoving fields ϕ 2 (u, x min , z), ψ v2 (u, x min , z) on the left boundary T xmin .The solution for which we next provide energy estimates consists of ϕ 2 (u f , x, z) on a surface N u f with u f > u i , ϕ 2 (u, x min , z) on the left boundary T xmin and ϕ 2 (u, x max , z), ψ v2 (u, x max , z) on the right boundary T xmax .
To obtain energy estimates in the characteristic domain, we can treat independently the left-and rightmoving fields, due to the fact that the right-moving fields do not have an explicit ∂ u time derivative.Let us first consider the left-moving field and start from the energy After acting with ∂ u and replacing with the rhs of Eq. (10c) we obtain By evaluating the total derivatives and integrating in u ∈ [u i , u f ] we arrive at where we have rearranged terms such that the lhs contains only the solution and the rhs the given data, as well as multiplied everything by a factor of 2, to cancel the characteristic speed in the worldtube integrals.This estimate allows us to completely control the solution for the left-moving data by given data.If there are source terms in Eq. (10c) that include the right-moving fields, then the above calculation does not work.In that case one needs to perform the calculation for the whole system at once.Potentially, one also needs to consider a slightly different computational domain such that there is an integrand involving both left-and right-moving variables and offer the chance to control them using Grönwall's inequality, as before.This however is beyond the scope of the current paper.Next, we consider the right-moving fields, focusing first on the SYMH case.Let us start by considering the energy By acting with ∂ x and replacing with the rhs of Eq. ( 10a), (10b) we obtain where the integrand on the right is a total ∂ z term, and hence vanishes by periodicity in z.Then, integrating in x ∈ [x min , x ′ ] for some x ′ ∈ (x min , x max ], we arrive at the estimate where the lhs is the solution and the rhs the given data.
As was also done in a similar calculation in [33], we can take the supremum of ||u 2 || 2 right(T x ′ ) within (x min , x max ] to obtain a better estimate, since ||u 2 || 2 right(T x ′ ) is not necessarily monotonically increasing with increasing x ′ .Then, since Eq. ( 12) is true for all x ′ ∈ (x min , x max ] we obtain If we consider the WH model ( 10), then we should consider the energy in analogy to the q norm (6).Following the same steps as before we can arrive at where we have also used that 2ψ 2 .Using again the Grönwall inequality and noticing that ||u 2 || 2 q−right(Tx min ) involves only given (boundary) data, and hence is a constant function of x or x f (so nondecreasing), we obtain where we again took the supremum of x ′ to obtain a better estimate and the lhs contains only the solution and the rhs the given data and multiplied overall with e xi−x f for later convenience in the continuum analysis.We can then proceed to add the estimates ( 11) and ( 12) to obtain the estimate for the homogeneous, symmetric hyperbolic system (10) or add (11) and ( 14) to obtain an estimate for the WH characteristic model Subsequently, we drop the factor of 2 in the lhs of Eq. ( 15) and ( 16) for convenience, but we still maintain the hypersurface integrals which provide the control of the solution on that region.
In analogy to the H 1 norm (8) for the SYMH IBVP, we also consider estimates in an H 1 norm adapted to the CIBVP, again by enlarging the SYMH characteristic system (10) with ∂ x u 2 and ∂ z u 2 .In the characteristic H 1 norm one can show the following estimate for the enlarged system with and by simple commutation of partial derivatives.

C. CCE and CCM
Using the energy estimates from the previous subsections, we can now discuss energy estimates for CCE and CCM.For CCE the estimate follows straightforwardly from the previous results, and if both the IBVP and CIBVP are individually well posed in some norm, then CCE is too.The only detail that needs attention for this conclusion to be drawn is to choose the boundary data on T 0 such that they are controlled in the appropriate norm for a WH CIBVP.Since these data are solutions of the IBVP, this may amount to having enough control over the derivatives of the solution on T 0 , so as to be controlled not only in the L 2 but also in the H 1 norm, as for instance shown in estimate (9).However, in this case we "lose" derivatives, in the sense that we control derivatives of more variables for the boundary data than for the solution to the WH CIBVP.
For the composite CCM problem to be well posed the situation is more complicated, as the norms in which the IBVP and the CIBVP are well posed need to be in some sense compatible.As shown in [33] the reason is that otherwise, there is a remaining integral over the interface T 0 as shown in Fig. 4. For completeness, we briefly repeat here this calculation.First we consider matching the two SYMH setups.By adding the estimates ( 5) and ( 15) and using that + The latter is the energy estimate for the composite CCM problem when matching the IBVP and CIBVP of the homogeneous symmetric hyperbolic Cauchy and characteristic systems (3) and (10), respectively.Notice that, as expected for the composite CCM problem, there is no worldtube integral over T 0 .A similar estimate can be obtained using the H 1 norms and adding the higher derivative estimates ( 9) and ( 17) for the SYMH IBVP and CIBVP, respectively.Let us next consider matching the IBVP of the SYMH Cauchy system (3) to the weakly well-posed CIBVP of the WH characteristic system (10).In this case we add the energy estimates ( 5) and ( 16), which yields where we have used that on T 0 , u 1 = u 2 in our models.The boxed term in the rhs of the latter is what prevents us from obtaining an energy estimate for the composite CCM in this setup and appears purely due to the incompatibility of the norms in which the IBVP and the CIBVP are well posed and weakly well posed, respectively.Since this term is not part of the given data, then the solution is not completely controlled by the intial and boundary data, and so the above is not a valid energy estimate for this composite CCM.If instead of using the L 2 estimate (5) for the IBVP, we use the H 1 estimate (9), the result would be similar, in that nonvanishing terms including an integral over the interface T 0 remain.Finally, we consider the matching of two WH setups.In this case, we combine the q norm estimates ( 7) and ( 16) for the IBVP and CIBVP, respectively, to obtain In this case we see that, due to the compatibility of the norms, the integral over T 0 vanishes and the lhs that contains only the solution, is completely controlled by the rhs, that includes only given data.

III. NUMERICAL EXPERIMENTS
We present numerical convergence tests for the individual IBVP and CIBVP, as well as for CCE and CCM, using discrete versions of the models analyzed in Sec.II.We work with the homogeneous models, unless explicitly stated otherwise.The implementation was performed using the Julia programming language [59] with the following details: 1.The numerical domain is defined by where u = t − ρ, so for ρ = 0, t f = u f .We denote the grid points for x as x i with i ∈ [1, . . ., N ] and the grid spacing as h x ≡ x 2 − x 1 , and so forth.
2. All the derivatives that appear on the rhs of the systems (3) and ( 10) are approximated by second order accurate finite difference operators.The derivative ∂ x is approximated with a centered stencil for the internal points of the x grid.The derivative on the first and last points is approximated with upwind stencils that match the truncation error of the centered one The same operators are used to approximate ∂ ρ , with h x → h ρ , whereas for ∂ z , the centered stencil is used for the first and last points as well, taking advantage of the periodicity in z.
3. All fields that involve explicit time derivatives, that is ϕ 1 , ψ v1 , ψ 1 , ψ 2 , are integrated in time using the explicit fourth order Runge-Kutta integration scheme.
4. The fields ϕ 2 , ψ v2 that satisfy equations intrinsic to the outgoing null hypersurfaces are integrated for every timestep from x = 0 to x = 1 using the trapezoidal rule, which is second order accurate.
10.For CCE, boundary data for the IBVP are given also on T ρ=0 for the left-moving field ϕ 1 .The values of the right-moving fields ϕ 1 , ψ v1 are propagated from the Cauchy to the characteristic domain via T ρ=x=0 , as in CCM.
12. For the tests presented in the paper, all the numerical boundary data are provided with pure injection, that is if g 1 (t, z) is the value of the grid function ϕ 1 for ρ = −1, then ϕ 1 (t, −1, z) = g 1 (t, z).We have also experimented with a different method for numerical boundary data, namely by providing the time derivative rather than the field itself.The results can be found in [60] and are qualitatively the same.
The main focus of this work is to examine the consequence of (a lack of) continuous dependence at the continuum level, as in (2), under numerical approximation.Detecting lack of convergence for WH systems can be subtle [61].Numerical convergence tests that include high frequency data are a good choice in order to achieve this [33,62,63].Low frequency data are not ideal to detect the effect of weak hyperbolicity in numerical simulations, since the solution may exhibit frequency dependent exponential or polynomial growth, that may not be evident for small frequencies in finite simulation time.Thus, the loss of convergence for an ill-posed problem may not be evident if the given data are not properly chosen.We run robust stability tests, which use high frequency given data, and monitor the convergence of the obtained numerical solutions in the norms presented in Sec.II.
Since we use uniform grids, we choose the high frequency data to be random noise of a certain amplitude.The random noise is added on top of an exact solution to the PDE problems, which in our tests is zero.We call these exact convergence tests, since we know the exact solution to the PDE problem.To fix ideas, next we present explicitly only the field ϕ 1 , but the discussion is the same for all fields ψ v1 , ψ 1 , ϕ 2 , ψ v2 , ψ 2 .If we denote by ϕ 1 the continuum solution and ϕ 1h the numerical at resolution h, and assume that the dominant error for our discretization is due to the second order finite difference operators, then we expect the following relation to be true: Given ϕ 1 = 0 in our setup, it follows that We denote by h c and h m the coarse and medium resolutions for which we solve the same PDE problem, and construct the convergence factor where ≃ here denotes equality up to terms of order O(h c ).Notice that the above relation is true only when the exact solution vanishes.Generically, when this is not true, one needs to replace ϕ 1c → ϕ 1 − ϕ 1c in the expression, and similarly for the medium resolution.Every time we increase resolution, we halve the grid spacing (for instance h m = h c /2), and therefore the expected convergence factor is Q = 4.We construct and monitor the following quantity during our simulations where || • || h is the discrete approximation to the continuum norm || • ||.The forms of the state vector u h and the norm || • || h depend on the specific discretized PDE problem solved, and so are clarified in the following subsections.Given that Q = 4, perfect second order convergence corresponds to C exact = 2 for our exact convergence tests.Motivated by the notation of the inequality (2), let us denote as f h the numerical initial and boundary data, which we specify as random noise of amplitude A h centered around zero.It is important to choose this amplitude appropriately such that C exact = 2 for all given data.Since f h ∼ A h , for a field that appears without any derivative in the norm under consideration, we find that when A hm = A hc /4, If however the norm includes a derivative of the field, then due to the finite difference operators used the exact convergence rate behaves as which for our setup is equal to 2 when we choose A hm = A hc /8.Consequently, our recipe to obtain given data that have C exact = 2 is the following: every time we double resolution, we drop the noise amplitude of the given data by a factor of 4 for any field with no derivative in the considered norm, and 8 for any field with at least one derivative in the considered norm.
We run tests with given data that exhibit C exact = 2 in the appropriate L 2 , q, or H 1 norm for the setup under consideration, and monitor the C exact of the solution in its L 2 , q, or H 1 norm.We consider the test as passed if the solution converges in the same norm as the given data, at the same order; that is its C exact tends to 2 with increasing resolution.In any other case we consider the test as failed.In our numerical experiments we see three different ways in which a test can fail: 1. C exact of the solution tends to a value different from 2 with increasing resolution, in the same norm as the given data.
2. C exact of the solution tends to 2 with increasing resolution, but in a norm different than the given data.
3. C exact does not tend to any fixed value with increasing resolution.
The first and second failing scenarios can actually happen for the same solution, if we just monitor its convergence in a different norm.
One argument to support our working definition of a passed and failed test is the following.Consider inequality (2) at the semidiscrete level, for instance at resolution h, namely with K, α real constants and K ≥ 1.For a passed test we have that where ≃ here denotes equality up to terms of order O(h), and therefore If this exact convergence rate is maintained for higher resolution runs, then we can imagine taking the limit of infinite resolution and recovering inequality (2) at the continuum.If, in contrast, during the tests we see the first type of failure, so that , with c some real positive constant, then after doubling the resolution n times we obtain In either case one does not recover the continuum version of inequality (2), which we interpret as failure of the solution to be controlled by the given data.In this argument we consider the theoretically expected value of C exact for both the solution and the given data, but in practice, we look for a C exact that tends to that value.
The second type of failure is an interesting scenario realized in our tests, where we obtain the same convergence rate for the solution and the given data, but only in different norms, where we have dropped the subscript ∼ h in the different discrete norms.In this case, one can at most obtain the following inequality: by considering again the limit of infinite resolution, This however is not the same as inequality (2), which is involved in the textbook definition of well-posedness.More importantly, since ||u|| q is smaller than ||u|| H1 (controls fewer derivative terms), it is perfectly possible to obtain a solution u with unbounded H 1 , but bounded q norm.In this case, the blowup would appear in derivatives of the solution that are included in H 1 but not in the q norm.We therefore consider this scenario also a failure.
The lowest resolution we use for our tests is denoted by h 0 and has N ρ = N x = 17 and N z = 16.The grid points and grid spacing for the different resolutions are set according to and h D labels the different resolutions, for instance D = 0 corresponds to the lowest resolution h 0 .We perform tests for D = 0, 1, 2, 3, 4, and compute C exact only for the timesteps that are common across all resolutions, in other words at each timestep with h 0 grid spacing.The code and parameter files for the convergence tests can be found in [60] and the data used to produce the following convergence plots in [64].In the following subsections, all the norms should be understood as evaluated on the numerical grid.Table II contains the summary of our numerical experiments.

A. IBVP
First we examine numerical convergence of the Cauchy setup.The discrete norms in which we test the conver-gence of the solution are In our tests, we choose given data that converge in either the L 2 , q, or H 1 norms.Based on the energy estimates given Sec.II A, these can be obtained by considering Eqs. ( 24)-( 26), but where the spacelike hypersurface sums are evaluated only for the initial data, and the worldtube sums on the left boundary for right-moving fields and on the right boundary for left-moving fields, i.e., the opposite of what is shown.We perform four different tests: 1. L 2 − L 2 test: convergence of the solution in the L 2 norm (24), for L 2 given data.
2. q − q test: convergence of the solution in the q norm (25), for q given data.
3. H 1 − H 1 test: convergence of the solution in the H 1 norm (26), for H 1 given data.
4. H 1 − q test: convergence of the solution in the q norm (25), for H 1 given data.
The L 2 − L 2 , q − q, and H 1 − H 1 tests examine, at the discrete level, continuous dependence of the solution on the given data as described by inequality (2).In contrast, the H 1 − q test examines the discrete version of inequality (23).If the latter is satisfied, it shows that a solution is controlled by given data that converge in a bigger norm.In our particular setup, the q norm (25) is smaller than the H 1 (26) in the sense that it allows control over fewer derivative terms.It then becomes clear that a blowup in a derivative of the solution can still happen, but fail to be detected by the q norm.We consider this test because it can explain why CCE involving a SYMH IBVP and a WH CIBVP can result in well-posed PDE problems that exhibit the expected numerical convergence.But the H 1 − q test does not address continuous dependence as understood in the standard sense for well-posedness as in inequality (2).The findings of our numerical tests are shown in Fig. 5.The top row of the figure shows C exact for a SYMH IBVP setup, whereas the bottom is for a WH one.Each column refers to one test, starting from the left and moving to the right with the L 2 − L 2 , q − q, H 1 − H 1 , and finally the H 1 − q test.The first three provide the expected convergence results, that is, for the SYMH case, good empirical convergence of the solution in L 2 and H 1 norms and loss of convergence in the q norm.In contrast in the WH case, we see C exact = 2 only for the q − q test.The H 1 − q test returns interesting and potentially confusing results: both solutions to the SYMH and WH setups exhibit good convergence in the q norm for given data controlled in their H 1 norm.This can be understood in the following terms: 1. SYMH: by monitoring the q norm, we take into account only the ∂ z ϕ 1 derivative of the solution.This part still converges appropriately, and due to its relative size in comparison to the nonderivative terms, dominates the convergence rate.This does not mean that the H 1 − q test provides numerical evidence for well-posedness of the SYMH IBVP in the q norm, which would be the case if the q − q test was passed.
2. WH: when we monitor the q norm, we ignore all the derivatives that appear in H 1 and not in q, which are responsible for the loss of convergence.We only monitor the ∂ z ϕ 1 term, which converges properly and due to its relative size, dominates the nonderivative terms of the q norm.Since we are considering the homogeneous case, the loss of convergence of the nonmonitored derivative terms does not propagate to ∂ z ϕ 1 due to the lack of a coupling among them.Finally, proper convergence of the given data in the H 1 norm implies the same in the q norm, since H 1 is bigger than q.
When we see loss of convergence in these tests we observe C exact = 1, and since this is not the expected rate of convergence for our numerical scheme, we consider these tests as failed.However, it is often the case in physically interesting simulations that the convergence order of the solution is lower than that of the numerical scheme used (for instance shockwave capturing or schemes that include interpolation).Here, we insist in considering these tests as failed, because we do not have additional elements in the scheme that can drop the order of convergence (the data are smooth plus the artificially exaggerated numerical error, no interpolations) and more importantly we try to make contact with the continuum energy estimates.As a final remark on this point, we should mention that it is not clear why C exact = 1 in these cases.A possibility is that it relates to the bulk integrals in the continuum analysis that prevent us from obtaining energy estimates for the scenarios of the failed tests.Obtaining however a concrete answer requires more testing and is beyond the scope of this work.

Synopsis
The SYMH IBVP exhibits second order convergence in the L 2 and H 1 norms for given data controlled in the same norms, and second order convergence in the q norm only when the given data are controlled in the H 1 norm.The WH IBVP shows second order convergence in the q norm for given data controlled in either their q and H 1 norms.These results are inline with the expectations from theory.

B. CIBVP
We repeat the tests for the CIBVP, monitoring the convergence of the solution in the following discrete norms   The left column presents the results for the L 2 − L 2 test, the second to left for the q − q test, and so on.The results of the L 2 − L 2 , q − q, and H1 − H1 tests are as expected by theory, that is convergence in the L 2 and H1 for the SYMH case and in the q norm for the WH one.The convergence in the q norm seen in the H1 − q test is merely a result of the fact that H1 (26) is a bigger norm than q (25), in the sense that it controls more derivative terms.We only monitor the ∂zϕ1 derivative term, which still converges appropriately for the SYMH case and due to its size in comparison to the nonderivate terms, dominates the q norm.We also see good convergence of the solution to the WH IBVP in the H1 − q test, because we ignore the derivative terms that do not exhibit convergence.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs. and

Synopsis
The results are qualitatively the same as for the IBVP, and are summarized in Fig. 6.We again find the expected convergence rate in L 2 and H 1 for the SYMH and in the q norm for the WH setup.The H 1 − q test provides apparent convergence in the q norm for both setups, which is understood exactly in the same way as in the IBVP.

C. CCE
The conclusions drawn from the previous tests for the IBVP and CIBVP shape our expectations for CCE.In this setup, convergence of the numerical solution is still examined for the individual IBVP and CIBVP parts of the CCE scheme.The important difference now is that part of the IBVP solution serves as boundary data (for the right-moving fields) for the CIBVP.When performing these tests for CCE we expect the following results: 1. CCE between a SYMH IBVP and a SYMH CIBVP: both IBVP and CIBVP solutions converge in L 2 and H 1 , but not in the q norm.The H 1 − q test has the same qualitative results as in Secs.III A and III B.

WH-WH CCE:
Both IBVP and CIBVP converge only in the q norm, but not in L 2 or H 1 .The H 1 −q test shows again that both solutions converge in their q norms for the H 1 given data.

SYMH-WH CCE:
The solution to the IBVP exhibits the same behavior as in Sec.III A. The solution to the CIBVP can converge properly only in the q norm.The boundary data to the CIBVP are part of the solution to the IBVP and so we need a solution to the IBVP that converges properly in the H 1 norm, such that the boundary data to the CIBVP also converge properly in the q norm (neglecting the derivatives that are in H 1 but not in q).The aforementioned setup is understood in the H 1 −q test.In addition, choosing characteristic initial data that converge in L 2 or H 1 does not affect the convergence of the CIBVP solution in the q norm.The reason is that the derivative term ∂ z ϕ 2 dominates the characteristic q norm (28).For our   Each column shows the results of each test.The L 2 − L 2 , q − q, and H1 − H1 tests exhibit the expected convergence in the L 2 and H1 for SYMH and in the q norm for WH.Convergence in the q norm seen in the H1 − q test is again a result of the fact that H1 (29) is bigger than the q (28) norm, as in the IBVP.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs.
tests we drop the amplitude of ϕ 2 (u = 0, x, z) by a factor of 8 every time the grid spacing is halved, that is compatible with H 1 .7. Cexact for all our tests, for the CIBVP part of a CCE setup between a SYMH IBVP and a WH CIBVP.As expected, only the H1 q test provides good convergence.In this case, the H1 given data provide an H1 controlled solution to the IBVP, part of which serves as boundary data for the WH CIBVP.That control is more than enough to have q controlled given data for the CIBVP, and hence, the solution to the CIBVP is controlled in the q norm.
Repeating these tests for CCE, our expectations were verified.The SYMH-WH scenario is the relevant one for current implementations of CCE in GR simulations, and therefore we present only these results in Fig. 7.We see that indeed the solution to the WH CIBVP setup exhibits the appropriate convergence rate only when H 1 data are given for the IBVP.The solution to that IBVP provides boundary data for the CIBVP that have H 1 control, which is more than the necessary control in the q norm.This result is compatible with that of the H 1 − q test in Sec.III B.

Synopsis
CCE between a SYMH IBVP and a weakly well-posed WH CIBVP, can only provide a characteristic solution with second order convergence in its q norm, if the solution of the SYMH IBVP is controlled in its H 1 norm.This is inline with our expectations from theory.

D. CCM
The discrete norms for which we run our tests for CCM are    The first column refers to the L 2 − L 2 test, the second to the q − q test, and so on.The L 2 − L 2 , q − q, and H1 − H1 tests show the expected convergence in the L 2 and H1 for SYMH-SYMH and in the q norm for WH-WH.The convergence in the q norm seen in the H1 − q test is again a result of the fact that H1 (32) is bigger than q (31), as in the IBVP.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs.
First, let use consider matching PDEs with the same degree of hyperbolicity, that is SYMH to SYMH or WH to WH.The convergence results of these setups are shown in Fig. 8, with SYMH-SYMH at the top and WH-WH the bottom row.As expected by the energy estimates of Sec.II C, we see good convergence for SYMH-SYMH in the L 2 − L 2 and H 1 − H 1 tests that is in norms L 2 (30) and H 1 (32) respectively, and in the q − q test for the WH-WH case, i.e. norm q (31).In the H 1 − q test we see convergence of the solution in the q norm (31) for H 1 given data, in line with the results of the same test for the IBVP and CIBVP.The explanation of this phenomenon is the same as earlier, that is for the SYMH-SYMH setup we see the convergence of the ∂ z ϕ 2 term that dominates the norm, whereas in the WH-WH setup we ignore the rest of the derivative terms that do not converge.
FIG. 9. Cexact for a SYMH-WH CCM setup.The L 2 − L 2 test is shown in top left, the q − q in top right, the H1 − H1 in bottom left, and the H1 − q test in bottom right.The first three tests fail to exhibit the appropriate convergence rate, as expected for this setup.The H1−q test manifests convergence in the q norm (31), again due to the fact that H1 (32) is a bigger norm.The fact that the behavior seen is similar to that of the CIBVP part of CCE, in Fig. 7, is a result of treating the homogeneous case.In Fig. 10 we see that the behavior is different for an inhomogeneous SYMH-WH CCM setup.
Next we consider the matching of a SYMH IBVP to a WH CIBVP, the actual case for CCM in GR as currently performed.The convergence plots are presented in Fig. 9.We see that the solution does not exhibit good convergence in any norm for the L 2 − L 2 , q − q, and H 1 − H 1 tests, which are the ones that reflect numerically continuous dependence on given data.Only the H 1 − q test shows good convergence of the solution in the q norm for H 1 given data.This test only provides numerical support for inequality (23) rather than (2).Comparing Fig. 7 with Fig. 9, we see the same qualitative behavior, namely good convergence only for the H 1 − q test.The latter refers to CCM and the former the CIBVP part of CCE, both for a SYMH IBVP and a WH CIBVP.The difference between CCM and CCE is only the additional propagation of information between the characteristic and Cauchy domain for the left-moving field for CCM.Since we are considering only the homogeneous case here, we see that there is no structure in the composite system that couples the right-moving variables of the characteristic domain that are responsible for weak hyperbolicity, to the left-moving one.Consequently, the qualitative behavior in the convergence tests for CCM and CCE is very similar.The important difference between CCE and CCM is that the H 1 − q test can actually provide good boundary data for the WH CIBVP in the CCE case, because then the IBVP and CIBVP are still treated individually regarding well-posedness.In contrast, for CCM the whole problem has to be treated simultaneously and so the H 1 − q test is not an indicator for well-posedness in this case.To clarify the difference between CCE and CCM further, we consider next the inhomogeneous setup (33) and (34).

An inhomogeneous example
We now consider an inhomogeneous IBVP and CIBVP, in which sources are added in such a way that rightand left-moving variables are coupled when CCM is performed.The inhomogeneous IBVP (labeled B1) is: where the boxed term controls the hyperbolicity of the system as before.The inhomogeneous CIBVP (labeled B2) is When the boxed term is included the systems are SYMH, and when omitted they are only WH.In comparison to the earlier homogeneous case, we made an addition of source terms that is minimal in the following sense: in the characteristic setup, the left-moving part is affected by the right-moving one through the sources, whereas in the Cauchy part, the opposite is true, that is the rightmoving part is affected by the left-moving one through the sources.Notice that the characteristic setup has still a nested structure as is desirable in a model for GR.The important difference in comparison to the homogeneous SYMH-WH CCM case is that now the right-moving characteristic variables that form a nontrivial Jordan block in the angular principal part, can affect the left-moving characteristic variable, which then influences the solution to the SYMH IBVP.The latter eventually provides inappropriate boundary data for the right-moving variables of the CIBVP, enhancing the effect of weak hyperbolicity.This coupling provides the aforementioned interaction loop between the left-and right-moving variables when CCM is performed, but not for CCE, since then boundary data for the IBVP are not given by the solution to the CIBVP.
In Fig. 10 we present all our tests for CCM performed with the SYMH B1 IBVP and the WH B2 CIBVP.In the top row we see the exact convergence rate C exact , and at the bottom the respective norms that are included in the calculation of C exact .In this case we do not see convergence in any of the tests, which can be explained by the fact that the respective norms increase in time, with a growth rate that increases with resolution.Therefore, the ratio between two norms calculated at different resolutions cannot be the desired constant in time.Notice that the way these tests fail is rather different from the way some earlier tests failed.Here we see complete loss of convergence, which is often called an instability-even though the simulation does not stop abruptly-whereas earlier we observed convergence at an order lower than the theoretically expected, or in a smaller norm than that of the given data.
For comparison, we perform the H 1 − H 1 test for CCM between SYMH B1 and SYMH B2, the q−q test for CCM with WH B1-WH B2, as well as the H 1 − q test for CCE for SYMH B1-WH B2.We chose the specific tests and combinations of PDEs, because these provided good convergence results earlier when monitoring derivative terms in the norms.The results are shown in Fig. 11.More specifically, CCM between the inhomogeneous SYMH setups exhibits the expected convergence in the H 1 norm, whereas for WH setups we see good convergence in the q norm only for a brief period.The latter should not be surprising, since standard numerical methods as the ones utilized here are developed having in mind PDE problems that are strongly well posed, such as strongly or symmetric hyperbolic.These results indicate that the main source of this rapid loss of convergence is the different degree of hyperbolicity between the IBVP and CIBVP in CCM.In addition, we see good convergence in q norm for the CIBVP part of CCE with SYMH B1 IBVP and WH B2 CIBVP, in line with earlier CCE results.Based   10, we present here the following (from left to right): the H1 − H1 test for SYMH B1-SYMH B2, the q − q test for WH B1-WH B2, and the H1 − q test for the CCE setup SYMH B1-WH B2 (plotting the CIBVP norms).
on these findings, we deduce that having the influence of the bad part of the WH CIBVP solution affecting the SYMH IBVP is exactly what causes the severe loss of convergence and the rapid growth of the respective norms.When this influence is prevented by doing only CCE, this exponential growth of the norms is not present.For completeness, we have performed convergence tests for CCM with SYMH B1 IBVP and WH B2 CIBVP, with smooth given data.In complete contrast to the random noise data, the smooth data tests exhibit convergence in the L 2 norm.Note that these are not exact, but self convergence tests, since we do not know the exact solution of the tested smooth setup.

Synopsis
CCM provides the expected second order convergence for the solution of the composite problem, when it is performed between an IBVP and CIBVP with the same de-gree of hyperbolicity, that is SYMH or WH.In the case of SYMH IBVP matched with WH CIBVP, we considered two cases.In the homogeneous case, most convergence tests failed with the convergence rate reduced from second to first order; whereas in the inhomogeneous case all tests failed with the solution growing at an exponential rate that increases with increased resolution.

IV. SUMMARY AND DISCUSSION
Modeling gravitational waves at future null infinity accurately is a challenge that has to be met to obtain waveforms of high fidelity.One possibility is to successfully develop a CCM scheme for GR.Such a scheme would combine the ability to capture in detail the physics of the strong gravity regime in the Cauchy setup, with the accurate propagation of gravitational waves due to the foliation with null hypersurfaces of the characteristic setup.CCM does not suffer from artificial boundary conditions on the outer boundary of the Cauchy domain, as for instance a CCE scheme.The characteristic setup is not the only way to obtain accurate gravitational waveform models at future null infinity.Alternative strategies include hyperboloidal compactification [65][66][67][68][69][70][71][72] and the use of the conformal EFEs [73][74][75].
In this paper we focused on mathematical properties of CCE and CCM schemes and their impact on numerical simulations.Our motivation comes from recent progress in understanding the mathematical structure of Einstein's field equations in the characteristic setup, in which the gauge freedom of GR is typically fixed by choosing Bondi-like coordinates which are built upon outgoing (in the case of CCE and CCM) null geodesics.In [33], certain characteristic PDE systems of GR were shown to be only weakly hyperbolic in some Bondi-like coordinates, and subsequently in [34] this structure was attributed to the choice of a Bondi-like gauge.Based on these results, it becomes clear that current CCE and CCM schemes of GR involve a strongly or symmetric hyperbolic system for the IBVP and the WH Bondi-like one in the CIBVP.Since PDE problems based on WH systems are generically not well posed, it is natural to ask whether these composite CCE and CCM setups are well posed, and to what extent should a solution based on them be trusted in making accurate predictions for gravitational waveforms.
To address this question we constructed toy models and utilized them as models for GR CCE and CCM.To address the effect of weak hyperbolicity we studied weakly and symmetric hyperbolic models for both the IBVP and the CIBVP, comparing different combinations of IBVP and CIBVP in CCE and CCM.The weakly hyperbolic structure of the CIBVP mimics that of the EFE in a Bondi-like gauge.Importantly, the WH system is constructed in such a way that it has a weakly well-posed IBVP and CIBVP.This assumption is the best possible scenario for the Bondi-like CIBVP of the EFE, involving up to second order metric derivatives.Given that there are symmetric hyperbolic characteristic formulations of GR in a Bondi-like gauge that include third order metric derivatives, one might hope to show weak well-posedness for the characteristic setups commonly used in numerical relativity.To the best of our knowledge, such a result is yet to be obtained.
For our toy models, we first studied well-posedness at the continuum level, where we considered the homogeneous case and provided energy estimates for the IBVP, CIBVP, CCE and CCM.We discussed energy estimates in norms with (q and H 1 ) and without (L 2 ) derivative terms, and adapted to the IBVP, CIBVP, CCE and CCM setups.A similar analysis was performed in [33] only for the L 2 norm.As in that earlier work, we found that CCM cannot provide a well-posed PDE problem when matching a WH with a SYMH PDE problem.On the contrary, CCE can be weakly well posed if the WH problem is weakly well posed.The essential difference is in the communication of the solution between the Cauchy and the characteristic domains.In terms of the energy estimates, this is translated into a worldtube integral over the interface between the Cauchy and the characteristic domains, that is not controlled by the given data.
Based on the continuum analysis, we performed several numerical convergence tests, a second order accurate implementation.Our main goal was to address continuous dependence of the solution on the given data, in discrete approximations to the L 2 , q, and H 1 norms previously constructed.Our first step was to verify the expected norm convergence for the individual IBVP and CIBVP.We saw the expected convergence for the SYMH setups in the L 2 and H 1 norms and for the WH in the q norm, for both the IBVP and CIBVP.An interesting scenario was realized when the given data were controlled in the H 1 norm and the convergence of the solution was tested in the q norm, which is appropriate only for the WH setup.In this scenario, we observed second order convergence for both the SYMH and WH cases.Our explanation for this phenomenon is that the H 1 norm provides more control over derivatives than q, so the given data are appropriate for both SYMH and WH setups.Furthermore, when monitoring only the q norm of the solution for the WH case we ignore the derivative terms that are responsible for the loss of convergence we see in H 1 .For the SYMH case, all derivative terms converge well.In the q norm, we only monitor the convergence of one of them, which dominates the norm.This scenario is crucial in understanding the reason that CCE between a SYMH IBVP and a WH CIBVP is numerically well behaved.The key point is that for H 1 controlled given data, the solution to the SYMH IBVP is also controlled in the same norm.This solution then provides boundary data for the WH CIBVP with enough control over derivative terms and so the solution to that PDE problem converges well in the q norm.
The expected convergence was also seen for CCM between IBVP and CIBVP with the same degree of hyper-bolicity, that is when matching SYMH to SYMH or WH to WH.In the first case, we could obtain good convergence for the solution in the L 2 and H 1 norms, for given data controlled in the same norms, respectively, whereas in the second, the same was true in the q norm.When we analyzed the convergence of the solution in the q norm for H 1 given data, we observed the same scenario as earlier, namely second order convergence in the q norm for both SYMH-SYMH and WH-WH setups.Our understanding of this behavior is exactly the same as for the individual IBVP and CIBVP, with the only difference that in this case the various norms considered were adapted to the composite CCM PDE problem.
When matching a SYMH IBVP to a WH CIBVP we did not recover good second order convergence in any of the L 2 , q, or H 1 norms for the solution, when the given data were controlled in the same norm.This result is in line with our expectation from the preceding continuum analysis.It is worth mentioning that as earlier, we were able to see good second order convergence in the q norm when the given data were controlled in the H 1 .Even though this behavior is qualitatively the same as for SYMH-WH CCE, the interpretation in terms of continuous dependence of the solution on the given data is different.In CCM, the composite IBVP-CIBVP PDE problem is treated as one and therefore the solution is controlled in a smaller norm than the given data (some derivative terms are not controlled).In CCE, this was not an issue, because we could verify convergence of the IBVP solution in the H 1 norm, for H 1 given data, and for the CIBVP solution in the q norm for H 1 given data, where the latter include also given data controlled in the q norm.
Our homogeneous models exhibit the same qualitative behavior for SYMH-WH CCE and CCM due to the lack of source terms.To demonstrate this, we considered an inhomogeneous example in which for the CIBVP we added a source term such that the right-moving variables which are responsible for weak hyperbolicity affect the left-moving variable, and in IBVP the left-moving variable couples to the first right-moving one.In this way, there is an interaction loop between left-and rightmoving variables that is closed only when we consider CCM but not CCE.Our understanding is that when this loop is closed, it allows the part of the CIBVP solution that is influenced by the WH structure to affect in turn the IBVP solution, which consequently provides bad boundary data to the CIBVP and drives the numerical solution to explode.In our numerical experiments we see severe loss of convergence for the inhomogeneous SYMH-WH CCM setup in all norms and for all given data, even for the combination of H 1 given data and q norm of the solution.In fact, we observe that the norms involved in the calculation of the convergence rate grow exponentially with time and faster for increasing resolution, which is a clear sign for loss of convergence.For comparison, we tested the inhomogeneous case for SYMH-SYMH and WH-WH CCM.In the first case we found good sec-ond order convergence in the CCM L 2 and H 1 norms for the solution, whereas in the second we saw second order convergence in the q norm only for a brief initial period, the cause of which is presently unclear to us.Finally, to assess that indeed the two-way communication between the Cauchy and the characteristic domain is the crucial point in this setup, we tested the inhomogeneous SYMH-WH CCE setup for H 1 given data.In this case, we find second order convergence for the inhomogeneous WH CIBVP, exactly as in the inhomogeneous case.We interpret these results as support of our understanding of the mechanism that drives this loss of convergence, as given above.
We should highlight that we used high frequency given data in all our tests.More specifically, we used random noise with amplitude that dropped appropriately with increasing resolution.This type of data is useful in detecting WH in numerical experiments.In fact, we were able to see second order convergence in the L 2 norm for our inhomogeneous SYMH-WH CCM model, for smooth given data (see [60]).We understand that random noise high frequency data may not be appropriate for nonuniform grids, as for instance in spectral methods.In that case, one has to construct a different type of high frequency given data.A possible suggestion could be a sinusoidal of some high but resolved frequency, that increases with increasing resolution, and an amplitude which is dropped with increasing resolution, mimicking the expected behavior of numerical error.

V. CONCLUSIONS
The take home message from this analysis for current CCE setups in GR is that if the solution to the Cauchy setup is sufficiently controlled (in terms of derivatives) and provided that the WH CIBVP is weakly well posed (a strong assumption that should be firmly established) then, depending on the discretization, the numerical solution of this CCE should converge to the true solution in the limit of infinite resolution.This is another way to say that error estimates on such a solution can be trusted to an arbitrary small level (given the necessary computational resources) to provide highly accurate gravitational waveforms.Ultimately any CCE scheme is still subject to artificial boundary conditions, for which a detailed consideration is beyond our present scope.
Another lesson we take from this study is that CCM for the EFE as currently performed between a SYMH IBVP and WH CIBVP cannot be guaranteed to provide numerical approximations that converge to the true solution of the problem in the limit of infinite resolution.We do not expect the results to be changed if symmetric is replaced with strong hyperbolicity for the IBVP.Developing a strongly, if not symmetric hyperbolic characteristic formulation of the EFE seems to be necessary to move toward an implementation of GR CCM that is (strongly) well posed and properly convergent.

4 .
FIG. 4. The CCM domain D3 = D1 + D2, with boundaries the spacelike Σt i , Σt f , the null Nu i , Nu f and the timelike Tρ min , Tx max .Initial data for all fields are provided on Σt i , as well as on Nu i only for left-moving fields.Boundary data are given on Tρ min for right-moving, and on Tx max for left-moving fields.Data are communicated between the domains D1 and D2 via the worldtube T0 = Tρ max = Tx min .The z direction is compact.

3 FIG. 5 .
FIG.5.The exact convergence rate Cexact is shown for the IBVP, for all the tests.The top row refers to the SYMH setup, and the bottom to the WH one.The left column presents the results for the L 2 − L 2 test, the second to left for the q − q test, and so on.The results of the L 2 − L 2 , q − q, and H1 − H1 tests are as expected by theory, that is convergence in the L 2 and H1 for the SYMH case and in the q norm for the WH one.The convergence in the q norm seen in the H1 − q test is merely a result of the fact that H1 (26) is a bigger norm than q (25), in the sense that it controls more derivative terms.We only monitor the ∂zϕ1 derivative term, which still converges appropriately for the SYMH case and due to its size in comparison to the nonderivate terms, dominates the q norm.We also see good convergence of the solution to the WH IBVP in the H1 − q test, because we ignore the derivative terms that do not exhibit convergence.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs.

3 FIG. 6 .
FIG.6.Exact convergence rate Cexact for the CIBVP solution, for our tests in the different norms.The SYMH setup is shown in top row and the WH one in the bottom.Each column shows the results of each test.The L 2 − L 2 , q − q, and H1 − H1 tests exhibit the expected convergence in the L 2 and H1 for SYMH and in the q norm for WH.Convergence in the q norm seen in the H1 − q test is again a result of the fact that H1 (29) is bigger than the q (28) norm, as in the IBVP.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs.

3 FIG. 8 .
FIG.8.The exact convergence rate of the CCM solution, for SYMH-SYMH (top) and WH-WH (bottom).The first column refers to the L 2 − L 2 test, the second to the q − q test, and so on.The L 2 − L 2 , q − q, and H1 − H1 tests show the expected convergence in the L 2 and H1 for SYMH-SYMH and in the q norm for WH-WH.The convergence in the q norm seen in the H1 − q test is again a result of the fact that H1 (32) is bigger than q (31), as in the IBVP.The plots for the H1 − H1 and H1 − q tests use data from exactly the same runs.
Table I contains the list of all the acronyms used throughout the text.

TABLE I .
List of acronyms used throughout the text.
2. The Cauchy domain D1, with boundaries that are spacelike Σt i , Σt f and timelike Tρ min , Tρ max .Initial data are provided on Σt i and boundary data on Tρ min , Tρ max for the right-and left-moving variables, respectively.The z direction is compact.
The characteristic domain D2, with boundaries that are null Nu i , Nu f and timelike Tx min , Tx max .Initial data for left-moving fields are provided on Nu i and boundary data on Tρ min , Tρ max for the right-and left-moving variables, respectively.The z direction is compact.