Construction of initial data

4.2 Construction of initial data

In order to start the evolution one has to determine the initial data. So far, the numerical methods used for their construction are based entirely on Theorem 7. It is assumed that the extrinsic curvature of the hyperboloidal initial surface is pure trace. Then the Yamabe equation is solved for a conformal factor, which chooses from a given conformal class of Riemannian metrics the one that has constant negative curvature. Finally, the other initial data are obtained from Equations (35

, 36

There are two problems with this kind of procedure. First, the Yamabe equation degenerates on the boundary. This is a difficult problem if one tries to prove existence of solutions of this equation, because the loss of ellipticity on the boundary means that one cannot simply appeal to known theorems. However, numerically this has not been a serious obstacle. The Yamabe equation is solved in a more or less standard way by using a Richardson iteration scheme to reduce the solution of the non-linear equation to a series of inversions of a linear operator.

The degeneracy of the equation on the boundary forces the solution and its derivative by regularity to have certain well defined boundary values. This and the global nature of the elliptic equations suggests that we use pseudo-spectral (or collocation) methods. Such methods are well suited for problems for which it is known beforehand that the solution will be sufficiently regular. Then, the solution can be expanded into a series of certain basis functions that are globally analytic. Therefore, the regularity conditions on the boundary are already built into the method. We refer readers interested in these methods to a recent review paper by Bonazzola et al. [25 ].

The more difficult problem in the determination of the initial data is the fact that the curvature components are obtained by successive division by $Ω$ . Since $Ω$ vanishes on the boundary one has to compute a value which is of the form $0∕0$ . Although the theorem tells us that this is well defined (provided that some boundary conditions are satisfied), it still poses a numerical problem, because a straightforward implementation of l’Hôpital’s rule loses accuracy due to the truncation error in the calculation of the numerator.

There exist two methods to overcome this problem. The first one [47 , 51 ] is again based on pseudo-spectral methods. Essentially, in this approach the multiplication with the function $Ω$ is expressed as a linear operation between the expansion coefficients of $Ω$ and the other factor. The pseudo-inversion of this linear operator (in the sense of finding its Moore–Penrose inverse [122 ]) then corresponds to division by $Ω$ .

A different method for computing the quotient has been devised by Hübner [95 ]. Here, the problem of dividing by $Ω$ has been reformulated into a problem for solving a second-order PDE for the quotient. This PDE is similar to the Yamabe equation in the sense that it also degenerates on the boundary where $Ω$ vanishes, but it is linear.

Figure 12: An axi-symmetric solution of the Yamabe equation with $u = 1$ on the boundary (top) and its third differences in the $z$ -direction (bottom). The location of $ℐ$ is clearly visible.

Another consequence of the singularity of the Yamabe equation on $ℐ$ is the following. It is usually desirable to extend the solution of the Yamabe equation beyond $ℐ$ in such a way that it remains a solution of the equation. The reason is mainly that $ℐ$ is usually not aligned with the boundary of the computational domain so that it is impossible to restrict the calculation to the physical domain inside $ℐ$ . Furthermore, the numerical stencils used in the discretisation of the evolution equations usually extend beyond $ℐ$ into the unphysical part of the computational domain, thus picking up information from the outside. Therefore, if the constraints are not satisfied outside $ℐ$ then it might be possible that constraint violating modes are excited. While the solution of the Yamabe equation inside $ℐ$ is uniquely determined, this is not so on the outside, where the solution is influenced by the choice of boundary conditions on the boundary of the computational grid. This means that the solution cannot be smooth across $ℐ$ unless the boundary conditions are chosen exactly right. But in general it is impossible to know the correct boundary values beforehand. It is guaranteed that the solution will be at least $𝒞2$ across $ℐ$ , but the third derivative may be discontinuous.

Figure 13: The difference between two solutions of the Yamabe equation obtained with different boundary conditions outside. Only the values less than $0.005$ are shown.

This situation does in fact occur as illustrated in Figure 12 , which shows in the top picture a solution $u$ of the Yamabe equation in an axi-symmetric space-time. The $z$ -axis ( $r = 0$ ) is the axis of rotation, and the $r$ -axis ( $z = 0$ ) corresponds to the equatorial plane. The solution is assumed to be equatorially symmetric, so that it is obtained by specifying $u = 1$ on the boundary $r = 1$ and $z = 1$ while using the symmetry conditions on the other parts of the computational grid. On the bottom are shown the third differences of $u$ in the $z$ -direction. It is clearly visible that $ℐ$ is located at the quarter circle $r2 + z2 = 0.8$ . In Figure 13 we illustrate the “shielding” property of $ℐ$ for solutions of the Yamabe equation. Here is shown the difference of two solutions obtained with different boundary data on the outside. The solutions agree to a very high accuracy inside the physical domain. Some of these issues and possible ways to overcome them are discussed in further detail in [51 ].