A Cosmology Forecast Toolkit – CosmoLib

# A Cosmology Forecast Toolkit – CosmoLib

Zhiqi Huang CEA, Institut de Physique Théorique, 91191 Gif-sur-Yvette cédex, France
CNRS, URA 2306, F-91191 Gif-sur-Yvette, France
July 17, 2019
###### Abstract

The package CosmoLib is a combination of a cosmological Boltzmann code and a simulation toolkit to forecast the constraints on cosmological parameters from future observations. In this paper we describe the released linear-order part of the package. We discuss the stability and performance of the Boltzmann code. This is written in Newtonian gauge and including dark energy perturbations. In CosmoLib the integrator that computes the CMB angular power spectrum is optimized for a -by- brute-force integration, which is useful for studying inflationary models predicting sharp features in the primordial power spectrum of metric fluctuations. As an application, CosmoLib is used to study the axion monodromy inflation model that predicts cosine oscillations in the primordial power spectrum. In contrast to the previous studies by Aich et al and Meerburg et al, we found no detection or hint of the osicllations. We pointed out that the CAMB code modified by Aich et al does not have sufficient numerical accuracy. CosmoLib and its documentation are available at http://www.cita.utoronto.ca/~zqhuang/CosmoLib.

## I Introduction

The hot big bang model and the cosmological perturbation theory, where the physical metric is perturbed around the spatially homogeneous and isotropic Friedmann-Robertson-Walker (FRW) metric Peebles:1980 (); Peacock:1999 (); Mukhanov:2005 (); Weinberg:2008a (), have led to a remarkable success in interpreting the plethora of observational data of the last two decades Riess/etal:2011 (); Amanullah/etal:2010 (); Sullivan/etal:2011 (); Reid/etal:2010 (); Percival/etal:2010 (); Komatsu/etal:2011 (). Observations of the temperature anisotropy in the Cosmic Microwave Background (CMB) have been playing an essential role in building the standard cosmological model and measuring its parameters Komatsu/etal:2011 (); Fixsen/etal:1996 (). In order to maximize the usage of the observational data, one would like to compute the theoretical prediction on the CMB anisotropy for a given model as accurately as possible, with tolerable time consumption. Computation tools developed over the years such as CMBFAST Seljak/Zaldarriaga:1996 (), CAMB Lewis/etal:2000 (), CMBEASY Doran:2005 () and CLASS Lesgourgues:2011 (); Blas/etal:2011 () are capable of computing a CMB angular power spectrum to percent-level accuracy within a few seconds on a modern desktop personal computer.

The crucial technique used in all the fast CMB codes to date is the line-of-sight integration approach Seljak/Zaldarriaga:1996 (); Hu/White:1997 (); Hu/etal:1998 () and an assumption that the primordial power spectrum of metric perturbations is smooth. (Here for readability we focus on the comoving curvature perturbations and temperature anisotropies, although the same arguments can be as well applied to the tensor perturbations and CMB polarization.) A CMB code first computes the radiation transfer function by solving the linear-order Boltzmann equations and using the line-of-sight integration method, then convolves with the primordial power spectrum to obtain the CMB angular power spectrum . The smoothness assumption allow us to compute only a few tens of multipoles spanning from to a few thousands and interpolate the remaining ’s. Furthermore, since is assumed to a smooth function, sparse sampling of the radiation transfer function has been implemented for the integration of each .

A smooth primordial power spectrum is a prediction of the simplest single-field slow-roll inflation models Abbott/Wise:1984 (); Lucchin/Matarrese:1985 (); Lucchin/Matarrese:1985a (); Lyth/Stewart:1992 (); Stewart/Lyth:1993 (). However, local signatures in the primordial power spectrum that makes it deviate from smoothness can arise in various alternative models, for instance, when the inflaton potential has sharp features Starobinsky:1992 (); Chen/etal:2007 (), when there is a transition between different stages in the inflaton evolution Contaldi/etal:2003 (); Cline/etal:2003 (), when more than one field is present Polarski/Starobinsky:1992 (); Langlois/Vernizzi:2005 (), from particle production during inflation Barnaby/etal:2009a (); Barnaby/Huang:2009 (), modulated preheating  Chambers/Rajantie:2008 (); Bond/etal:2009 (), or, more recently, in models motivated by monodromy in the extra dimensions Silverstein/Westphal:2008 () (see also Bean/etal:2008 ()). These features represent an important window on new physics because they are often related to UV scale phenomena inaccessible to experiments in the laboratory. For these models, the CMB angular power spectrum is not necessarily smooth, and therefore needs to be computed at each multipole without interpolation. This increases the computing time by a factor of a few tens. Moreover, for the numerical-integration of each , the sampling frequency in the wavenumber often needs to be increased, again, by a factor of a few tens. The required sampling frequency in is model-dependent. It is determined by the larger between the minimum width of the features in the primordial power spectrum and the minimum width of the oscillations in the radiation transfer function.

To keep track of the features in the primordial power spectrum, one can modify standard CMB codes by naively doing an -by- brute-force calculation with increased integration sampling frequency in . However, in the case where the features in the primordial power spectrum are really sharp (), this naive modification increases the computing time by a factor of (a few tens in sampling times a few tens in sampling). Moreover, the memory that is required to store all the transfer functions and tables of spherical Bessel functions can be too large for most desktop personal computers. One of the purposes of this paper is to introduce a more optimized algorithm to treat these problems. In fact, apart from increasing the sampling frequency, that cannot be avoided, all the other problems can be significantly alleviated by using the recurrence relation of spherical Bessel functions. An optimized algorithm, which we detail in Section III, is times slower than the standard algorithm for the smooth- case. This new algorithm has been implemented in the CosmoLib package, a self-contained package that we developed to compute cosmological perturbations, CMB angular power spectra, and the forecast constraints on cosmological parameters from future cosmological surveys using Fisher matrix analysis and Monte Carlo Markov Chain (MCMC) calculation. In particular, the cosmological surveys that we consider are CMB, large scale structure (LSS) and supernovae (SN).

In addition to the enhanced CMB integrator, CosmoLib has a few other features that are complementary to the other publicly available Boltzmann/CMB/MCMC codes. For instance, the MCMC engine in CosmoLib has a modified rejection rule that allows the proposal density (the probability of random-walking to a new point in the parameter space) to periodically depend some parameters. This is useful when one considers a likelihood that depends on some periodic parameter. This happens, for instance, in the context of inflation from axion monodromy Silverstein/Westphal:2008 (); Chen/etal:2008 (); McAllister/etal:2010 (); Flauger/etal:2010 (), where the oscillations in the predicted power spectrum depend on a free phase. Moreover, CosmoLib treats the dark energy equation of state (EOS) and the primordial scalar and tensor power spectra and , as free functions, which can be either chosen from a list of build-in models or defined by the user. This makes CosmoLib a convenient tool to study non-standard parametrizations of dark energy EOS and primordial power spectra. Finally, CosmoLib is written in Newtonian gauge (also called Poisson gauge) Mukhanov/etal:1992 (); Ma/Bertschinger:1995 (); Bertschinger:1996 (), while many other codes are mainly developed in synchronous gauge (see e.g. Peacock:1999 ()). This is a plus-and-minus point. We found that our Newtonian-gauge Boltzmann code is slightly slower than the codes written in synchronous gauge. However, many theoretical works in the literature have derived equations in Newtonian gauge. For instance, second-order Boltzmann equations have been derived in this gauge Bruni/etal:1997 (); Bartolo/etal:2006 (); Bartolo/etal:2007 (); Nakamura:2007 (); Enqvist/etal:2007 (); Nakamura:2009 (); Malik/Wands:2009 (); Boubekeur/etal:2009 (); Pitrou:2009 (); Nitta/etal:2009 (); Senatore/etal:2009 (); Nakamura:2010 (); Beneke/Fidler:2010 (); Bernardeau/etal:2011 (). Implementing these equations in a code already in Newtonian gauge would be much easier. To conclude the discussion, we list the differences between CosmoLib and other publicly available CMB codes in Table 2.

As an application, CosmoLib is used to study the “hints” of cosine osicllations in the primordial power spectrum that was recently found in Refs. Aich/etal:2011 (); Meerburg/etal:2012 (). In an accompanying paper Huang/etal:2012 (), CosmoLib is applied to forecast the constraining power of future CMB and galaxy survey data on the primordial power spectrum from inflation, with an emphasis on models generating features in the power spectrum.

This paper is organized as follows. In Section II we introduce the Boltzmann code in Newtonian gauge and discuss its stability and performance. Section III details the algorithm used in the enhanced CMB integrator. In Section IV we introduce the forecast technique and parameter estimation methods. Section V concludes.

Throughout this paper, unless otherwise specified, repeated indices are summed over. Greek indices run from 0 to 3. Latin indices run from 1 to 3, that is only over spatial dimensions. We use natural units and the reduced Planck Mass .

## Ii CosmoLib in Newtonian Gauge

### ii.1 The Background Solutions

Let us start discussing the background solutions. We consider a flat FRW metric , where is the scale factor and is the conformal time. The normalization of is arbitrary. We normalize it such that today. CosmoLib uses the e-fold number as the time variable. The physical Hubble expansion rate is defined as . Its present value is denoted by .

We assume a universe with cold dark matter (labeled with a subscript ), dark energy (labeled with a subscript ), baryons (labeled with a subscriber ), radiation (labeled with a subscript ), and 3 species of massless neutrinos (labeled with a subscript ). For a component () the background density is denoted as , and the background pressure . The present-day fractional energy density is written as . Dark energy is assumed to be a perfect fluid with known equation of state (pressure to density ratio) and a constant sound speed in its rest frame. The users can either choose from a list of build-in models or define their own functions. The build-in models of include the cosmological constant model Einstein:1917 (), a constant EOS , a linear function Chevallier/Polarski:2001 (); Linder:2003 (), and a general three-parameter parametrization for the minimally coupled quintessence/phantom models Huang/etal:2011 ().

For a given the background solutions are

 a=eN,ρc=3H20M2pΩc0a−3,pc=0,ρb=3H20M2pΩb0a−3,pb=0,ργ=3H20M2pΩγ0a−4,pγ=13ργ,ρν=3H20M2pΩν0a−4,pν=13ρν,ρΛ=3H20M2pΩΛ0a−3exp[3∫0Nw(a)dN],pΛ=w(a)ρΛ,H=1Mp√ρc+ρb+ργ+ρν+ρΛ3. (1)

We will also use the derived quantities (), , and

 ϵ=−dlnHdN=32[1+pΛ+pγ+pνρc+ρb+ρΛ+ργ+ρν]. (2)

The conformal time can be related to the scale factor via

 τ=∫a0daHa2. (3)

The electron number density is obtained using RecFast version 1.5 Seager/etal:1999 (); Wong/etal:2008 (), which has been incorporated into CosmoLib. We denote the differential optical depth (increment of optical depth per ) as

 κN≡dκdN=neσTH, (4)

where is the Thomson scattering cross section. The baryon sound speed is obtained by solving the differential equations (68-69) in Ref. Ma/Bertschinger:1995 ().

With these background solutions in hand, now we can write down the governing equations for scalar perturbations.

### ii.2 Scalar Perturbations

The metric in the (generalized) Newtonian gauge can be written as

 ds2=a2(τ){−(1+2Φ)dτ2+ωidxidτ+[(1−2Ψ)δij+hij]dxidxj}, (5)

where , and . The vector perturbation decays in an expanding universe and hence it is set to zero in CosmoLib. The tensor perturbation is gauge-invariant and its governing equations are identical in all gauges. Thus, we will only focus on the scalar perturbation equations that in CosmoLib differ from those in many other Boltzmann codes.

The linear-order relative density perturbation of is denoted by , and the linear-order velocity . Unless otherwise specified, and are all defined in Fourier space, that are functions of and the wave vector .

The radiation relative temperature fluctuation from direction seen by an observer at position is expanded as Hu/White:1997 ()

 ΔTT(x,n,τ)=∫d3k(2π)3∞∑ℓ=02∑m=−2Θγ(ℓ,m)(−i)ℓ√π4(2ℓ+1)Ymℓ(n)eik⋅x, (6)

where are the spherical harmonic functions. Note that the moments are functions of the wavenumber and the conformal time . The energy density fluctuation and velocity of photons are related to the moments and via

 δγ=Θγ(0,0); υγ=14Θγ(1,0). (7)

The neutrino moments are defined in the same way, by replacing the subscript with .

For the polarization of radiation, the Stokes parameters are expanded using the spin- harmonics Hu/White:1997 ()

 (Q±iU)(x,n,τ)=∫d3k(2π)3∞∑ℓ=02∑m=−2[E(ℓ,m)±B(ℓ,m)](−i)ℓ√π4(2ℓ+1)[±2Ymℓ(n)]eik⋅x, (8)

where and are functions of the wave vector and conformal time.

The linear-order Fourier modes are decoupled. The Fourier-space variables to be evolved are , , , , , , , , ( = , , , …, ), ( = , , , …, ), ( = , …, ). The truncations , and are adjustable integers. In CosmoLib their default values are taken to be , , , respectively. Without loss of generality we choose the azimuthal direction (the -axis direction that is used to define ) to be parallel to .

The gravitational potential can be obtained from the Einstein equations Bond/Efstathiou:1984 (); Ma/Bertschinger:1995 (); Hu/White:1997 ()

 Φ=Ψ−35k2H[ΩγΘγ(2,0)+ΩνΘν(2,0)], (9)

where we have introduced the reduced wavenumber

 kH≡kaH. (10)

Note that are functions of time. We do not treat as an independent-variable. Instead we view it as a function of the variables , and .

The close set of first-order differential equations including all the truncation schemes is:

 dΨdN =ΨN, (11) dδcdN =−kHυc+3ΨN, (12) dυcdN =−υc+kHΦ, (13) dδbdN =−kHυb+3ΨN, (14) dυbdN =−υb+kH(Φ+c2s,bδb)−κNR[υb−14Θγ(1,0)], (15) dδΛdN =−3(c2s,Λ−w)δΛ−9[c2Λ,s−(w−dw/dN3(1+w))]θΛkH−kHθΛ+3(1+w)ΨN, (16) dθΛdN (17) dΘγ(0,0)dN =−13kHΘγ(1,0)+4ΨN, (18) dΘγ(1,0)dN =kH[Θγ(0,0)−25Θγ(2,0)+4Φ]+κN[4υb−Θγ(1,0)], (19) dΘγ(2,0)dN =kH[23Θγ(1,0)−37Θγ(3,0)]−κN[910Θγ(2,0)+√610E(2,0)], (20) dΘγ(ℓ,0)dN =kH[ℓ2ℓ−1Θγ(ℓ−1,0)−ℓ+12ℓ+3Θγ(ℓ+1,0)]−κNΘγ(ℓ,0)(2<ℓ<ℓmax,γ), (21) dΘγ(ℓmax,γ,0)dN =2ℓmax,γ+12ℓmax,γ−1kHΘγ(ℓmax,γ−1,0)−(κN+ℓmax,γ+1aHτ)Θγ(ℓmax,γ,0), (22) dΘν(0,0)dN =−13kHΘν(1,0)+4ΨN, (23) dΘν(1,0)dN =kH[Θν(0,0)−25Θν(2,0)+4Φ], (24) dΘν(ℓ,0)dN =kH[ℓ2ℓ−1Θν(ℓ−1,0)−ℓ+12ℓ+3Θν(ℓ+1,0)](2≤ℓ<ℓmax,ν), (25) dΘγ(ℓmax,ν,0)dN =2ℓmax,ν+12ℓmax,ν−1kHΘν(ℓmax,ν−1,0)−ℓmax,ν+1aHτΘν(ℓmax,ν,0), (26) dE(2,0)dN =−kHK3,0,27E(3,0)−κN[25E(2,0)+√610Θγ(2,0)], (27) dE(ℓ,0)dN =kH[Kℓ,0,22ℓ−1E(ℓ−1,0)−Kℓ+1,0,22ℓ+3E(ℓ+1,0)]−κNE(ℓ,0)(2<ℓ<ℓmax,E), (28) dE(ℓmax,E,0)dN =2ℓmax,E+12ℓmax,E−1kHE(ℓmax,E−1,0)−(κN+ℓmax,E+1aHτ)E(ℓmax,E,0), (29) dΨNdN =12{(3c2s,b−1)δbΩb−δcΩc+[(3c2s,Λ−1)δΛ+9(c2s,Λ−w+dw/dN3(1+w))θΛkH]ΩΛ}−2Ψ −2(1−ϵ)Φ−k2H3(2Ψ−Φ)−(5−ϵ)ΨN+35k2H(ΩγdΘγ(2,0)dN+ΩνdΘν(2,0)dN). (30)

In the radiation and neutrino hierarchy equations (18-29) we have used the Clebsch-Gordan coefficients , which are defined as Hu/White:1997 ()

 Kℓ,m,s=⎧⎪⎨⎪⎩√(ℓ2−m2)(ℓ2−s2)ℓ, if\ ℓ≥max{|m|,|s|,1};0, otherwise. (31)

These equations are already written in the form that can be directly implemented into a generic first-order ordinary-differential-equation (ODE) solver. For a derivation of these equations, see Refs. Bond/Efstathiou:1984 (); Ma/Bertschinger:1995 (); Hu/White:1997 (); Fang/etal:2008 (). (The change of time variable from to can be done straightforwardly using .) The initial conditions can be found in Ref. Ma/Bertschinger:1995 (). For the tight-coupling approximation we follow Ref. Doran:2005a (), where the obvious typos in Eqs. (15-16) has been fixed. Since these treatments are identical to the original source, we do not repeat the discussion here. The interested readers are referred to these references for the governing equations and their derivation.

CosmoLib allows the user-input to be a phantom-crossing function, that is a function crossing the line . In this case we force and to be zero around the proximity of the phantom crossing. This is an approximation. Exact treatment requires input of at least one more degree of freedom Vikman:2005 (); Hu:2005 (); Caldwell/Doran:2005 (), which cannot be implemented in a generic code. In Ref. Fang/etal:2008 () the reader can find an alternative treatment that works better for multiple scalar field models.

Equation (30) is the key equation that guarantees the numerical stability of the code (even for isocurvature initial conditions). It is obtained by subtracting the components of the perturbed Einstein equations (pressure perturbations) from the component (density perturbations). This particular combination of the Einstein equations has been applied in the numerical code CMBquick Pitrou:2009 (); Pitrou/etal:2010 (), which assumes that dark energy is a cosmological constant and ignores the baryon sound speed. Eq. (30) is a generalized version that includes dark energy perturbations and a nonzero baryon sound speed.

We can use the energy constraint ( component of the perturbed Einstein equations, that is ) and the momentum constraint (-component of the perturbed Einstein equations, that is ) to estimate the numerical error of the code. As shown in shown in Figure 1, the relative errors are for a wide range of scales and different initial conditions.

## Iii CMB Angular Power Spectra

### iii.1 Algorithm

Optionally CosmoLib can compute the CMB angular power spectrum for each multipole by brute force, i.e., without interpolation. The angular spectrum for the temperature anisotropies is given by

 Cℓ=∫|Δkℓ|2Ps(k)dlnk, (32)

where is the temperature transfer function given by the line-of-sight integration

 Δkℓ=∫τ00S(k,τ)jℓ[k(τ0−τ)]dτ, (33)

where is the spherical Bessel function and is at redshift zero. The source can be computed from the perturbations , , , (, , , , ), and Seljak/Zaldarriaga:1996 (); Hu/White:1997 (). In Ref. Hu/White:1997 () the line-of-sight integration involves the functions , and . As shown in Ref. Seljak/Zaldarriaga:1996 (), however, the dependence on and can be eliminated by integrating by part. (We have corrected the obvious typos in eq. (12b) in Ref. Seljak/Zaldarriaga:1996 ().)

Since is evaluated numerically and it typically oscillates quickly, its sampling is time consuming. Indeed, in modern fast CMB codes – such as CAMB, CLASS, CMBEASY – the integral (32) is computed by sampling using a step size in that can be typically much larger than the oscillation period in . For instance, in Fig. 2 we show an example of for a fixed . A typical sampling scheme is shown by the red solid triangles in the upper-right panel, which zooms-in part of the figure. According to Parseval’s theorem, if is a smooth function, such sparse sampling of is enough.

However, when has local sharp features, the minimum sampling distance should be determined by the larger between the typical relative width (the width measured in ) of the oscillations in and , the typical relative width of the features in . The former is about , while the latter is model-dependent. For instance, if our goal is to sample features with width , the required sampling frequency is typically to times higher than that used for a smooth . Furthermore, as we wish to compute the ’s for each rather than interpolating it over few tens of ’s, the total time consumption will be again multiplied by a factor of . The naively estimated total time consumption is hence times more than that in the smooth- case. A final complication is due to the fact that, if all the transfer functions and the precomputed tables are to be stored, one has also to face a memory barrier that cannot be easily bypassed. For these reasons, simply increasing the and resolution in standard codes such as CAMB, CLASS or CMBEASY, will not be efficient enough for the purpose of scanning the whole parameter space.

The algorithm can be significantly improved, however, if we notice that the output from the Boltzmann code is a 2D matrix in - space. If is also a precomputed 2D matrix with the same structure, the integration (33) can be obtained by taking the inner product of the two matrices. Modern Fortran90 compilers can optimize such operation and make the computation much faster. The difficulty, however, is that the matrices for all ’s will occupy too much memory (can be up to a few tens of Giga bytes in the worst scenario). Our solution is then to only store the matrices for two neighboring ’s and update them using the recurrence relation of spherical Bessel functions.

Let us describe our strategy. We first compute two neighboring ’s by brute force. Two matrices of spherical Bessel functions and are stored in the memory for each indices. Then we compute . To do that, we update the matrix to the matrix using the recurrence relation

 jℓ−1(x)=2ℓ+1xjℓ(x)−jℓ+1(x) . (34)

Again, using and we then calculate and hence . This downward iteration is very stable for a few tens of steps, after which we need to recompute another couple of neighboring ’s and iterate downward again.

The initial neighboring ’s are calculated using precomputed 25-th order Chebyshev fitting formulas. (For the rapidly oscillating part at , the phase and amplitude of oscillations are fitted using Chebyshev polynomials.) Chebyshev fitting is slightly slower than the cubic-spline fitting used in other publicly available CMB codes, but it is more memory-efficient and more accurate – it has an accuracy of – and allows more downward iterative steps. Finally, note that the algorithm proposed here is more efficient both CPU-wise and memory-wise, enhancing the speed of -by- computation of ’s by a factor of to .

For CMB lensing we use the power spectrum approach as described in Refs. Seljak:1996 (); Zaldarriaga/Seljak:1998 ().

### iii.2 Testing the Code

The trivial comparison between CosmoLib and CAMB for smooth- models can be found in the online documentation at http://www.cita.utoronto.ca/~zqhuang/CosmoLib. Here we focus on the enhanced CMB integrator that does not assume the smoothness in . Since this feature is not available in other CMB codes, direct numerical comparison is not possible when there is very sharp features in . Thus, we need to study a model in which we have some theoretical insights. An ideal candidate is the axion monodromy inflation model, where the primordial power spectrum displays sinusoidal oscillations superimposed to a smooth power spectrum. It can be written as Flauger/etal:2010 ()

 Ps(k)=As(kk∗)ns−1[1+δnscos(ln(k/k∗)δlnk+φ)], (35)

where and are the amplitude and tilt, respectively. The parameter describes the width of the oscillations in , while gives their relative amplitude. The pivot scale is chosen to be in our computation.

We compute the CMB temperature power spectrum using the enhanced CMB integrator in CosmoLib and compare the results to the smooth- case. The relative difference between the non-smooth (for a series of ) and the smooth model is shown in Figure 3. For and we compare the results to CAMB output (both with lensing) and find good agreement. The CAMB outputs are obtained by a straightforward modification of CAMB, i.e., increasing the sampling frequency in the input file and increasing the sampling frequency in the source code. For the simple modification of CAMB fails due to insufficient memory to store the transfer functions.

For , the amplitude of oscillations in the CMB angular power spectrum (right-hand panels) is smaller than that in (left-hand panels). This suppression is generic when a 3D power spectrum is projected to a 2D one, even though in the CMB case it is further complicated by the finite duration of recombination and the recombination physics Adshead/etal:2011 (). As shown in Huang/etal:2012 (), when the frequency of oscillations is constant in , such as in eq. (35), the relative suppression is given by , as confirmed by the examples shown in Figure 3. Moreover, for , in addition to the projection effect, CMB lensing also significantly smears out the oscillations in at high . While for , the lensing smearing effect is almost negligible. See Adshead/etal:2011 (); Lewis/Challinor:2006 () for more detailed discussions about the lensing smearing effect. Finally, note that, although the oscillations in are damped, they maintain the same relative width of those of the left-hand panels, i.e., where . At low where the oscillations in space disappear in space due to the discreteness of .

This discussion shows that the enhanced CMB integrator can accurately compute the oscillations in CMB to . This does not mean, however, that the total is accurate to . The power spectrum can be systematically biased at subpercent level due to, e.g., recombination uncertainties Rubino-Mart/etal:2010 (). Understanding and eliminating these theoretical errors is important if we want to extract generic features in with accuracy. On the other hand, if we are only interested in a model predicting a specific feature in that cannot be mimicked by other effects, we can focus only on the relative difference in .

## Iv The forecast techniques

CosmoLib uses Fisher matrix analysis and MCMC method to forecast the constraints on cosmological parameters for future CMB, LSS and SN experiments. In this section we discuss the modeling of the likelihoods and the parameter estimation methods.

### iv.1 The likelihoods

#### iv.1.1 CMB simulation

Given a likelihood function , we define . For a nearly full-sky CMB experiment can be approximated by Verde/etal:2006 (); Baumann/etal:2009 ()

 χ2= ℓmax∑ℓ=ℓmin(2ℓ+1)fsky⎡⎢⎣^CBBℓCBBℓ−3+ln(CBBℓ^CBBℓ)+^CTTℓCEEℓ+^CEEℓCTTℓ−2^CTEℓCTEℓCTTℓCEEℓ−(CTEℓ)2+ln⎛⎝CTTℓCEEℓ−(CTEℓ)2^CTTℓ^CEEℓ−(^CTEℓ)2⎞⎠⎤⎥⎦ , (36)

where and are suitable cutoffs that are determined by the observed fraction of sky and the beam resolution of the experiment. In this formula, are the model-dependent theoretical angular power spectra (including the noise contributions) for the temperature, and polarizations and their cross-correlations, with . We compute the noise contribution assuming Gaussian beams. The mock data are for the fiducial model.

We use the model introduced in Verde/etal:2006 () (and later updated in Baumann/etal:2009 ()) to propagate the effect of polarization foreground residuals into the estimated uncertainties on the cosmological parameters. For simplicity, in our simulation we consider only the dominant components in the frequency bands that we are using, i.e., the synchrotron and dust signals. We assume that foreground subtraction can be done correctly down to a level of 5%. (This parameter is adjustable by the user.)

#### iv.1.2 SN simulation

For the SN simulation, we use the model given by the Dark Energy Task Force (DETF) forecast Albrecht/etal:2006 (). In this case

 χ2=∑i(mi−^miδmi)2, (37)

with going over the SN samples. More specifically, here and are the theoretical expectation and observed magnitude of the -th supernova, respectively. The uncertainty is computed by quadratically adding a peculiar velocity (a user-defined constant) to the intrinsic uncertainty in the supernova absolute magnitude (another user-specified constant).

The apparent magnitude of SN is modeled as

 m =M−μLz−μQz2+5log10(dLMpc)+25−μSδnear . (38)

The first three terms model the redshift evolution of the absolute magnitude of the supernova peak luminosity. In particular, is a free parameter with a flat prior over ; for and , Gaussian priors are applied. The widths of the Gaussian priors are user-input parameters. Finally, given that the nearby samples are likely to be a collection from many other experiments, an offset , where is unity for the nearby samples () and zero otherwise, is added to account for the systematics. For we also apply a Gaussian prior with a user-specified width. The threshold redshift is also user-defined. In conclusion, in this model there are four nuisance parameters , which we marginalized analytically.

#### iv.1.3 LSS likelihood

We model the galaxy power spectrum in redshift space as (e.g., Kaiser:1987 (); Peacock:1992 (); Peacock/Dodds:1994 ())

 Pg(k,μ;z)=(b+fμ2)2D2(z)Pm(k)exp(−k2μ2σ2r), (39)

where is the cosine of the angle between the wave vector and the line of sight, is the linear growth factor, is the linear growth rate, is the matter power spectrum today (at ) and parameterizes the effect of small scales velocity dispersion and redshift errors as explained below. The matter power spectrum is computed using Poisson’s equation, that is, .

The term accounts for the redshift distortions due to the large-scale peculiar velocity field Kaiser:1987 (), which is correlated with the matter density field. The exponential factor on the right-hand side accounts for the radial smearing due to the redshift distortions that are uncorrelated with the LSS. In particular, we consider two contributions. The first is due to the redshift uncertainty of the spectroscopic galaxy samples which is estimated to be Laureijs:2009 (). (In CosmoLib the user is allowed to change this value.) The second comes from the Doppler shift due to the virialized motion of galaxies within clusters, which typically has a pairwise velocity dispersion of the order of few hundred km/s. This can be parameterized as Peacock/Dodds:1994 (). The two contributions are quadratically added together

 σ2r=(1+z)2H2(z)(σ2z+σ2g/2), (40)

where is the Hubble parameter.

Practically, neither the redshift measurement nor the virialized motion of galaxies can be precisely modeled. In particular, the radial smearing due to peculiar velocity is not necessarily close to Gaussian. Thus, eq. (39) should not be used for wavenumbers , where the radial smearing effect is important. We introduce a UV cutoff as the smallest value between and , where is chosen such that the r.m.s. linear density fluctuation of the matter field in a sphere with radius is .

The survey volume is split into redshift bins from to , with all these parameters to be specified by the user. The number density of galaxies that can be used is , where is the fraction of galaxies with measured redshift to be specified by the user. Due to the high accuracy of the spectroscopic redshift and the width of the bins, we ignore the bin-to-bin correlations and write as

 χ2=∑k,μ,z bins(Pg,model−Pg,fiducialΔPg,fiducial)2 . (41)

As on large scales the matter density field has, to a very good approximation, Gaussian statistics and uncorrelated Fourier modes, the band-power uncertainty is given by Tegmark/etal:1998 ()

 ΔPg=[2(2π)3(2πk2dkdμ)(4πr2fskydr)]1/2(Pg+1¯n), (42)

where is the comoving distance given, for a flat FRW universe, by . The second term in the parenthesis is due to shot noise, under the assumption that the positions of the observed galaxies are generated by a random Poisson point process. In practice is not known a priori and is calibrated by galaxies themselves. The imperfect knowledge of can bias on the scale of the survey Tegmark/etal:1998 (). This is taken into account by using an IR cutoff . This is chosen such that , where is the comoving volume of the -th () redshift slice. Finally, the user has to specify the binning scheme for and . For we allow uniform binning in or in . For only uniform binning in is allowed.

In the special case where has sharp features, we must consider the smearing effect due to the fact that we are only observing a finite volume. This effect is approximated by replacing in (39) with its convolution with a Gaussian window, where the width of the Gaussian window has been chosen to be

 σW=√2ln22π(4π3)1/3kmin≃0.302kmin. (43)

In such a way, the real-space representation of the window, if cut off at its half-height, contains the same volume as that of the redshift bin. The fact that is smaller than allows us to neglect the overlap between window functions centered around neighboring values of .

### iv.2 Parameter Estimation

#### iv.2.1 Fisher Matrix Analysis

In general, the likelihood can be written as

 lnL(p;pfid)=−12[d(p)−d(pfid)]TC−1(p;pfid)[d(p)−d(pfid)], (44)

where is the data vector, the fiducial parameter vector, the parameter vector for which one wants to evaluate the likelihood, and the covariance matrix.

The fisher matrix for , (two components of ) is then

 Fij≡−∂2lnL∂pi∂pj∣∣∣p=pfid=∂d(p)∂piC−1(p;pfid)∂d(p)∂pj∣∣∣p=pfid, (45)

where the partial derivatives can be evaluated numerically:

 ∂d(p)∂pi=12Δpi[d(p1,p2,…,pi+Δpi,…,pn)−d(p1,p2,…,pi−Δpi,…,pn)]. (46)

The stepsize is initially supplied by the user, and then adjusted by the software in such a way that the variation of is of when is varied by . By doing this, we have assumed that the likelihood is approximately Gaussian in the proximity of where the variation of is . If the likelihood is highly non-Gaussian, Fisher matrix analysis does not give reliable estimations of the error bars of parameters. In this case, one should use the MCMC method to fully explore the structure of the likelihood.

#### iv.2.2 MCMC method

CosmoLib has an independent MCMC engine using the Metropolis-Hastings algorithm. The traditional approach is to define the proposal density (the probability of walking from to in the parameter space) using a roughly estimated covariance matrix

 Q(x;x′)∝exp[−12(x−x′)TC−1e(x−x′)]. (47)

Convergence can be achieved quickly if is close the posterior covariance matrix of .

However, sometimes we need to treat models where the likelihood periodically depends on some phase parameters. Here we take the axion monodromy inflation model for example. The likelihood is a periodic function of the axion phase parameter ,

 L(P,φ)=L(P,φ+2π), (48)

where we have used to represent the collection of other parameters. If is not well constrained, we will obtain multi-branches in the posterior, i.e., for a fixed value of and a chosen confidence level, the allowed values of locate in well separated regions in the parameter space.

Intuitively the separated regions can be more efficiently explored by restricting the range of to one period and proposing with wrap-around or, in a more rigorous language, by using a periodic proposal density. For and , we use

 Q(P,φ;P′,φ′)∝∞∑n=−∞exp[−12(x−x′n)TC−1e(x−x′n)], (49)

where . The estimation of covariance matrix, , is practically computed with a trial run that is terminated before the multi-branches of the posterior are explored by the random walk.

We find that the periodic proposal density (49) leads to significant improvement of the convergence. For the axion monodromy model, it takes about 5-10 times longer to achieve convergence using (47) than using (49).

The output chains in CosmoLib have the same format as those in CosmoMC Lewis/Bridle:2002 (). The chains can hence be directly analyzed using the GetDist tool in CosmoMC. For completeness, an independent postprocessing tool is supplied in CosmoLib to analyze and visualize the marginalized posterior of parameters. In the online documentation the reader can find the instructions on how to use this tool.

#### iv.2.3 Oscillations in the Current CMB Data?

Recently a hint of the axion monodromy cosine oscillations (see eq. (35)) in WMAP-7yr Komatsu/etal:2011 (); Larson/etal:2011 () and ACT CMB data Dunkley/etal:2011 () has been claimed in Ref. Aich/etal:2011 (). Ref. Meerburg/etal:2012 () confirms the finding that can be significantly improved in some regions of parameter space where oscillations in the primordial power spectrum are assumed. In this section we use CosmoLib to constrain the axion monodromy model with the same data sets. We find that when the CMB power spectrum is accurately computed and rigorous statistical method is used, there is no detectable axion monodromy oscillations in the CMB data.

In Refs. Aich/etal:2011 () the authors used their modified CAMB to compute the CMB power spectrum. As discussed in previous sections, such a modification is not trivial for . Since the best-fit found in Ref. Aich/etal:2011 () is small – (derived from Table III of Ref. Aich/etal:2011 () and equation 51 in Ref. Huang/etal:2012 ()), it is necessary to exam the numerical accuracy of the modified CAMB used in Aich/etal:2011 (). For , the modulation period in is . In the CMB power spectrum one should see same modulation period in , i.e., . Thus, from to there should be about oscillations in . However, in Figure 5 of Ref. Aich/etal:2011 () the number of oscillations in between and are much more than . This implies that the “modulations” in shown in Ref. Aich/etal:2011 () may just be numerical noises. In Figure 4 we show the CMB temperature angular power spectrum computed with CosmoLib, where the parameters are chosen to be close to the ones used in Figure 5 of Ref. Aich/etal:2011 (). Qualitative difference can be seen between the two figures. The spectrum computed using CosmoLib presents clear modulations that agrees with the value, while the modified CAMB used in Ref. Aich/etal:2011 () failed to produce the expected modulations.