Spherical symmetry as a test case for unconstrained hyperboloidal evolution II: gauge conditions

TL;DR

This paper investigates hyperboloidal evolution of a self-gravitating scalar field in spherical symmetry, exploring various gauge conditions and conformal gauges, including null infinity placement, using BSSN and Z4 formulations.

Contribution

It introduces numerical solutions with flexible gauge conditions and null infinity placement in hyperboloidal evolution, extending previous work by including the option of grid points at null infinity.

Findings

Robust evolutions achieved for various gauge choices.

Including null infinity in the grid is feasible and effective.

Different conformal gauges impact the singularity handling.

Abstract

We present numerical solutions of the hyperboloidal initial value problem for a self-gravitating scalar field in spherical symmetry, using a variety of standard hyperbolic slicing and shift conditions that we adapt to our hyperboloidal setup. We work in the framework of conformal compactification, and study both evolutions that employ the preferred conformal gauge, which simplifies the formal singularities of our equations at null infinity, and evolutions without this simplification. In previous work we have used a staggered grid, which excludes null infinity, while now we include the option of placing a gridpoint directly at null infinity. We use both the generalized BSSN and conformal Z4 formulations of the Einstein equations, study the effect of different gauge conditions, and show that robust evolutions are possible for a range of choices.

Figures

Equations

$${\bar{g}}_{\mu \nu }\equiv {\Omega }^2{\tilde{g}}_{\mu \nu }.$$

$$G_{\mu \nu }[\bar{g}]=8\pi T_{\mu \nu }-\frac{2}{\Omega }\left({\bar{\nabla }}_{\mu }{\bar{\nabla }}_{\nu }\Omega -{\bar{g}}_{\mu \nu }{\bar{\nabla }}^{\gamma }{\bar{\nabla }}_{\gamma }\Omega \right)-\frac{3}{{\Omega }^2}{\bar{g}}_{\mu \nu }({\bar{\nabla }}_{\gamma }\Omega ){\bar{\nabla }}^{\gamma }\Omega .$$

$${\bar{\square }\Omega |}_{\text{mathrsfs}{I}^{+}}=0.$$

$${{\bar{\nabla }}_{\mu }{\bar{\nabla }}_{\nu }\Omega |}_{\text{mathrsfs}{I}^{+}}=0,$$

$$\Omega (r)=\left(-K_{CMC}\right)\frac{r_{\text{mathrsfs}I}^2-r^2}{6r_{\text{mathrsfs}I}},$$

$${-{\bar{g}}_{tt}|}_{\text{mathrsfs}{I}^{+}}={\left({\alpha }^2-{\chi }^{-1}{\gamma }_{rr}{{\beta }^r}^2\right)|}_{\text{mathrsfs}{I}^{+}}=0,$$

$${\hat{\square }\hat{\Omega }|}_{\text{mathrsfs}{I}^{+}}=0.$$

$${\left[\left({\bar{\nabla }}^a\Omega \right)\left({\bar{\nabla }}_a\ln \omega \right)\right]|}_{\text{mathrsfs}{I}^{+}}={-\frac{\bar{\square }\Omega }{4}|}_{\text{mathrsfs}{I}^{+}},$$

$${\frac{\dot{\omega }}{\omega }|}_{\text{mathrsfs}{I}^{+}}={\frac{{\beta }^r{\bar{g}}_{tt}^{\prime }}{2{\alpha }^2}|}_{\text{mathrsfs}{I}^{+}},$$

$$d{t}_B=\frac{{\alpha }^2\omega }{{\beta }^r{\Omega }^{\prime }}dt.$$

$${{\alpha }^2|}_{\text{mathrsfs}{I}^{+}}={{\beta }^r{\Omega }^{\prime }|}_{\text{mathrsfs}{I}^{+}}.$$

$$\dot{\alpha }={\beta }^r{\alpha }^{\prime }-{\alpha }^{2}f(\alpha )\left(K-K_0\right)+L_0.$$

$$\dot{\alpha }={\beta }^r{\alpha }^{\prime }-{\alpha }^2\left(K-K_0\right)+L_0,$$

$$\dot{\alpha }={\beta }^r{\alpha }^{\prime }-n_{1+log}\alpha \left(K-K_0\right)+L_0.$$

$$n_{1+log}={\alpha |}_{\text{mathrsfs}{I}^{+}}=-\frac{K_{CMC}{r}_{\text{mathrsfs}I}}{3}.$$

$$\dot{\alpha }={\beta }^r{\alpha }^{\prime }-n_{cK}\left(K-K_0\right)+L_0.$$

$$n_{cK}={{\alpha }^2|}_{\text{mathrsfs}{I}^{+}}={\left(-\frac{K_{CMC}{r}_{\text{mathrsfs}I}}{3}\right)}^2.$$

$$\dot{\tilde{\alpha }}=$$

$$\hat{\tilde{\chi }}=\widehat{{\gamma }_{\theta \theta }}=1,\widehat{{\gamma }_{\tilde{r}\tilde{r}}}=\frac{1}{{\hat{\tilde{\alpha }}}^2},\hat{\tilde{\alpha }}=\sqrt{1+{\left(\frac{K_{CMC}\tilde{r}}{3}\right)}^2}\text{and}\widehat{{\beta }^{\tilde{r}}}=\frac{K_{CMC}\tilde{r}}{3}\sqrt{1+{\left(\frac{K_{CMC}\tilde{r}}{3}\right)}^2}.$$

$$\dot{\alpha }=$$

$$\hat{\chi }=\widehat{{\gamma }_{rr}}=\widehat{{\gamma }_{\theta \theta }}=1,\hat{\alpha }=\sqrt{{\Omega }^2+{\left(\frac{K_{CMC}r}{3}\right)}^2}\text{and}{\hat{\beta }}^r=\frac{K_{CMC}r}{3},$$

$$\dot{\alpha }=$$

$${\partial }_{\perp }{\beta }^a=\frac{3}{4}{B}^a,{\partial }_{\perp }{B}^a=\lambda {\partial }_{\perp }{\Gamma }^a-\eta B^a\text{with}{\partial }_{\perp }={\partial }_t-{\beta }^i{\partial }_i,$$

$${\Lambda }^a={\Gamma }^a-{\gamma }^{bc}{\hat{\Gamma }}_{bc}^a(+2{\gamma }^{ab}{Z}_b\text{in the Z4 formulation}),$$

$$\dot{{\beta }^r}=$$

$$\dot{B^r}=$$

$$\dot{{\beta }^r}=$$

$${\dot{\beta }}^r={\beta }^r{{\beta }^r}^{\prime }+{\alpha }^2\chi {\Lambda }^r+\frac{{\alpha }^2{\chi }^{\prime }}{2{\gamma }_{rr}}-\frac{\alpha \chi {\alpha }^{\prime }}{{\gamma }_{rr}}-\frac{2{\alpha }^2\chi }{r{\gamma }_{\theta \theta }}+L_0-\frac{{\xi }_{{\beta }^r}}{\Omega }{\beta }^r.$$

$${\dot{\beta }}^r={\beta }^r{{\beta }^r}^{\prime }-\widehat{{\beta }^r}{\widehat{{\beta }^r}}^{\prime }+{\alpha }^2\chi {\Lambda }^r+\frac{{\alpha }^2{\chi }^{\prime }}{2{\gamma }_{rr}}-\frac{\alpha \chi {\alpha }^{\prime }}{{\gamma }_{rr}}+\hat{\alpha }{\hat{\alpha }}^{\prime }+{\xi }_{{\beta }^r}\left(\frac{\widehat{{\beta }^r}}{\Omega }-\frac{{\beta }^r}{\Omega }\right).$$

$$f_{match}=\frac{3r_{\text{mathrsfs}I}}{(-K_{CMC})}\Omega =(r_{\text{mathrsfs}I}^2-r^2)$$

$$\dot{\alpha }:$$

Full text

1 Introduction

The energy loss and radiation of isolated systems are in general only well defined at future null infinity ($\mathrsfs{I}^{+}$). This is where we can consider observers of astrophysical events to be located [1, 2, 3], and where we would ideally want to extract radiation signals from numerically constructed spacetimes. In order to treat a physical system of infinite extent, we follow Penrose’s approach [4, 5], and conformally compactify spacetime: the physical metric $\tilde{g}_{\mu\nu}$, which diverges at infinity, is rescaled by a conformal factor $\Omega$ that is chosen to vanish at $\mathrsfs{I}^{+}$ to construct a finite metric $\bar{g}_{\mu\nu}$:

[TABLE]

The Einstein field equations written in terms of the rescaled metric $\bar{g}_{\mu\nu}$ formally diverge at $\mathrsfs{I}^{+}$:

[TABLE]

In our work we will solve an initial value problem where spacetime is foliated by smooth spacelike slices that reach $\mathrsfs{I}^{+}$, known as hyperboloidal slices. The hyperboloidal initial value problem was pioneered by Friedrich, who developed a reformulation of the conformally rescaled Einstein equations (2) in terms of a regular system of equations, the “Conformal Field Equations” [6, 7, 8]. Due to general relativity’s character as a gauge theory, approximation schemes that solve the Einstein equations as an initial value problem tend to suffer from instabilities in the continuum equations. As a consequence, small numerical errors can trigger growing modes in the gauge or constraint degrees of freedom, see e.g. [9, 10, 11].

In recent years, numerical simulations of mergers of compact objects have provided crucial information about the gravitational wave signals of such events, which was essential to identify the sources the signals discovered since late 2015 by the LIGO and Virgo detectors [12, 13, 14, 15]. In these simulations, Cauchy slices are evolved up to some finite distance from the source, and extrapolation methods or characteristic extraction [16, 17] are used to estimate the signal at future null infinity. Such simulations have only been possible after the breakthroughs of 2005 [18, 19, 20], which have established gauge conditions that are appropriate in particular for the strong field region in the vicinity of black holes, and avoid instabilities that had plagued the field of numerical relativity for many years. Our goal here is to develop a simple implementation of the hyperboloidal initial value problem: starting with a formulation close to those routinely used to study astrophysical situations like the coalescence of compact objects (such as BSSN and Z4 [21, 22, 23]), we then modify it to resolve instabilities rooted in the continuum formulation. Since the equations are formally singular, numerical schemes to solve the conformally rescaled equations have to be constructed to appropriately regularize the singular terms. A key idea in simplifying such regularization procedures is an appropriate choice of gauge, which is adapted to the inherent simplicity of null infinity. In order to achieve this goal, we face problems and instabilities similar to those that previously plagued the numerics of the strong field region. For a discussion of how such problems are resolved by current gauge conditions on Cauchy evolutions, see e.g. [24, 25]. The problems we face here are mainly due to the compactification of null infinity. As for the problems associated with black holes, the essential problems can already be captured in spherical symmetry, which motivates using spherical symmetry in our present work. While in the full 3-dimensional case new problems may appear, we do expect that resolution of the problems we address in the spherical symmetric context will carry over to the full 3-dimensional case. The extraction of gravitational wave signals at future null infinity is one of the main final applications of our work. As gravitational waves cannot arise in spherical symmetry, we use a massless scalar field to introduce physical dynamics and mimick gravitational wave effects, so that we can extract and analyse signals in our spherically symmetric implementation. The coupling of a scalar field with the Einstein equations has been widely studied in relation with its critical collapse [26], using asymptotically flat instead of hyperboloidal slices, where the scalar field does not introduce particular difficulties apart from requiring high resolution to resolve critical phenomena.

A natural class of simple choices for the coordinates and conformal factor at null infinity is known as the “preferred conformal gauge” [27, 28, 29], where the conformal factor $\Omega$ satisfies the wave equation with respect to the unphysical rescaled metric $\bar{g}_{\mu\nu}$,

[TABLE]

With this choice one can show that also

[TABLE]

so that the three divergent terms in (2) attain a finite limit at $\mathrsfs{I}^{+}$ independently, and consequently the problem to find a regularized and stable numerical scheme for the hyperboloidal initial value problem simplifies considerably. The preferred conformal gauge has been used to derive boundary conditions at null infinity for elliptic-hyperbolic systems, see [30, 31, 32], including constrained evolution schemes in a tetrad formalism [33, 34], and alternatively to define gauge conditions for a hyperbolic free evolution scheme in the approach by Zenginoğlu [35, 36, 37]. The preferred conformal gauge ensures that the null generators of $\mathrsfs{I}^{+}$ are geodesic and expansion-free (they are already shear-free in a conformally invariant way). A relevant consequence of the expansion-free null generators of $\mathrsfs{I}^{+}$ is that their affine parameter can be identified with the time coordinate at $\mathrsfs{I}^{+}$, which is then commonly called Bondi time and corresponds to the proper time measured by an inertial observer moving along $\mathrsfs{I}^{+}$ [7]. The Bondi time parameterization can be chosen separately from the preferred conformal gauge, and such a choice allows to extract signals at $\mathrsfs{I}^{+}$ without any deformations caused by non-inertial coordinate effects.

Pioneering work toward a free evolution scheme for the hyperboloidal initial value problem was performed by Zenginoğlu [38], who carried out evolutions of Schwarzschild initial data in a spherically symmetric code, using the generalized harmonic form of the conformally rescaled Einstein equations (2), and in particular the generalized harmonic coordinate gauge on the conformally compactified manifold. The latter included source functions that satisfied the preferred conformal gauge (3), although it was not clear how to impose the Bondi time parameterization at $\mathrsfs{I}^{+}$. The conformal factor $\Omega$ was a fixed, time-independent function of the coordinates, while the numerical grid did not include $\mathrsfs{I}^{+}$. No numerical instabilities directly associated with null infinity were detected in his work, however continuum instabilities originating in the bulk of the spacetime obstructed long-time evolutions.

Inspired by Zenginoğlu’s work, we have previously presented stable unconstrained evolutions in spherical symmetry for different types of regular initial data [39]. Instead of the generalized harmonic formulation of the Einstein equations used by Zenginoğlu, we used the BSSN [21, 22] and Z4 [23] systems, which have been used successfully for black hole spacetimes in the context of singularity-avoiding slicing conditions [40]. In this work we did not use the preferred conformal gauge, but a modified harmonic slicing condition with source terms designed to account for the non-trivial background geometry of hyperboloidal slices (with respect to standard Cauchy slices) and to avoid instabilities in the continuum equations; in particular formally singular damping terms were used to drive the lapse to a constant value at $\mathrsfs{I}^{+}$ throughout the evolution. The shift vector was not evolved, but chosen fixed to its value in flat spacetime and compatible with our requirements at $\mathrsfs{I}^{+}$ (see subsection 2.1).

In the present paper we implement and test a wider class of gauge conditions both for the lapse and the shift. Our lapse conditions are generalisations of the Bona-Massó family of slicing conditions [41], while (apart from a fixed shift) we consider modifications of the Gamma-driver shift condition [42] and the generalized harmonic shift conditions [43]. As chosen by Zenginoğlu, we formulate the gauge conditions in terms of the compactified metric instead of the physical one. For work that discusses time-independent hyperboloidal slices using the Bona-Massó family of slicing conditions formulated in terms of the physical (unrescaled) metric, see [44].

We achieve stable evolutions for various combinations of these gauge conditions, as well as for similar gauge conditions that satisfy the preferred conformal gauge (see section 4), and we describe the techniques we have used to suppress instabilities and choose appropriate parameters in our setup to achieve long-term stable evolutions. While our class of gauge conditions also includes the case of Bondi time parameterization (see sections 2.2, 4 and 6.2.1), in this particular case we have not yet been able to suppress code instabilities. We will describe how the solution can in general be reparameterized “a posteriori” to obtain undistorted signals at $\mathrsfs{I}^{+}$. Essentially, we find that while work may be required to make a particular class of gauge conditions work, and in some cases, such as the Bondi time parameterization, further work is still required, a robust treatment of $\mathrsfs{I}^{+}$ is indeed possible with respect to different gauge choices. We will present extensions of our work with black hole spacetimes in future work [45]. For a recent study of hyperboloidal black hole initial data see [46]. Another innovative approach to the hyperboloidal problem that uses the dual foliation formalism [47] is [48].

The implementation of stable gauge conditions that satisfy the preferred conformal gauge also simplified the inclusion of $\mathrsfs{I}^{+}$ in the numerical grid. Previously [39], we used a staggered grid that did not include $\mathrsfs{I}^{+}$ (where the evolution equations are formally singular) as a grid point, but we could show that extrapolation to $\mathrsfs{I}^{+}$ converged at the expected order, consistent with the convergence order in the interior of the grid. With the new non-staggered grid we can evaluate all evolved quantities exactly on future null infinity, without any kind of extrapolation. This approach has in particular two types of technical ingredients (see subsection 6.2): the identification and imposition of the regularity conditions at $\mathrsfs{I}^{+}$ and the transformation and rewriting of the divergent terms in the equations so that they attain a finite limit at future null infinity.

The paper is organized as follows: a detailed description of the preferred conformal gauge and its relation with the Bondi time is presented in section 2. How we adapted several standard slicing and shift conditions to the hyperboloidal problem is explained in section 3. The implementation of the preferred conformal gauge as a system of hyperbolic gauge conditions is described in section 4. The performance of some gauge conditions is illustrated in section 5, where the deformation of the scalar field signal at $\mathrsfs{I}^{+}$ is compared for several combinations of slicing and shift conditions. The implementation of the non-staggered grid that includes a gridpoint on $\mathrsfs{I}^{+}$, as well as some results that compare its behaviour to that of the staggered grid, are presented in section 6. Our notation, details on the formulation of the Einstein equations and the evolution variables we use is summarized in A. Explicit expressions of the gauge conditions that satisfy the preferred conformal gauge are presented in B, while previous versions of slicing conditions that were tested in our code are included in C.

2 Preferred conformal gauge

Here we will describe our basic coordinate setup (choice of conformal factor and fixed coordinate location of $\mathrsfs{I}^{+}$) as the framework where we consider the preferred conformal gauge (3), and relate both to the Bondi time parameterization at $\mathrsfs{I}^{+}$. How we implemented the preferred conformal gauge as hyperbolic gauge conditions will be explained in section 4, after introducing several modifications of standard general hyperbolic gauge conditions for the hyperboloidal initial value problem in section 3.

2.1 Our basic setup: time-independent conformal factor and scri-fixing

We start by choosing a conformal factor $\Omega$ as a fixed function of the compactified radial coordinate $r$, so that $\Omega$ is in particular time independent. Furthermore, we assume the following functional form such that $\Omega(r)$ is a regular function which vanishes at $\mathrsfs{I}^{+}$ with non-vanishing derivative (compare e.g. [9, 49]),

[TABLE]

where $r_{\!\!\mathrsfs{I}}$ is the coordinate location of future null infinity (which we set to $r_{\!\!\mathrsfs{I}}=1$ without restricting generality). The previous expression is obtained in spherical symmetry by imposing an explicitly conformally flat spatial metric in our coordinates (setting (20b) in [39], with $a=-3/K_{CMC}$, equal to unity) on a constant-mean-curvature (CMC) slice, where the trace of the physical extrinsic curvature is constant and corresponds to the negative parameter $K_{CMC}$. Note that even if the initial data in a simulation correspond to a CMC slice, this is very likely to change during the evolution, but the value of $K_{CMC}$ in the definition of the conformal factor $\Omega$ (5) is kept constant throughout the evolution. An alternative simple choice for a time-independent conformal factor is $\tilde{\Omega}=r_{\!\!\mathrsfs{I}}-r$, with $r_{\!\!\mathrsfs{I}}=1$, used by Zenginoğlu in [50, 35, 37, 51]. This choice provides a simpler relation between the compactified radial coordinate and the conformal factor ($r=r_{\!\!\mathrsfs{I}}-\tilde{\Omega}$, in order to use the latter as a coordinate at $\mathrsfs{I}^{+}$), but it may pose problems at the origin, where the conformal factor is not differentiable. The leading order with which the conformal factor goes to zero at $\mathrsfs{I}^{+}$ is the same in both cases, as our $\Omega\propto(r_{\!\!\mathrsfs{I}}+r)(r_{\!\!\mathrsfs{I}}-r)$.

As a further simplification we choose to set $\mathrsfs{I}^{+}$ to a fixed coordinate location in our numerical grid, which is known as scri-fixing [52, 35]. In our spherically symmetric hyperboloidal setup with a time-independent conformal factor we implement it by making the time coordinate flow along $\mathrsfs{I}^{+}$. This requires that the following relation is satisfied at all times:

[TABLE]

see also our discussion and figure 1 in [39].

2.2 Preferred conformal gauge and Bondi time at $\mathrsfs{I}^{+}$

For a metric and conformal factor that satisfy the preferred conformal gauge condition (3) we will use the notation $\hat{g}_{ab}=\hat{\Omega}^{2}\tilde{g}_{ab}$,

[TABLE]

A given conformal factor (5) can always be rescaled to satisfy the preferred conformal gauge condition, by choosing an appropriate $\omega>0$ such that $\hat{\Omega}=\omega\,\Omega$ and (7) is satisfied. From now on we will only refer to (3) as the preferred conformal gauge if our setup is such that (5) satisfies it. An in-depth description of the preferred conformal gauge can be found in chapter 11 of [53] and in [36], and in the same notation as here in section 4.2 in [54].

A real null vector $l^{a}$ in the Newman-Penrose tetrad [55] can be constructed tangent to the null geodesic generators of $\mathrsfs{I}^{+}$, see figure 1 in [39]. Written as $\hat{l}^{a}=\hat{g}^{ab}\hat{\nabla}_{b}\hat{\Omega}$, in the preferred conformal gauge it will satisfy the geodesic equation in the affine parameterization $\hat{l}^{a}\hat{\nabla}_{a}\hat{l}^{b}=0$ on $\mathrsfs{I}^{+}$ by virtue of (4). Its affine parameter is called Bondi time coordinate ($t_{B}$) at $\mathrsfs{I}^{+}$, scaled as $(\partial/\partial t_{B})^{a}\hat{\nabla}_{a}t_{B}=1$. World lines of increasingly distant geodesic observers converge to null geodesic generators of future null infinity, and proper time converges to Bondi time [7], thus the Bondi time parameter can be identified with the proper time of inertial observers at large distances from the source. Note that the preferred conformal gauge does not imply that the time coordinate is indeed Bondi time at $\mathrsfs{I}^{+}$; this is a separate condition, see the discussion and (11) below, and subsections 2.4 in [38].

When the conformal factor (5) chosen for a numerical construction of the spacetime does not satisfy the preferred conformal gauge, it is still possible to determine a conformal factor that does satisfy (7) “a posteriori”. To do so we take advantage of the conformal freedom introduced by $\omega$ and rewrite (7) in terms of $\omega$ and our fixed $\Omega$

[TABLE]

and solve this equation for $\omega$ using the existing numerical data. In our spherically symmetric reduction the previous relation takes the form

[TABLE]

where a dot denotes a time derivative and a prime a spatial derivative with respect to the compactified radial coordinate. If we substitute (6) and use it to rearrange the terms in (9), we will obtain (28) in [39]. In the spherically symmetric case we only need to solve a single ordinary differential equation (9), but in the 3-dimensional case we would expect a system of equations along the 2-dimensional spherical shell that corresponds to $\mathrsfs{I}^{+}$, which is a much more complicated procedure, and thus evolving with the preferred conformal gauge constitutes an important simplification.

The general relation between our numerical time coordinate $t$, whose scaling is determined by the behaviour of the lapse and shift during evolution, and the Bondi time at $\mathrsfs{I}^{+}$ is given by [39]

[TABLE]

This is obtained by identifying the Bondi time with the inertial time coordinate of flat Minkowski spacetime at $\mathrsfs{I}^{+}$. If the preferred conformal gauge is satisfied during evolution, the auxiliary conformal factor is trivial, $\omega=1$. The parameterization of our numerical time coordinate $t$ will depend on the chosen gauge conditions, as they determine the behaviour of $\alpha$ and $\beta^{r}$ during the evolution. From expression (10) we obtain the relation between lapse and shift that has to be satisfied at $\mathrsfs{I}^{+}$ for the code time $t$ to coincide with $t_{B}$ when the preferred conformal gauge (3) holds:

[TABLE]

Figure 6 in [39] shows a comparison of a signal on some non-affinely parameterized time (and where (3) does not hold) and with a Bondi rescaled time. A similar, more complete comparison is included in the present paper in figure 2.

As mentioned before, the preferred conformal gauge (3) also introduces several simplifications in the equations in our setup, among them: i) the independent regular limit of the divergent terms at $\mathrsfs{I}^{+}$ in (2); ii) the $tt$ component of the rescaled metric goes to zero as $\left.-\bar{g}_{tt}\right|_{\mathrsfs{I}^{+}}\propto\Omega^{2}$, faster than in the general case where it vanishes as $\left.-\bar{g}_{tt}\right|_{\mathrsfs{I}^{+}}\propto\Omega^{1}$; iii) the compactified radial coordinate $r$ corresponds to the areal radius at $\mathrsfs{I}^{+}$ for all times during the evolution, which in our notation is expressed as $\dot{\bar{g}_{\theta\theta}}=const$.

The derivation of hyperbolic gauge conditions that satisfy the preferred conformal gauge (3), is included in section 4.

3 Adapting standard hyperbolic gauge conditions to our hyperboloidal approach

The adaptation to the hyperboloidal background geometry is especially important for the equation of motion of the lapse, of which we consider several cases of the Bona-Massó family of slicing conditions, while some additions to the shift evolution equation also improve the simulations. We will explain in detail the slicing and shift conditions used in our approach, and also describe the use of mixed conditions (making some choice of lapse and shift conditions in the interior of the domain and another one near $\mathrsfs{I}^{+}$ - basically to ensure no incoming gauge characteristic speeds appear at $\mathrsfs{I}^{+}$) and how the matching between both was performed.

3.1 Outline of our approach to construct appropriate gauge conditions

Here we will follow and extend our approach in [39]. Taking advantage of the broad existing knowledge about hyperbolic gauge conditions commonly used in current numerical relativity codes, we adapt them to our hyperboloidal implementation, by adding extra damping terms and source functions, and tuning the introduced parameters experimentally to avoid instabilities. By “tune” we mean the process of finding a range of values for the free parameters (such as the damping strength denoted as $\xi$ in sections 3 and 4 and in the appendices) that provide a stable numerical evolution of the complete system.

Using a time-independent conformal factor, in particular the simple form we set in (5), simplifies the numerical implementation of the hyperboloidal initial value problem significantly: vanishing time derivatives of $\Omega$ simplify the complex and formally singular evolution equations (2), and the scri-fixing condition (6) prevents the location of $\mathrsfs{I}^{+}$ from moving in the numerical integration domain.

Our approach for adapting existing hyperbolic slicing conditions is the following: we first set the principal part [56, 57, 58] (the terms in the equations that determine hyperbolicity and the characteristic behaviour of the system), then add damping terms with adjustable strength and source functions, and adapt them according to experimental results until a stable numerical evolution with a well-behaved stationary end state is achieved. The lapse equation of motion can be constructed in the physical spacetime (as was done in [44]) and then translated to the conformal domain for its implementation, or be directly developed in the conformal domain. Both options are essentially equivalent thanks to the simple translation of the quantities between the two domains (see section 2.2.4 in [54] and subsection 3.2.4 below), but some of the non-principal part terms may have different coefficients; for instance, compare the 4th term in (20)’s and (22)’s right-hand-side (RHS). As previously Zenginoğlu, we have decided to directly work in the conformal domain, which appeared simpler, but it would be interesting to explore alternatives in future work. The detailed construction process for the slicing and shift conditions we have tested is presented in the following subsections.

3.2 Slicing conditions

The slicing condition that we use has the same basic form as the one in [39], a generalized Bona-Massó equation

[TABLE]

In the same way as in [39] we add two analytic functions of the radial coordinate to $\dot{\alpha}$’s RHS: $K_{0}(r)$ and $L_{0}(r)$. For hyperboloidal slices the extrinsic curvature $K$ approaches a finite negative value on future null infinity in our convention ($K$ approaches zero at spacelike infinity for spatially asymptotically flat slices). Setting the usual choice for $f(\alpha)=n\alpha^{(m-2)}$ with $n>0$ in (12) without $K_{0}$ yields $\dot{\alpha}=\beta^{r}\alpha^{\prime}-n\alpha^{m}K+L_{0}$, and it is clear that for $m=1,2$ the coefficient in front of $\alpha^{m}$ will be positive, causing an exponential growth that will render the numerical evolution unstable. The quantity $K_{0}$ is thus included to counteract the effect of the hyperboloidal background geometry. The source function $L_{0}$ is calculated from flat data on the hyperboloidal slice and its main purpose is to ensure that the conformal lapse $\alpha$ remains finite at $\mathrsfs{I}^{+}$. In the following subsections we will describe the three different choices for the function $f(\alpha)$ we have tested: $f(\alpha)=1$, $n_{1+log}/\alpha$ and $n_{cK}/\alpha^{2}$. The slicing condition (12) can be defined in both the physical or conformal domains, that is, taking $\alpha$ and $K$ to be physical or conformal variables. In all six cases (three choices of $f(\alpha)$ with the evolution equation set in the physical or conformal domains), the effect of the hyperboloidal geometry is encoded in the source terms (and an extra rescaling by a power of $\Omega$ in (50) and (51)), so that the hyperbolicity properties of the lapse’s evolution equation (assuming reasonable data) are qualitatively the same as those for the Cauchy case - with the exception of $\mathrsfs{I}^{+}$, where the equations are singular. However, all combinations of the slicing condition considered here form a strongly hyperbolic system of equations in the interior of the domain, and as $\mathrsfs{I}^{+}$ is an ingoing null hypersurface, no physical information can enter the domain. In subsection 3.2.4 we will describe a possible way of calculating the source functions for a derivation started in the physical domain, as well as giving a similar expression defined in the conformal one. The precise slicing conditions used in the code, expressed in terms of the actual evolution variables, are also included. Previous versions of the lapse equation of motion are summarized in C.

3.2.1 Harmonic condition:

The harmonic slicing condition ($f(\alpha)=1$) has been widely studied in the literature and it has already been tested in the hyperboloidal initial value approach in [38], together with the generalized harmonic formulation of the Einstein equations. The question of stationary hyperboloidal slices was already raised in [44], where the harmonic condition was considered (without source function $L_{0}$) and stationary hyperboloidal slices were found. For these reasons we expected that the harmonic lapse condition would provide successful results in our implementation, as was the case dealing with regular spacetimes described in [39]. In our formulation, the generalized harmonic slicing condition takes the form

[TABLE]

and the characteristic gauge speeds provided by this equation of motion are in perfect agreement with the requirements at $\mathrsfs{I}^{+}$, as they coincide with the speed of light. The effect that the tilting of the causal cone along the hyperboloidal slice has on the speed of light is shown in figure 1, where the characteristic speeds for flat spacetime data are displayed.

3.2.2 1+log condition:

The 1+log slicing condition is obtained by substituting $f(\alpha)={\textstyle\frac{n_{1+log}}{\alpha}}$ (with $n_{1+log}$ a real positive number) in (12),

[TABLE]

The common choice of $n_{1+log}=2$ has proven to be well-behaved in numerical Cauchy evolutions of spacetimes with strong gravitational fields [59, 60, 42]. It has not been directly set in the previous expression because it does not provide physical characteristic speeds at future null infinity for any value of $K_{CMC}$ - a property of the 1+log slicing condition is that gauge speeds can easily become superluminal. If a characteristic speed at $\mathrsfs{I}^{+}$ becomes larger than the lightspeed there, a negative characteristic speed (incoming mode) will also appear at $\mathrsfs{I}^{+}$. The behaviour of the lightspeeds $c_{\pm}$ in the flat hyperboloidal CMC geometry is illustrated by figure 1. As $\mathrsfs{I}^{+}$ is an ingoing null hypersurface, no physical information from the outside can cross it. Gauge information is however allowed to enter the domain, but then boundary conditions have to be implemented to treat the incoming gauge modes. Given that boundary conditions are a difficult problem and we do not know what a good choice for them at $\mathrsfs{I}^{+}$ in the hyperboloidal approach would be, we will for the moment use physical gauge propagation speeds. Thus, instead of $n_{1+log}=2$ we make the following default choice for $n_{1+log}$, which makes the 1+log expression compatible with the physical characteristic speeds at $\mathrsfs{I}^{+}$ in the stationary regime (but not in the dynamical one - see subsection 3.4):

[TABLE]

If (14) is defined in the physical domain and translated to conformal evolution quantities, an extra dividing factor $\Omega$ has to added to $n_{1+log}$ to prevent the slicing outgoing characteristic speed from vanishing at $\mathrsfs{I}^{+}$. This was the case for (53) in C.2.

The 1+log condition studied in [44] did not include a source function $L_{0}$ as we do here and the calculations were performed with the physical quantities on a non-compactified hyperboloidal slice. Therefore, [44]’s conclusions, which state that no stationary hyperboloidal slices could be found for the 1+log case with nonzero offset $K_{0}$, are not applicable to (14) in the way we have constructed it, as the 1+log conditions considered in [44] and the current paper are different. We have indeed performed successful hyperboloidal simulations with our adapted version of 1+log.

3.2.3 cK condition:

Motivated by our work to evolve black hole spacetimes, which we will report elsewhere [45] (in particular to increase characteristic speeds inside of the horizon to propagate the perturbations faster, which we found to aid the stability of simulations, see [61, 54]), we also implemented the Bona-Massó condition with $f(\alpha)={\textstyle\frac{n_{cK}}{\alpha^{2}}}$, where again $n_{cK}$ is a real positive number,

[TABLE]

We call this slicing condition the “cK” condition, for the form of the $-n_{cK}K$ term in its RHS: constant times K. Unlike in the harmonic and 1+log slicing conditions, here the proportionality factor of $K$ in the RHS will not vanish if $\alpha=0$. In a similar way as with the 1+log condition, the value of $n_{cK}$ is chosen such that the gauge propagation speeds at $\mathrsfs{I}^{+}$ correspond to the physical ones once the stationary state has been achieved, namely

[TABLE]

This choice of $f(\alpha)$ does not seem to have been studied in the literature. A related expression arose in the study of gauge shocks by Alcubierre in [62], where $f(\alpha)=1+k/\alpha^{2}$ with $k$ a vanishing or positive constant is the solution of the differential equation that guarantees that the characteristic fields will not generate shocks (as is the case for the harmonic slicing, with $k=0$). It has been further studied in [63, 64, 65], but we are not aware of a studies of the slicing condition (16) in the literature. For a vanishing shift, slicing conditions of the form (12) lead to a generalized wave equation of the lapse along the direction normal to the slices with wave speed $v_{g}=\alpha\sqrt{f(\alpha)\bar{\gamma}^{ii}}$ along the direction $x^{i}$ [63], where $\bar{\gamma}_{ij}$ is the spatial metric. Assuming spherical symmetry (and a flat spatial metric $\bar{\gamma}^{rr}=\chi/\gamma_{rr}=1$) [66], for this choice of $f(\alpha)$ the characteristic speed along the radial direction associated to the slicing condition is $v_{g}=\sqrt{n_{cK}}$, a constant. In the past, algebraic lapse conditions [67, 68] were considered for black hole spacetimes. Some of those algebraic slicings relate the lapse to the determinant of the spatial metric ($\bar{\gamma}$ in our notation) and are calculated from the Bona-Massó slicing conditions (12) assuming a zero shift and $K_{0}=L_{0}=0$. Using the relation $\partial_{t}\log\bar{\gamma}=-2\alpha K$, for the harmonic case we obtain $\alpha=f(x^{i})\sqrt{\bar{\gamma}}$, where $f(x^{i})$ denotes some chosen function of the spatial coordinates, and for the 1+log one, with coefficient $n_{1+log}$ as in (14), $\alpha=f(x^{i})+{\textstyle\frac{n_{1+log}}{2}}\log{\bar{\gamma}}$. For our cK condition the result is $\alpha=\sqrt{f(x^{i})+n_{cK}\log\bar{\gamma}}$.

3.2.4 Calculation of the source functions:

In the physical spacetime we define the following Bona-Massó-like slicing condition, where tildes denote physical quantities:

[TABLE]

The coefficient $f$ is a function of the lapse, the background lapse and the radial coordinate. The quantities $\hat{\tilde{\alpha}}$ and $\hat{\beta^{\tilde{r}}}$ are part of the background physical (uncompactified) metric components (corresponding to the flat metric on a hyperboloidal slice), which together with $\hat{\tilde{\chi}}$, $\hat{{\gamma_{\tilde{r}\tilde{r}}}}$ and $\hat{{\gamma_{\theta\theta}}}$ (the last one does not change under a conformal compactification) are given by

[TABLE]

The background trace of the extrinsic curvature $\hat{\tilde{K}}=K_{CMC}$ has already been included in (18) in the place of $K_{0}$ in (12). Note that all the previous expressions use the uncompactified radial coordinate $\tilde{r}$. The RHS in (18) has been designed to vanish when the variables $\tilde{\alpha}$, $\beta^{\tilde{r}}$ and $\tilde{K}$ are equal to their background values. The parameters $\xi_{1}$ and $\xi_{2}$ have been introduced to control the damping on the lapse. In most cases, enforcing the linear term controlled by $\xi_{1}$ is enough for stability. However, including the quadratic damping term with $\xi_{2}$ is useful when experimenting, as due to the non-linearities the behaviour of the complete system can be different.

The equation of motion for the conformal lapse $\alpha$ is obtained by substituting $\tilde{\alpha}=\alpha/\Omega$ and $\hat{\tilde{\alpha}}=\hat{\alpha}/\Omega$. The transformation of the advection terms is such that ${\beta^{\tilde{r}}}\partial_{\tilde{r}}\alpha(\tilde{r})={\beta^{r}}\partial_{r}\alpha(r)\equiv{\beta^{r}}\alpha^{\prime}$. For the relation between physical and conformally compactified quantities see subsection 2.2.4 in [54]. We maintain the physical $\tilde{K}$ (no mixing with the Z4 quantity $\Theta$ is included), but express it as the variation with respect to its background value $\Delta\tilde{K}=\tilde{K}-K_{CMC}$. The transformed slicing condition yields

[TABLE]

with $f=1$ for harmonic, $f=m_{1+log}\,\hat{\alpha}/\alpha$ for 1+log and $f=m_{cK}\,\hat{\alpha}^{2}/\alpha^{2}$ for the cK condition - the $\hat{\alpha}$ factors are included to make $f$ conformally invariant. Comparing with (14) and (16), the coefficients $n_{1+log}\equiv m_{1+log}\hat{\alpha}$ and $n_{cK}\equiv m_{cK}\hat{\alpha}^{2}$ are no longer constant: they are time-independent, but functions of the compactified radius $r$. The components of the background conformally compactified metric are

[TABLE]

corresponding to initial and stationary flat data and calculated from (20) in [39] with $a=-3/K_{CMC}$ and setting $\Omega^{\prime}$ from (5). The damping terms controlled by $\xi_{1}$ and $\xi_{2}$ appear now divided by $\Omega$, in the way that our numerical experiments so far seem to require.

A slicing condition similar to (20) can be derived from (12) in the conformal domain. For this, the conformal extrinsic curvature $K$ in (12) is transformed to its physical equivalent $\tilde{K}$ using $K={\textstyle\frac{\tilde{K}}{\Omega}}+{\textstyle\frac{3\beta^{r}\Omega^{\prime}}{\alpha\Omega}}$, and then the result is expressed in terms of $\Delta\tilde{K}$. After adding source functions and damping terms as in (20), the final expression yields

[TABLE]

The only difference with respect to (20) is the $3f$ coefficient in front of the 4th term in the RHS above. If $f=m_{cK}\,\hat{\alpha}^{2}/\alpha^{2}$ is set and at some point $\alpha=0$, which can happen for strong field data, the slicing condition will have a diverging term and (22) may render the simulation unstable. This can be solved by introducing the $f$ factor in front of $\alpha^{2}\tilde{K}$ after the transformation to the physical trace of the extrinsic curvature has been performed: then the coefficient in front of the 4th term in (22)’s RHS will be simply 3.

3.3 Shift conditions

3.3.1 Fixed shift:

For our first simulations on regular spacetimes, presented in [39], we used a time-independent shift equal to its initial (and also stationary) value, expressed in terms of our compactified radial coordinate. This is the simplest shift choice we could find that is compatible with the scri-fixing condition (6). Nevertheless, for simulations with richer physical content, such as those including a black hole or large perturbations of a scalar field, an evolved time-dependent shift is a more appropriate choice. Among the possible hyperbolic shift conditions, we tested the following ones.

3.3.2 Gamma-driver:

Our first test with a hyperbolic shift condition was based on the Gamma-driver shift [42], due to its common use in current numerical relativity simulations and its simple form,

[TABLE]

where the auxiliary variable $B^{a}$ is coupled to the contracted connection $\Gamma^{a}$. Given that the flat spacetime value of the three-dimensional contracted connection $\Gamma^{r}=\gamma^{ij}\Gamma^{r}_{ij}=-2/r$ is badly behaved at the origin in spherical coordinates, for numerical purposes it is much more convenient to use instead the related quantity $\Lambda^{a}$ (see [69] or A for more details),

[TABLE]

where $\hat{\Gamma}^{a}_{bc}$ is a time-independent background connection. Thus, we couple the evolution equation for $B^{a}$ (23) to $\Lambda^{a}$ instead of $\Gamma^{a}$. In a similar way as done with the slicing conditions in subsection 3.2.4, we add a source term to account for the stationary value of the advection term, and a damping term. The resulting generalized Gamma-driver shift condition yields

[TABLE]

The Gamma-driver shift can also be implemented in its integrated version [70], which in spherical symmetry takes the form

[TABLE]

The properties of the principal part of these equations are formally the same as for their non-hyperboloidal equivalent, because the added source terms do not change the principal part: versions (25) or (26) of the Gamma-driver together with the Einstein equations and any hyperbolic slicing condition from subsection 3.2 form a strongly hyperbolic system for regular data. The same analysis as in [71] applies here for the quantities evaluated on a hyperboloidal slice (except at $\mathrsfs{I}$, where the equations are singular). The gauge propagation speeds associated with the shift condition depend on the values of the parameters $\lambda$ and $\mu$. There will be no negative characteristic speeds (no incoming modes) at $\mathrsfs{I}^{+}$ if the product $\lambda\,\mu\leq\left({\textstyle\frac{r_{\!\!\mathrsfs{I}}}{3}}K_{CMC}\right)^{2}$ for the Gamma-driver condition with auxiliary variable $B^{r}$ (25) and if $\lambda\leq{\textstyle\frac{1}{12}}\left(r_{\!\!\mathrsfs{I}}K_{CMC}\right)^{2}$ for the integrated version (26). The reason for keeping the parameter $\mu$ separated from the parameter $\lambda$ in the Gamma-driver shift condition (25) is that even if the characteristic speeds depend only on $\lambda\,\mu$, different choices of $\lambda$ and $\mu$ may still lead to a different numerical behaviour. For simulations on spatially asymptotically flat slices, the standard choice for the parameter $\lambda$ in the Gamma-driver shift condition is $\lambda=3/4$ (with $\mu=1$), which implies physical characteristic speeds (i.e. coinciding with the speed of light). A detailed study on parameter choices that provide strong hyperbolicity in non-hyperboloidal black hole simulations can be found in [71]. In our hyperboloidal approach, setting $\lambda=3/4\,\alpha^{2}\chi$ in (26) gives exactly the physical characteristic speeds everywhere. Future research, specially when considering a full three-dimensional system, may show that choosing values for the gauge characteristic speeds different from the physical ones is a more convenient option. There is much flexibility in the three-dimensional case, so that better choices may be found, following e.g. [72]. However, at this stage of development of the spherically symmetric case, the choice of physical characteristic speeds for the gauge modes (at least asymptotically) seems like a reasonable option.

3.3.3 Doubly generalized harmonic shift:

The standard generalized harmonic shift condition [43, 73] is defined as the spatial part of the four-dimensional equation $\bar{g}^{ab}\bar{\Gamma}^{c}_{ab}=F^{c}$, where $F^{c}$ is usually chosen to be a function of the coordinates and evolution variables, and such that it does not affect the principal part of the Einstein equations. The time component gives the usual harmonic slicing condition (with dynamical part like (13)), while in spherical symmetry the radial component gives an evolution equation for the radial shift. The spatial three-dimensional $\Lambda^{a}$ (24) (as $\Lambda^{r}$) is introduced in the shift condition to substitute the first spatial derivatives of the spatial metric components, common practice [73] required in order to make the system strongly hyperbolic. We set the source term to $F^{r}={\textstyle\frac{\xi_{\beta^{r}}}{\alpha^{2}\Omega}}\beta^{r}-{\textstyle\frac{L_{0}}{\alpha^{2}}}-\beta^{r}F^{t}$, so that it provides a source function $L_{0}$ (in a similar way as done with the slicing condition) and a damping term (proportional to $-\beta^{r}/\Omega$) to ensure that the value of $\beta^{r}$ stays fixed at $\mathrsfs{I}^{+}$, a behaviour similar to the time-independent shift we used in [39]. Our doubly generalized harmonic shift condition then takes the form:

[TABLE]

This equation includes a term with $r$ in its denominator that formally diverges at the origin. In numerical tests we find that indeed this term results in an instability at the center $r=0$. Note also that this term would not have appeared if (27) had been derived from $\bar{g}^{ab}\bar{\Lambda}^{c}_{ab}=F^{c}$ instead of $\bar{g}^{ab}\bar{\Gamma}^{c}_{ab}=F^{c}$. As this diverging term does not belong to the principal part, we simply absorb it into $L_{0}$, and thus drop it in equation (27). The final form of our adapted harmonic shift condition is obtained after substituting $L_{0}=\hat{\alpha}\,\hat{\alpha}^{\prime}-\hat{\beta^{r}}\hat{\beta^{r}}^{\prime}+\xi_{\beta^{r}}{\textstyle\frac{\hat{\beta^{r}}}{\Omega}}$, which ensures that flat hyperboloidal data are a stationary solution of the equation:

[TABLE]

Note that this evolution equation provides physical characteristic speeds. Optionally, a damping term of the form $+\eta(\hat{\beta^{r}}-{\beta^{r}})$, like the one in the Gamma-driver condition, can be added. This damping term has been useful preventing instabilities at the origin for several configurations described in section 5.

3.4 Matching

The numerical results we obtained using the integrated Gamma-driver with $\lambda\neq 3/4\,\alpha^{2}\chi$, the 1+log (14) or the cK (16) slicing conditions, with choices (15) and (17) respectively, showed, at some point during the evolution, some fluctuations arising at $\mathrsfs{I}^{+}$ that propagated inwards. An appropriate choice of parameters and enough dissipation can keep these incoming fluctuations under control, but the loss in parameter freedom restricts the numerical experiments considerably.

The most plausible explanation for these fluctuations, as after treating this issue they did not appear anymore, is that at some point at least one of the characteristic speeds at $\mathrsfs{I}^{+}$ becomes negative and some information enters from outside of the domain. This cannot happen to the physical speeds, but it may happen for the gauge speeds: the speeds related to the shift and slicing conditions. The incoming physical characteristic speed is proportional to the expression in the scri-fixing condition (6) (see figure 1) and this ensures that no information will travel inwards across $\mathrsfs{I}^{+}$. Nevertheless, the fixed values that are set as coefficients in some principal part terms of the gauge conditions, like (15) and (17), and the mentioned ranges for $\lambda$ in (26), may produce an ingoing speed that is no longer proportional to (6), so that non-vanishing incoming gauge speeds can arise when the value of the metric variables at $\mathrsfs{I}^{+}$ differs from the stationary one. For instance, if $-\left.\sqrt{\frac{\gamma_{rr}}{\chi}}\beta^{r}\right|_{\mathrsfs{I}^{+}}$ is different from $\sqrt{n_{cK}}\equiv\alpha|_{\mathrsfs{I}^{+}}$ (given by (17)) at some point during the evolution, the speed of the slicing incoming mode at $\mathrsfs{I}^{+}$ does not vanish and uncontrolled information enters the domain.

Considering only physical gauge speeds at $\mathrsfs{I}^{+}$, this problem is solved by making the principal part of the slicing and shift conditions coincide exactly with those of the harmonic slicing and shift conditions respectively at $\mathrsfs{I}^{+}$ at all times during the evolution. However, we may want to have different characteristic speeds in the interior of the domain, especially when dealing with strong fields near the origin. To achieve this we can use a smooth transition between the required expressions at $\mathrsfs{I}^{+}$ and our choices in the interior.

A possibility is to match using the following function that vanishes at $\mathrsfs{I}^{+}$

[TABLE]

and construct the coefficient of the $K$ or the $\Lambda^{r}$ term in the corresponding slicing or shift condition respectively as:

[TABLE]

where $n_{cK}$, $\lambda$, $n_{\alpha}$ and $n_{\beta^{r}}$ are free parameters. In this way, the gauge characteristic speeds at future null infinity will always coincide with the physical speeds and no incoming modes will appear.

The region of matching by superposition in (30) extends from $r=0$ to $r=r_{\!\!\mathrsfs{I}}$. In order to control the extension and location of the matching region, we have experimented with a new matching that takes place between two given values of the compactified radial coordinate $r_{-}$ and $r_{+}$, with $r_{-},r_{+}\in[0,r_{\!\!\mathrsfs{I}}]$ and $r_{-}<r_{+}$. We use a different function,

[TABLE]

designed to have a smooth profile everywhere in the domain; assuming the natural extensions from (32) ($f_{match}=0$ for $r<r_{-}$ and $f_{match}=1$ for $r>r_{+}$), this matching function is $C^{3}$ at the joining points $r_{-}$ and $r_{+}$. The matching joins $G_{int}$ in the interior part and $G_{ext}$ in the exterior including $\mathrsfs{I}^{+}$,

[TABLE]

where $G$ can be identified with the coefficients $n_{coeff}$ or $\lambda_{coeff}$ in (30) or the complete RHS of $\dot{\alpha}$ or $\dot{\beta}^{r}$. In this way, for $r<r_{-}$ we can have larger characteristic speeds, and between $r_{+}$ and $r_{\!\!\mathrsfs{I}}$ the physical characteristic propagation speeds will ensure the correct treatment of $\mathrsfs{I}^{+}$.

4 Implementation of the preferred conformal gauge as hyperbolic gauge conditions

We will now construct source terms for our generalized harmonic gauge discussed in section 3, which will imply the preferred conformal gauge (3). The complete derivation and preliminary testing of these gauge conditions is included in subsections 4.5.2 and 7.3.10 in [54] respectively.

The main motivation behind the design of these gauge conditions is to try and find a more systematic way of expressing the source functions in the gauge conditions (see section 3 for how the source functions were set in the adapted gauge conditions) and also to enforce the preferred conformal gauge (3) at all times by means of the gauge conditions, as suggested by [36]. We start with the following expression

[TABLE]

where on the physical 4-dimensional spacetime we set the contraction of the difference between the connection $\tilde{\Gamma}^{c}_{ab}$ of the physical metric $\tilde{g}_{ab}$ and that $\hat{\tilde{\Gamma}}{}^{c}_{ab}$ of a time-independent physical background metric $\hat{\tilde{g}}_{ab}$ equal to some source function $\tilde{F}^{a}$ – note that these are all 4-dimensional quantities. This construction basically corresponds to the generalized harmonic gauge condition (see subsections 3.2.1 and 3.3.3) on the physical spacetime, but with the extra subtraction of background terms (especially convenient in spherical symmetry and on a hyperboloidal background), which was inspired by the definition of the 3-dimensional conformal spatial quantity $\Lambda^{a}$ (see A and [69]). The background terms do not change the principal part of the system (the same as for generalized harmonic with respect to the physical metric), because only non-principal part terms including inverse metric components are introduced. An equivalent construction can be performed in the conformal domain, that is, using the conformal metric $\bar{g}_{ab}$ and the conformal connections $\bar{\Gamma}^{c}_{ab}$ and $\hat{\bar{\Gamma}}{}^{c}_{ab}$ in (33). However, the obtained numerical results were not satisfactory: stable evolutions were achieved, but at the price of not satisfying the preferred conformal gauge. Thus, instead of including the details here, we point the interested reader to subsections 4.5.1 and 7.3.9 in [54], where the derivation and numerical experiments of the conformal version of (33) are described in detail.

Using (1) and the corresponding transformation from physical to conformal Christoffel symbols (see (2.4) in [54], for instance), the above condition is transformed to conformal quantities:

[TABLE]

where $\hat{\bar{g}}_{ab}$ is the time-independent conformal background metric and $\hat{\bar{\Gamma}}{}^{c}_{ab}$ its associated connection. The preferred conformal gauge (3) can be expressed as

[TABLE]

and substituting $\bar{\Gamma}^{c}_{ab}$ above from the gauge conditions (34) yields

[TABLE]

The quantity $\left.\Omega^{\prime}\right|_{\mathrsfs{I}^{+}}\neq 0$ by definition and in principle also $\left.\Omega^{\prime\prime}\right|_{\mathrsfs{I}^{+}}\neq 0$. The metric component $\bar{g}^{rr}$ in our variables is $\bar{g}^{rr}=\frac{\chi}{\gamma_{rr}}-\frac{{\beta^{r}}^{2}}{\alpha^{2}}$, which vanishes at $\mathrsfs{I}^{+}$ by virtue of the scri-fixing condition, see (6). The background Christoffel symbol $\hat{\bar{\Gamma}}{}^{r}_{ab}$ is calculated from flat spacetime data and it is such that at $\mathrsfs{I}^{+}$ the only nonzero component is $\left.\hat{\bar{\Gamma}}{}^{r}_{rr}\right|_{\mathrsfs{I}^{+}}=\left.{\textstyle\frac{1}{r}}\right|_{\mathrsfs{I}^{+}}={\textstyle\frac{1}{r_{\!\!\mathrsfs{I}}}}$. Thus, $\bar{g}^{ab}\hat{\bar{\Gamma}}{}^{r}_{ab}\equiv\bar{g}^{rr}\hat{\bar{\Gamma}}{}^{r}_{rr}$ will always vanish at $\mathrsfs{I}^{+}$. The third term also becomes zero at null infinity, because $\left.\bar{g}^{rr}\right|_{\mathrsfs{I}^{+}}\equiv\left.\hat{\bar{g}}{}^{rr}\right|_{\mathrsfs{I}^{+}}=\mathcal{O}(\Omega^{2})$ and $\bar{g}^{ab}\hat{\bar{g}}_{ab}$, a combination of conformally rescaled metrics, is expected to attain a finite value. The only remaining term will also vanish at future null infinity if $\left.\tilde{F}^{r}\right|_{\mathrsfs{I}^{+}}\propto\Omega^{q}$ with $q>2$.

The actual evolved gauge conditions for the lapse and shift are obtained by isolating $\dot{\alpha}$ and $\dot{\beta}^{r}$ from the time and radial components of (34) respectively and are presented in B. The choice of source function $\tilde{F}^{a}$ compatible with the preferred conformal gauge that provided a stable numerical behaviour in our experiments was $\tilde{F}^{r}=0$ and $\tilde{F}^{t}=\xi_{\alpha}\Omega(\alpha-\hat{\alpha})$, where $\xi_{\alpha}$ is a parameter to be tuned experimentally.

The relation between lapse and shift at $\mathrsfs{I}^{+}$ for the previously derived gauge conditions is given by (38). A vanishing $\xi_{\alpha}$ would make the code time $t$ and the Bondi time $t_{B}$ in (10) exactly equal, as (38) with $\xi_{\alpha}=0$ reduces to (11). This means that we would be able to have the Bondi time parameterization at $\mathrsfs{I}^{+}$ during the numerical evolutions. Unfortunately, the choice $\xi_{\alpha}=0$ is not yet possible in our implementation, because it does not provide stable numerical simulations. Only a range of non-negative values is allowed, as $\xi_{\alpha}$ is introduced as the strength of a damping term that would probably create an exponential growth if the wrong sign was used. The range $\xi_{\alpha}\in[2,4]$ works fine for most configurations. The conclusion is that despite satisfying the preferred conformal gauge (3), the time coordinate at $\mathrsfs{I}^{+}$ in our simulations will not be affinely parameterized at all times (but it can be easily rescaled using (10)).

If a match of type (30) is attempted with the gauge conditions derived here and included in B, the exponents $n_{\alpha}$ and $n_{\beta^{r}}$ have to be equal to or larger than 3, or else the preferred conformal gauge will not hold at $\mathrsfs{I}^{+}$.

When used together with our spherically symmetric reduction of the Z4c formulation of the Einstein equations, the previously derived gauge conditions provide an evolution that satisfies the preferred conformal gauge at $\mathrsfs{I}^{+}$, in the sense that the scheme converges when increasing the numerical resolution. However, evolving the same gauge conditions with our implementation of the generalised BSSN formulation does not satisfy (3) yet. We are currently working towards obtaining a stable evolution of BSSN together with the preferred conformal gauge.

5 Gauge condition results

The gauge conditions described in the sections 3 and 4 provide many possible combinations: slicing and shift conditions that can optionally be matched to others in the region close to $\mathrsfs{I}^{+}$ with a variety of matching functions. Among the parameters to set we have: the strength of the damping in the slicing and the shift conditions; the coupling in the gauge characteristic speeds for 1+log, cK and Gamma-driver conditions; the location of the possible matchings. It is beyond the scope of this work to describe all stable configurations with our current spherically symmetric implementation of the Einstein equations (see A for details), for the whole system is interrelated and changing one parameter can affect the values of many others. We will provide a few examples of gauge choices that appear convenient for simulations of regular data and briefly comment on their observed behaviour.

Our implementation uses the Method of Lines, with 4th order finite differences for the spatial discretization and a 4th order Runge-Kutta for the time integration. The simulations described here were performed with a staggered grid (right diagram in figure 4), 200 spatial gridpoints and a time step of $dt=0.0005$, which provides a Courant number of 0.1 (we can evolve the staggered grid with a Courant number of 0.2, but the simulations with a point on $\mathrsfs{I}^{+}$ presented in section 6 require 0.1 for stability). One point off-centered stencils were implemented in the derivatives in the advection terms, as is common practice in numerical relativity simulations [74]. Kreiss-Oliger dissipation terms [75] as used in [76] are added to the RHSs, with a strength parameter of $0.25$ for these simulations. We evolved our Z4c system (dropping the non-principal part terms labelled with $C_{Z4c}$), with parameters $\kappa_{1}=1.5$ and $\kappa_{2}=0$. The dynamics in the simulations was introduced by a Gaussian-like scalar field of the form $\bar{\Phi}_{0}=\tilde{\Phi}_{0}/\Omega=A/\Omega\,\rme^{-\frac{(r^{2}-c^{2})^{2}}{4\sigma^{4}}}$ (like the one presented in figures 5 and 6 in [39]) that perturbed flat initial data, with parameter choices $A=0.058$, $\sigma=0.1$ and $c=0.25$, and with the constant trace of the physical extrinsic curvature in the initial data set to $K_{CMC}=-3$. The coordinate location of future null infinity is set to $r_{\!\!\mathrsfs{I}}=1$.

The gauge conditions that satisfy the preferred conformal gauge (48) provide stable results with the parameter choice $\xi_{\alpha}=2$. Both the physical (20) and conformal (22) slicing conditions give good results with $\xi_{1}=2$ and $\xi_{2}=0$ or with $\xi_{1}=0$ and $\xi_{2}=2$. Choosing $\xi_{1}=1$ and $\xi_{2}=0$ together with a fixed shift does not provide enough damping to suppress the small fluctuations that appear in the variable $\Lambda^{r}$ right next to $\mathrsfs{I}^{+}$ (in a staggered grid the outermost gridpoint is half a spatial step away from $\mathrsfs{I}^{+}$, see figure 4) and propagate inwards while growing in amplitude. Nevertheless, $\xi_{1}=1$ is enough for long-term stability when both the harmonic lapse and shift are used near $\mathrsfs{I}^{+}$.

Even in regular spacetimes, instabilities at the origin (in the form of high frequency fluctuations that grow very fast) have appeared for the hyperbolic shift conditions where the coefficient in front of $\Lambda^{r}$ is proportional to $\alpha^{2}\chi$, that is, for the harmonic shift and the integrated Gamma-driver with $\lambda=3/4\,\alpha^{2}\chi$. These instabilities can be effectively suppressed by using a superposition of the form (30b) with $\lambda=0.1$ and $n_{\beta^{r}}=2$. As a Gamma-driver-type shift condition with constant coefficient $\lambda$ is a better choice for the interior part where $r=0$ or where a black hole would be located, we have not performed an in-depth study of this instability. More interesting is the performance of the same shift conditions (harmonic shift and integrated Gamma-driver with $\lambda=3/4\,\alpha^{2}\chi$) at $\mathrsfs{I}^{+}$: the first one is stable at long term and presents smooth features at future null infinity, while the second one requires a damping of at least $\xi_{\beta^{r}}=2$ (a non-zero value of $\eta$ is not enough) to suppress ingoing fluctuations in $\Lambda^{r}$ like those described in the previous paragraph. Even if the shift characteristic speeds are the same for the two shift conditions considered, the difference in behaviour can be due to their different characteristic fields or their different coupling to $\Lambda^{r}$.

The matching of the form (30) can be applied at the level of the coefficient in front of $K$ or $\Lambda^{r}$, but also to the whole gauge condition, as $G=G_{int}f_{match}^{n}+G_{ext}$ (using the same notation as in (32)). Implementing the latter for the cK condition and the integrated Gamma-driver as $G_{int}$ with $n=3$ and the preferred conformal gauge conditions (48) as $G_{ext}$ provided a stable evolution for the regular spacetime case considered. However, the relaxation to the final stationary state corresponding to flat spacetime was slower ($t_{relax}\sim 30$) than when using only the unmatched $G_{int}$ ($t_{relax}\sim 6$), for instance. The other matching considered (32) works quite nicely if $r_{-}\leq 0.4$, although the relaxation to the final state is also quite slow ($t_{relax}\sim 30$). If $r_{-}>0.4$, the evolution variables do not move towards their final state monotonically, but oscillate around it (this is especially visible in $\Lambda^{r}$) either with decreasing amplitude, so that the final state is reached in the end, or with an increasing one, which renders the simulation unstable.

5.1 Deformation of the scalar field signals at $\mathrsfs{I}^{+}$

The simulations whose data are presented in figures 2 and 3 used a staggered grid, which means that the curves shown correspond to data extrapolated to $\mathrsfs{I}^{+}$ for half a spatial step. The basic parameters were given at the beginning of this section, while the specific gauge parameters for each simulation will be indicated in the following paragraph. As mentioned above, the simulations included here are only examples of the vast range of possible configurations.

Figure 2 shows the signal of the rescaled scalar field $\bar{\Phi}=\tilde{\Phi}/\Omega$ at future null infinity as raw output from the code. The exception is the black solid line: it corresponds to the signal rescaled by $\hat{\Omega}$, the conformal factor that satisfies the preferred conformal gauge, and is plotted against the Bondi time coordinate. The following 3 signals listed in the legend were produced in simulations with fixed (time-independent) shift and i) harmonic lapse, ii) cK condition with $m_{cK}=1$, and iii) cK condition with $m_{cK}=1$ superposed as in (30a) (with $n_{\alpha}=2$) on the harmonic condition. The 5th listed signal corresponds to a simulation with 1+log ($n_{1+log}=1$) and integrated Gamma-driver ($\lambda=3/4$ and $\xi_{\beta^{r}}=2$, so that the shift is forced to remain fixed at $\mathrsfs{I}^{+}$, see right plot in figure 3), without any matching. The 6th curve was obtained from a simulation with harmonic lapse and shift, with only a small superposition (30b) (with $n_{\beta^{r}}=2$ and $\lambda=0.1$) to keep the shift condition stable at the origin without adding extra damping. The final two signals correspond to simulations where the preferred conformal gauge is satisfied at $\mathrsfs{I}^{+}$. Both include a matching of the form (31) in (32) with the cK condition ($n_{cK}=0.4$) and the integrated Gamma-driver ($\lambda=1$): in the first one the matching is performed practically at the origin ($r_{-}=0$ and $r_{+}=0.05$, this is why it is labelled simply as preferred conformal gauge), while the second one uses $r_{-}=0.4$ and $r_{+}=0.9$. For these two simulations $\left.\Omega\right|_{\mathrsfs{I}^{+}}=\left.\hat{\Omega}\right|_{\mathrsfs{I}^{+}}$ holds (see how $\left.\omega\right|_{\mathrsfs{I}^{+}}=1$ at all times on the left in figure 3) and thus the amplitude of the scalar field signal is the same as for the black solid line.

Other relevant quantities are plotted at $\mathrsfs{I}^{+}$ in figure 3: the auxiliary conformal factor $\omega=\hat{\Omega}/\Omega$ and the variations of the lapse and shift, $\Delta\alpha=\alpha-\hat{\alpha}$ and $\Delta\beta^{r}=\beta^{r}-\hat{\beta^{r}}$. The values of $\left.\omega\right|_{\mathrsfs{I}^{+}}$ differ from unity for the first five simulations, where the preferred conformal gauge does not hold. The lapse and shift stay fixed at $\mathrsfs{I}^{+}$ in the first four simulations (the lapse is “damped” in all of them, while the shift is fixed in the first three ones and “damped” in the fourth one), and so ensure that no incoming modes appear at $\mathrsfs{I}^{+}$. They are however allowed to move in the simulations with harmonic-like behaviour at $\mathrsfs{I}^{+}$: harmonic lapse and shift, and the two conditions (the one matched near the origin and the one matched in the middle of the domain) that satisfy the preferred conformal gauge. Snapshots at constant time of $\alpha(r)$ and $\Delta\beta^{r}(r)$ for similar gauge conditions are presented in figure 8.9 in [54].

It is evident from figure 2 that an evolved shift plays a very important role: its effect in controlling the behaviour of $\Lambda^{r}$ causes less distortion in the spacetime and the evolved quantities reach their final states earlier. The left plot in figure 3 also shows that $\left.\omega\right|_{\mathrsfs{I}^{+}}$’s deviation from unity is much smaller than in the cases with time-independent shift. Of the two cases that satisfy the preferred conformal gauge at $\mathrsfs{I}^{+}$, the one that is matched to a larger interior region with cK slicing condition and Gamma-driver (last listed case in legend) shows smaller deviations in lapse and shift at $\mathrsfs{I}^{+}$ and a scalar field signal at $\mathrsfs{I}^{+}$ closer to the black solid line. It is possible that larger characteristic speeds in the interior or stronger couplings with the $K$ and $\Lambda^{r}$ variables decrease the spacetime distortion caused by the perturbation more effectively. The signal from the simulation with unmatched 1+log and integrated Gamma-driver, not one of the best options due to its potential non-vanishing incoming modes, is surprisingly very close to the undistorted signal measured with the Bondi time coordinate - at least around the minimum in the signal, which corresponds to the part of the scalar field perturbation that was reflected at the origin (see figure 5 in [39]). This is a combination of the “damped” lapse and shift at $\mathrsfs{I}^{+}$ and the fact that the auxiliary conformal factor $\omega$ at $\mathrsfs{I}^{+}$ is practically one again when the reflected signal leaves the domain through $\mathrsfs{I}^{+}$ (the similarity with the undistorted “Bondi” signal is probably not as high in the first peak). This behaviour is not completely by chance, as the same gauge conditions with slightly different parameter choices give a similar result. However, changing the initial perturbation may provide different results.

6 Simulations with a non-staggered grid: point on $\mathrsfs{I}^{+}$

Here we will describe the implementation of the relevant technical improvement that allows us to include $\mathrsfs{I}^{+}$ on a gridpoint in our numerical domain.

6.1 Staggered vs. non-staggered grid

The Einstein equations expressed in terms of the conformally rescaled metric $\bar{g}_{ab}$ and reduced to spherical symmetry formally diverge at two points of the compactified radial coordinate. The first one is the coordinate origin $r=0$, that is just a coordinate singularity for regular data but in presence of a black hole corresponds to a compactified asymptotically flat end for wormhole data (expressed as a puncture [77]) or to the infinitely far end of a trumpet [24] for trumpet data. The other diverging point is the location of future null infinity $r=r_{\!\!\mathrsfs{I}}$, where the conformal factor $\Omega$ vanishes.

The results presented in [39, 61, 78, 54] were obtained using a staggered grid in our numerical implementation. In a staggered grid, the gridpoints where the variables and equations are evaluated during evolution do not coincide with any of the locations where the equations diverge ($r=0$ and $r=r_{\!\!\mathrsfs{I}}$), as illustrated on the right diagram in figure 4. This simplifies the implementation considerably, because no special treatment of the RHSs on the divergent points is required. However, the value of the quantities at $\mathrsfs{I}^{+}$ is only obtained by extrapolation (for less than a spatial step, half a step in our case).

The non-staggered grid in our setup includes the radial locations $r=0$ and $r=r_{\!\!\mathrsfs{I}}$ as actual gridpoints (see diagram on the left of figure 4) and allows us to directly calculate our quantities at $\mathrsfs{I}^{+}$, but its implementation is considerably more difficult: the formally diverging RHSs of the spherically symmetric evolution equations (in Appendix C of [39] and in Chapter 2 of [54]) have to be transformed so that a finite numerical value is obtained at the diverging points. The treatment of $r=0$ is an important issue: a non-staggered grid at the origin can be suitable for the evolution of regular spacetimes (using the parity behaviour of the evolution variables to determine their values at $r=0$111Odd variables will necessarily vanish there, while even ones have a vanishing first derivative at $r=0$, which can be used to calculate their value there. More details are included in subsection 6.3.2 in [54].), but a staggered origin is a better choice for evolving a black hole. For instance, in the moving puncture approach [79, 80] the location of the puncture that represents each of the black holes must never coincide with any of the gridpoints. Therefore, when including a black hole in the evolved spacetime [45], we will use a mixed grid, where the origin $r=0$ is staggered and $\mathrsfs{I}^{+}$ is located on an actual gridpoint.

6.2 Implementation of the equations for the non-staggered grid at $\mathrsfs{I}^{+}$

The modification of our spherically symmetric reduction of the equations of motion for their evaluation at $\mathrsfs{I}^{+}$ is performed in two steps. First, note that the terms divided by $\Omega$ (and $\Omega^{2}$) that appear in the RHSs are expected to cancel and provide finite RHSs that attain a finite value at $r=r_{\!\!\mathrsfs{I}}$ - Friedrich showed that the Einstein equations expressed as the Conformal Field Equations attain a regular limit at $\mathrsfs{I}^{+}$ [6, 8] for the type of regular initial data we are interested in. This means that for these kind of solutions the sum of the numerators of the diverging terms in our equations (in Appendix C of [39] and (2.82) in Chapter 2 of [54]) must necessarily vanish at $\mathrsfs{I}^{+}$ at the appropriate order. These relations (numerators = 0) is what we call the regularity conditions at $\mathrsfs{I}^{+}$, which for our formulation and reduction of the Einstein equations are given by (37), and they will be used to calculate the regular limit of the RHSs at $\mathrsfs{I}^{+}$.

6.2.1 Impose the regularity conditions on the variables at $\mathrsfs{I}^{+}$:

The regularity conditions for the conformally rescaled Einstein equations that we use in our code (included in Chapter 2 of [54] - note that the metric coefficient $\gamma_{\theta\theta}\equiv\gamma_{rr}^{-1/2}$ is eliminated in the equations implemented in our code) assuming that the preferred conformal gauge holds at future null infinity are the following:

[TABLE]

as well as two much lengthier expressions for $\left.A_{rr}\right|_{\mathrsfs{I}^{+}}$ and $\left.\Lambda^{r}\right|_{\mathrsfs{I}^{+}}$, obtained by solving for these quantities the numerators of the $\Omega^{-1}$ terms in $\dot{A}_{rr}$ and $\dot{\Lambda^{r}}$. The equation of motion of the massless scalar field does not introduce any extra regularity conditions.

The regularity condition that arises from the divergent terms in the gauge conditions that satisfy the preferred conformal gauge (described in section 2) is

[TABLE]

Substituting the value of $\left.\beta^{r}\right|_{\mathrsfs{I}^{+}}$ according to (38) into (10), using also $\left.\Omega^{\prime}\right|_{\mathrsfs{I}^{+}}=K_{CMC}{\textstyle\frac{r_{\!\!\mathrsfs{I}}}{3}}$ and $\xi_{\alpha}=0$ yields that the code time $t$ at future infinity equals the Bondi time $t_{B}$. Unfortunately, at the current state our simulations are not stable with a vanishing $\xi_{\alpha}$.

We implemented the regularity conditions in the code in the following way: as there is only one degree of freedom available at $\mathrsfs{I}^{+}$ for our spherically symmetric reduction of the Einstein equations (8 evolution variables - $\chi$, $\gamma_{rr}$, $\alpha$, $\beta^{r}$, $A_{rr}$, $\Delta\tilde{K}$, $\Lambda^{r}$ and $\tilde{\Theta}$ - and only 7 regularity conditions), we choose to identify it with the free value of $\gamma_{rr}$. Using relations (37) and (38), we express all regularity conditions in terms of $\gamma_{rr}$ (this may not be a requirement for stability) and the derivatives of the other variables. Note that the regularity conditions for lapse and shift as expressed above are circular: in order to obtain the expressions to impose in the code we substitute $\chi$ in (37c) in terms of $\gamma_{rr}$ as in (37a), and then solve (37c) and (38) for $\left.\alpha\right|_{\mathrsfs{I}^{+}}$ and $\left.\beta^{r}\right|_{\mathrsfs{I}^{+}}$. The correct result for $\xi_{\alpha}\neq 0$ is given by

[TABLE]

while for $\xi_{\alpha}=0$ we obtain the simple relations

[TABLE]

The final expressions are used to impose the values of the evolution variables at the gridpoint corresponding to $\mathrsfs{I}^{+}$ (except for $\gamma_{rr}$, which is determined by its equation of motion) after each time step.

6.2.2 Apply l’Hôpital rule to the divergent RHSs at $\mathrsfs{I}^{+}$:

Due to their considerable length, we are not including the equations of motion that are implemented numerically at $\mathrsfs{I}^{+}$ in our code, but we describe how to calculate them. The following equation is taken as an example: variable X’s RHS has a regular part at $\mathrsfs{I}^{+}$, denoted by A (it depends on other evolution variables and on $r$), as well as two formally singular terms at $\mathrsfs{I}^{+}$ (those divided by a power of $\Omega$):

[TABLE]

The l’Hôpital rule is applied on the diverging terms as

[TABLE]

where the last expression is regular at future null infinity. This limit is substituted in the RHS to be evaluated at $\mathrsfs{I}^{+}$ as

[TABLE]

This procedure is followed for each of the diverging RHSs. In the presence of terms that diverge as $\Omega^{-2}$ the l’Hôpital rule is applied twice:

[TABLE]

The regularity conditions (37) (as well as those for $\left.A_{rr}\right|_{\mathrsfs{I}^{+}}$ and $\left.\Lambda^{r}\right|_{\mathrsfs{I}^{+}}$) and (38) are substituted into the “l’Hôpitalized” RHSs at the analytical level.

6.3 Results and comparison to staggered grid

Comparison of the performance of the staggered and non-staggered grids at $\mathrsfs{I}^{+}$ is shown in figure 5 as a convergence test. The corresponding simulations used the same initial data as described in section 5. The Z4c formulation as well as the gauge conditions satisfying the preferred conformal gauge (section 4 and B) were used. The curves correspond to the differences between the signals at $\mathrsfs{I}^{+}$ in simulations with 1600 (low resolution), 2400 (medium) and 3600 (high) points. The non-staggered data (in blue/gray) has been shifted 0.75 forward in time to be able to compare both sets of data in a clearer way. The expected 4th order convergence is very good in both cases, even if it looks slightly better in the non-staggered one.

7 Conclusions

In [39] we presented our first approach to the hyperboloidal initial value problem in spherical symmetry, describing a very basic setup in terms of gauge conditions: harmonic slicing and fixed shift. The preferred conformal gauge was not satisfied, which made an affine reparameterization of the time coordinate at $\mathrsfs{I}^{+}$ more cumbersome, and future null infinity was not included in the numerical grid.

As part of the development of a robust framework to solve the hyperboloidal initial value problem, here we have successfully implemented the preferred conformal gauge in the form of hyperbolic gauge conditions (47) that enforce (3), without altering our time-independent conformal factor $\Omega$ (5) and the simplifications that the latter adds to our setup. These gauge conditions were also designed to systematically include appropriate source functions, required for the hyperboloidal evolution due to its non-trivial background. In our current implementation the time coordinate at $\mathrsfs{I}^{+}$ is related (10) in a non-trivial way to the Bondi time. The condition (11) that makes them equal is satisfied with the parameter choice $\xi_{\alpha}=0$, but this is not yet achievable in our simulations for stability reasons. Thus, for the moment our code does not directly provide an affine time coordinate at $\mathrsfs{I}^{+}$ at all times. However, the preferred conformal gauge makes the “a posteriori” calculation of the Bondi time much easier, as compared to the simpler setup in [39].

The study of the different slicing and shift conditions presented in this paper has increased our understanding of how to adapt current common options of hyperbolic gauge conditions (like the Bona-Massó family or the Gamma-driver) to the free evolution of the hyperboloidal initial value problem. The main ingredients are still those presented in [39]: the addition of damping terms and source functions. Rather than in the region near $\mathrsfs{I}^{+}$, the tested slicing and shift conditions show the most representative differences in behaviour in the interior part of the domain, towards the origin. This effect will be more relevant in the treatment of spacetimes that include a black hole or very large perturbations of regular spacetimes, and will be described in [45]. Our current setup, where we choose not to have incoming gauge modes at future null infinity, requires physical characteristic speeds at $\mathrsfs{I}^{+}$, so that the matching (using (30) or (32)) of at least some coefficients in the gauge conditions is necessary. Where in the integration domain and how this matching is performed is still an open question. At the current stage it is difficult to determine which gauge conditions are preferable, as a change in some of the input parameters, the type of initial data (regular or including a black hole) or the used formulation of the Einstein equations can have an important impact on the performance of different slicing and shift conditions. To exemplify the influence that gauge conditions have on the signals extracted at future null infinity, we compared the values of the rescaled scalar field at $\mathrsfs{I}^{+}$ for several combinations of slicing and shift conditions and included the undistorted signal expressed in terms of the Bondi time in figure 2. These examples served to illustrate and understand the differences in behaviour of the presented gauge choices: the main conclusions are that both an evolved shift, as opposed to a fixed time-independent shift, and larger characteristic speeds in the inner part of the domain (towards the origin) cause considerably less distortion at $\mathrsfs{I}^{+}$. We hope that the different lapse and shift equations compiled here will start to conform a testbed of hyperbolic gauge conditions to study the hyperboloidal problem, and that the insights gained will be useful to future attempts at solving the problem in 3 spatial dimensions.

We have also added a relevant technical improvement to our implementation: the inclusion of a point on $\mathrsfs{I}^{+}$ in the integration domain. This has greatly benefited from the cancellations introduced by the preferred conformal gauge (3), although a gridpoint on $\mathrsfs{I}^{+}$ should also be feasible without this choice. The point on $\mathrsfs{I}^{+}$ represents a step forward in both our understanding of the behaviour of our formulation of the Einstein equations at future null infinity and the degrees of freedom available there: only one of the evolution variables can be freely specified at $\mathrsfs{I}^{+}$. Evaluating the equations on $\mathrsfs{I}^{+}$ requires the adaptation of the divergent terms in the RHSs of the equations of motion as described in subsection 6.2.2, but extrapolation is no longer needed to evaluate the evolved quantities at future null infinity. We have obtained consistent confirmation (see figure 5) that our previous results with a staggered grid were already very good, although the non-staggered grid performs even better. Having a gridpoint on $\mathrsfs{I}^{+}$ will allow us to experiment with different boundary conditions for super-luminal gauge speeds and even to implement a reflecting-like boundary to perform simulations in asymptotically Anti-deSitter spacetimes.

We thank David Hilditch for valuable comments on the manuscript. AV was supported by the FPU-grant AP2010-1697 and partially by the European Research Council Consolidator Grant 647839. AV and SH acknowledge the support of the European Union FEDER founds, the MINECO grants FPA2010-16495, FPA2013-41042-P and CSD2009-00064 of the Spanish Ministry of Economy and Competitiveness and the Conselleria d’Economia i Competitivitat of the Govern de les Illes Balears. Most of the algebraic derivations were performed using the Mathematica package xAct [81] and some of the numerical simulations took advantage of computational resources provided by the STFC grant ST/I006285/1.

Appendix A Conformally compactified equations, notation and variables

In our simulations we use either the Generalized BSSN formulation [21, 22, 82] or a similar conformal version of the Z4 formulation [23, 83, 84], the Z4c equations [85, 86], in their spherically symmetric reduction. The derivation of the equations for the conformally rescaled metric is described in [39] and in Chapter 2 of [54] and the actual equations used in the simulations can be found in Appendix C of [39] and again in Chapter 2 of [54]. Our notation for the metrics is the same one as used in [39]: the 4-dimensional physical metric is denoted as $\tilde{g}$, the 4d-conformal metric as $\bar{g}$, the 3-dimensional conformal metric (induced by $\bar{g}$) as $\bar{\gamma}$, the 3d twice conformal metric as $\gamma$ and the 3d twice conformal background metric as $\hat{\gamma}$. The massless scalar field $\tilde{\Phi}$ satisfies the wave equation in the physical spacetime, and with respect to the conformally rescaled metric it satisfies

[TABLE]

The dynamical variables are the conformally rescaled spatial metric $\gamma_{ab}=\chi\bar{\gamma}_{ab}$, where $\bar{\gamma}_{ab}$ is the spatial metric induced from $\bar{g}_{ab}$, and the spatial conformal factor $\chi$. From the decomposition of the conformal extrinsic curvature tensor $\bar{K}_{ab}$ we evolve its conformal trace-free part $A_{ab}=\chi\bar{K}_{ab}-\frac{1}{3}\gamma_{ab}\bar{K}$ (with $\bar{K}=\bar{K}_{ab}\bar{\gamma}^{ab}\equiv K_{ab}\gamma^{ab}$) and its physical trace, mixed with the physical Z4 variable $\tilde{\Theta}$, $\tilde{K}=\Omega\bar{K}-{\textstyle\frac{3\beta^{a}\partial_{a}\Omega}{\alpha}}-2\tilde{\Theta}$. The actual evolution variable in the code is a variation with respect to its initial value $\Delta\tilde{K}=\tilde{K}-\tilde{K}_{0}=\tilde{K}-K_{CMC}$ (this last parameter corresponds to the value chosen for $\tilde{K}$ on the initial hyperboloidal constant-mean-curvature (CMC) slice and in our notation it will be a negative value). The Z4 quantity $\tilde{\Theta}$ is evolved as well. The evolution of the Z4 variable $Z_{a}$ is included into that of the vector $\Lambda^{a}=\gamma^{bc}\left(\Gamma^{a}_{bc}-\hat{\Gamma}^{a}_{bc}\right)+2\gamma^{ab}Z_{b}$, where $\Gamma^{a}_{bc}$ are the Christoffel symbols calculated from $\gamma_{ab}$ and $\hat{\Gamma}^{a}_{bc}$ the ones built from a time-independent background metric $\hat{\gamma}_{ab}$. The evolved gauge variables are the conformal lapse $\alpha$, the shift $\beta^{r}$ and, if required, its auxiliary variable $B^{a}$. The scalar field is expressed as $\bar{\Phi}=\tilde{\Phi}/\Omega$, together with its auxiliary variable $\bar{\Pi}=\tilde{\Pi}/\Omega=(\partial_{t}\tilde{\Phi})/\Omega$, mainly because these quantities do not vanish at $\mathrsfs{I}^{+}$ and so are more convenient for visualization purposes.

Our spherically symmetric ansatz for the line element is

[TABLE]

The freedom introduced by the spatial conformal factor $\chi$ is fixed by eliminating $\gamma_{\theta\theta}=\gamma_{rr}^{-1/2}$ in the equations implemented in the code. After explicitly imposing the trace-freeness condition of $A_{ab}$ and substituting (46), the only remaining independent component of $A_{ab}$ is $A_{rr}$. The only non-vanishing component of the quantities $\Lambda^{a}$, $\beta^{a}$, $B^{a}$ and $Z_{a}$ is the radial one, denoted respectively by $\Lambda^{r}$, $\beta^{r}$, $B^{r}$ and $Z_{r}$. Primes denote derivatives with respect to the compactified radial coordinate $r$ and dots denote time derivatives.

Appendix B Gauge conditions for the preferred conformal gauge

The decomposition of (33) into lapse and shift evolution equations is

[TABLE]

The values of the background metric components $\hat{\alpha}$, $\hat{\beta^{r}}$, $\hat{\chi}$, $\hat{\gamma_{rr}}$ and $\hat{\gamma_{\theta\theta}}$ are given by (21). In the following simplified version of the previous equations, $\hat{\gamma_{rr}}$ and $\hat{\gamma_{\theta\theta}}$ are set to unity, the evolution equation $\dot{\alpha}$ is substituted in $\dot{\beta}^{r}$ and $\Lambda^{r}$ (determined by setting the $Z_{r}$ constraint, (C.2c) in [39], to zero) is used to eliminate the remaining $\gamma_{rr}^{\prime}$ and $\gamma_{\theta\theta}^{\prime}$ terms in $\dot{\beta}^{r}$ to make the system hyperbolic:

[TABLE]

In the code we set the free specifiable source functions to $\tilde{F}^{r}=0$ and $\tilde{F}^{t}=\xi_{\alpha}\Omega(\alpha-\hat{\alpha})$. The equation of motion for $\dot{\chi}$ ($\gamma_{\theta\theta}$ is eliminated in terms of $\gamma_{rr}$, which makes the $\dot{\gamma}_{rr}$ and $\dot{\gamma}_{\theta\theta}$ terms cancel) is to be substituted above.

Appendix C Previous slicing conditions

Slicing condition (20) makes it easy to implement the harmonic, 1+log and cK conditions (or a matching between them), by providing source functions and damping terms that are valid for any of them. Here we present some of the previous equations of motion for the lapse that were used in our implementation before the simple equation (20) was found. In general, $\xi=1,2$ is a good choice for the $\xi_{harm}$, $\xi_{1+log}$ and $\xi_{cK}$ parameters.

C.1 Defined in the physical domain:

The harmonic slicing takes the following form:

[TABLE]

The $\Delta\tilde{K}$ term in the 1+log condition that follows has been divided by $\Omega$ in order to prevent the lapse gauge speeds from vanishing at $\mathrsfs{I}^{+}$ (this is an effect of the transformation from the physical to the conformal domain):

[TABLE]

In the cK case, it is by $\Omega^{2}$ that the $\Delta\tilde{K}$ terms needs to be divided:

[TABLE]

C.2 Defined in the conformal domain:

The difference with the previous set of equations is that the last damping term has a different coefficient (3 instead of 1). The harmonic gauge is implemented as:

[TABLE]

Again, the damping is controlled by the parameter $\xi_{harm}$. In the 1+log condition that follows, the equivalent damping term formally diverges at $\mathrsfs{I}^{+}$:

[TABLE]

The same is valid for the cK condition, whose difference with the 1+log one is the coefficient in front of $\Delta\tilde{K}$:

[TABLE]