Evolution equations

4.3 Evolution equations

The initial data constructed in one of the ways described above must be evolved with the propagation equations extracted from the conformal field equations. As we remarked upon above, there are several ways to write the evolution equations, and therefore there are also several ways to implement them numerically. To date, there exist two different implementations. One is based on the ADM-like formulation of the conformal field equations, while the other uses a tetrad formalism. In the ADM-like formulation developed by Hübner [95

] and Husa [98

, 100 ], one introduces the coefficients of the intrinsic metric on the space-like slices and their extrinsic curvature as variables, while in the other formulation by Frauendiener [46

], a tetrad and connection coefficients with respect to that tetrad are the basic variables. Concerning efficiency and usefulness, not much difference between these formulations is to be expected, provided the numerical schemes are the same. One minor advantage of the tetrad formulation might be the additional freedom in the orientation of the tetrad, which can be used on $ℐ$ to make the computational tetrad agree with the preferred tetrad on $ℐ$ (see below).

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 [21 , 22 ]. It allows one to compute the entire future of an initial hyperboloidal data set within that class up to time-like infinity $+ i$ .

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 $Δt$ 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.

Our next point is concerned with the extraction of the radiative information from the numerically generated data. This is the part of the entire numerical process where the superiority of the conformal approach becomes apparent. How does one determine the radiative field? First of all, one needs to find $ℐ$ on the current time-slice. Since $ℐ$ is the surface on which the conformal factor $Ω$ vanishes and since $Ω$ is explicitly known during the evolution, the location of $ℐ$ is a simple task. The next problem is concerned with the orientation of the tetrad on $ℐ$ . The asymptotic quantities are defined with respect to a specific geometrically characterized tetrad, a Bondi frame. But, in general, this tetrad is completely unrelated to the “computational” tetrad used for the evolution. Therefore, one needs to find the transformation from one to the other at each point of $ℐ$ . Without going into too many details (see [46 , 49 , 48 ]), we remark that most asymptotic quantities, in particular the radiation field, are of a local character, so they can be read off without constructing a Bondi frame. This is rather fortunate because there are global issues involved in the transformation from the computational tetrad to the Bondi frame. These have implications for the determination of global quantities like the Bondi energy-momentum four-vector, but they have no effect on the radiation field, which is defined as that (complex) component of $Kabcd$ that is entirely intrinsic to $ℐ$ :

ψ4 = Kabcd na ¯mbnc ¯md.

Here, $na$ is a null-vector tangent to the generators of $ℐ$ , i.e. $na ∝ ∇a Ω$ , and $ma$ is any complex space-like null-vector that is orthogonal to $a n$ . It is useful to require the space-like vector to be tangent to the intersection of $ℐ$ with the current time-slice. Augmenting these two vectors by a further real null-vector $la$ yields a null-tetrad which is fixed up to rotations in the (two-dimensional) tangent space of that intersection and boosts in the plane orthogonal to it. The behaviour of $ψ 4$ under these transformations is that of a GHP-weighted quantity [80 , 133 ] with boost weight $− 2$ and spin weight $− 2$ . This corresponds to the quadrupole-like character of the gravitational radiation field. However, $ψ4$ really depends only on the null-vector $a n$ . Suppose we perform a null-rotation around $na$ , then $ma$ transforms into $ma + αna$ for some complex valued function $α$ on $ℐ$ . But $ψ4$ is invariant under this transformation. So in order to find $ψ 4$ it is only necessary to transform from the given computational tetrad to the tetrad specified above, which is rather straightforward. In fact, the computation of $ψ4$ involves only the combination of certain components of the gravitational field with powers of $∇a Ω$ .

The final step in the correct determination of the radiation is to find the correct time parameter. Suppose we follow a specific null generator of $ℐ$ crossing through successive time-slices. On each slice we compute $ψ4$ on that generator. Then we obtain the radiation emitted by the source into the direction specified by the generator as a function of our computational time parameter. Since the time coordinate is rather arbitrary, this means that the waveform determined so far has no physical meaning. The problem is already present in Maxwell’s theory: Suppose we have an emitter that sends out a pure sine wave. A detector far away from the source cannot determine the absolute frequency of the signal because the relative velocity of emitter and receiver might be non-zero, but the detector should also find a pure sine signal. However, this will be true only if the detector records the signal as a function of proper time. Any other time parameter along the detector’s world-line will not produce a pure sine.

What one needs to do in the general case is to select among all parameters along the generator a specific, geometrically distinguished one, namely a Bondi parameter. A generator and such a parameter along it can be understood as a certain limit of freely falling observers with proper time clocks as they move towards infinity [49 ]. Bondi parameters are obtained as solutions of an ordinary, linear, second order differential equation, which is conformally invariant.

The computation of the Bondi energy-momentum is a global procedure, i.e. it depends on properties of the entire cut of $ℐ$ with the current time-slice. There are two steps involved in this procedure. First, one needs to obtain the asymptotic translation group (see e.g. [134 ]) on each cut. This provides four functions on the cut which are then, in a second step, integrated against the “mass aspect” which is another function obtained from the “Coulomb” part $ψ2$ of the gravitational field, and the “news function,” which is a combination of components of the Ricci tensor and connection coefficients. The first step, the determination of the translation group, is the global step because it involves solving a second order elliptic equation on the cut. These issues are discussed in more detail in [48 ].

All these procedures for finding the relevant data on $ℐ$ have been worked out analytically and they have also been tested (at least in part) numerically [46 ]. The tests have been performed under the assumption that null-infinity admits toroidal cuts, which has the advantage that one can actually compare the numerical results with analytical expressions because a whole class of exact solutions [151 , 93 ] is known to exist. Admittedly, such space-times are rather unphysical, but since most of the extraction procedures are local there is no doubt that they will also work in more realistic cases.