The initial condition saga

If reading the previous posts on N-body simulations you have though that initial conditions are a little and easy task, you were wrong! And me with you!
The first things I have understood banging against them were:

• Gadget requires initial conditions (ICs) generate by (for example) N-GenIC

• N-GenIC requires an initial power spectrum from CMBFast that is no longer used

• CMBeasy should substitute CMBFast but it doesn’t, it only works for CMB anisotropies

• CAMB should substitute CMBFast

• the CAMB interface is terrible and not very well documentated

CAMB

Let’s start with the main CAMB interface. If you, like me, need an initial power spectrum as input for your ICs generator you are interested in just few options of the interface. First, you should select “Transfer functions” from the “Actions to Perform” section. You can leave the default selection on “Scalar Cl’s” and “Linear”. After that, check the “Cosmological Parameters” section if it’s ok for you and maybe, leave the default “Initial Scalar Perturbation Mode” that is “Adiabatic”. The last things you should be interested in are the maximum $k$, for me it is $10^4$, “k per logint” (it should be something like the k-sampling, for me, 50), the number of redshift (1) and the (transfer) redshift (0). Now select between “Interpolated Grid” or “Calculated Values” (this parameter switch between and interpolated regular grid in log k or array at actual computed values that are better for later re-interpolation, according with the CAMB README) and choose if you want high precision computation. When you have finished you can click on “Go!”. What you obtain is a page with some links to download:

• camb_*.ini: the configuration file to run the standalone CAMB code on your own
• camb_*.log: the calculation log
• camb_*_scalcls.dat
• camb_*_scalcls.fits
• camb_*_transfer_out_z0.dat: this is the file containing the transfer functions for CDM, baryon, photon, massless neutrino, massive neutrinos, and total (massive) respectively as function of $k$
• camb_*_matterpower_z0.dat: it contains the conventionally normalized matter power spectrum (for baryons+cdm+massive neutrinos), in h/Mpc units

where * is a number that identifies the computation and z0 can be different if you choose to use a different redshift.
You can use camb_*_matterpower_z0.dat as normalized input power spectrum or you can calculate (also as a check) it on you own using the first and the last column of camb_*_transfer_out_z0.dat. You should use the last column because DM simulations represents with DM particles all the mass, included that of the baryons, so the initial power spectrum should be the total power spectrum.
If you want to understand better how CAMB works you can try to read the CAMB README.

Initial power spectrum theory

Here some theory if you want to understand how to calculate the power spectrum on your own. Maybe it is well-known and trivial, but for me it wasn’t and, like many of the trivial and well-known things in (not only) cosmology, it’s not easy to find anywhere.
So, the spectrum is defined as $P(k) = A k^n T(k)^2$, where k $n$ is the “primordial spectral index” and is taken near the unity. This is a “scale-fee” spectrum. $T(k)$ is the transfer function that give a synthetic and parametric description of how the initial spectrum survive to the microphysic. $A=\left[\frac{D(z_{\rm fin})}{D(z_{\rm in})}\right]^2$ is the amplitude.
The normalization of the spectrum is given by the value of $\sigma_8$, that is the mean square amplitude of the density field filtered on the scale of 8 Mpc/h. This values comes from the ‘80s, when Peebles and others (Davis & Peebles 1983) measured $\sigma_{\rm galaxies}$ and $\sigma_{\rm gal}(R=8)\sim1$ so they took that values as reference.
$\sigma_8$ is defined by $\sigma^2(R) = \frac{1}{(2\pi)^3}\int \mathrm{d}^3kP(k)\tilde W(kR)$ with $R=8\mathrm{Mpc/h}$ and $\tilde W(kR)$ the Fourier transform of the window (filter) function, usually a top-hat in the positions space.
The last thing we need to know to obtain the spectrum is the value of the amplitude, and it can be find by imposing its value so that $\sigma_8$ has a certain (observed) value.
Because
$$\sigma^2(R) = \frac{1}{(2\pi)^3}\int \mathrm{d}^3kP(k)\tilde W(kR) = \int \mathrm{d}^3k Ak^nT^2(k))\tilde W(kR)$$
we have
$$A_0 = \frac{s^2_R(R=8) }{ \int\mathrm{d}^3k\, k^n T^2(k) \tilde W(8*k)}$$.
Usually $A_0$ is calculated for $\sigma_8=1$ and then scaled with $A = A_0\sigma^2_{\rm 8;obs}$ where $\sigma_{\rm 8;obs}$ is the observed values for $\sigma_8$. With the last observations we have $\sigma_8 = 0.8118405$.

Just to have an idea you can do these computation with a simple python code. The code below compare the computation done in Python with the values from camb_*_transfer_out_z0.dat with the normalized power spectrum from camb_*_matterpower_z0.dat.

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt

# Load data
transfer = np.genfromtxt('camb_88704620_transfer_out_z0.dat', usecols = (0,6))
matterpower = np.genfromtxt('camb_88704620_matterpower_z0.dat')

# Python calculations
camb_k = transfer[:,0]
camb_tf = transfer[:,1] 
R = 8 #Mpc/h
s_8 = 0.8118405
sp_ind = 1

# Define the FT of the window function
def FTW(R, k):
  return 3*(np.sin(k*R)-k*R*np.cos(k*R)) / (k*R)**3

# camb_k**(2+sp_ind) that is k^(2+n) because d^3k=4\pi k^2dk
amp_integrand = camb_k**(2+sp_ind)*camb_tf**2 * FTW(8, camb_k)**2
amp_integral = integrate.trapz(amp_integrand, camb_k)
amp_0 = 2*np.pi**2/amp_integral
amp = amp_0*s_8**2
spectrum = camb_k**sp_ind*camb_tf**2 * amp

ax = fig.add_subplot(111)
ax.set_title('Fortran/Python CDM initial power spectrum')
ax.set_xlabel('k')
ax.set_ylabel('P(k)')
ax.set_xscale('log')
ax.set_yscale('log')

ax.plot(transfer[:,0], spectrum[:], color = "magenta", 
           linestyle = '-', marker = '', label = "* python amp") 
ax.plot(matterpower[:,0], matterpower[:,1], color = "black", 
            linestyle = '--', marker = '', label = "matterpower")
ax.legend(loc='best')
ax.grid(True)

# Adjust figure size and save
fig.set_size_inches(20, 20)
plt.savefig('camb_f90_py_check')

The result is this

If you want you can check the power spectrum we have obtained by integrating it to find $\sigma_8$:

sigma_integrand = camb_k**2 * spectrum * FTW(8, camb_k)**2
s_8_check = pow(integrate.trapz(sigma_integrand, camb_k)/(2*np.pi**2), 0.5)
print "s_8_calculated", s_8_check
print "s_8 observed", s_8

obtaining

ziosi@uno:~/Code/spettro_CMB$ ./CAMB_check_plot.py
s_8_calculated 0.8118405
s_8 observed 0.8118405

N-GenIC

Now we have the initial power spectrum ready for the ICs generator. After N-GenIC have been compiled (try to read this if you have compilation problems related to the parallel double precision FFTW libraries or if you want to know how to customize the Makefile) we should have a look at the configuration file.
We are interested in:

• Nmesh: the size (=the number of nodes) of the FFT grid used to compute the displacement field, should be Nmesh $\geq$ Nsample
• Nsample: this is the maximum k that the code uses, i.e. this effectively determines the Nyquist frequency that the code assumes, $k_{\rm Nyquist} = 2\cdot \pi/{\rm Box} \cdot {\rm Nsample}/2$ Normally, one chooses Nsample such that ${\rm Ntot} = {\rm Nsample}^3$, where Ntot is the total number of particles. Because the grid sample the particles quantities, Nmesh sets the Nyquist frequency of Nsample, so it’s good if ${\rm Nmesh} = 2\cdot {\rm Nsample}$
• ReNormalizeInputSpectrum: set it to 0 because we are going to use the previous spectrum that is already normalized, if you don’t remember this the code will have integration problems
• TileFac: it represents how many times you need to tile the glass file (for each dimension) to cover the number of particles you want to use. The glass file contains 4097 particles. Glass particles positions will be automatically stretched to cover the box dimension. When you download N-GenIC you find ${\rm Nsample} = 128$ and ${\rm TileFac} = 8$, this is because the total number of particles is ${\rm Ntot} = {\rm Nsample}^3 = 128^3=2097125$ and the number of glass particles is ${\rm TileFac}^3\cdot {\rm Nglass} = 8^3\cdot 4096 = 2097125$. In practice, if you want to know what TileFac should be, and you have Ntot particles in you simulation, TileFac will be $\frac{{\rm Ntot}^{1/3}}{4096} = \frac{{\rm Nsample}}{4096}$.
• WhichSpectrum: let you choose if you want to use an internal spectrum (calculated with a function) or the spectrum from CAMB
• FileWithInputSpectrum: it’s, obviously, the name of the file containing the spectrum

The other parameters are cosmological parameters or are related to the folders, the name and the number of files, the parallelization and to the internal measure units.
Other options are (more or less) documented with comments in the code or in the README.
We can now start N-GenIC with

mpiexec -np 2  ./N-GenIC  ics.param

where -np let you set the number of processors you want to use in parallel.
Before using the output files with gadget we should open those and calculate the power spectrum to check that this realization of it is a good one. This problem arise because of the sampling of the k-space where few modes are available, but I will deepen on those matter in a future post.
There is also an improved version of N-GenIC, 2LPTIC, but it need a different installation of the FFTW so I didn’t try it.

Many of these things can be found here.