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 CCM C and the numerical determination of ¯Ω.
The critical value ofCCM C calculated from the values ofM,QandKCM 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 CCM 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 CCM C is only determined once in the numerical simulation.
Once the trumpet value of CCM 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 ofR0, the double root of (3.41) that corresponds to the choice of critical CCM 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ˆ − R0 r
Ω.¯ (6.16)
Now the differential equation for ˆΩ is given by Ωˆ0 = 1
3rR02 h
9CCM C2 Ωˆ6−54CCM C2 Ωˆ5−6 ˆΩ3 30CCM C2 +CCM CKCM CR30 −3M R03 +135CCM C2 Ωˆ4+ 9 ˆΩ2 15CCM C2 + 2CCM CKCM CR30+R02 R20−6M R0
−18 ˆΩ 3CCM C2 +CCM CKCM CR30+R02 R20 −3M R0
+9CCM C2 + 6CCM CKCM CR30+R02 KCM C2 R40 −18M R0+ 9R201/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) = ∆r2A 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 =rI) = 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 =rI) = 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 |KCM 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 ofKCM C, the larger the numerical precision required andKCM 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, KCM C =−1 and corresponding critical CCM C. Future null infinity is located at rI = 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 Arr should therefore take initially its flat spacetime value there (Arr0|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 ˜Φ00 is zero, more specifically at r = 0 and r = rI. 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ψAis 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 KCM 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 CCM 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 KCM 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 10300, 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.