Criticality of inhomogeneous Nariai-like cosmological models

TL;DR

This paper constructs and analyzes inhomogeneous vacuum solutions with positive cosmological constant, revealing oscillatory dynamics at the boundary between collapse and expansion, using both analytical and numerical methods.

Contribution

It introduces new inhomogeneous solutions generalizing the Nariai model and explores their oscillatory behavior through combined analytical and numerical approaches.

Findings

Discovery of oscillatory dynamics in inhomogeneous Nariai-like models

Analytical and numerical evidence supporting the oscillatory behavior

Solutions at the threshold between collapse and exponential expansion

Abstract

In this paper, we construct and study solutions of Einstein's equations in vacuum with a positive cosmological constant which can be considered as inhomogeneous generalizations of the Nariai cosmological model. Similar to this Nariai spacetime, our solutions are at the borderline between gravitational collapse and de-Sitter-like exponential expansion. Our studies focus in particular on the intriguing oscillatory dynamics which we discover. Our investigations are carried out both analytically (using heuristic mode analysis arguments) and numerically (using the numerical infrastructure recently introduced by us).

Tables

Table 1. Table 1: Critical solutions for ϵ = − 10 − 4 italic-ϵ superscript 10 4 \epsilon=-10^{-4} and various ℓ ℓ \ell : Oscillatory behavior of δ ∗ , ℓ subscript 𝛿 ℓ \delta_{*,\ell} . Here, T L subscript 𝑇 𝐿 T_{L} is the prediction from Eq. ( 3.18 ). T N 1 subscript 𝑇 subscript 𝑁 1 T_{N_{1}} is the actual period length of the first full oscillation after t = 0 𝑡 0 t=0 , and T N 2 subscript 𝑇 subscript 𝑁 2 T_{N_{2}} and T N 3 subscript 𝑇 subscript 𝑁 3 T_{N_{3}} of the second and third one. Moreover, t L subscript 𝑡 𝐿 t_{L} is the predicted time of the first oscillation maximum according to Eq. ( 3.15 ) and t N subscript 𝑡 𝑁 t_{N} the actual numerical value.

Data	$t_L$	$t_N$	$T_L$	$T_{N_1}$	$T_{N_2}$	$T_{N_3}$
$\ell =2$	$0.723221$	$0.69999$	$3.14159$	$3.19999$	$3.5499997$	…
$\ell =3$	$0.47418$	$0.46662$	$1.98692$	$2.09979$	$2.3997598$	$2.3997598$
$\ell =4$	$0.359534$	$0.333$	$1.48096$	$1.5651$	$1.7649001$	$2.1644998$

Table 2. Table 2: Critical solutions for ϵ = − 10 − 4 italic-ϵ superscript 10 4 \epsilon=-10^{-4} and various ℓ ℓ \ell : Oscillatory behavior of δ ∗ , 0 subscript 𝛿 0 \delta_{*,0} . Here, T L subscript 𝑇 𝐿 T_{L} is the prediction from Eq. ( 3.20 ). T N 1 subscript 𝑇 subscript 𝑁 1 T_{N_{1}} is the actual period length of the first full oscillation after t = 0 𝑡 0 t=0 , and T N 2 subscript 𝑇 subscript 𝑁 2 T_{N_{2}} and T N 3 subscript 𝑇 subscript 𝑁 3 T_{N_{3}} of the second and third one.

Data	$T_L$	$T_{N_1}$	$T_{N_2}$	$T_{N_3}$
$\ell =2$	$1.5708$	$1.65$	$1.7500002$	$1.75$
$\ell =3$	$0.993459$	$1.1665499$	$1.1105$	$0.9665699$
$\ell =4$	$0.74048$	$0.7992$	$0.86580002$	$0.73259997$

Equations

$$G_{ab}+\Lambda g_{ab}=0,$$

$$g_{ab}=\frac{1}{\Lambda }(-\mathrm{d}{t}^2+{\cosh }^{2}t\mathrm{d}{\rho }^2+g_{{\mathbb{S}}^2}).$$

$${\nabla }_a{\nabla }^a\psi$$

$${\nabla }_a{\nabla }^a\omega$$

$$R_{ab}$$

$$\psi :=g_{ab}{\xi }^a{\xi }^b,$$

$${\nabla }_a\omega ={\Omega }_a:={ϵ}_{abcd}{\xi }^b{\mathfrak{D}}^c{\xi }^d.$$

$${\hat{h}}_{ab}:=g_{ab}-\frac{1}{\psi }{\xi }_a{\xi }_b,$$

$$h_{ab}:=\psi {\hat{h}}_{ab}.$$

$${\hat{R}}_{ab}:=R_{ab}+{\nabla }_{(a}{\mathcal{D}}_{b)},$$

$${\mathcal{D}}^{\alpha }:=(-{\Gamma }^{\alpha }{}_{\beta \gamma }+{\bar{\Gamma }}^{\alpha }{}_{\beta \gamma })h^{\beta \gamma }+f^{\alpha }.$$

$${\nabla }_a{\nabla }^a\psi$$

$${\nabla }_a{\nabla }^a\omega$$

$${\hat{R}}_{ab}$$

$${\nabla }_b{\nabla }^b{\mathcal{D}}_a-{\mathcal{D}}^b{R}_{ab}=0.$$

$$0={\mathcal{D}}^{\nu }=h^{\rho \sigma }({\bar{\Gamma }}^{\nu }{}_{\rho \sigma }-{\Gamma }^{\nu }{}_{\rho \sigma })+f^{\nu },$$

$${\nabla }_{\mu }{\mathcal{D}}_{\nu }=0$$

$$T^a:={\partial }_t^a,m^a:=\frac{1}{\sqrt{2}}\left({\partial }_{\theta }^a-\frac{i}{\sin \theta }{\partial }_{\phi }^a\right).$$

$$({\partial }_0^a,{\partial }_1^a,{\partial }_2^a)=(T^a,m^a,{\overline{m}}^a)$$

$${\omega }_a^0={\nabla }_{a}t,{\omega }_a^1=\frac{1}{\sqrt{2}}\left({\nabla }_a\theta +i\sin \theta {\nabla }_a\phi \right),{\omega }_a^2={\overline{\omega }}_a^1.$$

$$h_{ab}=\lambda {\omega }_a^0{\omega }_b^0+2{\omega }_{(a}^0\left(\beta {\omega }_{b)}^1+\bar{\beta }{\omega }_{b)}^2\right)+2\delta {\omega }_{(a}^1{\omega }_{b)}^2+\varphi {\omega }_a^1{\omega }_b^1+\bar{\varphi }{\omega }_a^2{\omega }_b^2.$$

$$f_0=f_a{T}^a,f_1=f_a{m}^a,f_2=f_a{\bar{m}}^a$$

$${\bar{h}}_{ab}=-{\omega }_a^0{\omega }_b^0+2{\omega }_{(a}^1{\omega }_{b)}^2$$

$$n^a=\frac{1}{\alpha }\left(T^a-\beta (m^a+{\overline{m}}^a)\right),n_a=-\alpha {\omega }_a^0;$$

$$\alpha =\sqrt{{\beta }^2-\lambda }.$$

$${\gamma }_{ab}=h_{ab}+n_a{n}_b={\beta }^2{\omega }_a^0{\omega }_b^0+2\beta {\omega }_{(a}^0\left({\omega }_{b)}^1+{\omega }_{b)}^2\right)+2\delta {\omega }_{(a}^1{\omega }_{b)}^2+\varphi \left({\omega }_a^1{\omega }_b^1+{\omega }_a^2{\omega }_b^2\right),$$

$$\alpha =\delta ,\beta =0,$$

$$\lambda =-\delta ,\beta =0.$$

$${\partial }_{tt}\delta +a{\partial }_{\theta \theta }\delta +b{\partial }_{\theta }{f}_1=...,$$

$${\partial }_{tt}\varphi +a{\partial }_{\theta \theta }\varphi +b{\partial }_{\theta }{f}_1=...,$$

$${\partial }_{tt}\psi +a{\partial }_{\theta \theta }\psi =...,$$

$${\partial }_{tt}\omega +a{\partial }_{\theta \theta }\omega =...,$$

$${\partial }_{tt}\lambda =-{\partial }_{tt}\delta$$

$${\partial }_{tt}\beta =0$$

Full text

1 Introduction

The Nariai spacetime was discovered by Nariai in 1950 (see the reprints of the original works in [34, 35]). It is the solution of Einstein’s vacuum equations444We use the signature convention $(-,+,+,+)$ for spacetime metrics and the sign convention for curvature tensors in [45]. In this convention the de-Sitter spacetime is a solution of Eq. (1.1) with $\Lambda>0$.

[TABLE]

where555In this whole paper, we consider $\Lambda$ either as any positive constant or set $\Lambda=1$. $\Lambda>0$ is the cosmological constant and $G_{ab}$ is the Einstein tensor associated with the metric $g_{ab}$, with spatial topology $\mathbb{S}^{1}\times\mathbb{S}^{2}$ and

[TABLE]

Here, $\rho$ is the standard coordinate along the spatial $\mathbb{S}^{1}$-factor and $g_{\mathbb{S}^{2}}$ is the metric of the standard round unit two-sphere. The time coordinate is $t\in(-\infty,\infty)$. In the last years, the Nariai spacetime has become an object of special interest since Ginsparg and Perry [24] found that it can be interpreted as a de-Sitter universe containing a black hole of “maximal size”. Thanks to its geometrical properties, the Nariai spacetime has been used to model several situations. One of the most remarkable applications was carried out by Bousso and Hawking [10, 11, 12, 9] who used this spacetime to study the quantum pair creation of black holes during inflation. These cosmological models, at the borderline between inflation and gravitational collapse, restricted to spherically symmetric perturbations of the Nariai spacetime. It was found that under certain conditions those models asymptotically approach the de Sitter universe in agreement with the cosmic no-hair conjecture [23, 27]. It states that this behavior is generic for inhomogeneous and anisotropic expanding solutions of Eq. (1.1). Although there is some mathematical evidence [29, 40, 44, 1] that supports the validity of this conjecture, the general case still remains unclear. A particular property of the Nariai spacetime itself is its peculiar time dependence which is not consistent with this. While the spatial $\mathbb{S}^{1}$-factor expands exponentially for large $t$, the geometry of the spatial $\mathbb{S}^{2}$-factor remains constant. Thus, the expansion of this solution is very anisotropic, and, in fact, is inconsistent with the cosmic no-hair paradigm. In more geometric terms, this corresponds to the fact, which was proven for the first time in [3] (an alternative proof was given in [20]), that the Nariai spacetime does not possess even a piece of a smooth conformal boundary. If the cosmic no-hair conjecture really holds, the Nariai spacetime must be therefore “very special”, and hence in particular be unstable under “generic” perturbations.

In [3] one of us has initiated the study of homogeneous (but fully nonlinear) perturbations of the Nariai spacetime. The Nariai solution is a member of the class of Kantowski-Sachs spatially homogeneous (but anisotropic) solutions [43] of Einstein’s vacuum equation with a positive cosmological constant. The spatial topology of all of these models is $\mathbb{S}^{1}\times\mathbb{S}^{2}$. It was found there that the Nariai solution is critical in this family in the following sense. For all Kantowski-Sachs models, except for the Nariai solution, we can choose the time orientation such that the solution either collapses to the future (a big crunch characterized by the formation of a curvature singularity and all future inextendible causal curves are incomplete) or expands eternally to the future in consistency with the cosmic no-hair picture (existence of a smooth future conformal boundary in consistency with the future asymptotics of de-Sitter space so that all future inextendible causal curves are complete). The Nariai solution is exactly at the borderline between these two extremes as the curvature is bounded everywhere and all inextendible causal curves are both future and past complete, but it nevertheless does not agree with cosmic no-hair.

The first rigorous work in [3] on this topic has recently been extended in [15]. The numerical studies in [4] of Gowdy-symmetric [25, 14] (see Section 2.1 for more details on Gowdy symmetry) inhomogeneous fully nonlinear perturbations of the Nariai solution have revealed evidence that the analogous critical phenomenon also exists in much larger classes of spacetimes. In particular, it was found that all solutions, which are obtained from initial data not too far away from the Nariai solutions, always either globally collapse or expand in the same manner as in the spatially homogeneous case — with the exception of critical solutions which are exactly at the borderline between these two cases. In particular, it was found that in contrast to the spherically symmetric models considered by Bousso et al. above [9], Gowdy symmetric models never locally collapse or expand, and the formation of cosmological black holes in this class was therefore ruled out. Because the perturbations considered in [3] were small in some sense, the question remained open whether this may be different for larger Gowdy symmetric perturbations. One of the finding in our paper here now suggests that Gowdy symmetric models indeed never form cosmological black holes. For future work, it will be interesting to pose this question again within more general classes of spacetimes, for example $U(1)$-symmetric spacetimes, and study whether cosmological black holes may be created by perturbations of the Nariai spacetime.

Before we continue, let us remind the reader about the heuristic idea of the study of the criticality of the cosmological models in [4]. There we worked with Gowdy symmetric initial data (which satisfy the constraint equation implied by Eq. (1.1)) given by two real parameters $\mu$ and $\Sigma_{\times}^{(1)}$ whose precise definition is irrelevant now (cf. [4] for the details). The special choice $\mu=\Sigma_{\times}^{(1)}=0$ corresponds to Nariai initial data while $\Sigma_{\times}^{(1)}=0$ and $\mu\in\mathbb{R}$ yields a class of spatially homogeneous models. The larger the value of $|\Sigma_{\times}^{(1)}|$ is, however, the “more spatially inhomogeneous” the initial data are. The idea was to fix some non-zero value of the “inhomogeneity parameter” $\Sigma_{\times}^{(1)}$ and then to study a sequence of (fully nonlinear) cosmological models given by a sequence of values of $\mu$. On the one hand, it was found that if $\mu$ is sufficiently large, the corresponding model expands globally to the future; in fact, the solution develops a smooth conformal boundary to the future in this case and is hence fully consistent with the cosmic no-hair conjecture666Notice that the studies in [4] made use of Friedrich’s conformal field equations [19] and therefore allowed us to calculate the full conformally compactified solutions, including the conformal boundary if it exists. Here, we will not make use of such compactification techniques.. If $\mu$ is sufficiently small on the other hand, the model collapses globally to the future and eventually forms a curvature singularity. At the borderline between these two regimes corresponding to a critical value for $\mu$, the corresponding model neither collapses nor expands to the future. However, no further information about such critical models was extracted in [4].

The purpose of our present paper is manifold. Again, we restrict to the class of Gowdy-symmetric models with a positive cosmological constant and we revisit the same situation, but tackle it with a more advanced approach. To this end, we use a different class of initial data which now depends on three parameters $\epsilon$ (which has a similar meaning as the “inhomogeneity parameter” $\Sigma_{\times}^{(1)}$ above), $C$ (which has a similar meaning as the parameter $\mu$ above) and an additional parameter $\ell$ which essentially controls the wave number of the inhomogeneous perturbation (the initial data in [4] restricted to the case $\ell=2$). The details are discussed in Section 2.3. On the one hand, we confirm and strengthen the findings in [4] by performing a similar numerical analysis. On the other hand, however, we shall focus in great detail on the critical solutions here and thereby reveal an interesting new oscillatory phenomenon which could potentially be interpreted as gravitational waves. The main finding of our paper are now summarized as three main results. The purpose of this paper is to provide the details and give justifications.

In all of what follows we shall assume777In this paper, we shall not be interested in the homogeneous perturbations associated with the case $\ell=0$. The case $\ell=1$ is a special borderline case which turns out to be not well described by our analytic method discussed in Section 3.1. We therefore completely disregard the case $\ell=1$ in this paper. $\ell\geq 2$. Notice that the well-understood [3] homogeneous case of our models corresponds to $\epsilon=0$. The Nariai solution is determined by $\epsilon=C=0$. One can easily check that if $\epsilon\not=0$ or $C\not=0$, the corresponding solution of Eq. (1.1) is not isometric to the Nariai solution by comparing the Kretschmann scalar with the particular globally constant value for the Nariai solution. Our first main finding is summarized as follows.

Result 1

Pick any real value $\epsilon$ and integer $\ell\geq 2$. Then there is a constant $C_{crit}$ such that the solution of Eq. (1.1), determined by initial data given by the parameter $\epsilon$, $\ell$ and any real value $C$ as in Section 2.3, globally collapses and forms a curvature singularity if $C>C_{crit}$ and globally expands in consistency with the cosmic no-hair conjecture if $C<C_{crit}$.

For small values of $\epsilon$, this result is in full consistency with the findings in [4]. Here we claim now that this also holds for large values of $\epsilon$. As mentioned earlier, this rules out in particular the possibility of “local” collapse and hence there are generically indeed only two kinds cosmological models in this class. In this paper here we provide some more refined numerical evidence complemented by a heuristic mode analysis (see Section 3.1). We call any of our models critical if $C=C_{crit}$ for any given $\epsilon$ and $\ell$, and almost critical or close-to-critical if $C\not=C_{crit}$, but $|C-C_{crit}|\ll 1$.

The second main finding of our work which significantly goes beyond the results in [4] is summarized as follows.

Result 2

For any non-zero value of $\epsilon$, the critical and close-to-critical solutions asserted in Result 1 are oscillatory.

Based on the before-mentioned heuristic analysis in Section 3.1, we are able to derive formulas for oscillation frequencies, amplitudes and phases and how these depend on the initial data parameters. The only non-oscillatory solutions correspond to the spatially homogeneous case $\epsilon=0$ in which the critical solution is known to coincide with the Nariai solution (see [3]) and therefore has the peculiar asymptotics discussed above. Our numerical work here suggests that all the critical models, also the inhomogeneous ones, are Nariai-like in the following sense.

Result 3

The critical solutions behave asymptotically as follows. While the spatial $\mathbb{S}^{1}$-factor expands exponentially, the spatial $\mathbb{S}^{2}$-factor geometry oscillates around the round unit $2$-sphere geometry and is therefore in particular bounded. All these models therefore violate the cosmic no-hair picture by these highly anisotropic asymptotics.

We emphasize that our work here is not actually concerned with the instability of the Nariai solution; this issue is addressed elsewhere [3, 15]. The point of our work here is now to identify and describe inhomogeneous critical models and their Nariai-like asymptotics all of which violate the cosmic no-hair paradigm. We remark that all our numerical studies were conducted with slight generalizations of the numerical code presented in [8]. More details and references regarding our numerical infrastructure are given in Section 2.2.

The paper is organized as follows. In Section 2, we discuss the general setup, i.e., the formulation of Einstein’s equation in the presence of symmetries via Geroch’s symmetry reduction, our extraction of evolution equations with a well-posed initial value problem and of constraint equations from this, our numerical implementation and our particular family of initial data. Section 3 is devoted to our analytical and numerical studies. First we discuss our heuristic mode analysis which is the basis for all of what follows. Then we provide numerical evidence for all of the main results above.

2 Setup: Formulation and implementation of Einstein’s equations

2.1 Einstein’s vacuum equations for Gowdy symmetry and spatial $\mathbb{S}^{1}\times\mathbb{S}^{2}$-topology

Geroch’s symmetry reduction

As mentioned earlier, we shall focus here on Gowdy symmetric spacetimes with spatial topology $\mathbb{S}^{1}\times\mathbb{S}^{2}$. We start by equipping the spatial manifold $\mathbb{S}^{1}\times\mathbb{S}^{2}$ with coordinates $(\rho,\theta,\varphi)$ where $\rho\in(0,2\pi]$ is the standard parameter on the $\mathbb{S}^{1}$-factor and $(\theta,\varphi)$ are standard polar coordinates on the $\mathbb{S}^{2}$-factor. With respect to these coordinates it turns out that a spacetime with spatial $\mathbb{S}^{1}\times\mathbb{S}^{2}$-topology is Gowdy (or $\mathrm{U(1)}\times\mathrm{U(1)}$-)symmetric if the metric is invariant both under translations along the $\mathbb{S}^{1}$-factor (i.e., is independent of $\rho$), and, under rotations of the $\mathbb{S}^{2}$-factor around the polar coordinate axis (i.e., is independent of $\varphi$). One can show that Einstein’s equations for this class of spacetimes can be formulated in almost exactly the same way as for the class of Gowdy symmetric spacetimes with spatial $\mathbb{S}^{3}$-topology which was discussed in detail in [8]. Because of the close similarity we shall refer to that paper for all the details and only give a quick summary of the necessary results now.

When Geroch’s symmetry reduction [21] is performed with respect to the Gowdy Killing vector fields $\xi^{a}=\partial_{\rho}^{a}$ for any $3+1$-dimensional Gowdy-symmetric metric $g_{ab}$ in a spacetime $M=\mathbb{R}\times\mathbb{S}^{1}\times\mathbb{S}^{2}$, one finds that Einstein’s vacuum field equations with cosmological constant (1.1) imply the system

[TABLE]

on the $2+1$-manifold $S=\mathbb{R}\times\mathbb{S}^{2}$ where $h_{ab}$ is a metric on $S$ with signature $(-,+,+)$, $\nabla_{a}$ is its covariant derivative and $R_{ab}$ is the corresponding Ricci tensor. The scalar field $\psi$ is defined on $M$ as

[TABLE]

for $\xi^{a}=\partial_{\rho}^{a}$ which is then projected888In order to simplify the notation of [8] slightly, we shall not distinguish here between quantities on $S$ and their counterparts on $M$ which are obtained by a pullback along the projection map $\pi$. In fact, all quantities, which carry a $\tilde{\,\,}$ in [8], shall be written without a $\tilde{\,\,}$ here. to $S$. The other scalar field $\omega$ is the well-defined global potential of the twist $1$-form $\Omega_{a}$ on $M$ defined by

[TABLE]

Here, $\mathfrak{D}_{c}$ is any derivative operator (for instance the covariant derivative associated with $g_{ab}$) and $\epsilon_{abcd}$ is a volume form associated with $g_{ab}$. We shall often refer to $\psi$ as the norm and to $\omega$ as the twist of $\xi^{a}$ respectively. Eqs. (2.1) can be interpreted as the $2+1$-dimensional Einstein equations999Strictly speaking, this is only the case when $\Lambda=0$. When $\Lambda\not=0$, the second term on the right-hand side of the third equations differs from a standard cosmological constant term. This will however not play any role in our discussions here. coupled to two scalar fields $\psi$ and $\omega$.

Without going into the details, see for example [8], let us mention that once a solution $(h_{ab},\psi,\omega)$ of Eqs. (2.1) has been found on $S$, one can construct the corresponding physical spacetime metric $g_{ab}$ on $M$ which then solves Eq. (1.1). It is important to notice that $g_{ab}$ and $h_{ab}$ are related as follows

[TABLE]

and

[TABLE]

The metric $h_{ab}$ in Eqs. (2.1) is therefore not the physical $2+1$-metric, but is related by the conformal transformation Eq. (2.5) to the physical $2+1$-metric $\hat{h}_{ab}$ defined by Eq. (2.4). Both metrics $h_{ab}$ and $\hat{h}_{ab}$ will play a role in our discussion later.

The system of equations (2.1) takes care of only one of the two Gowdy Killing vector fields $\xi^{a}=\partial_{\rho}^{a}$ (i.e., translations along the spatial $\mathbb{S}^{1}$-factor) so far. It turns out that the second Killing field $\partial_{\varphi}^{a}$ prevails on $S$, i.e., all quantities $h_{ab}$, $\psi$ and $\omega$ on $S$ defined above are axi-symmetric and therefore invariant under the action of $\partial_{\varphi}^{a}$ (i.e., under rotations around the polar coordinate axis of $\mathbb{S}^{2}$). We remark that one could perform Geroch’s symmetry reduction a second time, but now with respect to this Killing field (see also [22]). This however leads to explicit singularities at the poles of the two-sphere. As in [8], we shall therefore work in all of what follows with axially symmetric (i.e., $\varphi$-invariant) solutions of Eqs. (2.1) without any further symmetry reductions.

The (generalized) wave map formalism

The next task is to extract suitable evolution and constraint equations from Eqs. (2.1) in order to obtain a well-posed initial value problem. The first two equations of (2.1) are scalar wave equations. It therefore remains to deal with the third equation of (2.1). Since this is the $2+1$-Einstein equation for the metric $h_{ab}$ with a (as one can check) divergence free energy momentum tensor of the “matter source”, we can apply all kinds of standard techniques which were developed for $3+1$-Einstein’s equations. Because of its geometric nature, which is particularly useful for dealing with the spatial topology $\mathbb{S}^{2}$ of $S$, we work with the generalized wave map formalism [18] here. Again, all details are worked out in [8] and we just give a quick summary here.

The point is that the third equation in Eqs. (2.1) is a-priori not a system of wave equations for the components of $h_{ab}$ (with respect to any frame) and hence the initial value problem is in general not well-posed. This problem is overcome when we replace $R_{ab}$ in that equation by the new tensor field

[TABLE]

where the components $\mathcal{D}^{\alpha}$ of the vector field $\mathcal{D}^{b}$ with respect to any smooth local frame101010The components of tensor fields with respect to any such frame on $S=\mathbb{R}\times\mathbb{S}^{2}$ are denoted by greek indices. are given as

[TABLE]

The vector field $f^{a}$ can be specified freely and is referred to as a gauge source field. Its components $f^{\alpha}$ with respect to any frame are often called gauge source functions. The connection coefficients of the covariant derivative associated with $h_{ab}$ with respect to this frame are denoted above by $\Gamma^{\alpha}{}_{\beta\gamma}$, while $\bar{\Gamma}^{\alpha}{}_{\beta\gamma}$ are the corresponding connection coefficients associated with any freely specifiable reference metric $\bar{h}_{ab}$ on $S$. In total this produces a (complicated) system of quasilinear wave equations for the components of the metric and the fields $\psi$ and $\omega$:

[TABLE]

In particular, the initial value problem of these evolution equations is well-posed for suitable initial data.

Suppose now that a solution $(h_{ab},\psi,\omega)$ of the initial value problem of the evolution equations has been found on $S$. It is clear that this is a solution of the original system (2.1) if $\mathcal{D}^{a}$ vanishes and hence $\hat{R}_{ab}=R_{ab}$ on $S$. In this case we say that $h_{ab}$ is in generalized wave map gauge. We show now that $\mathcal{D}^{a}$ vanishes only if the initial data for the evolution equations satisfies certain constraints. As discussed for example in [41], the evolution equations, the fact that the energy momentum tensor of the matter source in Eqs. (2.1) is divergence free, and the contracted Bianchi identities together imply

[TABLE]

Since the metric $h_{ab}$ (and hence $R_{ab}$) is considered as known at this stage, this is a linear homogeneous system of wave equations for the unknown $\mathcal{D}_{a}$. It follows that $\mathcal{D}^{a}$ vanishes everywhere on $S$ if and only if $\mathcal{D}^{a}=0$ and $\nabla_{a}\mathcal{D}^{b}=0$ on the initial hypersurface; these two conditions therefore constitute constraints. The first constraint takes the form (with respect to any smooth local frame)

[TABLE]

which can be satisfied for any initial data $h_{ab}$, $\psi$ and $\omega$ on the initial hypersurface by a suitable choice of the free gauge source quantities $f^{a}$ and $\bar{h}_{ab}$. Eq. (2.9) is therefore referred to as the gauge constraint. Once we know that the gauge constraint is satisfied, it turns out that the second constraint above

[TABLE]

is equivalent to the standard Hamiltonian and Momentum constraints which we discuss in more detail in Section 2.3. We emphasize the surprising fact that Eq. (2.10) is not another restriction on the gauge source quantities $f^{a}$ and $\bar{h}_{ab}$; these are only constrained by the gauge constraint. Eq. (2.10) therefore represents the “physical constraint” imposed on the initial data for $h_{ab}$, $\psi$ and $\omega$.

2.2 Formulation and implementation of the evolution equations

Choice of frame and parametrization of the metric

If $t$ is any time function on $S$ and $(t,\theta,\varphi)$ are coordinates as before, we set

[TABLE]

Then we choose

[TABLE]

as our local frame which is defined almost everywhere on $S=\mathbb{R}\times\mathbb{S}^{2}$ (excluding the poles of the two-sphere). Notice, that this frame is in general not an orthonormal frame with respect to $h_{ab}$. It is merely a particular linear combination of the coordinate frame which is motivated by the spin-weight formalism below. In the following we shall express all tensor fields on $\mathbb{S}^{2}$ with respect to this frame and its dual frame $(\omega^{0}_{a},\omega^{1}_{a},\omega^{2}_{a})$ which is given by

[TABLE]

It terms of this frame, we can write

[TABLE]

For the spin-weight formalism [37, 8, 6, 5, 7] we assume that the fields $T^{a}$, $m^{a}$ and $\overline{m}^{a}$ have spin-weights [math], $+1$ and $-1$, respectively, which implies that the spin-weights of $\omega^{0}_{a}$, $\omega^{1}_{a}$, $\omega^{2}_{a}$ are [math], $-1$ and $+1$, respectively. The quantities $\lambda,\delta,\beta$ and $\phi$ in Eq. (2.14) therefore have spin-weights [math], [math], $+1$, $+2$, respectively, and the complex conjugates carry the corresponding negative spin-weights. It is of fundamental importance for all of what follows that once the gauge freedom in terms of the smooth quantities $f^{a}$ and $\bar{h}_{ab}$ has been fixed, the whole system of evolution equations can be written as a quasilinear coupled system of six complex wave equations for the six complex unknowns $\lambda$, $\delta$, $\beta$, $\phi$, $\psi$ and $\omega$. Moreover, once all directional derivatives along $m^{a}$ and $\bar{m}^{a}$ have been replaced by the so-called $\eth$- and $\bar{\eth}$-operators via Eq. (A.4), each term in each equation of this system has a consistent well-defined spin-weight and is explicitly regular at the poles $\theta=0$ and $\theta=\pi$ of the $2$-sphere. Indeed, this explicit regularization of the “pole problem” is the main advantage of the spin-weight formalism.

We recall that Gowdy symmetry implies that $\partial_{\varphi}^{a}$ is a Killing vector field on $S$. Since $\partial_{\varphi}^{a}$ commutes with each of the fields in Eqs. (2.11) and (2.13), it follows that $h_{ab}$, as given by Eq. (2.14), is invariant under the action of $\partial_{\varphi}^{a}$ if and only if all the quantities $\lambda,\delta,\beta$ and $\phi$ are functions of $t$ and $\theta$ only. This means in particular that all these functions can be expanded in terms of axi-symmetric spin-weighted spherical harmonics, see Section A. We also know that all quantities with spin-weight [math] must be real, while other quantities could in principle be complex. However, in the particular representation in $(\theta,\varphi)$-coordinates used exclusively in this whole paper, one can show that if all unknown quantities (and their time derivatives) in the evolution equations are real at the initial time, if all gauge source quantities

[TABLE]

are real for all times, and, if the background metric $\bar{h}_{ab}$ is once and for all chosen as

[TABLE]

for all times, then all unknown quantities are real for all times. Below we see that this restriction to real quantities is purely a gauge restriction. In summary, without further notice we shall assume in all of what follows that all quantities $\lambda$, $\delta$, $\beta$, $\phi$, $\psi$, $\omega$, $f_{0}$, $f_{1}=f_{2}$ are real and only depend on $t$ and $\theta$.

Gauge drivers and conformal time gauge

Instead of fixing the gauge freedom by choosing the gauge source quantities $f_{0}$ and $f_{1}=f_{2}$ as outlined in the previous section, it may sometimes be advantageous numerically to fix the gauge by choosing the “lapse” and “shift” of the $2+1$-dimensional metric $h_{ab}$. The equations which then determine $f_{0}$ and $f_{1}$ are often called gauge drivers. Even though this approach has been used successfully in some situations (see for instance [38, 39]), it may cause numerical instabilities. The reason lies in the fact that the resulting total system of evolution equations (including the gauge drivers) may not have a well-posed initial value problem despite the fact that the original evolution equations in the wave gauge formalism do. Some general proposals for gauge drivers which do not suffer from this problem can be found in [30, 31]. In this work here now, we will construct particular gauge drivers now and then show that the total system of evolution equations is strongly hyperbolic.

To start with, we consider the unit normal vector to the $t=const$-surfaces (recall that $\beta$ is assumed to be real)

[TABLE]

recall Eqs. (2.11), (2.13) and (2.14), where

[TABLE]

We can therefore interpret $\beta$ as the shift111111The shift vector is $\beta^{a}=\beta(m^{a}+\overline{m}^{a})$. and $\alpha$ as the lapse. The induced metric on the $t=\mathrm{const}$-hypersurfaces is therefore (recall that $\phi$ is assumed to be real)

[TABLE]

cf. Eq. (2.14). Now, according to our discussion of gauge drivers above, let us attempt to fix the gauge by picking

[TABLE]

during the whole evolution. This corresponds to

[TABLE]

Heuristically, the idea is that the lapse is proportional to the area $\delta$ of the spatial $2$-sphere. From the point of view of any Eulerian observer, the coordinate clock will therefore tick faster or slower depending on whether the $2+1$-spacetime is expanding or collapsing. An important consequence is that the foliation tends to “freeze” in the collapsing case. This gauge therefore avoids singularities. This sort of gauge is commonly known in the standard cosmology literature as conformal time gauge, and is used frequently in the linear theory of cosmological perturbations [33].

In this work, we wish to implement gauge drivers which preserve this gauge during the evolution. To do so, we use Eq. (2.18) to express the evolution equations for $\lambda$ and $\beta$ in the wave gauge formalism as evolution equations for the gauge source functions $f_{0}$ and $f_{1}$. Hence, from now on, $\lambda$ and $\beta$ will not be considered as unknown variables anymore, but $f_{0}$ and $f_{1}$ will. The question is whether the resulting evolution system is hyperbolic, and, if yes, in which sense. In what follows we consider this question in detail.

We continue to assume that all unknown fields in the evolution equations are real functions, and all the partial derivatives with respect to the coordinate $\varphi$ vanish. Then, expanding the covariant derivatives and expressing the frame vectors $m^{a}$ and $\overline{m}^{a}$ in terms of the coordinate vector $\partial_{\theta}$, we obtain evolution equations for $\delta$, $\phi$, $\psi$ and $\omega$ of the form

[TABLE]

where $a=h^{11}/h^{00}$ and $b=\sqrt{2}/h^{00}$. Note that we have used $f_{2}=f_{1}$ in the evolution equation for $\delta$. The ellipses in the right-hand side of the equations denote lower order terms which are irrelevant for this analysis. Setting $\lambda=-\delta$ and $\beta=0$ we obtain evolution equations for the gauge source functions $f_{0}$ and $f_{1}$, respectively as

[TABLE]

Naturally, these evolution equations, which we call gauge drivers, control the behavior of the generalized gauge source functions such that the conformal time gauge is preserved during the evolution. Next, in order to analyze the hyperbolicity of the resulting system of evolution equations Eqs. (2.19)–(2.24), we rewrite it in first-order form as

[TABLE]

where we have defined the vector

[TABLE]

and where $\Pi$ is a $14\times 14$ (non-symmetric) matrix121212Because of the size of $\Pi$ we do not write it here explicitly. with eigenvalues

[TABLE]

After a straightforward calculation we show that all eigenvectors are linearly independent. Note that $h^{00},-h^{11}>0$ provided that $h_{ab}$ is a Lorentzian metric with signature $(-,+,+)$, thus $a<0$ and hence all the eigenvalues of $\Pi$ are real. This implies that $\Pi$ has a complete set of eigenvectors with real eigenvalues. Our system Eq. (2.25) is therefore strongly hyperbolic. This implies the well-posedness of the initial value problem.

Constraint damping terms

In order to deal with the widely known problem of the growth of constraint violations during numerical evolutions, Brodbeck et al. [13] have suggested a general approach such that the constraint surface is an attractor. Later, following this idea, Gundlach et al. [26] introduced so called constraint damping terms into Einstein’s equations by adding to the right side of Eq. (2.6) the term

[TABLE]

with $\eta_{a}$ being a timelike vector and $\kappa$ a constant. With this new term, the subsidiary equation Eq. (2.8) takes the form

[TABLE]

They showed by means of perturbations of the Minkowski spacetime that all the “short wave length” modes in the solutions of the subsidiary system Eq. (2.27) are damped if $\kappa>0$ at either the rate $e^{-\kappa t}$ or $e^{-\kappa t/2}$. In the last years, a good amount of numerical simulations have been successfully conducted using this approach (see for instance [31, 38]), which confirms its effectiveness for several situations.

Nevertheless, a complete understanding of how the “long wave length modes” solutions are damped (or not) for generic spacetimes is still missing. Due to the expanding (or collapsing) behavior of most cosmological spacetimes, the “long wave length modes” are expected to be dominant during the evolution. For our particular interest here it is therefore a-priori not clear whether these constraint damping terms really improve the evolution of constraint violations. In order to address this question, we simplify our analysis now by using the assumption that the violation covector is approximately only $t$-dependent

[TABLE]

However, since $\mathcal{D}_{\mu}$ is a covector, the projection of its spatial components $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ to the frame $(m^{a},\overline{m}^{a})$ must have spin-weight $1$ and $-1$, which directly implies that $\mathcal{D}_{1}=\mathcal{D}_{2}=0$ under the above assumption because only a function with spin-weight [math] can have a mode of $l=0$ (which is spatially independent). If the perturbed metrics are “close” to the Nariai metric during the initial part of the evolution, we can use it as the background metric for writing the subsidiary equation that rules the evolution of $\mathcal{D}_{0}$. Thus, replacing the Nariai metric in Eq. (2.27) with $\eta_{a}=(-1,0,0)$, we obtain the evolution equation that rules the behavior of $\mathcal{D}_{0}(t)$ as

[TABLE]

for the same constant $\kappa$ as in Eq. (2.26). Evidently, the solution of this equation is

[TABLE]

where $A$ and $B$ are constants. The constraint violation $\mathcal{D}_{0}(t)$ does therefore not grow if $\kappa>0$.

In order to check numerically the validity of the above statements let us consider the constraint error as131313We have excluded $\mathcal{D}_{2}(t)$ from this definition because $\mathcal{D}_{1}(t,\theta)=\overline{\mathcal{D}_{2}(t,\theta)}=\mathcal{D}_{2}(t,\theta)$ in Gowdy symmetry.

[TABLE]

where the norm $\|\>.\>\|_{L^{2}(\mathbb{S}^{2})}$ is numerically computed by using Eq. (2.42). In Figs. 2 and 2, we plot this error obtained for different values of $\kappa$. Here, we have used the initial data family which we will describe in Section 2.3 and pick $\epsilon=-10^{-4}$, $C=0$ and $\ell=2$. In order to calculate $E(t)$ in these figures, we calculate the numerical solution of the evolution equations with constraint damping terms corresponding to these data using the pseudo-spectral method described in [8] with the Runge-Kutta method as time integrator and the axial symmetric spin weighted transform for computing the spatial derivatives (more details on the numerical infrastructure are given below).

As expected, the error for $\kappa=-8$ grows exponentially whereas for $\kappa=8$ it is bounded. Thus, from now on, we will keep this value for all the following numerical calculations in this paper.

Numerical infrastructure

The main difficulty for the numerical treatment of tensorial equations on manifolds with spherical topology is the fact that these cannot be globally covered by a single regular coordinate patch. In the literature this problem is commonly known as the pole problem because in standard polar coordinates for $\mathbb{S}^{2}$ these issues appear at the poles. Based on the previous works [6, 5], we introduced in [8] a pseudo-spectral infrastructure to overcome this issue numerically. It consists in using the spin-weight formalism for expressing tensor components in terms of spin-weighted spherical harmonics, which are a generalization of the well known spherical harmonics [37]; see Appendix A. This allows us to work with polar coordinate representations of Eqs. (2.1) that do not suffer from any polar singularity. This becomes manifest when all the spatial derivatives are replaced by the eth-operators in Eqs. (A.3).

As mentioned earlier, we shall exclusively restrict to Gowdy symmetric models which implies axial symmetry for all the fields in Eqs. (2.1). For the numerical treatment of such fields, we have introduced the one-dimensional variant of the spin-weighted transform introduced by Huffenberger and Wandelt [28], which we call axially symmetric spin-weighted transform in [8]. Our numerical infrastructure is therefore a pseudo-spectral scheme based on the method of lines where the temporal integration is carried out by certain Runge-Kutta integrators.

2.3 Constraints and initial data

Formulation of the constraints and choice of free data

As explained in Section 2.1, initial data for the evolution equations of Einstein’s equations must satisfy, first, the gauge constraint

[TABLE]

recall Eq. (2.9). Second, we must respect the Hamiltonian and Momentum constraints associated with Eq. (2.1) which take the form

[TABLE]

In this subsection, we use abstract indices $i,j,k,\ldots$ to represent two-dimensional purely spatial fields. Notice that in this subsection only, we use all the symbols $\delta$, $\phi$ etc., which we had introduced for fields on $S$ before, now to denote the restriction of these quantities to any $t=\mathrm{const}$-surface. The values of their time derivatives are denoted as $\dot{\delta}$, $\dot{\phi}$ etc. The quantity $\gamma_{ik}$ above is the induced $2$-metric (see Eq. (2.16)) and $D_{k}$ the corresponding covariant derivative. ${}^{(2)}R$ is the scalar curvature associated with $\gamma_{ik}$ and $\mathcal{K}_{ik}$ represents the extrinsic curvature with $\mathcal{K}=\mathcal{K}^{i}{}_{i}$. If $T_{ab}$ is the energy-momentum tensor of the matter source in the $2+1$-Einstein equations in Eqs. (2.1), then

[TABLE]

and hence

[TABLE]

Here the vector $j^{i}$ has been expressed in terms of the spatial frame $(\partial_{1}^{i},\partial_{2}^{i})=(m^{i},\overline{m}^{i})$; recall Eqs. (2.11) and (2.12). The corresponding spatial dual frame $(\omega^{1}_{i},\omega^{2}_{i})$ is defined as in Eq. (2.13). Notice that Eqs. (2.16) and (2.18) yield

[TABLE]

Because the shift $\beta$ vanishes as a consequence of Eq. (2.18), the extrinsic curvature $\mathcal{K}_{ik}$ is proportional to the time derivative of $\gamma_{ik}$ and is therefore determined by $\dot{\delta}$ and $\dot{\phi}$. The first step of finding a complete set of initial data on the $t=0$-surface is to find the quantities $\delta,\dot{\delta},\phi,\dot{\phi},\psi,\dot{\psi},\omega,\dot{\omega}$ as solutions of the Hamiltonian and Momentum constraints. Once this is done we find initial data for $f_{0},f_{1}$ as a solution of the gauge constraint in a second step. Recall that all these quantities are assumed to be real and only depend on $\theta$.

We shall construct solutions of the Hamiltonian and Momentum constraints using the York-Lichnerowicz conformal decomposition; see [2] and references therein. We shall not describe the general procedure here (which, in two dimensions, is slightly different from the standard $3$-dimensional case), but restrict to the simple time symmetric case

[TABLE]

in all of what follows. According to Eq. (2.33) and the choice of vanishing shift, this is the case if and only if

[TABLE]

The momentum constraint is therefore satisfied if $j^{1}=j^{2}=0$. According to Eq. (2.32), this is in particular the case if we pick

[TABLE]

The Momentum constraint is now satisfied and hence we attempt to solve the Hamiltonian constraint next. Because the topology of the spatial slices is $\mathbb{S}^{2}$, we can assume that $\gamma_{ik}$ is initially conformal to the standard round unit two-sphere metric. According to Eq. (2.33), we can therefore pick

[TABLE]

which yields

[TABLE]

where $\mathring{\gamma}_{ik}$ represents the metric for the unit round two-sphere

[TABLE]

The quantity $\delta$ can therefore be considered as the conformal factor in the standard conformal decomposition of the Hamiltonian constraint. We express the two-dimensional Ricci scalar in the Hamiltonian constraint as

[TABLE]

where $\mathring{R}=2$ is Ricci scalar of the unit round two-sphere. Replacing frame derivatives by eth-operators by means of Eq. (A.4), a straightforward calculation recasts the Hamiltonian to the form

[TABLE]

See Eq. (A.8) for our definition of the Laplace operator $\Delta_{\mathbb{S}^{2}}$ on the $2$-sphere. Using Eqs. (2.31) and (2.35), we get

[TABLE]

If we now pick

[TABLE]

the only remaining free function is $\dot{\psi}$ in terms of which the Hamiltonian constraint becomes

[TABLE]

Note that for $\dot{\psi}=0$ and $\delta=\text{const}$, the trivial (but not the only) solution of Eq. (2.39) is $\delta=1$ which yields the initial data of the Nariai metric. Hence, in order to obtain perturbations of the Nariai spacetime, we only have to provide non-zero functions $\dot{\psi}$ and solve numerically Eq. (2.39). In the whole paper, we choose

[TABLE]

where $\epsilon$ and $C$ are free real parameters and $\ell$ is any fixed positive integer. We list the necessary information about the functions $Y_{\ell}$ in Section A. Observe that $Y_{0}(\theta)={1}/(2\sqrt{\pi})$.

It only remains to provide initial data for $f_{0}$ and $f_{1}$ as solutions of the gauge constraint. Given the choices above, it turns out that the gauge constraint is satisfied if and only if

[TABLE]

Once all this is done, so, in particular, once we have solved Eq. (2.39) with Eq. (2.40), our initial data set is complete and satisfies all the required constraints: the Hamiltonian constraint, the Momentum constraint, and, the gauge constraint.

Numerical method to solve the Hamiltonian constraint

Now, we describe the basic idea for using a spectral implementation based on the spin-weighted spherical harmonics in the axi-symmetric case (no $\varphi$-dependence) for solving Eq. (2.39) with Eq. (2.40). We follow the approach in [17] for solving non-linear elliptic equations. For more information about these kind of methods, the interested reader is referred to [16, 36] and references therein.

Let us start by writing the right-hand side of Eq. (2.39) as a non-linear function $f(\delta,\eth\delta)$ with spin-weight [math]. The idea is then to construct a sequence of linearized problems whose solutions hopefully converge to the solution of the non-linear problem. For the Richardson’s iteration procedure this sequence of solutions $(\delta_{n})$ is constructed by solving

[TABLE]

for each $n=0,1,2,\ldots$ for some initial guess $\delta_{0}$ and then to set

[TABLE]

We call $\zeta$ the correction factor. The right-hand side of this equation is known as the residual $r_{n}$ at the step $n$, that measures how well $\delta_{n}$ satisfies the equation at the step $n$.

In our pseudo-spectral approach, we introduce suitable collocation points $\theta_{1},\ldots,\theta_{N}$, and impose Eq. (2.41) at those. Using the properties of the eth-operators listed in Eqs. (A), this yields an algebraic linear system of equations for $N$ spectral coefficients of $\zeta$ when written in the spin-weighted spherical harmonics basis. We shall not discuss the details of the solvability of this linear system here. However, it is guaranteed to have a unique solution in each step if the coefficients $(\partial f/\partial\delta)_{n}$ and $(\partial f/\partial\>\eth\delta)_{n}$ satisfy certain algebraic conditions in each step. If this is the case, then the iteration converges quickly and thereby allows us to construct accurate approximations of solutions of the nonlinear equation.

In Fig. 3, choosing $\epsilon=-10^{-4}$, $C=0$ and $\ell=2$, we show the behavior of the norm $\|r_{n}\|_{L^{2}(\mathbb{S}^{2})}$ as a function of $n$ which is numerically approximated by

[TABLE]

In this figure we observe that the norm of $r_{n}$ decays rapidly until it reaches a satisfactory order of $\sim 10^{-14}$.

Approximate analytic solutions of the Hamiltonian constraint

The heuristic analytical approach in Section 3.1 below relies on the following analytic approximations of solutions which is meaningful at least when the parameters $\epsilon$ and $C$ are small.

To this end, we shall now assume that the family of solutions $\delta$ of Eq. (2.39) with Eq. (2.40) depends smoothly on the parameters $\epsilon$ and $C$ in a neighborhood of $\epsilon=C=0$. Then we express $\delta$ approximately as

[TABLE]

for some so far unknown functions $\delta_{(1)},\ldots,\delta_{(5)}$ which are assumed to be independent of $C$ and $\epsilon$. With this ansatz, we find that Eqs. (2.39) and (2.40) are satisfied up to cubic order in the parameters, if

[TABLE]

which is implied by the linear orders in $\epsilon$ and $C$ and which yields that

[TABLE]

and, if

[TABLE]

where Eq. (2.44) and $Y_{0}={1}/(2\sqrt{\pi})$ have been used to simplify these equations. It is well known that for any PDE of the form

[TABLE]

defined on $\mathbb{S}^{2}$ given by any smooth source term function $f$ and any non-negative integer $p$, the uniquely determined solution is

[TABLE]

Regarding Eqs. (2.45), this implies that

[TABLE]

The coefficients $a_{\ell,k}$ here are defined implicitly by the equation

[TABLE]

which can therefore be calculated explicitly from the well-known Clebsch-Gordon coefficients [42]. When we combine all this with Eq. (2.43) we find

[TABLE]

This is an approximation of solutions of Eq. (2.39) with Eq. (2.40) which is expected to be valid for small values of the parameters $C$ and $\epsilon$.

In Fig. 4 we provide numerical evidence which supports the claim that Eq. (2.47) is a good approximation of solutions of the Hamiltonian constraint in many of the cases of interest.

3 Analysis and results

3.1 Heuristic mode analysis

Our interpretation of our numerical results and conclusions below are based on a heuristic mode analysis technique which we shall discuss first now. Recall that the unknowns of our dynamical equations are $\delta$, $\phi$, $\psi$, $\omega$, $f_{0}$ and $f_{1}$ which are all real quantities of spin-weight [math], $2$, [math], [math], [math] and $1$, respectively, depending only on $t$ and $\theta$. In order to facilitate the following analysis we define

[TABLE]

which are related to the physical $2+1$-metric $\hat{h}_{ab}$ as

[TABLE]

in the gauge Eq. (2.18); cf. Eqs. (2.4) and (2.5). Moreover, we define

[TABLE]

and set

[TABLE]

All the components of $u_{*}$ – and hence in particular the quantity $\delta_{*}$ which we shall mostly focus on in the following – can be decomposed using spin-weighted spherical harmonics of appropriate spin-weight (see Section A), for example

[TABLE]

Consistently with this, we write the collection of the $l$-th coefficients of all those components of $u_{*}$ for which these are defined schematically as $u_{*,l}$. We shall refer to this as the mode decomposition of $u_{*}$ and $\delta_{*}$, respectively. The $l=0$-mode $u_{*,0}$, so, in particular, $\delta_{*,0}$ will often be called fundamental mode. We shall also often write

[TABLE]

in order to emphasize that these are the relevant modes in the spatially homogeneous case. For the Nariai spacetime, we write $u_{*}=u_{*N}$ and $\delta_{*}=\delta_{*N}$ with

[TABLE]

Using this, the evolution equation for $\delta_{*}$ can be written schematically as

[TABLE]

We may rewrite this for each $l=0,1,\ldots$ as

[TABLE]

where

[TABLE]

We emphasize that Eq. (3.6) just an algebraic manipulation of Eq. (3.5). Also, it should be clear that similar decompositions can be performed for any of the other components of $u_{*}$. In any case, a lengthy calculation now reveals that

[TABLE]

for all $l=0,1,\ldots$, and hence that

[TABLE]

Now, suppose we are in a regime where $|F_{l}|$ is negligible in comparison to the other terms in Eq. (3.7) and that the dynamics is therefore dominated by the left-hand side. Then, this equation together with (3.4) implies

[TABLE]

where $l\geq 2$. Hence, in this regime, the $l=0$ mode is in general unstable (in fact, this is the heuristic explanation for the before-mentioned instability of the Nariai solution in the class of homogeneous spacetimes) while the $l\geq 2$-modes are all oscillatory. The $l=1$-mode is “somewhere in between”.

Before we continue, we wish to emphasize that the way the approximation Eqs. (3.8) – (3.10) is not a complete linearization of the evolution equations around the Nariai spacetime. It is therefore questionable whether Eqs. (3.8) – (3.10) are useful in any sense. In any case, our numerical experiments below show that the rather simplistic description above turns out to be sufficient as a basic for our main results.

Recall now that Eq. (2.47) is an approximation of the solution $\left.\delta\right|_{t=0}$ of the Hamiltonian constraint for our particular family of initial data. This approximation is expected to be valid for small parameter values $C$ and $\epsilon$. Let us now use Eq. (2.47) to express Eqs. (3.8) – (3.10) in terms of the initial data parameters $\ell$, $\epsilon$ and $C$. First, we see that Eqs. (3.1), (2.38) and (2.34) yield

[TABLE]

Eqs. (2.47) and (2.40) therefore give us the following result

[TABLE]

Now it turns out that in our applications, $C$ is typically much smaller than $\epsilon$. In fact, $C$ is often of the order $\epsilon^{2}$ (as justified below). When we combine Eqs. (3.8) – (3.10) with Eqs. (3.11) and (3.12) and only keep terms of order $\epsilon$, $\epsilon^{2}$ and $C$, we get

[TABLE]

for all $l=2,\ldots,2\ell$ provided $\ell\geq 1$. Here we use the notation

[TABLE]

The approximate description Eqs. (3.13) – (3.15) of the dynamics of the quantity $\delta_{*}$ can of course be expected to hold only for small values of the initial data parameters $\epsilon$ and $C$ (i.e., close to the exact Nariai spacetime given by $\epsilon=C=0$) and only for short times $t$ close to the initial time $t=0$. As long as this approximation holds, it suggests that the criticality of the cosmological models, i.e., the borderline between collapse and expansion globally in space, is mainly governed by the fundamental mode $\delta_{*,0}$ because all other modes are bounded if $\ell\geq 2$. Moreover, for any choice of $\epsilon\in\mathbb{R}$, the critical value of $C$, i.e., the value when $\delta_{*0}$ is exactly at the borderline in this approximation according to Eq. (3.13), should be close to

[TABLE]

where we use that $a_{\ell,0}=1/(2\sqrt{\pi})$ for all $\ell\geq 0$; recall the definition of $a_{\ell,l}$ by Eq. (2.46). In our applications, we typically set $\epsilon=-10^{-\kappa}$ for some positive integer $\kappa$, which therefore yields

[TABLE]

Later we provide numerical evidence that the actually critical value of $C$ for solutions of the fully nonlinear equations is indeed somewhat close to Eq. (3.16).

If we now choose any $\ell\geq 2$ and pick initial data parameters $\epsilon$ and $C$ close to the actual critical values of the fully nonlinear problem for which all modes are expected to be bounded, the oscillatory nature of the modes $\delta_{*,l}$ with $l=2,\ldots,2\ell$ suggested by Eq. (3.15) should dominate the dynamics. According to Eq. (3.15) the oscillation period is independent of $l$ (so long as $l$ is between $2$ and $2\ell$) and is given by

[TABLE]

The phase and the amplitudes of the oscillations however depend significantly on $l$. If $l=\ell$ the amplitude is proportional to $\epsilon$ (in leading order), while it is proportional to $\epsilon^{2}$ for $l\not=\ell$; moreover there is a phase shift of approximately $\pi/2$ between these two kinds of modes. In the next subsections, we present numerical evidence which supports all claims in this subsection.

3.2 Numerical evidence for Result 1: Existence of critical models

Based on the heuristic analysis in Section 3.1 we now present our numerical findings which support Result 1 in Section 1. Recall from our previous discussion that the critical behavior is mainly governed by the $l=0$-mode of $\delta_{*}$. The dynamics of this mode is approximated by Eq. (3.13) which suggests that, for initial data given by any fixed values $\epsilon$ and $\ell\geq 2$, the corresponding solution of the (fully nonlinear) evolution equation should be eventually expanding globally in space if $C<C_{crit}$, where $C_{crit}$ is given by Eq. (3.16), and be eventually collapsing globally in space if $C>C_{crit}$. The critical case should therefore be $C=C_{crit}$. Our numerical results now indeed confirm this, but with a slightly different value $C_{crit}$ than the value in Eq. (3.16). This suggests that there are nonlinear effects in the evolution equations, in particular effects of order $\epsilon^{2}$, which are not taken into account by our mode analysis (which was based on setting $F_{l}$ in Eq. (3.7) to zero). In any case, we find that the actual value is proportional to $\epsilon^{2}$ in leading order in consistency with Eq. (3.16); cf. Fig. 5.

In practice, we use the following algorithm to determine the actual value of $C_{crit}$ for any choice $\epsilon\in\mathbb{R}$ and $\ell\geq 2$ which is suggested by Eqs. (3.13) – (3.15):

Construct the full set of initial data as outlined in Section 2.3 for the given values of $\epsilon$ and $\ell$, and for the value $C$ given by Eq. (3.16). 2. 2.

Evolve the initial data to the future using the fully nonlinear evolution equations and gauges in Section 2.2. Determine whether the solution collapses (i) or expands (ii) globally in space. 3. 3.

Construct new ID in the same way as before for the same value of $\epsilon$ and $\ell$, but with some slightly decreased value of $C$ if (i) in Step 2, or, with some slightly increased value of $C$ if (ii); cf. Eq. (3.13). 4. 4.

Go back to Step 2 and repeat this process until a sufficiently good approximation of the critical solution has been obtained.

This algorithm is now used in Fig. 6 to approximate the actual fully nonlinear critical solution for $\epsilon=-10^{-4}$ and $\ell=2$. It demonstrates that the time period for which $\delta_{*h}$ is bounded (and oscillatory, see below) is longer the closer $C$ is to the critical value.

In Figs. 8 and 8 we apply the algorithm now to various values of $\epsilon$ and fixed $\ell$, only plotting our best numerical approximation of the critical solution obtained by our algorithm.

All these plots Figs. 6, 8 and 8 therefore confirm Result 1 in Section 1. In the next subsection, we shall study the oscillations. Before we get to this, let us emphasize that none of our numerical solutions is (as a matter of principle) exactly critical. Eventually during the evolution, the solutions “all make a decision whether to expand or to collapse”. Let us discuss how this happens in terms of the following alternative decomposition of the evolution equation of the $l=0$-mode

[TABLE]

The first term on the right-hand side captures all terms (also all nonlinear ones) in the equation which are present in the spatially homogeneous case in which we fully understand the criticality of the Nariai solution [3]. The second term can then be considered as “inhomogeneous corrections” to the equations close to the homogeneous case. Now, Fig. 9 suggests that the expected unstable behavior is triggered once the homogeneous term $G_{0}^{(H)}(t,u_{*,0},\dot{u}_{*,0})$ dominates the right-hand side of Eq. (3.19) and the evolution therefore displays the well-known Nariai-like instability. The solution plotted there is again our best numerical approximation of the critical solution in Fig. 6. We find that whatever the value of $\delta_{*,0}$ is at the time when the homogeneous term takes over determines whether the solution eventually expands or collapses globally in space.

3.3 Numerical evidence for Result 2: Nonlinear oscillatory dynamics

As discussed before, if the solution is critical (or almost critical) and hence $\delta_{*}$ is bounded for an extended period of time, the oscillatory nature of the modes $\delta_{*,l}$ for $l\geq 2$ suggested by Eq. (3.15) should dominate the dynamics. We shall discuss the dynamics of the modes $l\geq 2$ first now, before we explain the oscillatory behavior of the $l=0$-mode in Figs. 6, 8 and 8 (which is clearly not explained by Eq. (3.13)).

As before we assume that $\ell\geq 2$. Recall that the oscillation period of the modes $\delta_{*,l}$ for $l\geq 2$ is expected to be $T_{L}$ given by Eq. (3.18). In Table 1 we confirm the validity and accuracy of this heuristic prediction.

As expected the agreement with Eq. (3.15) is better for smaller $t$ and indeed gets worse for the second and third oscillation when non-linear effects clearly become significant.

Recall from Eq. (3.15) that there should be a phase difference of $\pi/2$ between the oscillations of the modes with $l\not=\ell$ ($l=2,\ldots,2\ell$), and the mode $l=\ell$. Moreover, the amplitudes of all the former modes should be proportional to $\epsilon^{2}$ while the amplitude of the latter is proportional to $\epsilon$. All this is confirmed in Figs. 11 and 11.

We also remark here that Fig. 8 confirms, in consistency with the heuristic predictions, that the oscillation period only depends on $\ell$ but not on $\epsilon$ or $C$.

Now, while our heuristic analysis explains that the fundamental mode is bounded for the (almost) critical solutions and all modes $l\geq 2$ are oscillatory, it misses the oscillatory behavior of the fundamental mode which is obvious in Figs. 6, 8 and 8. The basic assumption for the results in Section 3.1 was that the term $F_{l}$ in Eq. (3.7) is negligible. For the fundamental mode, this is clearly not the case and terms which are $O((\delta_{*,\ell})^{2}+(\partial_{t}\delta_{*,\ell})\delta_{*,\ell}+(\partial_{t}\delta_{*,\ell})^{2})$ are expected to change the dynamics significantly; recall that the amplitudes of all modes $\delta_{*,l}$ with $l\not=\ell$ are of order $\epsilon^{2}$ and hence of higher order than $\delta_{*,\ell}$. Because of the quadratic coupling of the fundamental mode and the $l=\ell$-mode, whose amplitude is proportional to $\epsilon$ and whose oscillation period length is approximately given by Eq. (3.18), we expect that the amplitude of the oscillation of the fundamental mode is proportional to $\epsilon^{2}$ and the oscillation period length is half of Eq. (3.18), i.e.,

[TABLE]

The statement about the amplitude is indeed confirmed by Fig. 8. The accuracy of the prediction about the oscillation period length is studied in Table 2.

In the rest of this subsection we consider the question whether these oscillations are a real physical effect of our models rather than just a gauge effect. To this end, we define a physical time function as the solution $\tau$ of the Eikonal equation

[TABLE]

with zero initial data on any of our models. As explained in [8], the value of $\tau$ represents the proper time along the congruence of unit timelike geodesics which start perpendicularly to the initial hypersurface. In Figs. 13 and 13, we plot this function $\tau$ in one of our cosmological models, mainly to demonstrate that our numerical solutions cover a significant part of physical time.

Let us now describe our oscillations in terms of the physical time before. In order to simplify the discussion a little, we exploit the fact that for $3+1$-Gowdy symmetric spacetimes and hence for axi-symmetric $2+1$-spacetimes, the poles of the spatial two-sphere is a geometrically distinguished point at all times. It is therefore geometrically (and physically) meaningful to look at the Kretschmann scalar of the $3+1$-metric as a function of $\tau$ at the pole $\theta=0$ of the spatial $2$-spheres only. This is done in Fig. 14 for various critical solutions. The oscillations are evident in this representation and hence are a real physical phenomenon

All this confirms Result 2 in Section 1.

3.4 Numerical evidence for Result 3: Late time behavior

We have now used numerical evidence to support our claim that for any $\epsilon\in\mathbb{R}$ and $\ell\geq 2$ it is possible to keep the quantity $\delta_{*}$ bounded for as long as we like by picking $C$ sufficiently close to some critical value. Fig. 15 now shows that also the quantity $\phi_{*}$ (see the definition in Eq. (3.1)) is bounded and, in fact, oscillatory. Since $\delta_{*}$ and $\phi_{*}$ determine the physical geometry of the spatial two-spheres (recall Eq. (3.2)), it follows that the spatial $\mathbb{S}^{2}$-factor of critical $3+1$-models do not deviate much from the geometry of the standard round unit $2$-sphere for an arbitrary long time, on the one hand.

On the other hand, the geometry of the spatial $\mathbb{S}^{1}$-factor, which is described by the quantity $\psi$, changes exponentially, as suggested by Fig. 17, and its growth is almost unaffected by whether $C$ is larger or smaller than the critical value. See also Fig. 17 which shows the deviation of this quantity from the corresponding Nariai values. Recall that $\psi_{*}$ is defined in Eq. (3.3).

This supports the claim that the long-term behavior of our inhomogeneous critical solutions is very similar to that of the exact Nariai solution. In particular, as for the Nariai solution, the highly anisotropic timelike future is expected to be inconsistent with the cosmic no hair picture. This supports Result 3 in Section 1. As explained earlier, we have convincing evidence now that whenever the solutions are non-critical, they either expand or collapse to the future eventually. In the expanding case, we expect the solutions to behave in accordance with the cosmic no-hair conjecture. This is indeed confirmed by our numerical results. For instance, Figs. 19 and 19 show (for a clearly non-critical solution; the value $C=-10^{-4}$ is far in the expanding regime) that, while the $3+1$-Kretschmann scalar starts off close to the Nariai value, it eventually approaches the expected de-Sitter value after all the oscillations have died out.

Acknowledgments

J. F. would like to thank the Department of Mathematics at the University of Oslo for hospitality. Part of this research was supported by the European Research Council through the FP7-IDEAS-ERC Starting Grant scheme, Project No. 278011 STUCCOFIELDS. The author L. E. was partly funded by the University of Otago Research Grant “Dynamical dark energy in the young universe and its consequences for the present and future history” in 2016. We would like to thank Dr Chris Stevens for helpful discussions.

Appendix A Spin-weighted spherical harmonics

Let $(\theta,\varphi)$ be the standard polar coordinates in $\mathbb{S}^{2}$. A function $f$ on $\mathbb{S}^{2}$ has spin-weight $s$ if it transforms under a local rotation by an angle $\tau$ in the tangent plane at every point $(\theta,\varphi)\in\mathbb{S}^{2}$ as $f\to e^{is\tau}f$. In this case, $f$ can be written as

[TABLE]

where ${}_{s}Y_{lm}(\theta,\varphi)$ are the spin-weighted spherical harmonics [37] and $a_{lm}$ are complex numbers. These functions are normalized as

[TABLE]

For any function $f$ of spin-weight $s$, this identity can be used to calculate the complex coefficients $a_{lm}$ in Eq. (A.1).

The eth operators $\eth$ and $\bar{\eth}$ are defined by

[TABLE]

cf. Eq. (2.11), for any function $f$ on $\mathbb{S}^{2}$ with spin-weight $s$. Using Eqs. (A.3), we can therefore express the frame vectors $(m^{a},\overline{m}^{a})$ in Eq. (2.11) in terms of the eth-operators as

[TABLE]

The properties of raising and lowering spin are

[TABLE]

In fact, we can use the relations (A) to define spin-weighted spherical harmonics ${}_{s}Y_{lm}$ with any integer spin-weight $s$ from the standard spherical harmonics

[TABLE]

It is easy to check that from any function $f$ with spin-weight $s$, we can obtain a function with either spin $s+1$ from $\eth(f)$ or spin $s-1$ from $\bar{\eth}(f)$. Thus, they are also known in the literature as the raising and lowering operators [32]. We can also check that

[TABLE]

The Laplace operator (in our sign convention) of the two-sphere can be written in terms of eth-operators as

[TABLE]

Further, using the commutation relation Eq. (A.7) we obtain the useful expressions

[TABLE]

Finally, any function $f$ with spin-weight $s$ can be expanded in terms of axi-symmetric spin-weighted spherical harmonics

[TABLE]

since the latter is independent of $\varphi$, i.e., Eq. (A.1) becomes

[TABLE]

for complex numbers $a_{l}$. In analogy to Eq. (A.6), we shall often write

[TABLE]