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)|| ≤Ke^{at}||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:

∂_{t}u =Pu with u=
u

v

, (5.2)

with the matrix P given by P =

A^{i}∂_{i}+B C
D^{ij}∂_{i}∂_{j} +E^{i}∂_{i}+F G^{i}∂_{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 M^{n} ≡
M^{i}n_{i}. 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 ˆP^{0}, 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 ˆP_{R} that are multiplied
byiω:

Pˆ_{R}^{0} =

0 0 0

0 iωA^{n} iωC
0 iωD^{nn} iωG^{n}

. (5.7)

Now writing E = ˆP_{R}^{0}/(iω) to simplify the notation, the principal part of the system is
formally solved as:

∂_{t}uˆ_{R}= ˆP_{R}^{0}uˆ_{R}=iωEuˆ_{R} ⇒ uˆ_{R}(t) =e^{iωEt}uˆ_{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: A_{rr}, K, Λ^{r}, Θ, B^{r} (if present) and ˜Π.

The eigenfields of the system will include the spatial derivatives of theuvariables, denoted
byχ^{0}, γ_{rr}^{0} , γ_{θθ}^{0} , α^{0}, β^{r}^{0} and ˜Φ^{0}.

The lightspeeds of the system are calculated from the line element by imposingd¯s^{2} = 0
and solving for c = ^{dr}_{dt}. The result is c± = −β^{r} ±αq _{χ}

γrr. There is an extra speed due
to the change in the coordinates given by c_{0} = −β^{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|K_{CM 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*+

*c*0
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*+

*c*0
Schwarzschild

Figure 5.1: Zero speed c_{0} = −β^{r} and lightspeeds c_{±} = −β^{r} ± αq _{χ}

γrr plotted for flat
spacetime (left) and a Schwarzschild BH withM = 1 (right), usingK_{CM C} =−1 for both.

Again, r_{I} = 1. The vertical line in the right plot indicates the radial position of the
horizon, where c_{0}, c−<0 and c_{+} = 0. AtI^{+} we have c_{0}, 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 c_{0} =−β^{r} and c± =−β^{r}±αq _{χ}

γrr.

System Eigenfields Eigenspeeds at I^{+}

GBSSN / Z4c ^{γ}_{γ}^{rr}^{0}

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
λ µ≤ ^{r}^{I}_{3}K_{CM C}2

of the Gamma-driver condition with auxiliary variableB^{r} (4.29) and if
λ≤ _{12}^{1} (r_{I}K_{CM 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.