The code used to run the simulations presented in this thesis is a one-dimensional spher-ically symmetric code. The time integration by the MoL is a 4th order RK, as described in (6.7). Finite differences of orders 2nd, 4th, 6th and 8th are implemented, but for the results shown here 4th order differences were used. Kreiss-Oliger dissipation is added to the discretized RHSs as described in section 6.4.

6.6. Numerical calculation of initial data 101

### 6.6.1 Calculation of the compactification factor Ω ¯

In flat spacetime ¯Ω ≡ Ω, where there is a simple analytical expression. However,
trum-pet initial data require the calculation of the critical value of C_{CM C} and the numerical
determination of ¯Ω.

The critical value ofC_{CM C} calculated from the values ofM,QandK_{CM C} by isolating
it from 3.42 is given by an analytical expression only for Q = 0,1. This analytical
expression can acquire an imaginary part due to round-off errors, so that in practice it
is more convenient to calculate the critical value of C_{CM C} numerically (for any value of
Q≤M). In the code this is done by means of a simple bisection method. The bisection
method is robust but relatively slow, although the second aspect is not a problem because
the C_{CM C} is only determined once in the numerical simulation.

Once the trumpet value of C_{CM C} is determined, we proceed to construct numerically
the profile of the compactification factor ¯Ω. As a first step it is convenient to calculate
the value ofR_{0}, the double root of (3.41) that corresponds to the choice of critical C_{CM C}.
This is also done with help of the bisection method.

The differential equation for ¯Ω is given by (3.28) substituting the appropriate
expres-sion for A(_{Ω}^{r}_{¯}). We now transform the compactification factor to a new ˆΩ given by

Ω = 1ˆ − R_{0}
r

Ω.¯ (6.16)

Now the differential equation for ˆΩ is given by
Ωˆ^{0} = 1

3rR_{0}^{2}
h

9C_{CM C}^{2} Ωˆ^{6}−54C_{CM C}^{2} Ωˆ^{5}−6 ˆΩ^{3} 30C_{CM C}^{2} +CCM CKCM CR^{3}_{0} −3M R_{0}^{3}
+135C_{CM C}^{2} Ωˆ^{4}+ 9 ˆΩ^{2} 15C_{CM C}^{2} + 2C_{CM C}K_{CM C}R^{3}_{0}+R_{0}^{2} R^{2}_{0}−6M R_{0}

−18 ˆΩ 3C_{CM C}^{2} +CCM CKCM CR^{3}_{0}+R_{0}^{2} R^{2}_{0} −3M R0

+9C_{CM C}^{2} + 6C_{CM C}K_{CM C}R^{3}_{0}+R_{0}^{2} K_{CM C}^{2} R^{4}_{0} −18M R_{0}+ 9R^{2}_{0}1/2

(6.17)
Note that there is no dependence on r inside of the square root. The transformation
(6.16) is such that ˆΩ and ˆΩ^{0} vanish atr = 0 and ˆΩ goes to unity atI^{+}. The transformed
compactification factor ˆΩ corresponding to the ¯Ω that was already presented in figure 3.6
is shown here in figure 6.5

The quantity ˆΩ is determined via a shooting and matching approach. First an initial
guess for the value of a quantity A is set at the origin as ˆΩ(r = ∆r) = ∆r^{2}A and it is
integrated along the radial coordinate with a RK of 4th order until I^{+} is reached. The
value obtained there is compared to ˆΩ(r =r_{I}) = 1. Using the bisection method a better
initial guess for A is found and the process is iterated until the discrepancy between the
calculated value at I^{+} and ˆΩ(r =r_{I}) = 1 is smaller than a chosen threshold. Now only
transforming the obtained profile back to ¯Ω using (6.16) is left.

There is a limit to the maximum value of |K_{CM C}| that could be set in the BH
simu-lations. The reason is that the just described calculation of the compactification factor ¯Ω
requires more thandouble precision to be calculated numerically. The larger the absolute
value ofK_{CM C}, the larger the numerical precision required andK_{CM C} <−1.5 needs more
than quad precision, the maximum available by the compiler used.

0.0 0.2 0.4 0.6 0.8 1.0 0.0

0.2 0.4 0.6 0.8 1.0

r

W`

Figure 6.5: Transformed compactification factor ˆΩ for M = 1, Q= 0, K_{CM C} =−1 and
corresponding critical C_{CM C}. Future null infinity is located at r_{I} = 1.

### 6.6.2 Numerical solution of the spherically symmetric constraint equations

The momentum constraint (3.48b) is independent of the quantity ψ, so that it can be
solved separately. The spacetime at I^{+} is asymptotically flat, which means that the
scalar field ˜Φ tends to zero when approaching I^{+} ( ¯Φ = ^{Φ}_{Ω}^{˜} is finite there). The quantity
A_{rr} should therefore take initially its flat spacetime value there (A_{rr0}|_{I}+ = 0). This
implies that the value of ψ_{A} defined in (3.46) has to satisfy ψ_{A}|_{I}+ = 0.

We will choose initial data of compact support for ˜Φ as described in subsection 3.5.1.

The parameters will be set in such a way that the perturbation vanishes at the origin and
at null infinity. This means that the source term in (3.48b) corresponding to the scalar
field vanishes everywhere where ˜Φ^{0}_{0} is zero, more specifically at r = 0 and r = r_{I}. Thus
the solution forψ_{A} will be zero there.

Knowing this we can set ψ_{A}|_{r=0} = 0 as initial value and integrate towards I^{+} using a
RK of 4th order. This method requires calculating the solution at half spatial steps. One
possibility to go around this is to integrate first the even gridpoints and then the odd ones,
which does not pose any problems given the simplicity of the differential equation. At the
end of the integration procedure we obtain that the value of ψ_{A} in the neighborhood of
I^{+} vanishes, as expected. Two example solutions for ψ_{A} are shown in figure 6.6, where
sign = 1 was was chosen to obtain a mostly outgoing scalar field perturbation. The
corresponding ingoing versions (sign=−1) would have the exact same form below zero,
sosign→ −sign⇒ψ_{A}→ −ψ_{A}. The parameters are the same ones used in the plots in
subsection 3.5.2.

Once the value ofψ_{A}is known, it is substituted into the Hamiltonian constraint (3.48a)
and the solution for ψ, defined in (3.46), can be obtained. This is done via a shooting
and matching method. An initial guess is made for the value of ψ at the origin and then

6.6. Numerical calculation of initial data 103

0.0 0.2 0.4 0.6 0.8 1.0

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

*r*
Ψ*A*

Regular-data-based

0.0 0.2 0.4 0.6 0.8 1.0

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

*r*
Ψ*A*

Schwarzschild-based

Figure 6.6: Examples of solution of ψ_{A} with K_{CM C} = −1 and sign = 1: on the left the
solution for regular data and on the right the one in presence of a BH with M = 1, Q= 0
and critical C_{CM C}. If sign <0, thenψ_{A} would take negative values.

the differential equation is integrated using a 4th order RK until I^{+} is reached, where
by comparing the result with ψ|_{I}+ = 0 (asymptotic flatness condition) a new guess is
produced. The procedure is repeated until the error in the ψ|_{I}+ obtained is below the
chosen threshold. Examples of ψ corresponding to perturbations of regular initial data
and a Schwarzschild BH, calculated using the ψ_{A} profiles in figure 6.6, are presented in
figure 6.7.

0.0 0.2 0.4 0.6 0.8 1.0

1.000 1.005 1.010 1.015 1.020 1.025 1.030 1.035

*r*

Ψ

Regular-data-based

0.0 0.2 0.4 0.6 0.8 1.0

1.000 1.002 1.004 1.006 1.008 1.010 1.012 1.014

*r*

Ψ

Schwarzschild-based

Figure 6.7: Examples of solution of ψ with K_{CM C} = −1 and sign = 1: on the left the
flat-spacetime-based solution and on the right the one in presence of a BH with M = 1,
Q= 0 and critical CCM C.

## Chapter 7

## Numerical experiments

### 7.1 Results of a straightforward implementation

The first numerical test I performed with the Einstein equations expressed in terms of the conformal metric was with the GBSSN evolution system in (2.76), with a fixed shift, the slicing condition (4.26) and with flat spacetime initial data as shown in figure 3.15.

The spatial grid was staggered, so that the location of I^{+} was avoided. The numerical
evolution was unstable and it crashed after a few time-steps. With “crash” I mean that
the variables (at least one of them) acquire a very large absolute value, of the order of
10^{300}, or they become indeterminate, with output “NaN” (not-a-number), in some part
of the integration domain. These effects are the result of either a continuum instability
of the equations or a numerical instability of the implementation.

The instability observed in this first numerical test presented high-frequency
fluctua-tions close to I^{+} that would appear after very few time-steps. At some point the values
of the variables become very large or NaN at one of the nearest gridpoints to I^{+} and
the simulation was considered to have crashed. The “crash-time” was not affected by
changes in the timestep, but it became smaller with an increase in spatial resolution
-meaning that the outermost gridpoint was closer toI^{+}. This indicates that the
instabil-ity is originated by a property of the continuum equations and not by a numerical artifact
of the implementation. Knowing that the system of equations is complicated and that
the CFEs are prone to continuum instabilities [103, 100], such a behaviour was expected.

Nevertheless, the problem can be analyzed and fixed.