Tutorial: Power Spectra - galacticusorg/galacticus GitHub Wiki

In this tutorial we'll use Galacticus to compute and output the cosmological linear theory power spectrum.

If you haven't already installed Galacticus, you can find the installation instructions here.

Running the calculation

Galacticus works by reading a parameter file that describes what you want it to do, and gives the numerical values for all relevant parameters (e.g. cosmological density parameters). For this tutorial we'll use the parameters/tutorials/powerSpectrum.xml parameter file. To run Galacticus on this parameter file issue the following command from the source directory where you installed Galacticus (note that the initial $ represents your command line prompt - don't type it in!):

$ ./Galacticus.exe parameters/tutorials/powerSpectrum.xml

If everything is working you should see output which looks something like:

              ##                                     
   ####        #                  #                  
  #   #        #             #                       
 #       ###   #  ###   ### ###  ##   ### ## ##   ## 
 #       #  #  #  #  # #  #  #    #  #  #  #  #  #   
 #   ###  ###  #   ### #     #    #  #     #  #   #  
  #   #  #  #  #  #  # #     #    #  #     #  #    # 
   ####  #### ### ####  ###   ## ###  ###   #### ##  

 © 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016,
   2017, 2018, 2019, 2020
   - Andrew Benson

M: Memory:         code +      nodes +       misc =      total
M:            22.246Mib +   1.000  b +   9.000  b =  22.246Mib
M: -> Begin task: power spectrum
M: <- Done task: power spectrum

and a file powerSpectrum.hdf5 will have been created.

Understanding the input

You can look at the entire parameter file for this tutorial here. Below we'll explore each section of the file to understand what it does. Galacticus parameter files use XML.

  <!-- Specify tasks to perform -->
  <task value="powerSpectra"/>

This block tells Galacticus what task we want it to perform - here we're just asking it to compute the power spectrum.

  <!-- Cosmological parameters -->
  <cosmologyFunctions  value="matterLambda"/>
  <cosmologyParameters value="simple">
    <HubbleConstant  value="70.20000"/>
    <OmegaMatter     value=" 0.27250"/>
    <OmegaDarkEnergy value=" 0.72750"/>
    <OmegaBaryon     value=" 0.04550"/>
    <temperatureCMB  value=" 2.72548"/>
  </cosmologyParameters>

In this block we're setting up a cosmological model. The first option specifies that cosmological functions (e.g. the expansion factor as a function of time) are computed assuming a universe containing matter and a cosmological constant ("lambda"). The next option specifies that we want to use the "simple" cosmological parameters type which just allows us to directly specify the values of all cosmological parameters. Within this option we specify the actual values for the Hubble constant, density parameters, and CMB temperature.

  <!-- Power spectrum options -->
  <transferFunction value="eisensteinHu1999">
    <neutrinoNumberEffective value="3.046"/>
    <neutrinoMassSummed      value="0.000"/>
  </transferFunction>
  <powerSpectrumPrimordial value="powerLaw">
    <index               value="0.961"/>
    <wavenumberReference value="1.000"/>
    <running             value="0.000"/>
  </powerSpectrumPrimordial>
  <powerSpectrumPrimordialTransferred value="simple"/>
  <cosmologicalMassVariance value="filteredPower">
    <sigma_8   value="0.807" />
    <tolerance value="1.0e-3"/>
  </cosmologicalMassVariance>

Next we set up the linear theory power spectrum. We choose to use the Eisenstein & Hu (1999) CDM transfer function, a power-law primordial power spectrum (with no running of the index). The "powerSpectrumPrimordialTransferred" options specifies how to combine these into the linear theory power spectrum. The "simple" option just multiplies the primordial power spectrum by the transfer function squared. Finally, we specify a normalization for the power spectrum by choosing an option for computing the variance of the mas density field ("filteredPower" just means that we integrate the power spectrum under some window function - a top-hat function by default - to compute the variance), and then giving a value for the sigma_8 parameter. We also specify a tolerance for the integration of the power spectrum here - this determines how accurately the variance of the mass density field will be computed.

  <!-- Structure formation options -->
  <linearGrowth value="collisionlessMatter"/>

This section deals with structure growth. We choose to model linear growth assuming collisionless matter.

  <!-- Output options -->
  <galacticusOutputFileName value="haloMassFunction.hdf5"/>
  <outputTimes value="list">
    <redshifts value="0.0 1.0"/>
  </outputTimes>

Finally we specify a file to output the results to, and a set of redshifts at which we want the power spectrum to be computed and output.

Understanding the output

The output file powerSpectrum.hdf5 is an HDF5 file - it contains datasets organized into a hierarchical tree of "groups" (like files in a directories in a filesystem). Most programming languages have tools to interact with HDF5 files. Here we'll explore the output file using the command-line tools h5ls and h5dump.

To see the content of the file use:

$ h5ls powerSpectrum.hdf5
Build                    Group
Outputs                  Group
Parameters               Group
Version                  Group

The h5ls command shows the content of the top-level group of the output file. We'll ignore most of the groups in this tutorial, and focus just on Outputs which contains the main output quantities. We can look inside this group using:

$ h5ls powerSpectrum.hdf5/Outputs
Output1                  Group
Output2                  Group

There are two output groups - corresponding to the two output redshifts specified in our parameter file. Output groups are numbered from the earliest to latest times, so in this case Output1 is the earliest time (z=1) and Output2 is the latest time (z=0).

Let's look at the output for z=0:

$ h5ls powerSpectrum.hdf5/Outputs/Output2
alpha                    Dataset {61}
growthFactor             Dataset {61}
growthFactorLogDerivative Dataset {61}
mass                     Dataset {61}
powerSpectrum            Dataset {61}
sigma                    Dataset {61}
wavenumber               Dataset {61}

Finally we find the power spectrum information! We'll ignore most of these datasets in this short tutorial and just focus on wavenumber and powerSpectrum. You can see the content of these datasets using:

$ h5ls -d powerSpectrum.hdf5/Outputs/Output2/wavenumber
wavenumber               Dataset {61}
    Data:
        (0) 0.001, 0.00125892541179417, 0.00158489319246111, 0.00199526231496888, 0.00251188643150958, 0.00316227766016838, 0.00398107170553497, 0.00501187233627273, 0.00630957344480193, 0.00794328234724282, 0.01,
        (11) 0.0125892541179417, 0.0158489319246111, 0.0199526231496888, 0.0251188643150958, 0.0316227766016838, 0.0398107170553497, 0.0501187233627272, 0.0630957344480193, 0.0794328234724281, 0.1, 0.125892541179417,
        (22) 0.158489319246111, 0.199526231496888, 0.251188643150958, 0.316227766016838, 0.398107170553497, 0.501187233627272, 0.630957344480193, 0.794328234724282, 1, 1.25892541179417, 1.58489319246111,
        (33) 1.99526231496888, 2.51188643150958, 3.16227766016838, 3.98107170553497, 5.01187233627272, 6.30957344480193, 7.94328234724282, 9.99999999999999, 12.5892541179417, 15.8489319246111, 19.9526231496888,
        (44) 25.1188643150958, 31.6227766016837, 39.8107170553498, 50.1187233627273, 63.0957344480194, 79.4328234724281, 100, 125.892541179417, 158.489319246111, 199.526231496888, 251.188643150958, 316.227766016838,
        (56) 398.107170553497, 501.187233627272, 630.957344480193, 794.328234724282, 1000

and

$ h5ls -d powerSpectrum.hdf5/Outputs/Output2/powerSpectrum
powerSpectrum            Dataset {61}
    Data:
        (0) 19752.765289969, 24275.8564961514, 29641.8593656998, 35882.4927362934, 42954.3023662321, 50704.5590320977, 58839.9340448319, 66901.937276083, 74241.9978355227, 79967.821801863, 82844.7851128752,
        (11) 81397.937733453, 75005.5647253165, 65337.3202631851, 54893.6236162408, 44836.2992006317, 35472.7444185093, 27064.7575914213, 19889.6111811954, 14100.9335074837, 9675.51818378156, 6450.72067033403,
        (22) 4195.23476538419, 2670.86480989484, 1669.46464710895, 1026.96411022781, 622.851533931182, 372.984109335854, 220.787916727825, 129.318408458248, 75.008972039714, 43.1187493270437, 24.5824755625454,
        (33) 13.9083861955195, 7.81421081814523, 4.36211699619959, 2.42070001882387, 1.33605767868285, 0.733738111969742, 0.401109051059581, 0.21834693384002, 0.11839620191705, 0.0639681164280232, 0.034446182833649,
        (44) 0.0184916545920113, 0.00989835104433317, 0.00528431342684245, 0.00281404516782447, 0.00149506907802831, 0.000792580963481224, 0.00041931165533301, 0.000221408739155846, 0.000116698731224898,
        (53) 6.14039039448803e-05, 3.22571173377652e-05, 1.69197052567665e-05, 8.8620191155726e-06, 4.63529660381791e-06, 2.42134808751638e-06, 1.26327961359142e-06, 6.58309379694076e-07

The "-d" option to h5ls tells it to actually output the content of the dataset.

As you'd expect, wavenumber contains a list of wavenumbers (in units of Mpc-1) at which thepower spectrumn was tabulated, while powerSpectrum contains the power spectrum itself (in units of Mpc3).

You can also explore the content of the HDF5 file using the h5dump tool, which shows some extra, useful information (I'll cut out much of the very long output for brevity):

$ h5dump -A -g Outputs/Output2 powerSpectrum.hdf5
HDF5 "haloMassFunction.hdf5" {
   ATTRIBUTE "outputRedshift" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SCALAR
      DATA {
      (0): 0
      }
   }
   ATTRIBUTE "outputTime" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SCALAR
      DATA {
      (0): 13.7915
      }
   }
   .
   .
   .
   DATASET "powerSpectrum" {
   COMMENT "The power spectrum."
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 61 ) / ( 61 ) }
      ATTRIBUTE "unitsInSI" {
         DATATYPE  H5T_IEEE_F64LE
         DATASPACE  SCALAR
         DATA {
         (0): 2.938e+67
         }
      }
   }
   .
   .
   .
}
}

For example, you can see that the output group contains some "attributes" such as the redshift and time of the output, and that each dataset has a description and a "unitsInSI" attribute which is the value you would multiply this dataset by to convert to SI units (and from there you can convert it to whatever units you want).

The full list of output datasets is:

  • alpha: The logarithmic slope of σ(M): α = d ln σ / d ln M;
  • mass: The mass scale, M, corresponding to the given wavenumber, k, defined such that M=4 π ΩM ρcrit / 3 k3 (in units of M);
  • powerSpectrum: The linear theory power spectrum at z=0: P(k) in units of Mpc3;
  • sigma: The dimensionless linear theory mass fluctuation at z=0: σ(M);
  • wavenumber: The wavenumber in units of Mpc-1.
⚠️ **GitHub.com Fallback** ⚠️