1 Introduction

FLRW spacetime is a homogeneous, isotropic, expanding solution to Einstein’s equations. This solution is the basis for the current standard cosmological model; \(\Lambda \)CDM. The spatially-flat FLRW spacetime sourced by pure-dust with no cosmological constant is the Einstein-de Sitter (EdS) spacetime.

Here we provide a thorn to give initial conditions for cosmology using the Einstein Toolkit (ET). Currently, FLRWSolver only implements the EdS solution. We provide pure-EdS spacetime (no perturbations), and linearly-perturbed EdS spacetime for various kinds of perturbations. Namely, you can choose a single (3-dimensional) sine-wave perturbation, a powerspectrum of perturbations, or your own perturbations read in from files. Currently this thorn implements only a linearly-accurate solution of the constraint equations, so there is a second-order violation in the constraints in the solution from the beginning of the simulation.

Links to detailed tutorials on setting up, compiling, and running simple simulations using FLRWSolver with the ET can be found in doc/tutorials/.

In this documentation we have set \(G=c=1\).

2 Use of this thorn

To use this thorn to provide initial data for the ADMBase variables \(\alpha \), \(\beta \), \(g_{ij}\) and \(K_{ij}\), and HydroBase variables \(\rho \), \(v^i\) activate the thorn and set the following parameters:

• HydroBase::initial_hydro = “flrw”

• ADMBase::initial_data = “flrw”

• ADMBase::initial_lapse = “flrw”

• ADMBase::initial_shift = “flrw”

• ADMBase::initial_dtlapse = “flrw”

• ADMBase::initial_dtshift = “zero”

Template parameter files for a pure-FLRW, single-mode perturbation, and powerspectrum of perturbations to FLRW spacetime are provided in the par/ directory.

This thorn was first presented in [2], and used further in [3]. It is written for use with GRHydro (see Section 6) and EOS_Omni. At the time of writing, neither of these thorns can handle a pure dust description, i.e. \(P=0\). Instead we use a polytropic EOS and ensure \(P\ll \rho \) by setting the polytropic constant accordingly. In [2] this was shown to be sufficient to match a dust FLRW evolution.

3 Choosing initial conditions

3.1 FLRW spacetime & background for perturbations

The FLRW line element in conformal time, \(\eta \), is given by \begin {equation} \relax \expandafter \ifx \csname cur:th\endcsname \relax \expandafter \:label \else \expandafter \l:bel \fi {EinsteinInitialData_flrwsolver_eq:FLRWmetric} ds^2 = a^2(\eta ) \left ( - d\eta ^2 + \delta _{ij}dx^i dx^j \right ) \end {equation} where \(a(\eta )\) is the scale factor describing the size of the Universe at time \(\eta \). To initialise this spacetime, the user must specify FLRW_perturb = “no”.

The user chooses the initial value of \(\mathcal {H}L\), i.e. the ratio of the initial Hubble horizon \(d_H(z_{\rm init})\equiv c/\mathcal {H}(z_{\rm init})\) to the comoving length of the box \(L\) in physical units. Since this quantity is dimensionless, this is the same ratio of the Hubble parameter to the box length in code units. This is controlled using the parameter FLRW_init_HL. This choice of \(\mathcal {H}L\), along with the size of the box set in code units via CoordBase, sets the initial Hubble rate, \(\mathcal {H}_{\rm ini}\), and therefore the initial background density in code units via the Friedmann equation for a matter-dominated spacetime, namely \begin {equation} \mathcal {H}_{\rm init} = \sqrt {\frac {8\pi \rho _{\rm init} a_{\rm init}^2}{3}}, \end {equation} where the initial value for the scale factor, \(a_{\rm ini}\), is set via FLRW_init_a. The initial lapse is set via FLRW_lapse_value. The user should also ensure the initial coordinate time (Cactus::cctk_initial_time) is consistent, i.e. for EdS we have \(\eta _{\rm ini} = 2/\mathcal {H}_{\rm ini}\).

When setting pure FLRW or single-mode perturbations to FLRW (see Section 3.2.1 below), the physical side of the box (or, equivalently, the initial redshift \(z_{\rm init}\)) is not necessary for creating initial data. However, setting a powerspectrum of initial data does require the physical length of the box at the initial time, see Section 3.2.2.

3.2 Linear perturbations to FLRW spacetime

Including scalar only perturbations to the FLRW metric in the longitudinal gauge gives \begin {equation} \relax \expandafter \ifx \csname cur:th\endcsname \relax \expandafter \:label \else \expandafter \l:bel \fi {EinsteinInitialData_flrwsolver_eq:perturbed_metric} ds^2 = a^2(\eta ) \left [ - \left (1 + 2\psi \right ) d\eta ^2 + \left (1 - 2\phi \right ) \delta _{ij}dx^i dx^j \right ], \end {equation} where \(\phi ,\psi \) coincide with the Bardeen potentials [1] in this gauge. Assuming \(\phi ,\psi \ll 1\) allows us to solve Einstein’s equations using linear perturbation theory, giving the system of equations [2, 3]

and \(\phi =\psi \), and we have left \(G, c\neq 1\) in the above for clarity of physical units. The perturbed rest-mass density is \(\rho = \bar {\rho } \left (1 + \delta \right )\), with \(\bar {\rho }\) the background FLRW density. We have \(v^i = \delta v^i\), since \(\bar {v}^i = 0\) for FLRW. In the above, \(\mathcal {H}\equiv a'/a\) is the conformal Hubble parameter, where \('\equiv \partial /\partial \eta \) and \(\partial _i \equiv \partial /\partial x^i\). Solving the above system, choosing only the growing mode (see [3] for more details), gives

where we have the freedom to choose the form of \(f(x^i)\), so long as it has amplitude such that \(\phi \ll 1\). Equations (??) denote the standard form of linear perturbations implemented in this thorn, with \(a\rightarrow a_{\rm init}\) and \(\mathcal {H}\rightarrow \mathcal {H}_{\rm init}\). See [3] for more details on the derivation of these relations.

The user must still specify the FLRW background as in Section 3.1.

Below we outline several different choices for perturbations, and how to set these using the thorn parameters.

3.2.1 Single-mode perturbation

This initial condition sets \(\phi \) as a sine-wave function, and the corresponding density and velocity perturbations are set using (??). Choose FLRW_perturb_type=“single_mode”. The parameter FLRW_perturb_direction controls in which spatial dimension to apply the perturbation, and will set either \(\phi =f(x^1),f(x^2),f(x^3)\), or \(f(x^i)\) depending on the choice. For example, choosing FLRW_perturb_direction=“all” we have \begin {equation} \relax \expandafter \ifx \csname cur:th\endcsname \relax \expandafter \:label \else \expandafter \l:bel \fi {EinsteinInitialData_flrwsolver_eq:phi} \phi = \phi _{0} \sum _{i=1}^{3} \mathrm {sin}\left (\frac {2\pi x^{i}}{\lambda } - \theta \right ), \end {equation} where \(\lambda \) is the wavelength of the perturbation, \(\theta \) is some phase offset, and \(\phi _0\ll 1\). This gives the density and velocity perturbation as, respectively, [2] \begin {align} \delta &= - \left [ \left (\frac {2\pi }{\lambda }\right )^{2} \frac {a_{\mathrm {init}}}{4\pi \rho ^{*}} + 2\right ] \phi _{0} \sum _{i=1}^{3} \mathrm {sin}\left (\frac {2\pi x^{i}}{\lambda } - \theta \right ),\label {EinsteinInitialData_flrwsolver_eq:initial_delta}\\ v^{i} &= -\left (\frac {2\pi }{\lambda }\right )\frac {\mathcal {H}_{\rm init}}{4\pi \rho ^{*}}\, \phi _{0}\, \mathrm {cos}\left (\frac {2\pi x^{i}}{\lambda } - \theta \right ). \label {EinsteinInitialData_flrwsolver_eq:initial_deltav} \end {align}

The wavelength, \(\lambda \), of the perturbation is controlled by single_perturb_wavelength, given relative to the total box length. The phase offset \(\theta \) is controlled by phi_phase_offset, and the amplitude is set by phi_amplitude, which must be set such that \(\phi _0\ll 1\) to ensure that the corresponding density and velocity perturbations are also small enough to satisfy the linear approximation. The code will produce a warning if either of the (dimensionless) density or velocity perturbations have amplitude \(> 10^{-5}\).

3.2.2 Power spectrum of perturbations

We can instead choose the initial conditions to be a power spectrum of fluctuations, to better mimic the early state of the Universe (see [3]). These fluctuations are still assumed to be linear perturbations to an EdS background.

The initial data are generated from a user-specified 3-dimensional matter power spectrum, \(P(k)\) (in units of (Mpc/\(h)^3)\), in the synchronous gauge at the chosen initial redshift, where the wavenumber \(k=\sqrt {k_x^2+k_y^2+k_z^2}\) (in units of \(h\)/Mpc). Here, \(h\) is defined such that \(H_0 = 100 h\) km/s/Mpc. It is assumed the text file is the same format as the matter power spectrum output from CAMB¹ or CLASS² , i.e. a text file with two columns \(k, P(k)\) and no entry for \(k=0\) (this is added by the code).

FLRWSolver then generates a Gaussian random field and scales it to the given power spectrum. This field is centred around zero and represents the initial (dimensionless) density perturbation to the EdS background in synchronous gauge, \(\delta _s\). We choose to first calculate \(\delta _s\) because CAMB/CLASS by default output \(P(k)\) in the synchronous gauge. We calculate the gauge-invariant Bardeen potential from \(\delta _s\) with the Poisson equation \begin {equation} \relax \expandafter \ifx \csname cur:th\endcsname \relax \expandafter \:label \else \expandafter \l:bel \fi {EinsteinInitialData_flrwsolver_eq:poisson} \nabla ^2 \Phi = \frac {3\mathcal {H}^2}{2}\delta _s, \end {equation} and the metric perturbation in longitudinal gauge is \(\phi =\Phi \). We use equations (??) and (??) in Fourier space to find the corresponding density and velocity perturbations in the longitudinal gauge, which are then inverse Fourier transformed back to give the perturbations in real space.

This initial condition is chosen by setting FLRW_perturb_type=“powerspectrum”. Set the FLRW background according to Section 3.1, and set the path to (including the name of) the text file containing the matter power spectrum using FLRW_powerspectrum_file. There is an example power spectrum in FLRWSolver/powerspectra/ which was generated using CLASS³ (see the README in that directory for the exact parameters used), however, we recommend you run CLASS or CAMB⁴ with your desired cosmological parameters. The Gaussian random field will be drawn with seed from FLRW_random_seed, so keep this constant to draw the same realisation more than once. You also must specify the size of the domain in comoving Mpc\(/h\) using FLRW_boxlength, which is required to use the correct physical scales from the power spectrum.

3.2.3 Linear perturbations from user-specified files

There is also an option for you to set up your own initial conditions and read them in to FLRWSolver. These must be in ascii format, and the 3D data must be stored in a specific way. Namely, the data must be written in 2D \(N\times N\) slices of the \(x-y\) plane, with the \(z\) index increasing every \(N\) rows, where \(N=n+n_{\rm gh}\) is the resolution (\(n\)) including the total number of boundary zones in each dimension (\(n_{\rm gh}=2\times \texttt {driver::ghost\_size}\)). Therefore, the columns are the \(x\) indices, and the rows are the \(y\) indices, repeating every \(N\) rows. The files must include boundary zones, although these may be set to zero, since your chosen boundary conditions will be applied after FLRWSolver is called.

To use this initial condition set FLRW_perturb_type=“fileread”. The files you must supply are the linear perturbations around an EdS background, and it is up to you to ensure these satisfy the constraint equations. You must set several string parameters containing the path to each file. The parameter FLRW_deltafile specifies the fractional density perturbation \(\delta \), FLRW_phifile for the metric perturbation \(\phi /c^2\), and FLRW_velxfile, FLRW_velyfile, and FLRW_velzfile for the components of the contravariant fluid three-velocity with respect to the Eulerian observer, \(v^i/c\) (see HydroBase documentation for an explicit definition of \(v^i\)).

4 Tools

There are some useful tools located in FLRWSolver/tools/. Here we describe the purpose and use of each of these files.

• get_init_HL.py: This script is useful for ensuring all background parameters are set consistently. Given a numerical resolution, comoving box length in Gpc\(/h\) (i.e. the proper length at \(z=0\)), and initial redshift, it will give the corresponding initial \(\mathcal {H}L\) and conformal time in code units. The script will also output the final time corresponding to a particular (approximate; based on the EdS prediction) scale factor, and intermediate times corresponding to particular redshifts (again based on EdS predictions) to assist with restarting simulations for increased output frequency. First, please check the parameters in this script carefully to ensure they are consistent with your simulation parameters. The script assumes an EdS evolution for \(\mathcal {H}\) in scaling \(H_0=100 h\) km/s/Mpc back to an initial value of \(\mathcal {H}_{\rm ini}\).

• cut_powerspectrum.py: This script makes a minimum scale cut to a provided power spectrum to remove small-scale structure. We might be interested in doing this if we want to ensure our simulation initially samples all perturbation modes with a certain number of grid cells. This can help reduce the amount of constraint violation due to under-sampled structures. Follow the instructions inside the script and ensure you have set all parameters within the “user changes” section of the script.

• make_test_ics.py: This script’s purpose is to generate initial data for convergence testing. For numerical convergence studies of quantities calculated at the grid scale (i.e., not statistical convergence), we need physical gradients to be consistent between all resolutions. This script will generate three sets of power spectrum initial data (at specified resolutions) which have identical structure at the initial time. This script uses cut_powerspectrum.py to ensure the modes are large-scale and thus remain in the linear regime.

• split_HDF5_per_iteration2(3).py: This script takes the *.xyz.h5 files in your simulation directory output by the ET (using IO::out_mode = "onefile" and IO::out_unchunked = "yes", namely, one file per variable containing all time steps) and splits them into individual files for each time step containing all variables. The format output by this script can be read SPLASH (see tools/README.splash) for visualisation as well as into the Python notebook Plot_ET_data.ipynb.

• Plot_ET_data.ipynb: This notebook provides an example of how to read in ET data with Python, using the split HDF5 file format described above.

• README.tutorial: This file contains links to two documents containing step-by-step tutorials for downloading, compiling, and running the ET with FLRWSolver.

Please correctly cite FLRWSolver if you use any of these tools in your analysis, and do let us know of any improvements that can be made.

5 Gauge

We note here that a useful gauge is the “harmonic-type” gauge specified by \begin {equation} \partial _t \alpha = - f \alpha ^n K, \end {equation} where \(K\) is the trace of the extrinsic curvature. This is implemented in ML_BSSN, and the parameters \(f\) and \(n\) are controlled using ML_BSSN::harmonicF and ML_BSSN::harmonicN, respectively. The initial value of the lapse, \(\alpha \) is set using FLRW_lapse_value. In a perturbed FLRW spacetime (using single_mode or powerspectrum) this is the background value of the lapse, which is then perturbed according to (??). This thorn uses \(\beta ^i=0\).

Setting \(f=1/3\) and \(n=2\) in the above yields \(\partial _t \alpha = \partial _t a\) for EdS (where \(a\) is the scale factor), i.e. we have \(t=\eta \) and our coordinate time of the simulation is the conformal time in the metric (??). This is valid at all times for a pure-FLRW evolution, and within the linear regime for the case of linear initial perturbations. Once any perturbations grow nonlinearly this connection is no longer valid. See, e.g., [3] for a discussion about appropriate gauges for nonlinear structure formation.

6 Notes on using GRHydro for cosmology

The work done in [2] and [3] using this thorn used GRHydro for the hydrodynamic evolution. Below we offer a few notes and suggestions on using this thorn for cosmological evolutions.

6.1 Equation of state

In deriving the initial conditions above, we assume a dust fluid in an FLRW spacetime, i.e. \(P=0\). If using GRHydro for the hydrodynamical evolution, it should be noted that this thorn does not work with an identically zero pressure. To compensate this, it is recommended to use a polytropic EOS, \(P=K \rho ^{\gamma }\) such that \(P\ll \rho \). In [2] this is shown to be sufficient to match the evolution of a dust FLRW spacetime.

The polytropic constant \(K\) (set via EOS_Omni::poly_k) will need to be adjusted accordingly for different choices of initial \(\mathcal {H}\) (since this directly sets the initial density, and larger \(\mathcal {H}\) implies larger \(\bar {\rho }\) in code units) to ensure \(P\ll \rho \), otherwise the evolution will not be close enough to a dust FLRW evolution. Currently no other evolution thorns or EOS choices have been tested with this thorn, but we invite users to explore this and let us know of any success or failures.

6.2 Potential issues with the atmosphere

In GRHydro, the atmosphere is used to solve the problem that part of the computational domain is essentially vacuum (when simulating compact objects; see the relevant documentation). For cosmology, we do not need an atmosphere as the matter fluid is continuous across the whole domain (in the absence of shell-crossings). GRHydro decides whether the position on the grid coincides with the atmosphere by checking several conditions. In some cases we have found that regions of the domain are flagged as being in “the atmosphere” for our cosmological simulations when the value of density in code units is very small. This causes these regions in the domain to be overwritten to the value of GRHydro::rho_abs_min. A way around this is to comment out the line which reliably causes this behaviour (as found in our cases). This will have no effect on the evolution as the purpose of the atmosphere is for simulations of compact objects, and not cosmology.

If you run into this problem, and if you intend to use GRHydro exclusively for cosmology, in the code GRHydro/src/GRHydro_UpdateMask.F90 comment out line 76 that sets atmosphere_mask_real(i,j,k) = 1. If using GRHydro for any simulations of compact objects, be sure to change this back to the distributed form (and recompile), as in these cases the atmosphere is necessary.

References

7 Parameters

flrw_phifile	Scope: restricted	STRING
Description: Path to (including name of) text file specifying the initial metric perturbation, phi/c. This file is assumed to be in a 2D-stacked format, see documentation for more details.
Range		Default: (none)
	Any string containing the path to the metric perturbation file

flrw_velxfile	Scope: restricted	STRING
Description: Path to (including name of) text file specifying the initial velocity, v/c. This is the contravariant fluid three-velocity w.r.t the Eulerian observer (see HydroBase doc). This file is assumed to be in a 2D-stacked format, see documentation for more details.
Range		Default: (none)
	Any string containing the path to the velocity v/c

8 Interfaces

General

Implements:

flrwsolver

Inherits:

hydrobase

eos_omni

admbase

constants

staticconformal

coordbase

Uses header:

SpaceMask.h

9 Schedule

This section lists all the variables which are assigned storage by thorn EinsteinInitialData/FLRWSolver. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage

NONE

Scheduled Functions

HydroBase_Initial (conditional)

  flrw_noperturb

  sets up initial data for flrw spacetime

HydroBase_Initial (conditional)

  flrw_singlemode

  sets up initial data for a single-mode perturbation to flrw spacetime

HydroBase_Initial (conditional)

  flrw_powerspectrum

  sets up initial data for a power spectrum of perturbations to flrw spacetime

HydroBase_Initial (conditional)

  flrw_fileread

  sets up initial data for linear perturbations to flrw spacetime from user-specified ascii files

CCTK_PARAMCHECK

  flrw_paramcheck

  check some parameters are set correctly