H.  J. Macpherson <hayleyjmacpherson@gmail.com >,
P. D. Lasky,
D. J. Price

June 14, 2024


This thorn provides cosmological initial conditions based on a Friedmann-Lemaitre-Robertson-Walker (FLRW) spacetime, either with or without small (linear) perturbations.

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:

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 [23]

\begin {subequations} \label {EinsteinInitialData_flrwsolver_eqs:perturbed_einstein} \begin {align} \nabla ^{2}\phi - 3 \mathcal {H}\left (\frac {\phi '}{c^2} + \frac {\mathcal {H} \phi }{c^2}\right ) &= 4\pi G \bar {\rho }\,\delta a^{2}, \label {EinsteinInitialData_flrwsolver_eq:einstein_1} \\ \mathcal {H} \partial _{i}\phi + \partial _{i}\phi ' &= -4\pi G\bar {\rho } \,a^{2} \delta _{ij}v^{j}, \label {EinsteinInitialData_flrwsolver_eq:einstein_2} \\ \phi '' + 3\mathcal {H}\phi ' &=0, \label {EinsteinInitialData_flrwsolver_eq:einstein_3} \end {align} \end {subequations}

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

\begin {subequations} \label {EinsteinInitialData_flrwsolver_eqs:linear_solnsg0} \begin {align} \phi &= f(x^{i}), \label {EinsteinInitialData_flrwsolver_eq:linear_phi}\\ \delta &= \frac {2}{3\mathcal {H}^2} \nabla ^{2}f(x^{i}) - 2 f(x^i) \\ v^{i} &= -\frac {2}{3 a \mathcal {H}} \partial ^{i}f(x^{i}), \end {align} \end {subequations}

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 CAMB1 or CLASS2 , 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 CLASS3 (see the README in that directory for the exact parameters used), however, we recommend you run CLASS or CAMB4 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.

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.


[1]   J. M. Bardeen, Phys. Rev. D 22, 1882 (1980)

[2]   H. J. Macpherson, P. D. Lasky, and D. J. Price, Phys. Rev. D 95, 064028 (2017), arXiv:1611.05447

[3]   H. J. Macpherson, D. J. Price, and P. D. Lasky, Phys. Rev. D99, 063522 (2019), arXiv:1807.01711

7 Parameters

Scope: restricted  REAL

Description: Size of box used in comoving Mpc – only used in generating power spectrum perturbations

Range   Default: 1.0
Must be positive.

Scope: restricted  STRING

Description: Path to (including name of) text file specifying the initial fractional density perturbation, delta. 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 density perturbation file

Scope: restricted  REAL

Description: FLRW initial scale factor

Range   Default: 1.0
Must be positive and non-zero

Scope: restricted  REAL

Description: FLRW initial H_* L – ratio of box length to horizon size. This sets initial background density

Range   Default: 1.0
Must be positive

Scope: restricted  REAL

Description: Value of background initial lapse (alpha) – for perturbations this becomes alp=FLRW_lapse_value*sqrt(1 + 2 Phi)

Range   Default: 1.0

Scope: restricted  KEYWORD

Description: Add perturbation to matter and metric

Range   Default: no
Add perturbation to FLRW spacetime
Do not add perturbation, will run FLRW spacetime

Scope: restricted  KEYWORD

Description: Choose which direction to perturb (only used if FLRW_perturb_type = single_mode)

Range   Default: all
perturb x-direction only
perturb y-direction only
perturb z-direction only
perturb all dimensions

Scope: restricted  KEYWORD

Description: Add perturbation of the following type

Range   Default: single_mode
”Single mode perturbation (sine wave in density). Set FLRW_perturb_directi on to choose dimensions of perturbation and single_perturb_wavel ength.”
”A power spectrum of (linear) perturbations. Supply the path to your chosen power spectrum in FLRW_powerspectrum_f ile”
Linear perturbations to FLRW in longitudinal gauge – read in from supplied ascii files (specific format - see doc)

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

Scope: restricted  STRING

Description: Path to (and name of) text file specifying the matter P(k) (in synchronous gauge) to be used to generate scalar perturbations if perturb_type=powerspectrum is chosen.

Range   Default: (none)
Any string containing the path to the powerspectrum file

Scope: restricted  INT

Description: Random seed for generating the Gaussian random density perturbation from the power spectrum given in FLRW_powerspectrum_file.

Range   Default: 10
Must be positive

Scope: restricted  STRING

Description: Path to (including name of) text file specifying the initial velocity, vˆx  /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ˆx  /c

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

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

Scope: restricted  REAL

Description: Amplitude of perturbation in phi - implicitly sets amplitude of rho,vel perturbation

Range   Default: 1.e-8
Must be positive

Scope: restricted  REAL

Description: The phase offset, eta, in phi=A*sin(kx*x - eta) for single mode perturbation – applied same in all dimentions if FLRW_perturb_direction = all

Range   Default: 0.0

Scope: restricted  REAL

Description: Wavelength of single mode scalar perturbation relative to box size

Range   Default: 1.0
”Must be in (0,1]. Wavelength will be single_perturb_wavel ength * FLRW_boxlength”

Scope: shared from HYDROBASE  KEYWORD

Extends ranges:

FLRWSolver initial hydrobase variables

8 Interfaces











Uses header:


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.



Scheduled Functions

HydroBase_Initial (conditional)


  sets up initial data for flrw spacetime


 Type: function

HydroBase_Initial (conditional)


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


 Type: function

HydroBase_Initial (conditional)


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


 Type: function

HydroBase_Initial (conditional)


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


 Type: function



  check some parameters are set correctly


 Type: function