Regularity of the equations at null infinity

The regularity of the conformal factor terms was already mentioned in section 4.2, where the preferred conformal gauge was described. As a reminder, the divergent conformal factor terms arising from the transformation of the Ricci or the Einstein tensors due to the conformal metric rescaling, given by (2.3) and (2.6) respectively, are formally singular atI⁺. Certain relations between the variables have to be satisfied, so that the numerators cancel with the appropriate order of Ω and in the limit Ω→0 the divergent terms attain a finite value. These relations between the variables will be referred to as regularity conditions. We will assume that they are valid atI⁺.

The regularity conditions arising from (2.3) (or equivalently (2.6)) are summarized as:

∇¯^aΩ ¯∇aΩ

I⁺ = 0, (5.9a)

∇¯a∇¯bΩ− 1

4g¯ab¯ Ω

_I+

= 0, (5.9b)

Ω→0lim

2 ¯∇^aΩ ¯∇_aΩ

Ω −¯ Ω

I⁺

= 0. (5.9c)

The first one is obtained by multiplying (2.3) by Ω² and evaluating it at I⁺. For the other relations, (2.3) (or (2.6)) is multiplied by Ω and its trace-free part for (5.9b) and its trace for (5.9c) are evaluated at I⁺.

5.2.1 Regularity of the tensorial equations: approach

The regularity conditions of the tensorial equations (2.64) and (2.66) can be studied by separating the angular components from the radial one. For this, the variables and equations have to be decomposed in [2+1]+1 form.

The starting point is to write the line element as

ds² = −α²dt²+L²(dr+β^rdt)(dr+β^rdt)

+q_AB(dθ^A+b^Adr+β^Adt)(dθ^B+b^Bdr+β^Bdt), (5.10) where L⁻² = ¯γ^abD¯_arD¯_br equivalently as with the lapse function. The normal vector to the level sets of r is denoted by ¯s^a and is defined as ¯s^a = LD¯^ar. It is related to the two-metric ¯q_ab on the spheres orthogonal to these level sets by

¯

q_ab = ¯γ_ab−s¯_as¯_b. (5.11)

5.2. Regularity of the equations at null infinity 87 The conformal spacetime metric can be expressed in terms of the decompositions as

¯

g_ab =−¯n_an¯_b + ¯γ_ab =−¯n_an¯_b+ ¯s_as¯_b+ ¯q_ab. (5.12) The coordinate light-speeds of light-rays traveling in the ±sⁱ directions are,

c^r_± =−β^r±L⁻¹α, (5.13a)

c^A_± =−β^A∓b^AL⁻¹α, (5.13b)

radially and in the angular directions respectively. Obviously, c^r_± coincide with the light-speeds c± of the spherically symmetric reduction presented in subsection 5.1.4.

Expressing the regularity condition (5.9a) in the [2+1]+1 decomposition gives

−(L_n¯Ω)²+ (L_¯sΩ)²+q^ab∇¯_aΩ ¯∇_bΩ

I⁺ = 0. (5.14)

Assuming that the boundary Ω = 0 is also a level-set of the radial coordinate r=r_I, the previous expression becomes simply,

(L_¯nΩ)²

I⁺ = (L_¯sΩ)²

I⁺. (5.15)

One of the conditions for the scri-fixing described in 4.1 was that in our implementation Ω is time independent. Assuming the level set condition this translates to

∂_tΩ|_I+ = (αL_n¯Ω +Lβ^rL_¯sΩ)|_I+ = 0. (5.16) This is the basic picture where the regularity conditions can be derived. The calcula-tion itself is still ongoing work.

5.2.2 Regularity of the spherically symmetric equations

The regularity conditions of the spherically symmetric equations can be easily derived by looking at the numerator of terms divided by Ω. The regularity conditions for the evolution system (2.78), the one that uses the physical ˜K and is almost the same as the one tested numerically, are the following:

• There is only one term with Ω² in the denominator: ^2α_Ω2^ΘΩ˜γrr⁰, which appears in (2.78f), and regularity requires that the numerator vanishes at I⁺. This implies that

Θ˜

I⁺ = 0. (5.17a)

Using the l’Hˆopital rule, the mentioned term is substituted for the rest of the analysis by ^2α_Ωγ^Θ˜0

rr.

• The rest of the conditions will be given by the terms that are divided by Ω in the evolution equations. The condition that makes the numerator of Ω in ˙χ vanish is

K˜

I⁺ = −3β^rΩ⁰ α

I⁺

. (5.17b)

• Before obtaining the regularity conditions for the variables A_rr and Λ^r, it is conve-nient to do a separate calculation for the GBSSN and Z4c equations.

– GBSSN: The regularity condition for A_rr is computed from the quantities over Ω in ˙A_rr and then substituted into ˙Λ^r to obtain the expression for Λ^r in an equivalent way. The resulting conditions are:

Arr|^GBSSN_I+ = – Z4c: To calculate the regularity conditions corresponding to the Z4c equations,

first all Z_r have to be substituted using (2.79c). The A_rr and Λ^r variables involved in the regularity conditions appear in the singular terms in both ˙A_rr and ˙Λ^r, so that the system has to be solved together. The solution yields:

A_rr|^Z4c_I+ =

• The regularity condition required for equationsK˙˜ and Θ is the scri-fixing condition:˙˜

β^r2

• The only missing term in the evolution equation of the scalar field auxiliary variable Π is satisfied if˙˜

Π˜

I⁺ = 0. (5.17h)

5.2. Regularity of the equations at null infinity 89 Only if the relations (5.17) are satisfied, the formally singular terms in the equations will attain a regular limit at I⁺. The regularity conditions for the system (2.82) can be obtained substituting ˜K by ∆ ˜K using (2.80), for instance from (5.17b)

∆ ˜K

I⁺ =−K_{CM C} −3β^rΩ⁰ α

I⁺

. (5.18)

5.2.3 Regularity conditions arising from the gauge conditions

Extra regularity conditions may arise from the gauge conditions. To illustrate this, we will consider some examples from the actual gauge equations of motion that have been tested numerically and described in section 7.3.

• The gauge conditions where the lapse and shift are damped present the following behaviour:

– The harmonic slicing condition (4.26) with explicit source terms given by (7.8) and the 1+log slicing condition (4.24) with (4.25) and the source terms in (7.14) are such that the divergent terms in the equations can only cancel appropriately if, after using (5.18), the lapse satisfies

α|_I+ =−K_{CM C}r_I

3 . (5.19a)

– The same happens with the damping terms added in the Gamma-driver shift conditions (4.29) and (4.30) with source terms as in (7.12) and (7.13) respec-tively. If ξ_β^r 6= 0, then necessarily

β^r|_I+ = K_{CM C}r_I

3 . (5.19b)

The previous regularity conditions on the gauge quantities have the effect of fixing their values at I⁺, as will be described in chapters 7 and 8. This puts further constraints on the regularity conditions derived in the previous subsection, so that (5.17g) and (5.18) now reduce to

χ|_I+ = γ_rr|_I+ and ∆ ˜K

I⁺ = 0, (5.19c)

among other simplifications in the regularity conditions (5.17).

• The harmonic gauge conditions with source terms calculated from the background conformal metric components (4.38) have no terms divided by Ω. However, in subsection 7.3.9, I will describe the form of a source function that maintains the numerical simulation stable, namely ¯F_α^t = ^K_Ωα^{CM C}3 ( ˆα²−α²). The regularity condition induced by this term is the same as (5.19a). No regularity condition appears for the shift - it is allowed to move at I⁺, see figure 8.9.

• The harmonic gauge conditions with source terms including background physical metric components (4.43) do include terms that diverge at I⁺. After setting the

source functions as indicated in subsection 7.3.10 ( ˜F^r= 0 and ˜F^t=ξ_αΩ(α−α)) andˆ using the regularity conditions (5.17), the numerator of the Ω terms in ˙β^r vanishes, while the one in ˙α gives the relation

2β^rΩ⁰−ξ_αα³− 1

3α²ξ_αK_{CM C}r_I − 6β^r2γ_rrΩ⁰ K_{CM C}r_Iγ_θθ

I⁺

= 0. (5.19d)

In this case neither α nor β^r are fixed at I⁺, but their values have to satisfy this relation.

The regularity conditions have to be satisfied at I⁺ so that the formally divergent terms take finite limits there. The conditions also have to vanish with the appropriate power of Ω to avoid the appearance of stiff terms. If a staggered grid (see subsection 6.3.2) is used, a stable evolution will automatically satisfy the regularity conditions. Neverthe-less, in the case of a non-staggered grid, the conditions will have to be explicitly imposed in order to be able to evaluate finite values of the RHSs exactly onI⁺.

Chapter 6

Numerical methods

The main ingredients required for the implementation of the equations in the code are:

the discretization in space and time of the variables, the approximation of the spatial derivatives in the RHSs, the integration in time of the equations, boundary conditions, dissipation and the calculation of the initial data. The importance of the Courant factor and of convergence checks is also briefly discussed.

6.1 Spatial discretization

A continuum system of equations has to be discretized in order to be implemented and solved in a numerical code. There are several discretization procedures available for discretizing the spatial part. For smooth fields and simple geometries, excluding fluid dynamics, the pseudo-spectral methods and the finite difference approach are the most popular ones.

Pseudo-spectral methods, which are a variation of spectral methods, basically con-sist of writing the solutions to the differential equations as a sum of some given basis functions with time-dependent coefficients. The basis functions are selected depending on the properties of the problem that has to be solved and common options are Fourier series or Chebyshev polynomials. The number of terms in the sum is the same as the number of points in the discrete grid. The spectral method thus uses all the gridpoints in the interpolations, taking a global approach to the system. The coefficient equations are solved using common numerical techniques and then the solution is reconstructed.

Pseudo-spectral methods provide exponential convergence when the solution is smooth and the number of points they require is not very large. In spite of these advantages, pseudo-spectral methods are not optimal for our purposes, because their global character is likely to make them more difficult to stabilize, especially in presence of the divergent terms at I⁺ that appear in the equations we consider. For a simple description of the pseudo-spectral approach see [41].

The finite differences approach does not rely on a basis in the same way as the pseudo-spectral methods do (the solution is assumed to be represented by a polynomial around every gridpoint) and it only uses a certain number of gripoints in the interpolations (local approach). The convergence order is finite, but due to its robustness it will better fulfill our requirements.

91

6.1.1 Finite differences

The first step is to discretize the values of the variables to be integrated along a given number of points in the coordinates, or more generally, the variables they depend on.

Suppose a given functionu(t, x) is one of those integration variables. A grid containing a given number of points intandxis set and the functionu(t, x) is evaluated at those points.

Assuming that the chosen grid is equidistant (mesh refinement will not be considered here), the discretized quantities are

x_i = x₀+i∆x , i= 0,1, ..., I, (6.1a) t_n = t₀+n∆t , n= 0,1, ..., N, (6.1b)

uⁿ_i = u(t_n, x_i). (6.1c)

The spatial derivatives of the variables are approximated by a difference quotient. In the context of finite differences a difference quotient is given by a fractional expression that includes a linear combination of the values of the variable at a certain number of close points of the grid (called the stencil) in the numerator and a multiple of the grid-spacing to the same power as the order of the derivative in the denominator. In a more intuitive way, it is as if the derivatives were substituted by a discrete version of their definition.

For example, a first derivative is defined by

∂u

∂x = lim

∆x→0

u(x+ ∆x)−u(x)

∆x (6.2)

and the corresponding finite difference (first order and forward) with ∆x grid-spacing is

∂u

∂x i

= ui+1−ui

∆x +O(∆x)≈ ui+1−ui

∆x . (6.3)

Finite differences are approximations to the analytical derivatives, so that they have an associated error proportional to a certain power of the grid-spacing. The value of this power gives the convergence order of the finite differences. The higher it is, the faster the numerical solution will converge to the analytical one when the grid resolution is increased. On the other hand, the higher the convergence order, the more terms the stencil requires, which increases the calculation time.

The finite difference operators can be calculated using Taylor series, Lagrange poly-nomials or a simple algorithm derived by Fornberg [68]. As an example, the 4th order centered first and second derivatives are given by:

∂u

∂x i

= −u_i+2+ 8u_i+1−8u_i−1+u_i−2

12 ∆x +O(∆x⁴),

∂²u

∂x² i

= −u_i+2+ 16u_i+1−30u_i+ 16u_i−1−u_i−2

12 ∆x² +O(∆x⁴).

The stencils can be centered (taking the same number of points on both sides) or asymmetric. For the same order, the latter have larger errors. Depending on the fea-tures of the problem, one or the other may be more convenient. For instance, one-sided stencils may be useful at the boundaries, whereas using one-point off-centered stencils for

6.2. Method of Lines 93

In document Free evolution of the hyperboloidal initial value problem in spherical simmetry (sider 104-111)