xx FLRWSolver

FLRWSolver

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

May 16, 2024

Abstract

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; 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, FLRWSolverxonly 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 FLRWSolverxwith the ET can be found in doc/tutorials/.

In this documentation we have set .

2 Use of this thorn

To use this thorn to provide initial data for the ADMBase variables , , and , and HydroBase variables , activate the thorn and set the following parameters:

     HydroBase::initialxhydro = “flrw”

     ADMBase::initialxdata = “flrw”

     ADMBase::initialxlapse = “flrw”

     ADMBase::initialxshift = “flrw”

     ADMBase::initialxdtlapse = “flrw”

     ADMBase::initialxdtshift = “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 Sectionx6) and EOSxOmni. At the time of writing, neither of these thorns can handle a pure dust description, i.e. . Instead we use a polytropic EOS and ensure 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, , is given by where is the scale factor describing the size of the Universe at time . To initialise this spacetime, the user must specify FLRWxperturb = “no”.

The user chooses the initial value of , i.e. the ratio of the initial Hubble horizon to the comoving length of the box 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 FLRWxinitxHL. This choice of , along with the size of the box set in code units via CoordBase, sets the initial Hubble rate, , and therefore the initial background density in code units via the Friedmann equation for a matter-dominated spacetime, namely

where the initial value for the scale factor, , is set via FLRWxinitxa. The initial lapse is set via FLRWxlapsexvalue. The user should also ensure the initial coordinate time (Cactus::cctkxinitialxtime) is consistent, i.e. for EdS we have .

When setting pure FLRW or single-mode perturbations to FLRW (see Sectionx3.2.1 below), the physical side of the box (or, equivalently, the initial redshift ) 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 Sectionx3.2.2.

3.2 Linear perturbations to FLRW spacetime

Including scalar only perturbations to the FLRW metric in the longitudinal gauge gives where coincide with the Bardeen potentials [1] in this gauge. Assuming allows us to solve Einstein’s equations using linear perturbation theory, giving the system of equations [2,x3]

and , and we have left in the above for clarity of physical units. The perturbed rest-mass density is , with the background FLRW density. We have , since for FLRW. In the above, is the conformal Hubble parameter, where and . 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 , so long as it has amplitude such that . Equations (??) denote the standard form of linear perturbations implemented in this thorn, with and . See [3] for more details on the derivation of these relations.

The user must still specify the FLRW background as in Sectionx3.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 as a sine-wave function, and the corresponding density and velocity perturbations are set using (??). Choose FLRWxperturbxtype=“singlexmode”. The parameter FLRWxperturbxdirection controls in which spatial dimension to apply the perturbation, and will set either , or depending on the choice. For example, choosing FLRWxperturbxdirection=“all” we have where is the wavelength of the perturbation, is some phase offset, and . This gives the density and velocity perturbation as, respectively, [2]

The wavelength, , of the perturbation is controlled by singlexperturbxwavelength, given relative to the total box length. The phase offset is controlled by phixphasexoffset, and the amplitude is set by phixamplitude, which must be set such that 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 .

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, (in units of (Mpc/, in the synchronous gauge at the chosen initial redshift, where the wavenumber (in units of /Mpc). Here, is defined such that km/s/Mpc. It is assumed the text file is the same format as the matter power spectrum output from

CAMB1 ^{1https://camb.info}

or CLASS2 ^{2https://class-code.net}  ,

i.e. a text file with two columns and no entry for (this is added by the code).

FLRWSolverxthen 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, . We choose to first calculate because CAMB/CLASS by default output in the synchronous gauge. We calculate the gauge-invariant Bardeen potential from with the Poisson equation

and the metric perturbation in longitudinal gauge is . 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 FLRWxperturbxtype=“powerspectrum”. Set the FLRW background according to Sectionx3.1, and set the path to (including the name of) the text file containing the matter power spectrum using FLRWxpowerspectrumxfile. There is an example power spectrum in FLRWSolver/powerspectra/ which was generated using

CLASS3 ^{3https://class-code.net} (see the README in that directory for the exact parameters used), however, we recommend you run CLASS or

CAMB4 ^{4https://camb.info} with your desired cosmological parameters. The Gaussian random field will be drawn with seed from FLRWxrandomxseed, so keep this constant to draw the same realisation more than once. You also must specify the size of the domain in comoving Mpc using FLRWxboxlength, 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 slices of the plane, with the index increasing every rows, where is the resolution () including the total number of boundary zones in each dimension (). Therefore, the columns are the indices, and the rows are the indices, repeating every rows. The files must include boundary zones, although these may be set to zero, since your chosen boundary conditions will be applied after FLRWSolverxis called.

To use this initial condition set FLRWxperturbxtype=“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 FLRWxdeltafile specifies the fractional density perturbation , FLRWxphifile for the metric perturbation , and FLRWxvelxfile, FLRWxvelyfile, and FLRWxvelzfile for the components of the contravariant fluid three-velocity with respect to the Eulerian observer, (see HydroBase documentation for an explicit definition of ).

4 Tools

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

getxinitxHL.py: This script is useful for ensuring all background parameters are set consistently. Given a numerical resolution, comoving box length in Gpc (i.e. the proper length at ), and initial redshift, it will give the corresponding initial 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 in scaling km/s/Mpc back to an initial value of .

cutxpowerspectrum.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.

makextestxics.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 cutxpowerspectrum.py to ensure the modes are large-scale and thus remain in the linear regime.

splitxHDF5xperxiteration2(3).py: This script takes the *.xyz.h5 files in your simulation

     directory output by the ET (using IO::outxmode = "onefile" and IO::outxunchunked = "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 PlotxETxdata.ipynb.

PlotxETxdata.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 where is the trace of the extrinsic curvature. This is implemented in MLxBSSN, and the parameters and are controlled using MLxBSSN::harmonicF and MLxBSSN::harmonicN, respectively. The initial value of the lapse, is set using FLRWxlapsexvalue. In a perturbed FLRW spacetime (using singlexmode or powerspectrum) this is the background value of the lapse, which is then perturbed according to (??). This thorn uses .

Setting and in the above yields for EdS (where is the scale factor), i.e. we have 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. . 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, such that . In [2] this is shown to be sufficient to match the evolution of a dust FLRW spacetime.

The polytropic constant (set via EOSxOmni::polyxk) will need to be adjusted accordingly for different choices of initial (since this directly sets the initial density, and larger implies larger in code units) to ensure , 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::rhoxabsxmin. 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/GRHydroxUpdateMask.F90 comment out line 76 that sets atmospherexmaskxreal(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 [1]xxxJ. M. Bardeen, Phys. Rev. D 22, 1882 (1980) [2]xxxH. J. Macpherson, P. D. Lasky, and D. J. Price, Phys. Rev. D 95, 064028 (2017), arXiv:1611.05447 [3]xxxH. J. Macpherson, D. J. Price, and P. D. Lasky, Phys. Rev. D99, 063522 (2019), arXiv:1807.01711 7 Parameters

11flrwxboxlength 10Scope: restrictedx01 REAL

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

11  Range  10  x01  Default: 1.0
110.0:*  10Must be positive.

11

11flrwxdeltafile  10Scope: restrictedx01  STRING

11Description: 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.

11  Range  10  x01  Default: (none)
11                     10Any string containing the path to the density perturbation file

11

11flrwxinitxa  10Scope: restrictedx01  REAL

11Description: FLRW initial scale factor

11 Range 10 x01 Default: 1.0
110.0001:*  10Must be positive and non-zero

11

11flrwxinitxhl  10Scope: restrictedx01  REAL

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

11  Range  10  x01  Default: 1.0
110.0:*  10Must be positive

11

11flrwxlapsexvalue  10Scope: restrictedx01  REAL

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

11  Range  10  x01  Default: 1.0
11*:*  10Anything

11

11flrwxperturb  10Scope: restrictedx01  KEYWORD

11Description: Add perturbation to matter and metric

11 Range 10 x01 Default: no
11yes  10Add perturbation to FLRW spacetime
11no  10Do not add perturbation, will run FLRW spacetime

11

11flrwxperturbxdirection 10Scope: restrictedx01  KEYWORD

11Description: Choose which direction to perturb (only used if FLRWxperturbxtype = singlexmode)

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

11flrwxperturbxtype  10Scope: restrictedx01  KEYWORD

11Description: Add perturbation of the following type

11 Range 10 x01 Default: singlexmode
11singlexmode  10”Single mode perturbation (sine wave in density). Set FLRWxperturbxdirecti on to choose dimensions of perturbation and singlexperturbxwavel ength.”

11powerspectrum  10”A power spectrum of (linear) perturbations. Supply the path to your chosen power spectrum in FLRWxpowerspectrumxf ile”
11fileread  10Linear perturbations to FLRW in longitudinal gauge – read in from supplied ascii files (specific format - see doc)

11

11flrwxphifile  10Scope: restrictedx01  STRING

11Description: Path to (including name of) text file specifying the initial metric perturbation, phi/cˆ2 . This file is assumed to be in a 2D-stacked format, see documentation for more details.

11  Range  10  x01  Default: (none)
11                     10Any string containing the path to the metric perturbation file

11

11flrwxpowerspectrumxfile10Scope: restrictedx01  STRING

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

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

11

11flrwxrandomxseed  10Scope: restrictedx01  INT

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

11  Range  10  x01  Default: 10
110:*  10Must be positive

11

11flrwxvelxfile  10Scope: restrictedx01  STRING

11Description: Path to (including name of) text file specifying the initial velocity, vx̂ /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.

11  Range  10  x01  Default: (none)
11                     10Any string containing the path to the velocity vx̂ /c

11

11flrwxvelyfile  10Scope: restrictedx01  STRING

11Description: Path to (including name of) text file specifying the initial velocity, vx/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.

11  Range  10  x01  Default: (none)
11                     10Any string containing the path to the velocity vx/c

11

11flrwxvelzfile  10Scope: restrictedx01  STRING

11Description: Path to (including name of) text file specifying the initial velocity, vx/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.

11  Range  10  x01  Default: (none)
11                     10Any string containing the path to the velocity vx/c

11phixamplitude  10Scope: restrictedx01  REAL

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

11 Range 10 x01 Default: 1.e-8
110.0:*  10Must be positive

11

11phixphasexoffset  10Scope: restrictedx01  REAL

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

11  Range  10  x01  Default: 0.0
11*:*  10Anything

11

11singlexperturbxwavelength 10Scope: restrictedx01  REAL

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

11 Range 10 x01 Default: 1.0
110.0001:*  10”Must be in (0,1]. Wavelength will be singlexperturbxwavel ength * FLRWxboxlength”

11

11initialxhydro  10Scope: shared from HYDROBASEx01  KEYWORD

11Extends ranges:

11flrw  10FLRWSolver initial hydrobase variables

11

8 Interfaces

General Implements: flrwsolver

Inherits: hydrobase eosxomni 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

HydroBasexInitial (conditional) x flrwxnoperturb x sets up initial data for flrw spacetime x

11x11Language:11fortran
11 x 11Type: 11function

11

HydroBasexInitial (conditional) x flrwxsinglemode x sets up initial data for a single-mode perturbation to flrw spacetime x

11x11Language:11fortran
11 x 11Type: 11function

11

HydroBasexInitial (conditional) x flrwxpowerspectrum x sets up initial data for a power spectrum of perturbations to flrw spacetime x

11x11Language:11fortran
11 x 11Type: 11function

11

HydroBasexInitial (conditional) x flrwxfileread x sets up initial data for linear perturbations to flrw spacetime from user-specified ascii files

11x11Language:11fortran
11 x 11Type: 11function

11

CCTKxPARAMCHECK x flrwxparamcheck x check some parameters are set correctly x

11x11Language:11fortran
11 x 11Type: 11function