### 4.3 Evolution equations

So far, the discretization of the equations has been rather straightforward. One of the schemes that has been used is the higher dimensional Lax–Wendroff scheme, also called the rotated Richtmyer scheme, a discretization scheme with second order accuracy. It has been employed alone [46] or together with Strang splitting [95] to treat the principal part of the equations differently from the source part. Since a second order scheme requires much more computing resources compared to higher order methods to achieve the same accuracy, Hübner started to use the method of lines [86] with a fourth order scheme to compute the spatial derivatives and fourth order Runge–Kutta for the evolution in time. He reports [95] that the fourth order method is very much superior to the second order scheme in terms of efficiency. (The feasibility of the method of lines in relativity has been studied by Hungerbühler [97], using pseudo-spectral methods for the spatial derivatives and a combination of Adams–Bashforth and Adams–Moulton schemes for the time evolution.) A similar method has been employed by Frauendiener and Hein [52] in a code that evolves the conformal field equations under the assumption of axi-symmetry. The code makes use of the “cartoon method” [2] (see also [50]), which can be used to treat axi-symmetric systems using Cartesian coordinates. The code has been used to reproduce exact solutions such as the boost-rotation symmetric solutions of Bičák and Schmidt [2122]. It allows one to compute the entire future of an initial hyperboloidal data set within that class up to time-like infinity .

The conformal field equations and the propagation equations derived from them are quasi-linear. This implies that the characteristics of the system depend on the current solution, and this, in turn, means that one has to be able to change the time-step between successive time-slices in order to keep a stable evolution scheme. This is necessary because the schemes are explicit schemes and, therefore, subject to the Courant–Friedrichs–Lewy condition [37], which states that the numerical domain of dependence of a point should always include the analytical domain of dependence. This requirement already excludes the popular leap-frog scheme, which is nevertheless used sometimes also for evolving the Einstein equations. A general criterion for computing the maximal time-step allowed in each iteration in arbitrary dimension has been derived in [46].

Another important point in the development of evolution codes is the numerical treatment of the boundaries. As explained already above, it is one of the advantages in the conformal approach that the outer boundary is not as influential as it is in the conventional approach using the standard Cauchy problem. It was also pointed out that in this case it is enough to impose a boundary condition that results in a numerically stable code because the outer boundary is located in the unphysical region and, therefore, cannot influence the physical space-time.

The proper way to treat the boundary is to prescribe conditions that are compatible with the full conformal field equations, in particular with their restriction to the boundary manifold. This has not been done so far. Since the outer boundary is not important for the physical effects, other ways of dealing with the boundary have been devised. One way is to forget about the restrictions of the conformal field equations to the boundary and to analyze the possible boundary conditions for the propagation equations. To first order, one can define ingoing and outgoing fields on the boundary. Then a sufficiently general boundary condition will be obtained by specifying the ingoing fields in terms of the outgoing ones. Although this boundary treatment is not necessarily compatible with the restriction of the conformal field equations to the boundary, it is compatible with the evolution equations. This means that the evolution can remain stable, although the solution will not satisfy the constraints in the domain of influence of the boundary, which, however, is always in the unphysical part of space-time. This method has been used in [46] with satisfactory results. In particular, the boundary did not give rise to non-physical modes. These findings are in agreement with the analysis of numerical boundary conditions by Trefethen [158].

Another method for dealing with the boundary has been found by Hübner [94]. He realized that it is sufficient to solve the conformal field equations in the physical space-time only, and not necessary to solve them in the unphysical region as long as the characteristics remain such that the information created in the unphysical part of the computational domain cannot reach the physical part. Consequently, in his treatment the grid is divided into three zones: the inner zone, the outer zone, and a transition zone. The inner zone covers the physical space-time (flagged by a positive conformal factor) and some part of the adjacent unphysical region. On this part of the grid the conformal field equations are solved. In the outer zone, which is located in a neighbourhood of the grid boundary, one solves an advection equation that propagates outwards, off the grid. In the transition zone, a sufficiently smooth interpolation between these two systems of equations is solved. The effect is that the boundary condition that has to be imposed on the grid boundaries is very simple, and that the noise generated in the transition region is propagated away from the physical region outward towards the grid boundary.