The Einstein equations form a non-linear system, which makes the study of stability of their implementation especially difficult. If an instability arises in any of the equations, it may contaminate the evolution of the rest of variables very fast due to the non-linear cou-plings between the equations. This makes the task of finding the origin of the instability very difficult.
Stability does not only depend on the evolution equations themselves, but also on the background around which the evolution takes place. A stationary background (like flat or Schwarzschild spacetimes) is relevant for our purposes and also a simple choice.
The initial data will consist of perturbations of the given stationary solution. Choosing 105
a stationary background also allows to study a linearized version of the equations around it. To do so, all variables in the equations are expressed as a Taylor expansion to first order around the given background. An unstable linearized system is unlikely to become stable in its non-linear form, so first trying the linearized system can be a useful approach, because it is a simpler problem to solve.
When studying the equations in search of the origin of the instabilities, it is very helpful to divide the system of equations into smaller parts and analyze those. The smallest division is to consider the single equations, then subsystems of two or three equations will be constructed and analyzed.
7.2.1 Analysis of the single equations
The simplest subsystem that can be analyzed are the single equations. For each equation in (2.76) or (2.78), as well as the gauge equations, all quantities except the one that is being evolved are substituted by their stationary solutions. For simplicity in the analysis we will suppose that the stationary value of the variable under study, that will be denoted byA, is zero. This can be assumed without loss of generality, because with the appropriate variable transformation this can be achieved for any evolved quantity. A model equation forAincluding the relevant terms (it is actually a linearization of the complete equation) is
A˙ =aA0 + b+ c
Ω
A, (7.1)
where the coefficients a, b and c are time-independent functions of the radius. Higher order terms in A may appear in the single equations, but they do not have a relevant effect on stability. This equation serves to model the behaviour of the quantitiesArr, Λr, K or ˜K and possibly α andβr, which are the ones more likely to pose stability problems.
Regular equations of motion (without divergent terms), like the ones forγrr orγθθ, should not require any numerical tests.
The objective of the study is to detect exponential growths. They will arise when the coefficient b+Ωc
> 0, so that simple inspection of the values taken by b and c will give us useful information. The advection term with coefficient a only propagates information on the grid and does not play an important role for long-term solutions of the equations where all perturbations have already been radiated away, so it does not need to be considered in the stability analysis. If c > 0, the growth at I+ will be very fast and the simulation is very likely to crash after only a few time-steps.
If b >0 (or c > 0), one of the checks that should be performed on the equation is to see if there is a stable stationary solution of the equation, because A = 0 turns out not to be one. This can be easily done by plotting ˙A as a function of A (for specific values of the radius and dropping the advection term) keeping all of the non-linear terms. If there is another A1 6= 0 where ˙A(A1) = 0 and the slope of ˙A is negative there, A1 is a stable solution of the equation.
The mentioned stability conditions on the coefficients are a property of the continuum equations. At the numerical level there are also important stability conditions to satisfy.
One is the value of the Courant factor defined in (6.8), which has to satisfy the CFL limit condition. In presence of stiff terms it may even be smaller than expected. The initial perturbations to test the stability of the equations are chosen small enough so that
7.2. Steps towards stability 107 they do not “break” the coordinate system. It is useful to first select small perturbation amplitudes to avoid exciting the non-linear effects in the equations.
7.2.2 Analysis of subsystems of equations
Once the single equations present stable and appropriate numerical behaviour, a hierar-chical way of constructing a well-behaved complete system consists of first testing small groups of equations. Each single equation of (2.76), (2.78) or (2.82) by itself is hyperbolic, but not every combination of equations taken from (2.76), (2.78) or (2.82) is. With the purpose of reducing the number of variables as much as possible, it is useful to fix the determinant freedom of the conformal spatial metric and eliminate γθθ in terms of γrr using (2.71), to solve the vacuum equations (not evolving the scalar field), to use a fixed shift and to consider the GBSSN equations, not evolving the Z4 quantity Θ.
Another requirement of our setup is that the values of the eigenspeeds have to be larger or equal to zero atI+. If a BH is present in the simulated spacetime and excision is being used, then the eigenspeeds have to be negative or vanishing at the horizon.
For the possible subsystems that can be constructed from (2.76) (and appropriate lapse condition), the resulting characteristic speeds satisfy these conditions.
The hyperbolic subsystems of (2.76) or (2.78) can be divided into two groups:
• Subsystems that involve the trace-free part of the extrinsic curvature (Arr in our ansatz) and the difference of contracted connections (Λr): certain combinations of the metric component γrr (and γθθ if not eliminated) with Arr and/or Λr are hyperbolic. The conformal factor χ can also be added to any of them without breaking hyperbolicity. The possible hyperbolic combinations are listed in the first part of subsection 7.3.2.
• Subsystems that involve the trace of the extrinsic curvature (K): the equation of motion of K (or equivalently ˜K) only includes α00 as second derivatives of the metric components, so that in our setup it can only couple hyperbolically to the lapse evolution equation. The possible hyperbolic choices are αK and αχK. In the same way as with the single equations, the hyperbolic subsystems have to be numerically checked for exponential growths of any of the variables.
7.2.3 Tracking down instabilities
Possibly the more powerful tool to locate the origin of instabilities is by visual inspection of the numerical data. In order to simplify the visualization, it is useful to have very simple stationary values (like 1 or 0) for the variables. This can be easily achieved by simple variable transformations, see next subsection. Another aid in tracking down the origin of the instability is plotting the RHSs of the equations in time and search for the first exaggerate growth.
If an exponential growth has been detected, plotting the data in a logarithmic scale allows to estimate the slope of the resulting straight line. Knowing its value can help identify which term in the equation is causing the exponential growth. Another way of trying to find where the instability has its source or checking if an existing guess is correct
is to add damping terms of the form−λA or−ΩξAwith λ, ξ ≥0 to the corresponding ˙A.
Obviously, the solution of the resulting set of equations no longer satisfies the Einstein equations, but the terms added by hand serve as a trick to find and understand the instabilities.
Checking if the regularity conditions listed in subsections 5.2.2 and 5.2.3 are satisfied by the numerical data during the evolution can also provide hints about the origin of the instabilities.
7.2.4 Variable transformations
Certain well-chosen variable transformations can help solve some of the continuum in-stabilities present in the equations. Some transformations change the principal part of the system, while others do not. To see this, consider the u and v variables defined in subsection 5.1.2. Let us denote a given u variable as Au and a given v variable as Bv. The possible variable transformations can be classified as follows:
• The principal part remains the same after performing transformations of the form Bv → Bv +f(u), where f is the function that defines the variable transformation and u denotes any u-type variables involved in the transformation. The principal part is left unchanged because the terms arising from the mixing with theuvariables are all non-principal part terms.
• The principal part is likely to change under transformations of the form Au → f(Au, u),Bv →f(Bv, v) andBv →f(Bv, u0), because the same type of variables (in the sense ofu, v andw≡u0) are involved and thus the spatial derivatives that form the principal part will also change. However, as these are similarity transformations, even if the eigenvectors of the system may change, the eigenvalues and hyperbolicity properties of the system remain the same.
7.2.5 Constraint damping
At the continuum level the constraints are satisfied exactly, so that adding a multiple of them to any evolution equation gives a new system of equations that is exactly equivalent to the original one. In a numerical simulation the constraints are only satisfied to some extent, there is always some numerical error. However, if the value of the constraints con-verges to zero, this manipulation of the equations is also valid for the evolution equations implemented in the code. Such a manipulation is performed for instance in the BSSN formulation, where a multiple of H is added to ˙¯K’s RHS to eliminate the Ricci scalar, while in the Z4 formulation the same effect is achieved with the mixing of ¯K and Θ in (2.48).
Adding a multiple of a constraint can be useful to eliminate a problematic term in an equation of motion or to enforce the vanishing of the given constraint. The terms withκ1 in (2.1) are the constraint damping terms [83] of the Z4 formulation and tend to minimize the value of the constraint variables Θ andZa. It can also help eliminate terms that cause exponential growths. Nevertheless, the resulting system of equations has to be hyperbolic and must have the correct eigenspeeds, properties which may change with the new terms brought in by the constraint expression.