Periodic boundary conditions and G₂ cosmology

As noted above, in standard cosmology the spatial curvature is assumed to be constant and zero (or at least perturbatively small) ³ . However, there is no completely independent constraint that guarantees a value for the magnitude of the effective normalized spatial curvature $\Omega_k$ of less than about 1% . Moreover, any possible such small non-zero measurement of $\Omega_k$ can perhaps be taken to indicate that the assumptions in the standard model are not valid. We also note that spatial curvature is, in general, evolving in relativistic cosmological models [8].

In order to develop models to be utilized for cosmological predictions and comparison with observational data, it is necessary to make assumptions. However, it is crucial to check how the assumptions affect the results. More importantly, we can only validate the consistency of assumptions and we cannot rule out alternative explanations. The assumption of a FLRW background on cosmological scales presents a number of problems [1]. For example, the assumption that there exists a 1+3 spacetime split and a global time and a background inertial (Gaussian normal) coordinate system over a complete Hubble scale background patch in the standard model leads to the constraint that the spatial curvature (and the vorticity) must be sufficiently small. In addition, any assumptions of exact periodic boundary conditions (PBCs) (appropriate on scales comparable to the homogeneity scale) necessarily imply that the global spatial curvature is exactly zero [1] (since in the actual standard model the Universe is taken to be simply connected), and thus any suitable approximation will amount to $\Omega_k$ being less than the perturbation (e.g. LPT) scale. We shall return to this point below.

There are also assumptions related to the weak field approach, the actual applicability of perturbation theory and Gaussian initial conditions, that imply spatial curvature is neglected. It is also asserted that cosmological backreaction can be neglected, but in LPT the initial fluctuations are usually taken to be Gaussian, which implies that at the linear level all averages are zero by construction. In addition, we should be aware that any intuition based on Newtonian theory may be misplaced; e.g. the average spatial curvature of voids and clustered matter (with negative and positive curvature, respectively) is not necessarily zero within GR. Thus, within standard cosmology the spatial curvature is assumed to be zero, or at least very small and at most first order in terms of the perturbation approximation, in order for any subsequent analysis to be valid. Any prediction larger than this indicates an inconsistency in the approach. In this sense the standard model cannot be used to predict a small spatial curvature.

It has recently been re-emphasized that, in an era of precision cosmology, any physical approximations (including, for example, the use of LPT), in standard cosmological models should be analyzed to check if they are sufficiently accurate for comparison with observations used to draw conclusions about the actual properties of our Universe [9].

In the general inhomogeneous cosmological case, the formulation of the evolution equations consists of utilizing global Cartesian coordinates and metric variables, where initial conditions on the metric are specified. In the actual N-body simulations a cell of length L (a fraction of the Hubble scale) is taken and the system is integrated where the PBC at L is applied at all times (and is used at every step to compute spatial derivatives near the edges of the cell). Hence spatial curvature is held zero on the boundary, and in this sense the growth of spatial curvature is suppressed, relatively speaking. It is of importance to verify numerically that initial PBC constrain the spatial curvature to be zero at the boundary for all time, which is maintained at zero even below numerical errors that begin to grow.

In addition, if in N-body simulations the cell is restricted to a finite size L, then periodicity implies that longer wavelength inhomogeneities than L are effectively neglected. This will have an effect on the numerical simulations, which is a general concern, especially when the large-scale cutoff is indeed below the Hubble scale as will be necessary if small (${\sim}Mpc$) scales are to be resolved. (There will be a cutoff in practice in the power spectrum of the inhomogeneities immediately from the initial conditions.)

Finally, it is not absolutely clear whether periodic initial conditions guarantee periodic evolution for all times. There is no general mathematical theorem on the preservation of discrete symmetries (that periodicity implies; however, see [10]). But note that if it were true that periodic initial conditions are propagated, then it follows that the spatial curvature is indeed constrained at L for all times.

In the few current relativistic cosmological simulations, the spatial curvature can have large local fluctuations but remains small when averaged over large scales. It becomes extremely small, $\Omega_R \sim 10^{-8}$, when averaged over the whole periodic box of the simulation [11]. Work is in progress for a more in-depth investigation of average spatial curvature from numerical relativity, along with the constraints due to the PBC, in a more frame-independent formalism [12]. However, the preliminary evidence is that these numerical solutions have an exceptionally and unnaturally small spatial curvature in a neighborhood of the periodic boundary, which affects the average spatial curvature in the whole cell, and can perhaps be ascribed to the use of PBC.

In different formulations of the EFE (e.g. the G₂ case below), it is difficult to compare with the standard numerical results of the full 3D code and the usual N-body simulations. Indeed, since the EFE are hyperbolic, and depending on the particular formulation of the evolution equations, once the initial conditions are fully specified it is not clear whether the application of spatial PBC are even required (beyond the specification of the initial data that is).

The simplest spatially inhomogeneous cosmological models have two commuting Killing vector fields (i.e. models admitting a 2-parameter Abelian isometry group acting transitively on spacelike 2-surfaces), which consequently have a single degree of spatial inhomogeneity. This class of models, which are referred to briefly as G₂ cosmologies, are governed by the EFE which constitute partial differential equations (PDE) in two independent variables [14]. A lot of the research to date has focused on the so-called Gowdy vacuum spacetimes [15], because such models are much more tractable than non-vacuum ones and the existence of compact space sections facilitates numerical simulations since the problem of boundary conditions at spatial infinity is avoided. Recently, a plane-parallel model (or wall Universe) [16], in which light rays pass through a series of plane-symmetric perturbations around a FLRW background, was used to study the back-reaction from the small-scale inhomogeneities [17, 18].

Belinskii, Khalatnikov and Lifshitz (BKL) [19] have conjectured that within GR the approach to the generic spacelike singularity to the past is vacuum dominated (assuming $p \lt \rho$), local, and oscillatory (labeled 'Mixmaster'). Locality implies that the contribution of terms in the evolution equations with spatial derivatives can be neglected. Numerical studies of the asymptotics of the simplest inhomogeneous vacuum Gowdy models demonstrate that on approach to the singularity, spiky structures form [20]. These spikes become narrower as the singularity is approached. Studies of G₂ and more general models have produced numerical evidence that the BKL conjecture generally holds except possibly at isolated points where spiky structures form [21, 22].

We are particularly interested in late time evolution here. We shall concentrate on G₂ cosmologies with a perfect fluid matter content, and especially dust [23]. At any instant of time, the state of a G₂ cosmology is described by a finite-dimensional dynamical state vector of functions of the spatial coordinate. That is, the dynamical state space of G₂ cosmologies is a function space and is thus infinite-dimensional. The evolution of a G₂ cosmology is consequently described by an orbit in this infinite-dimensional dynamical state space. In a G₂ cosmology, at each point there exists a preferred timelike 2-space that is orthogonal to the G₂-orbits. Therefore, there is an infinite family of geometrically preferred timelike congruences and gauge freedom in the choice of the orthonormal frame (and associated with the choice of the local coordinates). Following [14, 24], we shall utilize the orthonormal frame formalism adapted to the G₂ orbits. Such a dynamical formulation has proved effective in studying G₂ cosmologies. Here we define scale-invariant dependent variables by normalization with the area expansion rate of the G₂-orbits, in order to obtain the evolution equations as a system of PDE in first-order symmetric hyperbolic (FOSH) format. This FOSH format also provides a natural framework for numerical studies.

The geometry of the general G₂ class of spacetimes is given in [22], where all of the metric functions depend only on the time coordinate t and the spatial coordinate x. We will display the G₂ evolution system in terms of the timelike area gauge in the Gowdy subcase $\Sigma_2 = 0$ (with unit lapse, thereby defining t time). We utilize β-normalized variables [14, 23] in the orthonormal frame formalism [24], where the t and x subscripts denote partial differentiation. We also assume dust (γ = 1) and a single non-zero tilt component, v. We are particularly interested in future dynamics close to 'FL'. We also assume that the cosmological constant is zero, so that we have an unconstrained evolution system in FOSH format [14]:

$$\begin{align} \partial_{t}E_{1}{}^{1} & = \left(q+3\Sigma_{+}\right)\,E_{1}{}^{1} \end{align} \tag{ 1 }$$

$$\begin{align} \partial_{t}\Sigma_{-} + E_{1}{}^{1}\,\partial_{x}N_{\times} & = \left(q+3\Sigma_{+}-2\right)\,\Sigma_{-} + 2\sqrt{3}\,\Sigma_{\times}^{2} - 2\sqrt{3}\,N_{-}^{2} \end{align} \tag{ 2 }$$

$$\begin{align} \partial_{t}N_{\times} + E_{1}{}^{1}\,\partial_{x}\Sigma_{-} & = \left(q+3\Sigma_{+}\right)\,N_{\times} \end{align} \tag{ 3 }$$

$$\begin{align} \partial_{t}\Sigma_{\times} - E_{1}{}^{1}\,\partial_{x}N_{-} & = \left(q+3\Sigma_{+}-2-2\sqrt{3}\Sigma_{-}\right)\,\Sigma_{\times} - 2\sqrt{3}\,N_{\times}\,N_{-} \end{align} \tag{ 4 }$$

$$\begin{align} \partial_{t}N_{-} - E_{1}{}^{1}\,\partial_{x}\Sigma_{\times} & = \left(q+3\Sigma_{+}+2\sqrt{3}\Sigma_{-}\right)\,N_{-} + 2\sqrt{3}\,\Sigma_{\times}\,N_{\times} \end{align} \tag{ 5 }$$

$$\begin{equation} \left(\partial_{t} + v\,E_{1}{}^{1}\,\partial_{x}\right)\,\Omega + \Omega E_{1}{}^{1}\,\partial_{x}v = 2\,\Omega \ \left[\left(q+1\right) - \frac{1}{2}\,\left(1-3\Sigma_{+}\right)\,\left(1+v^{2}\right) - 1\right] \end{equation} \tag{ 6 }$$

$$\begin{equation} \left(\partial_{t} + v\,E_{1}{}^{1}\,\partial_{x}\right)\,v = \left(1-v^{2}\right)\, \left[\ 3\left(N_{\times} \Sigma_{-} - N_{-} \Sigma_{\times}\right) + \frac{3}{2}\Omega v - \left(1-3\Sigma_{+}\right)\,v\right] \end{equation} \tag{ 7 }$$

where

$$\begin{equation} \left(q+3\Sigma_{+}\right) \equiv 2 - \frac{3}{2}\,\left(1-v^{2}\right)\,\Omega \ \end{equation} \tag{ 8 }$$

and

$$\begin{equation} \Sigma_{+} \equiv \frac{1}{2}\,\left(1-\Sigma_{-}^{2}-N_{\times}^{2}-\Sigma_{\times}^{2} -N_{-}^{2}-\Omega\right). \\ \end{equation} \tag{ 9 }$$

Note that in the present case we have from equations (8) and (9) that

$$\begin{equation} q = \frac12 + \frac32\left( \Sigma_{-}^2 + N_{\times}^2 + \Sigma_{\times}^2 + N_{-}^2\right) + \frac32v^2\Omega \unicode{x2A7E} \frac12, \end{equation} \tag{ 10 }$$

which, together with

$$\begin{equation} \partial_t \beta = -\left(q+1\right) \beta, \end{equation} \tag{ 11 }$$

guarantees that β is single signed and is monotone (for small ε—see below). The (negative) spatial curvature ($-\Omega_k$) is defined via

$$\begin{equation} \Omega_k \equiv N_{\times}^{2} + N_{-}^{2}. \end{equation} \tag{ 12 }$$

The relation of the variables above with metric components is given by (see (6.39) and (6.40) of [23], or appendix A.3 of [14]):

$$\begin{equation} ds^2 = \beta^{-2} \left[ -dt^2 + \left(E_1{}^1\right)^{-2} dx^2 \right] + e^{2t} \left[ e^{P\left(t,x\right)} \left(dy + Q\left(t,x\right) dz\right)^2 + e^{-P\left(t,x\right)} dz^2 \right] \end{equation} \tag{ 13 }$$

$$\begin{align} \left( \Sigma_{-} ,\ N_{\times} \right) & = \frac{1}{2\sqrt{3}} \left( \partial_t,\ - E_1{}^1 \partial_x \right) P \end{align} \tag{ 14 }$$

$$\begin{align} \left( \Sigma_{\times} ,\ N_{-} \right) & = \frac{1}{2\sqrt{3}} e^P \left( \partial_t, E_1{}^1 \partial_x \right) Q.\end{align} \tag{ 15 }$$

2.2.1. Initial data.

2.2.2. Linear regime.

For small ε, we integrate equation (1) to obtain:

$$\begin{equation} {\tilde E} = e^{\frac{t}{2}} + o\left(\epsilon^2\right), \end{equation} \tag{ 16 }$$

where we have normalized ${\tilde E}$ to be unity at t = 0. In the linear regime we have that $e^{\frac{t}{2}}$ is small (which for $\epsilon \sim 10^{-4}$ corresponds to t < 10). Equations (6) and (7) also constitute o(ε²) corrections to zeroth order evolution equations. In the linear regime, $\Omega_k \sim o(\epsilon^2)$.

2.2.3. Shear and curvature.

Equations (2)–(5) represent the first order evolution equations for the normalized shear and curvature variables with second order corrections. We immediately see the growing modes for the normalized curvature variables ${\sim}e^{\frac{t}{2}}$ and the decaying modes for the normalized shear variables ${\sim}e^{-\frac{3t}{2}}$, corresponding to the familiar eigenvalues $\{\frac{1}{2}, -\frac{3}{2} \}$ for the flat FL 'saddle point' solution [so that the growth of the shear is suppressed relative to that of the spatial curvature in the initial linear regime].

Solving equations (2)–(5) we obtain:

$$\begin{align} {\tilde \Sigma}_- & = e^{-\frac{3t}{2}}\left[\sigma_- + \epsilon{\bar \Sigma}_-\right], \end{align} \tag{ 17 }$$

$$\begin{align} {\tilde N}_\times & = e^{\frac{t}{2}}\left[\nu_\times + \epsilon{\bar N}_\times\right], \end{align} \tag{ 18 }$$

$$\begin{align} {\tilde \Sigma}_\times & = e^{-\frac{3t}{2}}\left[\sigma_\times + \epsilon{\bar \Sigma}_\times\right], \end{align} \tag{ 19 }$$

$$\begin{align} {\tilde N}_- & = e^{\frac{t}{2}}\left[\nu_- + \epsilon{\bar N}_-\right], \end{align} \tag{ 20 }$$

where $\{\sigma_-, \sigma_\times, \nu_-, \nu_\times \}$ are slowly varying (and, for example, $\partial_{x}\Sigma_{-} \sim (\epsilon e^{\frac{t}{2}}) {\bar \Sigma}_-{}^{^{\prime}}$, $\partial_{x}\Sigma_{\times} \sim (\epsilon e^{\frac{t}{2}}) {\bar \Sigma}_\times{}^{^{\prime}}$, where a prime denotes $\partial_{x}$). In the regime in which the shear is sub-dominant, these quantities are approximately constant.

2.2.4. Spatial curvature.

From equations (3) and (5), to $o(\epsilon)$ equation (12) becomes:

$$\begin{equation} \partial_{t}\Omega_k = \Omega_k, \end{equation} \tag{ 21 }$$

so that $\Omega_k$ remains small (and of second order initially). More precisely, to $o(\epsilon^4)$:

$$\begin{align} \partial_{t}\Omega_k & = \left(\epsilon e^{\frac{t}{2}}\right)^2[\left(\nu_-^{2} + \nu_\times ^{2}\right) + \epsilon \left(\nu_\times {\bar{N}_\times} + \nu_- {\bar{N}_-}\right) +4 \sqrt{3}\left(\epsilon e^{\frac{t}{2}}\right)e^{-2t} \left(\sigma_- \nu_- + \sigma_\times \nu_{\times}\right)\nonumber\\ &\quad + \left(\epsilon e^{\frac{t}{2}}\right)e^{-2t} \left({\nu_-}{\bar \Sigma}_\times{}^{^{\prime}} + \nu_\times{\bar \Sigma}_-{}^{^{\prime}}\right). \end{align} \tag{ 22 }$$

Clearly to second order we duplicate equation (21); the leading order correction to this equation comes from the term $\epsilon (\nu_- {\bar{N}_\times} + \nu_\times {\bar{N}_-})$, which corresponds to equation (22) to next order. The remaining terms are of order $o(\epsilon e^{-2t})$ and $o(\epsilon^{2})$. For early times $\partial_{t}\Omega_k \gt 0$, and so the magnitude of the spatial curvature grows. The sign of the next order terms is not necessarily positive.

2.2.5. Boundary conditions.

2.2.6. Numerical methods: zooming.