The gauge conditions described above have been presented from a theoretical point of view. Their implementation in the numerical code will require further tuning and ex-perimenting, changing the source functions and damping coefficients in order to obtain a stable numerical simulation. The tuning and changes performed in practice to the gauge conditions have to respect certain requirements:
• the gauge equations of motion have to be well-behaved as single equations: no exponential growths should be present;
• they have to provide an appropriate stationary solution, i.e. that is compatible with the rest of the system and ensures that the slices in the evolution continue to be hyperboloidal;
• they form a hyperbolic system (see section 5.1) with the Einstein equations;
• they have appropriate eigenspeeds - this has to be checked especially at I+, where no incoming speeds should appear, and at the horizon (if present and if excision is being used).
I have not presented here the exact form of the source functions and damping coeffi-cients because they are tuned according to the rest of equations involved in the simulations.
The actual gauge evolution equations used in the numerical tests are presented in section 7.3.
Chapter 5
Well-posedness and regularity
In this chapter we will treat two important properties of the continuum equations: the well-posedness of their initial value formulation, which ensures the existence of a unique solution, and their regularity, especially important at future null infinity, where the con-formal factor terms diverge. If the system is not well-posed or the regularity conditions are not satisfied, even with a perfect numerical implementation the simulations will be unstable and crash at some point.
5.1 Treatment of hyperbolic equations
In chapter 2 the Einstein equations were decomposed into a system of elliptic (constraints) and hyperbolic (equations of motion) differential equations.
Hyperbolic equations usually describe dynamical processes that involve the propa-gation of information in the form of waves from one place to another following a time evolution. The difficulty is that not any choice of variables or equations is appropriate;
they need to be formulated as a well-posed problem to avoid the instabilities that will appear otherwise.
5.1.1 Well-posedness
Well-posedness is a necessary requirement for successfully solving a system of equations iteratively in time. It implies that an unique solution to the equations exists and that this solution depends continuously on the initial data provided. This last condition can be expressed as (see e.g. [85])
||u(t)|| ≤Keat||u(0)||. (5.1) The growth of the solution u(t) with respect to the initial data u(0), calculated using an appropriate norm ||.||, is bounded by certain finite values of K (not to be confused with the trace of the extrinsic curvature) and a, so that it cannot be infinite.
Quasilinear systems of equations are systems which are linear in the highest derivatives.
This is the case of the Einstein equations considered here. The study of well-posedness can be performed on this type of systems using a simple standard method. In this method it is sufficient to consider the principal part of the equations. The spherically symmetric reduced equations that will be used here are expressed as a mixed first order in time and
81
second order in space system. The exact terms that belong to the principal part of such a system are specified in the following subsection, but the basic rough idea is that it involves only the highest derivatives of each type of variables. General theorems like theorem 4.3.2 in [85] state that lower order terms do not contribute to the well-posedness of the system.
The algebraic analysis performed on the principal part will determine the hyperbolicity of the system of equations, which allows to make a claim on the well-posedness of the corresponding initial value problem.
5.1.2 First order in time and second order in space systems
The equations obtained from the 3+1 decomposition contain first order in time derivatives and second order in space derivatives. Their hyperbolicity analysis is more complicated than in the first order in space reduction, whose features are described by plenty of textbooks, but an implementation of the equations as a second order in space system has a similar numerical behaviour and avoids introducing extra variables and constraints.
The analysis summarized here follows [55, 101].
A common first order in time and second order in space system of hyperbolic differ-ential equations has the following structure:
∂tu =Pu with u= u
v
, (5.2)
with the matrix P given by P =
Ai∂i+B C Dij∂i∂j +Ei∂i+F Gi∂i+J
. (5.3)
There are two different types of variables: u represents those which are derived twice in space, whereas v stands for the ones which are derived at most once.
A Fourier transformation is now performed on the system, usingωi ≡ω ni and Mn ≡ Mini. This gives the equivalent matrix to P in the Fourier space, ˆP. Its principal part is given by the higher order derivatives in each of the four blocks and is denoted as ˆP0, the second order principal symbol: The system will now be reduced to first order in the Fourier space using a pseudo-differential reduction. To do so, the new variable ˆw ≡ iωuˆ is introduced, giving rise to the following system (and the constraintC = ˆw−iωu):ˆ
5.1. Treatment of hyperbolic equations 83 The principal symbol of the first order is given by the elements in ˆPR that are multiplied byiω:
PˆR0 =
0 0 0
0 iωAn iωC 0 iωDnn iωGn
. (5.7)
Now writing E = ˆPR0/(iω) to simplify the notation, the principal part of the system is formally solved as:
∂tuˆR= ˆPR0uˆR=iωEuˆR ⇒ uˆR(t) =eiωEtuˆR(0). (5.8) IfE is diagonalizable and has real eigenvalues, the solution will be purely oscillatory and bounded by the initial data. In this case, the initial value problem is well-posed.
5.1.3 Hyperbolicity
The analysis of the principal part matrix E is performed using a Jordan decomposition.
It can be considered a generalization of the Eigen decomposition: apart from the ele-ments in the diagonal, the Jordan canonical form can also have non-zero values in the superdiagonal. These appear due to the existence of Jordan-blocks, which are placed in the diagonal of the Jordan canonical form. If a matrix is diagonalizable, its diagonal and Jordan canonical forms coincide.
Both complex eigenvalues or Jordan-blocks in the Jordan canonical form of the ma-trix E can cause growing modes that will give an ill-posed problem. The properties of the matrix E determine the hyperbolicity of the equations, and thus the well-posedness (provided that the similarity matrix is regular). The different possibilities are:
• If all eigenvalues (the characteristic speeds) of the principal part matrix are real, the system is called (weakly) hyperbolic, and it is well-posed in absence of lower order terms.
• If a complete set of (left) eigenvectors exists and no Jordan-blocks appear, the system is called strongly hyperbolic, and admits a well-posed initial value problem.
• If the system is strongly hyperbolic and also admits a (direction-independent) con-served energy it is called symmetric hyperbolic. In one dimension strongly hyper-bolic systems are always also symmetric hyperhyper-bolic, because they are not direction-dependent and a conserved energy is to be found without exception.
The BSSN-type formulations of the Einstein equations are known to be strongly hyperbolic [135, 117, 81, 29]. Indeed the systems of evolution equations (2.76), (2.78) and (2.82) together with the gauge conditions considered in chapter 4 have a diagonalizable principal part matrix with real eigenspeeds and a full set of eigenfields. Adding spherical symmetry to it means that the system is symmetric hyperbolic and thus the problem is well-posed.
5.1.4 Characteristic decomposition of the equations
The system of spherically symmetric GBSSN and Z4c equations presented in (2.76) can be analyzed with the tools presented in the previous subsection. As a first order in time and second order in space system, its variables are divided into
u variables, those which appear spatially derived twice: χ,γrr, γθθ, α, βr and ˜Φ.
v variables, with only first spatial derivatives: Arr, K, Λr, Θ, Br (if present) and ˜Π.
The eigenfields of the system will include the spatial derivatives of theuvariables, denoted byχ0, γrr0 , γθθ0 , α0, βr0 and ˜Φ0.
The lightspeeds of the system are calculated from the line element by imposingd¯s2 = 0 and solving for c = drdt. The result is c± = −βr ±αq χ
γrr. There is an extra speed due to the change in the coordinates given by c0 = −βr. These three velocities are plotted in figure 5.1 for flat spacetime and in presence of a Schwarzschild BH. The speeds are proportional to|KCM C|. As expected for asymptotic flatness, the behaviour of the speeds in the vicinity of I+ does not change if a compact object is present at the origin. There are no negative speeds atI+, which is a consequence of the fact that future null infinity is an ingoing null hypersurface in the conformal picture. Similarly, at the horizon of the Schwarzschild BH there are no positive speeds, because all radiation is infalling.
0.0 0.2 0.4 0.6 0.8 1.0
-0.2 0.0 0.2 0.4 0.6
r c
-c+
c0 Flat
0.0 0.2 0.4 0.6 0.8 1.0
-0.2 0.0 0.2 0.4 0.6
r c
-c+
c0 Schwarzschild
Figure 5.1: Zero speed c0 = −βr and lightspeeds c± = −βr ± αq χ
γrr plotted for flat spacetime (left) and a Schwarzschild BH withM = 1 (right), usingKCM C =−1 for both.
Again, rI = 1. The vertical line in the right plot indicates the radial position of the horizon, where c0, c−<0 and c+ = 0. AtI+ we have c0, c+ >0 and c−= 0.
The eigendecomposition (or characteristic decomposition) of (2.76) is presented in table 5.1, to which the eigenfields χ, γrr, γθθ, α, βr (if evolved) and ˜Φ have to be added with zero propagation speed. For generality we have determined the eigenfields and eigenspeeds keeping the angular metric component γθθ as an evolution variable, but it can be eliminated using the freedom of the determinant of the conformal metric. Indeed it is only after fixing this degree of freedom that we obtain the standard Z4c/CCZ4 formulation out of our equations. Ifγθθ is eliminated, the first eigenfield listed in table 5.1 vanishes and the substitution (2.71) has to be applied to the rest of eigenfields.
5.1. Treatment of hyperbolic equations 85
Table 5.1: Eigenfields and eigenspeeds, using c0 =−βr and c± =−βr±αq χ
γrr.
System Eigenfields Eigenspeeds at I+
GBSSN / Z4c γγrr0
As already shown in figure 5.1 and also indicated in the last column of table 5.1, the eigenspeeds at I+ are either zero or positive, so that no incoming radiation exists.
The first group of five eigenspeeds and eigenfields is given by the Einstein equations and the generalized harmonic lapse (chosen as (4.26) or (4.37a)) and is common to all combinations (except for the first one that disappears if γθθ is eliminated). In the second group the first row corresponds to an evolution without Θ variable, whereas the other two appear in case of the Z4c system. The choice of a fixed shift or the generalized harmonic shift condition (given by (4.32b) or (4.37b), as their principal part is the same) determines which eigenquantities to select from the third group. The fourth group accounts for the scalar field equations, whose principal part decouples from the rest of the system.
This eigendecomposition can be compared with the ones presented in [49, 67], making the substitution of Λr by Γr, rearranging some eigenvectors, omitting the scalar field part and substituting the harmonic lapse by the 1+log slicing condition and the harmonic shift by a Gamma-driver condition.
The eigendecomposition corresponding to the Gamma-driver shift conditions (de-scribed in subsection 4.4.2) is not shown here because their free parameters λ and µ make the eigenfield expressions very complicated and long. What is important is the fact that evaluating the eigenspeeds at I+ allows us to determine the allowed values of the parameters. There will be no negative eigenspeeds (no incoming modes) at I+ if λ µ≤ rI3KCM C2
of the Gamma-driver condition with auxiliary variableBr (4.29) and if λ≤ 121 (rIKCM C)2 for the integrated version (4.30). The reason for keeping the parameter µ separated from the parameterλ when tuning the Gamma-driver shift condition (4.29)
is that even if the eigenspeeds depend only on λ µ, different compatible choices of λ and µcan present a different numerical behaviour.
The change from conformal K to physical ˜K or to ∆ ˜K only involves a rescaling by Ω in K in the eigenfields. The same is valid when transforming from Θ to ˜Θ. The eigendecomposition for (2.78) is obtained by substituting K = KΩ˜ and Θ = ΘΩ˜, and the one for (2.82) by transforming K = ∆ ˜ΩK and Θ = ΘΩ˜. The other terms involved in the transformation (2.39e) become non-principal part terms and therefore do not appear in the eigendecomposition.