Einstein’s Universe: Cosmological structure formation in numerical relativity

To evolve a fully general relativistic cosmology we use the open-source Einstein Toolkit (Löffler et al., 2012), a collection of codes based on the Cactus framework (Goodale et al., 2003). Within this toolkit we use the ML_BSSN thorn (Brown et al., 2009a) for evolution of the spacetime variables using the BSSN formalism (Shibata and Nakamura, 1995; Baumgarte and Shapiro, 1999), and the GRHydro thorn for evolution of the hydrodynamics (Baiotti et al., 2005; Giacomazzo and Rezzolla, 2007; Mösta et al., 2014). In addition, we use our initial-condition thorn, FLRWSolver (Macpherson et al., 2017), to initialise linearly-perturbed FLRW spacetimes with perturbations of either single-mode or CMB-like distributions.

We assume a dust universe, implying pressure $P=0$, however GRHydro currently has no way to implement zero pressure for hydrodynamical evolution. Instead we set $P\ll\rho$, with a polytropic equation of state,

where $K_{\mathrm{poly}}$ is the polytropic constant, which we set $K_{\mathrm{poly}}=0.1$ in code units. We have found this to be sufficient to match the evolution of a homogeneous, isotropic, matter-dominated universe. Deviations from the exact solution for the scale factor evolution, at $80^{3}$ resolution, are within $10^{-6}$ (see Macpherson et al., 2017).

We perform a series of simulations with varying resolutions, $64^{3},128^{3}$, and $256^{3}$, and comoving physical domain sizes, $L=100$ Mpc, 500 Mpc, and 1 Gpc, to study different physical scales. We simulate all three domain sizes at $64^{3}$ and $128^{3}$ resolution, and only the $L=1$ Gpc domain size at $256^{3}$ resolution due to computational constraints. During the evolution we do not assume a cosmological background, and for convenience, since we have not yet implemented a cosmological constant in the Einstein Toolkit, we assume $\Lambda=0$.

Post-processing analysis is performed using the mescaline code, which we introduce and describe in Section V.2.

We choose the comoving length unit of our simulation domain to be 1 Mpc, implying a domain of $L=100$ in code units is equivalently $L=100$ Mpc. In geometric units $c=1$, and so we can relate our length unit, $l=1$ Mpc, and our time unit, $t_{c}$, via the speed of light (in physical units)

To find our background FLRW density we use $H(z=0)=H_{0}$, with units of s^-1. This implies

where $H_{0,\mathrm{code}}$ and $H_{0,\mathrm{phys}}=100\,h\,\mathrm{km\,s^{-1}\,Mpc^{-1}}$ are the Hubble parameter expressed in code units and physical units, respectively. We use (2) together with (3) and the Friedmann equation for a flat, matter-dominated model

where an overdot represents a derivative with respect to proper time, $\bar{\rho}$ is the homogeneous density, and $a$ is the FLRW scale factor. We find the background FLRW density, evaluated at $z=0$, in code units, to be

For computational reasons we adopt the initial FLRW scale factor $a_{\mathrm{init}}=a(z=1100)=1$, whilst the usual convention in cosmology is to set $a_{0}=a(z=0)=1$. The density (5) was calculated using the Hubble parameter $H_{0,\mathrm{phys}}$ evaluated with $a_{0}=1$. The comoving (constant) FLRW density is $\rho^{*}=\bar{\rho}\,a^{3}=\bar{\rho}_{0}\,a_{0}^{3}$, and so (5) is the comoving density $\rho^{*}$. We choose $h=0.704$, and our choice $a_{\mathrm{init}}=1$ implies our initial background density is the comoving FLRW density.

Simulations are initiated at $z=1100$ and evolve to $z=0$. We quote redshifts computed from the value of the FLRW scale factor at a particular conformal time,

where $z_{\mathrm{cmb}}=1100$. Since we set $a_{\mathrm{init}}=1$, we have $a_{0}=1101$. The evolution of the FLRW scale factor in conformal time is

where $\xi$ is the scaled conformal time defined in Section III.1. Importantly, the redshifts presented throughout this paper are indicative only of the amount of coordinate time that has passed, and are not necessarily indicative of redshifts measured by observers in an inhomogeneous universe.

We solve the linearly-perturbed Einstein equations to generate our initial conditions. Assuming only scalar perturbations, the linearly-perturbed FLRW metric in the longitudinal gauge is

In this gauge the metric perturbations $\phi$ and $\psi$ are the Bardeen potentials (Bardeen, 1980). These are related to perturbations in the matter distribution via the linearly perturbed Einstein equations

where an over-bar represents a background quantity, and $\delta X$ represents a small perturbation in the quantity $X$, with $\delta X\ll X$. A matter-dominated (dust) universe has stress-energy tensor

where $\rho$ is the rest-mass density, $u^{\mu}=dx^{\mu}/d\tau$ is the four-velocity of the fluid, and $\tau$ is the proper time. Assuming small perturbations to the matter we have

where the fractional density perturbation is $\delta\equiv\delta\rho/\bar{\rho}$, and $v^{i}=dx^{i}/d\eta$ is the three-velocity.

Solutions to (9) are found by taking the time-time, time-space, trace and trace-free components, given by

respectively, where we have assumed all perturbations are small such that second-order (and higher) terms can be neglected. Here, $\partial_{i}\equiv\partial/\partial x^{i}$, $\nabla^{2}=\partial^{i}\partial_{i}$, $\partial_{\langle i}\partial_{j\rangle}\equiv\partial_{i}\partial_{j}-1/3\,\delta_{ij}\nabla^{2}$, a ^′ represents a derivative with respect to conformal time, and $\mathcal{H}\equiv a^{\prime}/a$ is the conformal Hubble parameter. Solving these equations, we find

where $f,g$ are arbitrary functions of spatial position, we introduce the scaled conformal time coordinate

and we have defined

Equations (14) contain both a growing and decaying mode for the density and velocity perturbations. We choose $g=0$ to extract only the growing mode of the density perturbation, and our solutions become

implying $\phi^{\prime}=0$ in the linear regime.

We use (14) along with the Code for Anisotropies in the Microwave Background (CAMB) (Lewis and Bridle, 2002) to generate the matter power spectrum at $z=1100$, with parameters consistent with Planck Collaboration et al. (2016) as input. Figure 1 shows the matter power spectrum from CAMB (grey curve), as a function of wavenumber $|\textbf{k}|=\sqrt{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$. We use the Python module c2raytools ¹¹1https://github.com/hjens/c2raytools to generate a 3-dimensional Gaussian random field drawn from the CAMB power spectrum. This provides the initial density perturbation. The magenta curve in Figure 1 shows the region of the matter power spectrum sampled in our highest resolution ($256^{3}$), largest domain size ($L=1$ Gpc) simulation. The smallest k component sampled represents the largest wavelength of perturbations — approximately the length of the box, $L$ — and the largest k component sampled represents the smallest wavelength of perturbations — two grid points. To relate the initial density perturbation to the corresponding velocity and metric perturbations, we transform (14) into Fourier space. Initially, $\xi=1$ which gives a density perturbation of the form

where $\textbf{k}=(k_{x},k_{y},k_{z})$, so we can define an arbitrary function $\delta(\textbf{k})$, and construct the metric perturbation and velocity, respectively, using

where $i^{2}=-1$. With the Fourier transform of the Gaussian random field as $\delta(\textbf{k})$, we calculate the velocity and metric perturbations in Fourier space using (19), and then use an inverse Fourier transform to convert the perturbations to real space. The density perturbation $\delta$ is already dimensionless, and we normalise by the speed of light, $c$, to convert $v^{i}$ and $\phi$ to code units. Figure 2 shows initial conditions at $256^{3}$ resolution for box sizes $L=1$ Gpc, 500 Mpc, and 100 Mpc in the left to right columns, respectively. The top row shows the density perturbation $\delta$, the middle row shows the normalised metric perturbation $\phi/c^{2}$, and the bottom row shows the magnitude of the velocity perturbation normalised to the speed of light $|v|/c$. These initial conditions are sufficient to describe a linearly-perturbed FLRW spacetime in FLRWSolver.

We assume a flat FLRW cosmology for the initial instance only. Simulations begin with small perturbations at the CMB, and so the assumption of a linearly-perturbed FLRW spacetime is sufficiently accurate.

The dimensionless cosmological parameters describe the content of the Universe. From (28) we define

giving the Hamiltonian constraint in the form

We require this to be satisfied at all times. Here, $\Omega_{m}$ is the matter energy density, $\Omega_{R}$ is the curvature energy density, $\Omega_{Q}+\Omega_{L}$ is the energy density associated with the backreaction terms; a purely general relativistic effect. For a standard $\Lambda$CDM cosmology, these cosmological parameters are $\Omega_{m}=0.308\pm 0.012$, $|\Omega_{R}|=|\Omega_{k}|<0.005$, $\Omega_{Q}=0$, and $\Omega_{L}=0$ (Planck Collaboration et al., 2016).

The Universe is measured to be homogeneous and isotropic on scales larger than $\sim 80-100h^{-1}\mathrm{Mpc}$ Scrimgeour et al. (2012). Above these scales it is unclear whether the evolution of the average of our inhomogeneous Universe coincides with the FLRW (or $\Lambda$CDM) equivalent. In attempt to address this, we calculate averages over our entire simulation domain, but also over subdomains within the simulation to sample a variety of physical scales. We measure averages over spheres of varying radius $r_{\mathcal{D}}$ embedded in the total volume, from which we calculate the dimensionless cosmological parameters (32), the Hubble parameter (26), and consequently the effective matter expansion $a_{\mathcal{D}}$.

The spatial Ricci tensor $R_{ij}$ is the contraction of the Riemann tensor. We calculate this directly from the metric using

where the spatial connection coefficients are

We use our analysis code mescaline, written to analyse three-dimensional HDF5 data output from our simulations. The code reads in the spatial metric $\gamma_{ij}$, the lapse $\alpha$, the extrinsic curvature $K_{ij}$, the density $\rho$, and the velocity $v^{i}$ from the Einstein Toolkit three-dimensional output. From these quantities we calculate the spatial Ricci tensor $R_{ij}$ from the spatial metric, and hence the Ricci scalar via $R=\gamma^{ij}R_{ij}$. We take the trace of the extrinsic curvature $K=\gamma^{ij}K_{ij}$ and with the set of equations defined in Appendix A we calculate averages and the resulting backreaction terms. We also use mescaline to calculate the Hamiltonian and momentum constraint violation, discussed in Appendix C. We compute derivatives using centred finite difference operators, giving second order accuracy in both space and time, the same order as the Einstein Toolkit’s spatial discretisation.

Figure 6 shows the global evolution of the effective scale factor, $a_{\mathcal{D}}$. The blue curve shows $a_{\mathcal{D}}$ calculated over the whole $L=1$ Gpc, $256^{3}$ resolution domain with (31). The purple dashed curve in the top panel shows the corresponding FLRW solution for the scale factor, $a_{\mathrm{FLRW}}$, found by solving the Hamiltonian constraint for a flat, matter-dominated, homogeneous, isotropic Universe in the longitudinal gauge,

giving the solution (7). The bottom panel of Figure 6 shows the residual error between the two solutions, which remains below $10^{-3}$ for the evolution to $z=0$.

Analysing the cosmological parameters as an average over the entire simulation domain we find agreement with the corresponding FLRW model in our chosen gauge. Globally, at $z=0$, we find $\Omega_{m}\approx 1$, $\Omega_{R}\approx 10^{-8}$, and $\Omega_{Q}+\Omega_{L}\approx 10^{-9}$. Systematic errors on these values are discussed in Appendix D.

We find global cosmological parameters consistent with a matter-dominated, flat, homogeneous, isotropic universe, and therefore no global backreaction. The evolution of the effective scale factor $a_{\mathcal{D}}$, evaluated over the whole domain, coincides with the corresponding FLRW model, as shown in Figure 6. The $<10^{-3}$ discrepancy between the two solutions does not correlate with the onset of nonlinear structure formation, indicating that this difference is most likely computational error.

We find a globally flat geometry in our simulations with $\Omega_{R}\approx 10^{-8}$. This could be a result of our treatment of the matter as a fluid. We cannot create virialised objects and so any “clusters” will continue to collapse towards infinite density. In reality, a dark matter halo or galaxy cluster would form, be supported by velocity dispersion, and stop collapsing. The surrounding voids would continue to expand and potentially contribute to a globally negative curvature (see e.g. Bolejko, 2017, 2018). Without a particle description for dark matter alongside numerical relativity we cannot properly capture this effect.

Any contribution from backreaction, $\mathcal{Q_{D}}$ or $\mathcal{L_{D}}$, is due to variance in the expansion rate and shear. The left panel of Figure 5 shows the matter expansion rate $\theta$, where collapsing regions (yellow, orange, and red) balance the expanding regions (green) due to our treatment of matter. While we see spatial variance in $\theta$, there is no global contribution from backreaction under our assumptions. Therefore, in our chosen gauge and under the caveats described in Section VII.3 below, backreaction from structure formation is unlikely to explain dark energy.

We find strong positive curvature on scales below the homogeneity scale of the Universe. Variations in measured cosmological parameters are up to 31% based purely on location in an inhomogeneous matter distribution. Our result is similar to that of Bolejko (2017) on small scales, but with larger variance in $\Omega_{R}$ because of increased small-scale density fluctuations due to our fluid treatment of dark matter.

On the approximate homogeneity scale of the Universe we find mean cosmological parameters consistent with the corresponding FLRW model to $\sim 1\%$. Aside from this, we find the parameters can deviate from these mean values by 4-9% depending on physical location in the simulation domain. This implies that, although on average these coincide with a flat, homogeneous, isotropic Universe, an observers interpretation may differ by up to 9% based purely on her position in space.

As we approach larger averaging radii within a 1 Gpc³ volume, we begin to move away from independent spheres, and each sphere begins to overlap with others; effectively sampling the same volume. Due to this, the confidence intervals contract, and eventually at $r_{\mathcal{D}}\approx 400$ Mpc most spheres become indistinguishable from the mean. The beginning of this is evident in Figure 8 as we approach $r_{\mathcal{D}}=250$ Mpc. This transition appears to be due to overlapping spheres, although could in part be due to the statistical homogeneity of the matter distribution at these scales.

Local observations of type 1a supernovae generally probe scales of $75-450\,h^{-1}$Mpc (Wu and Huterer, 2017). Nearby objects are excluded from the data in an effort to minimise cosmic variance on the result (Riess et al., 2016, 2018b, 2018a). In this work, we cannot meaningfully sample scales above 250 Mpc because our maximum domain size is only 1 Gpc³. In order to sample all scales used in nearby SNe surveys, we would need a domain size of $L\gtrsim 10\,h^{-1}$Gpc, with a resolution up to $1024^{3}$. Current computational constraints, and the overhead of numerical relativity, currently restrict us to domain sizes and resolutions used in this work. To address scales as similar as possible to those used in local surveys, we consider $75<r_{\mathcal{D}}<180\,h^{-1}$Mpc. On these scales we find $\Omega_{m}=1.002\pm 0.06$, $\Omega_{R}=0.002\pm 0.04$, and $\Omega_{Q}+\Omega_{L}=0.001\pm 0.02$, where variances are the 68% confidence intervals due to local inhomogeneity. This implies based on an observers physical location, measured deviations from homogeneity on these scales could be up to 6%. We expect this variance to decrease when including the full range of observations; including radii up to $450\,h^{-1}$Mpc. We investigate this further, including extrapolation to larger scales, in our companion paper Macpherson et al. (2018).

While the global effective scale factor demonstrates pure FLRW evolution, we find inhomogeneous expansion within spheres of 100 Mpc radius. Figure 9 shows the expansion rate differs by $2-8\%$ depending on the relative density of the region sampled. These differences agree with linear perturbation theory, to within 1%, on $\gtrsim 80$ Mpc scales, with smaller scales showing differences of up to 10%. These differences are most likely due to the nonlinearity of the density field on these scales, although, in addition, could involve general relativistic corrections. To properly test this we would require an equivalent Newtonian cosmological simulation to compare this relation at nonlinear scales, which we leave to future work.

1. 1.

   Our treatment of dark matter as a fluid is the main limitation of this work. Under this assumption, we are unable to form bound structures supported from collapse by velocity dispersions. In cosmological N-body simulations, particle methods are adopted so as to capture the formation of galaxy haloes, and local groups of galaxies as bound structures. Adopting a fully general relativistic framework in addition to particle methods would allow us to adopt a proper treatment of dark matter in parallel with inhomogeneous expansion.

2. 2.

   We take averages over purely spatial volumes. In reality, an observer would measure her past light cone, and hence the evolving Universe. Our results can thus be considered an upper limit on the variance due to inhomogeneities, since any structures located in the past light cone will be more smoothed out.

3. 3.

   Our results are explicitly dependent on the chosen averaging hypersurface. The result of averaging across different hypersurfaces has been investigated (Adamek et al., 2017; Giblin et al., 2018), and the results can show significant differences. It is clear the physical choice of hypersurface can be important for quantifying the backreaction effect.

4. 4.

   We assume $\Lambda=0$, and begin our simulations assuming a flat, matter dominated background cosmology with small perturbations. Throughout the evolution, on a global scale, we find the average $\Omega_{m}\approx 1$; consistent with this model. It is extremely well constrained that our Universe is best described by a matter content $\Omega_{m}\approx 0.3$ (e.g. DES Collaboration et al., 2017; Bonvin et al., 2017; Planck Collaboration et al., 2016; Bennett et al., 2013). The growth rate of cosmological structures in our simulations will therefore be amplified relative to the $\Lambda$CDM Universe.

5. 5.

   Given our limited spatial resolution, we underestimate the amount of structure compared to the real Universe. In addition, we resolve structures down to scales of two grid points, which means these structures may be under resolved.