{'fraction_non_alphanumeric': 0.03011217083558275, 'fraction_numerical': 0.007618019314978105, 'mean_word_length': 4.520529801324503, '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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dealing with the center and boundary problems in 1D Numerical Relativity

Instabilities in finite difference codes due to the singularity of spherical coordinates at the center are studied. In typical Numerical Relativity applications, standard regularization techniques by themselves do not ensure long term stability. A proposal to remedy that problem is presented, which takes advantage of redundant quantities introduced in recent hyperbolic formulations of Einstein's evolution equations. The results are discussed through the example case of a boson star, where a significant improvement in the implementation of boundary conditions is also presented.


I. INTRODUCTION
Spherically symmetric systems have deserved a lot of interest since the early stages of Numerical Relativity. In the pioneering work of May and White [1], the evolution of a selfgravitating spherical star was modelled by coupling the hydrodynamical equations with two first integrals of Einstein field equations: the energy and momentum constraints. These are two ordinary differential equations, which in spherical symmetry are enough to determine the evolution of the whole metric. The validity of this approach, however, is restricted to the spherical case, because in more general situations these first integrals are not enough to determine the solution of Einstein field equations. Furthermore, even in the spherical case, the treatment of boundary conditions with this system is known to be rather difficult. The latter is especially dangerous in systems with outgoing matter or radiation, which need to be carefully dealt with.
The alternative approach is to couple the hydrodynamical equations with the full set of Einstein partial differential equations. In that way one must evolve both the hydrodynamical quantities and the spacetime metric on the same footing. The resulting formalism can then be extended in principle to axially symmetric (2D) or even completely general situations. In the spherically symmetric case, however, the major drawback of this approach is not just the higher number of metric quantities to compute. It is rather the fact that the spherical coordinate system is singular at the origin, even when the spacetime itself is regular. The region close to the center of symmetry poses serious problems to numerical codes, which are to be overcome by using special techniques like artificial viscosity or algorithm matching. In most cases, this is usually done by trial and error and, together with the treatment of the outer boundary, becomes the main limiting factor when constructing finite difference codes.
The purpose of this paper is to present a new way of writing the full set of Einstein equations in spherical symmetry (with a matter content), which takes advantage of some redundant quantities that appear in recent hyperbolic formulations of Numerical Relativity [2,3], in order to solve this problem. In that way, one will be able to treat the center region with exactly the same numerical algorithms as any other region of the numerical grid, with a significative reduction of the time devoted to code development. The same formalism also allows us to impose more accurate external boundary conditions, related to the physics of the problem. The combination of these two improvements (center and boundary treatment) helps in increasing the algorithm stability and maintaining the required accuracy. As far as tested, the resultant code does not show any stability problem, even after very long runs. We will work out the description of the method through the case of a complex massive scalar field describing a radiating boson star, where we easily reproduce previous results from other authors [4].

II. FORMALISM
The object of Numerical Relativity is to solve the Einstein equations in an approximate way. This is necessary when we deal with some of the most interesting systems in astrophysics: black holes, neutron stars, close binaries, supernovae, etc.
G µν is the Einstein tensor associated to a certain spacetime (represented by a metric g µν ), G is the gravity constant and T µν is the momentum-energy tensor describing matter in this spacetime.
The metric element can be written, in normal coordinates, as where γ ij = g ij is the space metric and α is the lapse, a gauge degree of freedom.
The Einstein system (1) is constituted by 10 coupled second order partial differential equations, when written in terms of the 4-dimensional metric g µν . Six of them are evolution equations and the four remaining are constraint equations [5][6][7]. We deal with spherical symmetry problems, so the system reduces to a one-dimensional (1D) problem. In spherical coordinates, (2) is written where dΩ = dθ 2 + sin 2 θ dφ 2 .
The Einstein system can be written in first order form by introducing new variables. For simplicity, we will use the vacuum equations to illustrate our procedures. Following [2], we will express it as it follows: C ∂tgrr grr and q θ θ ≡ 1 C ∂tg θθ g θθ . Also, allowing for the momentum constraint, V r = 2D rθ θ . Note that we have written the Einstein system as a first order balance law system. The variable V r is redundant: it was introduced with the purpose of making the system hyperbolic, which allows us to directly deal with diagonalized variables with a well-defined speed and sense of propagation (very useful when imposing boundary conditions, as we will see later).
The function f (α) determines the gauge. We have chosen the gauge degree of freedom α to follow the harmonic slicing [8], which is the simplest choice (f = 1) and good enough for singularity-free spacetimes.

III. GEOMETRICAL FACTORS AND NUMERICAL ERROR
Due to spherical symmetry g θθ has a strong geometrical factor: it is proportional to r 2 when r → 0. It follows that both D rθ θ and V r are proportional to 1/r , a singular behaviour at the origin. Simulations dealing directly with these equations, or equivalent ones, are extremely unstable. This is well-known and some alternatives are usually taken. In previous works dealing with black holes, the problem was avoided by putting internal boundary conditions at the apparent horizon. In this way the grid was cut at the apparent horizon, far enough from both the physical and coordinative singularity at r = 0. The problem is completely different when we deal with singularity-free systems (stars). There, there is no way of cutting the solution before r = 0 without loosing physically relevant information. The coordinative singularity at r = 0 turns up as a big danger for accuracy and stability near the origin when using the standard system (4-12). The first step for dealing with this is to analytically extract the geometrical factors, as far as possible, from the equations, so they only deal with the regular part. This is done by writing the line element as − α 2 dt 2 + g rr dr 2 + r 2g θθ dΩ (13) so that the system, in our particular case, is rewritten with g θθ = r 2g θθ , D rθ θ = 1 r +D rθ θ , trD = D rr r + 2D rθ θ . The new equations read This is yet an improvement because the evolved quantitiesD rθ θ andg θθ are now regular at the center. But further improvement should be achieved, provided we want to deal with a stable code, because rewriting (4-12) with (13) introduces explicit 1/r factors in the terms containing space derivatives (flux terms) and 1/r 2 factors show up in the right hand side (source terms) of equations (11) and (12). We have somehow explicitly segregated the geometrical term, but it is still not out of the system. Now we can appreciate the kind of error that is introduced: both fluxes and sources are singular due to 1/r and 1/r 2 terms and exact cross-cancellation is needed between them to keep the system regular. However, this is not achieved by the finite difference algorithm due to the truncation error when calculating the space derivatives of the fluxes. The truncation error when computing a given variable with a second order code is proportional to (∆r) 2 , where ∆r is the distance between neighbouring points in the grid. The terms 1/r in the fluxes give errors like (∆r) 2 /r 4 . Where r → 0 (r ∼ ∆r) the code does not converge. As a particular instance we tried to solve the boson star 1D problem [4] by using these modified quantities. The numerical truncation error caused an instability that prevented us from obtaining any meaningful result before the code crash. These are well-known regularization techniques to handle this problem. In the next section, however, we present a new approach which regularizes the equations in a deeper way.

IV. THE ROLE OF THE MOMENTUM CONSTRAINT
Our proposal goes beyond the standard regularization procedures and succeeds in removing the most dangerous remaining 1/r factors. We take advantage from the way in which the momentum constraint was used to build the system (4-12). Remember that the variable V r is redundant and it was introduced with the objective of making the system hyperbolic [2]. As said, in (4-12) the momentum constraint is satisfied if and only if V r = 2D rθ θ . We can monitorize the level of violation of the momentum constraint by testing this algebraic identity through the evolution. From the considerations of the preceding section, it will be natural to evolve instead of V r the regular quantityṼ r ≡ 2D rθ θ . But it turns out to be more convenient to define instead which is also regular (provided that g rr /g θθ − 1 → 0 as r 2 , which demands smooth initial data near the origin) and removes both the 1/r terms in the fluxes and the 1/r 2 terms in the sources. The singularities have now disappeared from all the evolved variables as well as the numerical error caused by the incomplete cross-cancellation between fluxes and sources. The modified equations are There is still a point to be taken into account. In eqs. (19-20) the sources contain terms like 1/r times trD, L r or other variables that are radial derivatives of the metric coefficients. These terms do not pose any problem as r → 0, as radial derivatives of any differentiable function vanish at the origin. However, in the finite difference approach, we cannot put a grid point at r = 0 due to the explicit 1/r factors. We will therefore place the first grid point at r = ∆r/2. With this new way for writing the Einstein equations in spherical symmetry we are able to evolve the whole grid, without needing to apply special techniques or different algorithms to the center region. The equations are now intrinsically stable and can be directly coded.

V. EXAMPLE: THE BOSON STAR
We have applied this idea (17) to a Boson Star model [4]. Mathematically, the system is constituted by the Einstein equations plus the Klein-Gordon ones. It has stable configurations known as Boson Stars, which determine a mass-radius curve, as in the neutron star case, with a local maximum for the star masses. The stars located at the left (S-Branch) from this maximum (small radius) are known to be stable against small perturbations. The stars at the right (U-Branch) are, on the other hand, unstable. An U-Branch star, when perturbed can migrate towards a stable configuration in the S-Branch. The migration produces strong periodical oscillations in the star radius around the final position, as well as scalar field radiation. We tested our equations against this phenomenon. We could not evolve the system with (4-12) far enough to measure any relevant physical quantity. With (5-7, 14-16, 18-20), we have been able to reproduce the migration and the complex behaviour of its stationary states, with run time just limited by the required accuracy and not by algorithm stability, which is not easily achieved with the alternative formulation. We have to point out that part of the success comes from the hyperbolic balance law formulation itself, by allowing more physically consistent external boundary conditions. In the next section we go further in that point. In figure 1 we show the migration from an unstable configuration to a stable one [4]. With the original form of (4-12) we were unable to see nor half an oscillation.

VI. THE BOUNDARY CONDITIONS
We illustrate the method for boundary conditions by describing the treatment of the Klein-Gordon equations in our example system. Let's take a look at the equation for the real part, φ 1 , of the Klein-Gordon field (g µν φ ;µν = m 2 φ, in tensor form): where ρ ≡ ∂ r φ, ψ ≡ 1 C ∂ t φ. These equations, like the rest of the system, are written in balance law form. The flux terms are those that would appear in a flat space wave equation (transport terms). The m 2 term in the right hand side is an oscillatory term which comes from the m 2 φ in the Klein-Gordon equation. The rest of the source terms are due to spherical coordinates and the coupling between the Klein-Gordon field and the spacetime metric. To put a condition over the transport terms would be quite easy, but the combination with the source terms makes the problem more difficult. What one would like to do is to treat separately each term of the equation. Fortunately, this is possible in a numerical code, by using the 'operator splitting' approach. With the Strang splitting (second order operator splitting) each step in the numerical evolution can be decomposed in the following way: The whole step, ǫ(∆t), is split into three sequential processes. The first one is S( ∆t 2 ). There we evolve during half a step only the 'source' terms in the equations. If we restrict to the Klein-Gordon part, we have: This is an ordinary differential equation system for each point in the net. In this part there is no connection between points. Each point evolves by itself. So we do not have to state any boundary condition in this substep.
In the second substep, R(∆t), we only evolve the transport part (during a whole time step), so we have The remaining terms here are the wave equation-like ones. We know that the solution is formed by two components (ρ + φ and ρ − φ), which propagate in opposite senses at speed C (in general we need the matter part to be written in hyperbolic form, so it can be diagonalized). Now here is where we need to put a boundary condition. What we want is: first, the outgoing component not to be touched; any condition over it would be superfluous or contradictory.
Simply let it go out. Second, as we suppose that there is no scalar field outside the net, nothing is coming in (we disregard nonlinear effects between outgoing radiation and external metric), so the boundary condition to impose is simply The third step is again S( ∆t 2 ). This is all. The only boundary condition we need to put is (28), which is very simple indeed. Remember that there is not gravitational radiation in spherical symmetry. The momentum constraint plus the coordinate gauge condition fully determine the remaining (metrical) part of the boundary conditions. Note that the relation of D rθ θ with V r holds when the momentum constraint is satisfied, so we can use the computed values of V r at the boundary to obtain the value of the incoming non-gauge field there.

VII. CONCLUSIONS
The stability and accuracy near the origin of spherically symmetric dynamical systems in Numerical Relativity is strongly improved by a new formalism which provides an automatic care-free treatment of the geometrical factors arising from spherical symmetry. These treatment involves making use of the redundant variable V r , previously introduced for other purposes (to ensure the hyperbolicity of the system). The same formalism also allows boundary treatment based on physical considerations and helps in increasing accuracy and stability. That is, the same characteristic of our system, i.e. the introduction of the new quantity V r , allows us to simultaneously improve the treatment of the two main problems in spherically symmetric star-like systems: the center (through a redefinition of V r , which brings the system to a regular form) and the boundary (through hyperbolicity achieved thanks to V r , which lets us distinguish between incoming and outcoming fields and treating them according to physical criteria).
We have started from equations (4-12), which are just the 1D version of a general 3D system [2,3] which is actually being applied in multidimensional Numerical Relativity. Our point is that this system, when rearranged to deal with the geometrical factors, works fine in the 1D case. We are presently considering the axially symmetric case (2D), where again the 'redundant' quantities V i could be helpful to get rid of the geometrical factors that arise from the coordinate singularity at the symmetry axis. But this is still an open question.