{'fraction_non_alphanumeric': 0.03248640190745846, 'fraction_numerical': 0.00804709038074659, 'mean_word_length': 4.534845360824742, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Flux Limiter Methods in 3D Numerical Relativity

New numerical methods have been applied in relativity to obtain a numerical evolution of Einstein equations much more robust and stable. Starting from 3+1 formalism and with the evolution equations written as a FOFCH (first-order flux conservative hyperbolic) system, advanced numerical methods from CFD (Computational Fluid Dynamics) have been successfully applied. A flux limiter mechanism has been implemented in order to deal with steep gradients like the ones usually associated with black hole spacetimes. As a test bed, the method has been applied to 3D metrics describing propagation of nonlinear gauge waves. Results are compared with the ones obtained with standard methods, showing a great increase in both robustness and stability of the numerical algorithm.


Introduction
From the very beginning, 3D numerical relativity has not been an easy domain. Difficulties arise either from the computational side (the large amount of variables to evolve, the large number of operations to perform, the stability of the evolution code) or from the physical side, like the complexity of the Einstein equations themselves, boundary conditions, singularity avoidant gauge choices, and so on. Sometimes there is a connection between both sides. For instance, the use of singularity avoidant slicings generates large gradients in the vicinity of black holes. Numerical instabilities can be produced by these steep gradients. The reason for this is that the standard evolution algorithms are unable to deal with sharp profiles. The instability shows up in the form of spurious oscillations which usually grow and break down the code.
Numerical advanced methods from CFD (Computational Fluid Dynamics) can be used to avoid this. Stable codes are obtained which evolve in a more robust way, without too much dissipation, so that the shape of the profiles of the evolved quantities is not lost. These advanced methods are then specially suited for the problem of shock propagation, but they apply only to strongly hyperbolic systems, where one is able get a full set of eigenfields which generates all the physical quantities to be evolved. In the 1D case, these methods usually fulfill the TVD (Total Variation Diminishing) condition when applied to transport equations. This ensures that no new local extreme appear in the profiles of the eigenfields, so that spurious oscillations are ruled out ab initio (monotonicity preserving condition). Unluckily, there is no general method with this property in the 3D case, mainly because the eigenfield basis depends on the direction of propagation. We will show how this can be achieved at least in some cases.
The specific methods we will use are known as flux limiter algorithms. We will consider plane waves in 3D as a first generalization of the 1D case, because the propagation direction is constant. This specific direction leads then to an specific eigenfield basis, so that the 1D numerical method can be easily generalized to the 3D case.
The algorithm will be checked with a "Minkowski waves" metric. It can be obtained by a coordinate transformation from Minkowski space-time. All the metric components are transported while preserving their initial profiles. The line element has the following form: where H(x − t) is any positive function. We can choose a periodic profile with sharp peaks so both the space and the time derivatives of H(x − t) will have discontinuous step-like profiles. If we can solve well this case (the most extreme), we can hope that the algorithm will work as well in more realistic cases where discontinuities do not appear. Minkowski waves are a nice test bed because the instabilities can arise only from the gauge (there are pure gauge after all!). This is a first step to deal with evolution instabilities in the Einstein equations by the use of flux limiter methods. This will allow us to keep all our gauge freedom available to deal with more physical problems, like going to a co-rotating frame or adapting to some special geometry. Advanced numerical methods take care of numerical problems so that 'physical' gauge choices can be used to take care of physics requirements.

The system of equations
We will use the well known 3+1 description of spacetime [1,2,3] which starts by decomposing the line element as follows: where γ ij is the metric induced on the three-dimensional slices and β i is the shift. For simplicity the case β i = 0 (normal coordinates) will be considered. The intrinsic curvature of the slices is then given by the three-dimensional Ricci tensor (3) R ij , whereas their extrinsic curvature K ij is given by: In what follows, all the geometrical operations (index raising, covariant derivations, etc) will be performed in the framework of the intrinsic three-dimensional geometry of every constant time slice. With the help of the quantities defined in (2,3), the ten fourdimensional field equations can be expressed as a set of six evolution equations: plus four constraint equations The evolution system (4) has been used by numerical relativists since the very beginning of the field (see for instance the seminal work of Eardley and Smarr [4]), both in spherically symmetric (1D) and axially symmetric (2D) spacetimes. By the turn of the century, the second order system (4) has been rewritten as a first-order flux conservative hyperbolic (FOFCH) system [5,6,7] in order to deal with the generic (3D) case, where no symmetries are present. But the second order system (4) is still being used in 3D numerical calculations [8], mainly when combined with the conformal decomposition of K ij as introduced by Shibata and Nakamura [9,10]. In the system (3,4) there is a degree of freedom to be fixed because the evolution equation for the lapse function α is not given. In the study of Black Holes, the slicing is usually chosen in order to avoid the singularity [11]: where: Three basic steps are needed to obtain a FOFCH system from the ADM system. First, one must introduce some new auxiliary variables to express the second order derivatives in space as first order. These new quantities correspond to the space derivatives: The evolution equations for these variables are: At the second step the system is expressed in a first order balance law form where the array u displays the set of independent variables to evolve and both "fluxes" F k and "sources" S are vector valued functions. At the third step another additional independent variable is introduced to obtain a strongly hyperbolic system [11]: and its evolution equation is obtained using the definition of K ij from (3) and switching space and time derivatives in the momentum constraint (6). The result is an independent evolution equation for V i while the previous definition (13) in terms of space derivatives can be instead be considered as a first integral of the extended system. The extended array u will then contain the following 37 functions u = (α, γ ij , K ij , A i , D kij , V i ).

The numerical algorithm
Due to the structure of the equations, the evolution(represented by the operator E(∆t)) described by (12) can be decomposed into two separate processes; the first one is a transport process and the second one is the contribution of the sources. The sources step (represented by the operator S(∆t)) does not involve space derivatives of the fields, so that it consists in a system of coupled non-linear ODE (Ordinary Differential Equations): The transport step (represented by the operator T (∆t)) contains the principal part and it is given by a set of quasi-linear transport equations: The numerical implementation of these separated processes is quite easy. Second order accuracy in ∆t can be obtained by using the well known Strang splitting.

E(∆t) = S(∆t/2) T (∆t) S(∆t/2)
According to (3,7) the lapse and the metric have no flux terms. It means that a reduced set of 30 quantities u = (K ij , A i , D kij , V i ) are transported in the second step over an inhomogeneous static background composed by (α, γ ij ). The equations for the transport step (15) are given by: where: and m is an arbitrary parameter.
To evolve the transport step, we will consider flux-conservative numeric algorithms [12], obtained by applying the balance to a single computational cell. In the 1D case the cell goes from n to n+1 in time (t = n · ∆t) and from j-1/2 to j+1/2 in space (x j = j · ∆x), so that we have: Interface fluxes can be calculated in many different ways, leading to different numerical methods. We will use here the well known MacCormack method. This flux-conservative standard algorithm works well for smooth profiles, as it can be appreciated in Figure 1. But this standard algorithm is not appropriate for step-like profiles because it produces spurious oscillations near the steep regions, as it can be appreciated in Figure 2.   Fig. 1 with the step-like initial data for K xx . Continuous line is the initial condition. Dashed line is after 10 iterations. Note the spurious oscillations around the corners More advanced numerical methods must be used to eliminate (or at least to reduce) these oscillations. These advanced methods use information about the eigenfields and the propagation direction, so the flux characteristic matrix along the propagation direction must be diagonalized.

Eigenfields
We will use a convenient method to compute the eigenfields. Let us study the propagation of a step-like discontinuity in the transported variables u which will move along a specific direction n with a given velocity v. Information about the corresponding eigenfields can be extracted from the well known Rankine-Hugoniot shock conditions : where [ ] represents the jump in the discontinuity. In our case where we must note that both the background metric coefficients and the propagation direction are supposed to be continuous, so they are transparent to the [ ] symbol. If we develop this expressions we arrive at the following conclusions, where S n = n r S r is the projection of the quantity S over n and S ⊥ = S k − S n n k are the transverse components: and [A n − f trD n ] propagate along n with speed v = 0. There are 18 such eigenfields. For the line element given by (1) n k is along the x axis and all these fields are actually zero.
2) [λ n ij − trλ n n i n j ] and [K ij − trKn i n j ] do generate eigenfields propagating along n with speed v = ±α (light cones). There are only 10 such eigenfields because all of them are traceless. For Minkowski waves, where there is only gauge, all these combinations are zero. This indicates that the correct way to get the traceless part of a given tensor S ij in this context is just to take S ij − trS n i n j , so that the contribution of gauge modes will disappear.
3) [A n ] and [trK] do generate eigenfields propagating along n with speed v = ± √ f α (gauge cones). There are 2 such eigenfields corresponding to the gauge sector. For Minkowski waves, there are the only non-zero components. We are left with: so that [A k ] is proportional to n k . Now we can get the gauge eigenfields: These eigenfields propagate along n according simple advection equations, a familiar situation in the 1D case. Although this decomposition and diagonalization is trivial in 1D, it is very useful in the multidimensional case for a generic direction n.

Flux limiter methods
The flux limiter methods [12] we will use can be decomposed into some basic steps. First of all the interface fluxes have to be calculated with any standard second order accurate method (MacCormack in our case). Then, the propagation direction n and the corresponding eigenfields can be properly identified at every cell interface. Two advection equations (one for every sense of propagation) are now available for the gauge eigenfluxes (30).
Let us choose for instance the eigenflux which propagates to the right (an equivalent process will be valid for the other eigenflux propagating to the left). This interface eigenflux F n+1/2 j+1/2 can be understood as the grid point flux F n j plus some increment ∆ j+1/2 = F n+1/2 j+1/2 − F n j . In general, the purpose of the limiter is to use of a mixture of the increments ∆ j+1/2 and ∆ j−1/2 to ensure monotonicity. In our case we are using a robust mixture which goes by applying the well known minmod rule to ∆ j+1/2 and 2∆ j−1/2 . In that way, the limiter acts only in steep regions, where the proportion between neighbouring increments exceeds a factor of two.
We can apply this method to the step-like initial data propagating along the x axis. We can see in Figure 3 that the result is much better than before. It can be (hardly) observed a small deviation from the TVD condition, which is produced by the artificial separation produced by the Strang splitting into transport and non-linear source steps.
This method can be applied, with the general decomposition described in section 4, to discontinuities which propagate along any constant direction, and not only to the trivial case, aligned with the x axis, that we have considered until now. To prove it, we have rotated the metric of Minkowski waves in the x-z plane to obtain a diagonal propagation of the profile. The line element in this case has the following form: We show the results in the Figure 4. We can also see in Figure 5 a z=constant section of the same results to allow a more detailed comparison with the initial data.