Observational Constraints on the Dark Energy Equation of State

TL;DR

This paper uses observational data and statistical methods to constrain the evolution of dark energy's equation of state, finding no significant deviations from the standard cosmological model, Lambda-CDM.

Contribution

It introduces a nearly model-independent approach combining multiple datasets to constrain dark energy evolution using Markov Chain Monte Carlo sampling.

Findings

No strong evidence for deviations from Lambda-CDM

Constraints consistent with a cosmological constant

Multi-dataset analysis enhances robustness

Abstract

Late-time cosmic acceleration is one of the most interesting unsolved puzzles in modern cosmology. The explanation most accepted nowadays, dark energy, raises questions about its own nature, e.g. what exactly is dark energy, and implications to the observations, e.g. how to handle fine tuning problem and coincidence problem. Hence, dark energy evolution through cosmic history, together with its equation of state, are subjects of research in many current experiments. In this dissertation, using Markov Chain Monte Carlo sampling, we try to constrain the evolution of the dark energy equation of state in a nearly model-independent approach by combining different datasets coming from observations of baryon acoustic oscillations, cosmic chronometers, cosmic microwave background anisotropies and type Ia supernovae. We found no strong evidence that could indicate deviations from $\Lambda$CDM model, which is the standard model in cosmology accepted today.

Figures

Tables

Table 1. Table 2.1: BAO data for galaxy surveys.

Survey	Parameter	Effective redshift	Measurement
6DF	$r_s/D_V$	0.106	$0.336\pm 0.015$
SDSS-MGS	$D_V/r_s$	0.15	$4.4656\pm 0.1681$
BOSS-LOWZ	$D_V/r_s$	0.32	$8.250\pm 0.170$
WiggleZ	$A$	0.44	$0.474\pm 0.034$
BOSS-CMASS	$D_V/r_s$	0.57	$13.773\pm 0.134$
WiggleZ	$A$	0.60	$0.442\pm 0.020$
WiggleZ	$A$	0.73	$0.424\pm 0.021$

Table 2. Table 2.2: Inverse covariance matrix for WiggleZ data. Measurements were performed in the overlapping redshifts quoted and, as the matrix is symmetric, only the upper diagonal is displayed. Extracted from [ 81 ] .

Redshift slice
	1040.3	-807.5	336.8
		3720.3	-1551.9
			2914.9

Table 3. Table 2.3: Cosmic chronometers dataset. Hubble parameter is in [ k ​ m ​ s − 1 ​ M ​ p ​ c − 1 ] delimited-[] 𝑘 𝑚 superscript 𝑠 1 𝑀 𝑝 superscript 𝑐 1 [km\,s^{-1}\,Mpc^{-1}] . Extracted from [ 100 ] .

$z$	$H\left(z\right)$	$z$	$H\left(z\right)$
$0.07$	$69.0\pm 19.6$	$0.4783$	$80.9\pm 9.0$
$0.09$	$69.0\pm 12.0$	$0.48$	$97.0\pm 62.0$
$0.12$	$68.6\pm 26.2$	$0.5929$	$104.0\pm 13.0$
$0.17$	$83.0\pm 8.0$	$0.6797$	$92.0\pm 8.0$
$0.1791$	$75.0\pm 4.0$	$0.7812$	$105.0\pm 12.0$
$0.1993$	$75.0\pm 5.0$	$0.8754$	$125.0\pm 17.0$
$0.2$	$72.9\pm 29.6$	$0.88$	$90.0\pm 40.0$
$0.27$	$77.0\pm 14.0$	$0.9$	$117.0\pm 23.0$
$0.28$	$88.8\pm 36.6$	$1.037$	$154.0\pm 20.0$
$0.3519$	$83.0\pm 14.0$	$1.3$	$168.0\pm 17.0$
$0.3802$	$83.0\pm 13.5$	$1.363$	$160.0\pm 33.6$
$0.4$	$95.0\pm 17.0$	$1.43$	$177.0\pm 18.0$
$0.4004$	$77.0\pm 10.2$	$1.53$	$140.0\pm 14.0$
$0.4247$	$87.1\pm 11.2$	$1.75$	$202.0\pm 40.0$
$0.4497$	$92.8\pm 12.9$	$1.965$	$186.5\pm 50.4$
$0.47$	$80\pm 50$

Table 4. Table 4.1: List of parameters and definitions used in this work. First block represents cosmological parameters, second block is a set of “nuisance” parameters and in the third block some derived parameters are displayed. In prior range column, square brackets mean flat priors while parentheses represent Gaussian priors. There is an additional Gaussian prior which is a linear combination of two parameters: A kSZ + 1.6 × A 143 tSZ = 9.5 ± 3.0 superscript 𝐴 kSZ 1.6 subscript superscript 𝐴 tSZ 143 plus-or-minus 9.5 3.0 A^{\mathrm{kSZ}}+1.6\times A^{\mathrm{tSZ}}_{143}=9.5\pm 3.0 . Central and width columns represent an estimated value (taken from CosmoMC) for each parameter close to the peak of the posterior and a range (relative to the central value) where initial positions for walkers are randomly sampled, respectively. A detailed description of “nuisance” parameters from Planck 2015 is given in reference [ 95 ] .

Parameter	Prior range	Central	Width	Definition
${\omega }_b\equiv {\Omega }_b{h}^2$	[0.005, 0.1]	0.0221	0.0001	Baryon density today
${\omega }_c\equiv {\Omega }_c{h}^2$	[0.001, 0.99]	0.120	0.001	Cold dark matter density today
$n_s$	[0.8, 1.2]	0.965	0.004	Scalar spectrum index
ln$\left({10}^{10}{A}_s\right)$	[2, 4]	3.100	0.001	Log power of scalar amplitude
$H_0$	[20, 100]	67.3	1.0	Current expansion rate in $\text{km}{\text{s}}^{-1}{\text{Mpc}}^{-1}$
$\tau$	[0.01, 0.8]	0.08	0.01	Thomson optical depth due to reionization
$w_1\left(z=0.00\right)$	[-10, 8]	-1.0	1.0	1st dark energy equation of state parameter
$w_2\left(z=0.25\right)$	[-10, 8]	-1.0	1.0	2nd dark energy equation of state parameter
$w_3\left(z=0.50\right)$	[-10, 8]	-1.0	1.0	3rd dark energy equation of state parameter
$w_4\left(z=0.85\right)$	[-10, 8]	-1.0	1.0	4th dark energy equation of state parameter
$w_5\left(z=1.25\right)$	[-10, 8]	-1.0	1.0	5th dark energy equation of state parameter
$w_6\left(z=2.00\right)$	[-10, 8]	-1.0	1.0	6th dark energy equation of state parameter
$A_{217}^{\mathrm{CIB}}$	[0, 200]	65	10	CIB contamination at $\ell =3000$ (217-GHz)
${\xi }^{\mathrm{tSZ}\times \mathrm{CIB}}$	[0, 1]	0	0.1	SZ$\times$CIB cross-correlation
$A_{143}^{\mathrm{tSZ}}$	[0, 10]	5	2	tSZ contamination at 143 GHz
$A_{100}^{\mathrm{PS}}$	[0, 400]	255	24	Point source contribution in 100$\times$100
$A_{143}^{\mathrm{PS}}$	[0, 400]	40	10	Point source contribution in 143$\times$143
$A_{143\times 217}^{\mathrm{PS}}$	[0, 400]	40	12	Point source contribution in 143$\times$217
$A_{217}^{\mathrm{PS}}$	[0, 400]	100	13	Point source contribution in 217$\times$217
$A^{\mathrm{kSZ}}$	[0, 10]	0	3	kSZ contamination
$A_{100}^{\mathrm{dustTT}}$	[0, 50]	7	2	Dust contamination at $\ell =200$ in 100$\times$100
	$\left(7\pm 2\right)$
$A_{143}^{\mathrm{dustTT}}$	[0, 50]	9	2	Dust contamination at $\ell =200$ in 143$\times$143
	$\left(9\pm 2\right)$
$A_{143\times 217}^{\mathrm{dustTT}}$	[0, 100]	17	4	Dust contamination at $\ell =200$ in 143$\times$217
	$\left(21.0\pm 8.5\right)$
$A_{217}^{\mathrm{dustTT}}$	[0, 400]	80	15	Dust contamination at $\ell =200$ in 217$\times$217
	$\left(80\pm 20\right)$
$c_{100}$	[0, 3]	0.999	0.001	Power spectrum calibration for the 100 GHz
	$\left(0.999\pm 0.001\right)$
$c_{217}$	[0, 3]	0.995	0.002	Power spectrum calibration for the 217 GHz
	$\left(0.995\pm 0.002\right)$
$y_{\mathrm{cal}}$	[0.9, 1.1]	1.000	0.002	Absolute map calibration for Planck
	$\left(1.0000\pm 0.0025\right)$
$\alpha$	[0.0, 0.3]	0.135	0.020	Time stretching coefficient for supernovae
$\beta$	[1.5, 4.0]	3.0	0.3	Supernovae color coefficient
$M_B^1$	[-25, -15]	-19.05	0.15	Supernovae absolute magnitude
${\Delta }_M$	[-0.3, 0.3]	-0.05	0.10	Correction for supernovae absolute magnitude
${\Omega }_m$	-	-	-	Matter density today
${\Omega }_{\Lambda }$	-	-	-	Dark energy density today
$\mathrm{bar}{\chi }^2$	-	-	-	Mean chi square

Table 5. Table 4.2: Marginalized estimates from BAO+CC+SN datasets. Best fit, 68 % percent 68 68\% and 95 % percent 95 95\% confidence levels are displayed. In these estimates, w i subscript 𝑤 𝑖 w_{i} ’s are correlated. For comparison, Planck constraints for Λ Λ \Lambda CDM model at 1 ​ σ 1 𝜎 1\sigma [ 49 ] are also displayed.

Parameter	Best fit	$68\%$ limits	$95\%$ limits	Planck (TT+lowP)
${\Omega }_b{h}^2$	0.0257	${0.0234}_{-0.0048}^{+0.0040}$	${0.0234}_{-0.0085}^{+0.0089}$	$0.02222\pm 0.00023$
${\Omega }_c{h}^2$	0.106	$0.125\pm 0.015$	${0.125}_{-0.030}^{+0.030}$	$0.1197\pm 0.0022$
$H_0$	68.2	$69.8\pm 2.3$	${69.8}_{-4.5}^{+4.6}$	$67.31\pm 0.96$
$w_1$	-1.80	$-1.81\pm 0.63$	$-{1.8}_{-1.2}^{+1.2}$	-
$w_2$	-0.80	$-0.88\pm 0.20$	$-{0.88}_{-0.40}^{+0.37}$	-
$w_3$	-1.63	$-{1.75}_{-0.39}^{+0.47}$	$-{1.75}_{-0.89}^{+0.85}$	-
$w_4$	1.11		-	-
$w_5$	-1.6		-	-
$w_6$	2.8	-	-	-
${\Omega }_m$	0.282	$0.305\pm 0.027$	${0.305}_{-0.051}^{+0.055}$	$0.315\pm 0.013$
${\Omega }_{\Lambda }$	0.718	${0.695}_{-0.026}^{+0.028}$	${0.695}_{-0.055}^{+0.051}$	$0.685\pm 0.013$
$\mathrm{bar}{\chi }^2$	486.2	${499.2}_{-6.2}^{+4.0}$	${499}_{-9.6}^{+11}$	-

Table 6. Table 4.3: Marginalized estimates from CMB dataset. Best fit, 68 % percent 68 68\% and 95 % percent 95 95\% confidence levels are displayed. In these estimates, PCA method has already been applied in order to obtain uncorrelated w i subscript 𝑤 𝑖 w_{i} ’s. For comparison, Planck constraints for Λ Λ \Lambda CDM model at 1 ​ σ 1 𝜎 1\sigma [ 49 ] are also displayed.

Parameter	Best fit	$68\%$ limits	$95\%$ limits	Planck (TT+lowP)
${\Omega }_b{h}^2$	0.02204	$0.02207\pm 0.00015$	${0.02207}_{-0.00028}^{+0.00028}$	$0.02222\pm 0.00023$
${\Omega }_c{h}^2$	0.1211	${0.1199}_{-0.0014}^{+0.0012}$	${0.1199}_{-0.0025}^{+0.0027}$	$0.1197\pm 0.0022$
$n_s$	0.9648	$0.9641\pm 0.0040$	${0.9641}_{-0.0082}^{+0.0078}$	$0.9655\pm 0.0062$
ln$\left({10}^{10}{A}_s\right)$	3.0998	${3.0989}_{-0.0018}^{+0.0023}$	${3.0989}_{-0.0058}^{+0.0048}$	$3.089\pm 0.036$
$H_0$	65.5	$66.6\pm 2.4$	${66.6}_{-5.2}^{+4.9}$	$67.31\pm 0.96$
$\tau$	0.0822	${0.0820}_{-0.0029}^{+0.0022}$	${0.0820}_{-0.0047}^{+0.0055}$	$0.078\pm 0.019$
$w_1$	-0.89	$-{0.97}_{-0.40}^{+0.32}$	$-{0.97}_{-0.72}^{+0.81}$	-
$w_2$	-1.05	$-{0.97}_{-0.13}^{+0.12}$	$-{0.97}_{-0.25}^{+0.27}$	-
$w_3$	-1.05	$-0.99\pm 0.26$	$-{0.99}_{-0.52}^{+0.51}$	-
$w_4$	-1.21	$-{1.17}_{-0.24}^{+0.37}$	$-{1.17}_{-0.63}^{+0.58}$	-
$w_5$	-0.94	$-{0.83}_{-0.34}^{+0.39}$	$-{0.83}_{-0.69}^{+0.67}$	-
$w_6$	-1.10	$-{1.21}_{-0.52}^{+0.72}$	$-{1.21}_{-1.30}^{+1.10}$	-
${\Omega }_m$	0.333	${0.321}_{-0.023}^{+0.019}$	${0.321}_{-0.041}^{+0.052}$	$0.315\pm 0.013$
${\Omega }_{\Lambda }$	0.665	${0.677}_{-0.019}^{+0.023}$	${0.677}_{-0.053}^{+0.041}$	$0.685\pm 0.013$
$\mathrm{bar}{\chi }^2$	11268.7	${11288.4}_{-8.1}^{+5.5}$	${11290}_{-13}^{+14}$	-

Table 7. Table 4.4: Marginalized estimates from BAO+CC+CMB+SN datasets. Best fit, 68 % percent 68 68\% and 95 % percent 95 95\% confidence levels are displayed. In these estimates, PCA method has already been applied in order to obtain uncorrelated w i subscript 𝑤 𝑖 w_{i} ’s. For comparison, Planck constraints for Λ Λ \Lambda CDM model at 1 ​ σ 1 𝜎 1\sigma [ 49 ] are also displayed.

Parameter	Best fit	$68\%$ limits	$95\%$ limits	Planck (TT+lowP)
${\Omega }_b{h}^2$	0.02215	$0.02219\pm 0.00020$	${0.02219}_{-0.00038}^{+0.00041}$	$0.02222\pm 0.00023$
${\Omega }_c{h}^2$	0.1199	$0.1199\pm 0.0019$	${0.1199}_{-0.0039}^{+0.0039}$	$0.1197\pm 0.0022$
$n_s$	0.9639	$0.9645\pm 0.0055$	${0.965}_{-0.011}^{+0.011}$	$0.9655\pm 0.0062$
ln$\left({10}^{10}{A}_s\right)$	3.091	${3.094}_{-0.018}^{+0.024}$	${3.094}_{-0.053}^{+0.044}$	$3.089\pm 0.036$
$H_0$	70.2	$69.0\pm 1.7$	${69.0}_{-3.3}^{+3.2}$	$67.31\pm 0.96$
$\tau$	0.079	${0.079}_{-0.010}^{+0.013}$	${0.079}_{-0.028}^{+0.024}$	$0.078\pm 0.019$
$w_1$	-0.94	$-1.12\pm 0.20$	$-{1.12}_{-0.40}^{+0.39}$	-
$w_2$	-0.958	$-0.993\pm 0.072$	$-{0.99}_{-0.14}^{+0.14}$	-
$w_3$	-1.42	$-1.33\pm 0.25$	$-{1.33}_{-0.48}^{+0.50}$	-
$w_4$	-0.75	$-{0.78}_{-0.17}^{+0.25}$	$-{0.78}_{-0.46}^{+0.41}$	-
$w_5$	-1.55	$-{1.55}_{-0.49}^{+0.71}$	$-{1.5}_{-1.2}^{+1.1}$	-
$w_6$	-1.13	$-{1.7}_{-1.0}^{+1.6}$	$-{1.7}_{-2.8}^{+2.3}$	-
${\Omega }_m$	0.288	$0.299\pm 0.015$	${0.299}_{-0.029}^{+0.030}$	$0.315\pm 0.013$
${\Omega }_{\Lambda }$	0.710	$0.699\pm 0.015$	${0.699}_{-0.030}^{+0.029}$	$0.685\pm 0.013$
$\mathrm{bar}{\chi }^2$	11756.3	${11779.1}_{-8.6}^{+6.1}$	${11780}_{-14}^{+15}$	-

Equations

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=-d{t}^2+a^2(t)d{\sigma }^2,$$

$$d{\sigma }^2={\gamma }_{ij}d{x}^{i}d{x}^j=\frac{d{r}^2}{1-K{r}^2}+r^2(d{\theta }^2+{\sin }^2\theta d{\varphi }^2).$$

$$\eta \equiv {\int }_0^t\frac{d{t}^{\prime }}{a(t^{\prime })},$$

$${\Gamma }_{\nu \lambda }^{\mu }=\frac{1}{2}{g}^{\mu \alpha }(g_{\alpha \nu ,\lambda }+g_{\alpha \lambda ,\nu }-g_{\nu \lambda ,\alpha }).$$

$$R_{\mu \nu }={\Gamma }_{\mu \nu ,\alpha }^{\alpha }-{\Gamma }_{\mu \alpha ,\nu }^{\alpha }+{\Gamma }_{\mu \nu }^{\alpha }{\Gamma }_{\alpha \beta }^{\beta }-{\Gamma }_{\mu \beta }^{\alpha }{\Gamma }_{\alpha \nu }^{\beta }$$

$$R=g^{\mu \nu }{R}_{\mu \nu }.$$

$$G_{\mu \nu }\equiv R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R,$$

$$G_{\nu }^{\mu }=8\pi G{T}_{\nu }^{\mu },$$

$${\Gamma }_{ij}^0=a^{2}H{\gamma }_{ij},{\Gamma }_{0j}^i={\Gamma }_{j0}^i=H{\delta }_j^i,$$

$${\Gamma }_{11}^1=\frac{Kr}{1-K{r}^2},{\Gamma }_{22}^1=-r(1-K{r}^2),{\Gamma }_{33}^1=-r(1-K{r}^2){\sin }^2\theta ,$$

$${\Gamma }_{33}^2=-\sin \theta \cos \theta ,{\Gamma }_{12}^2={\Gamma }_{21}^2={\Gamma }_{13}^3={\Gamma }_{31}^3=\frac{1}{r},{\Gamma }_{23}^3={\Gamma }_{32}^3=\cot \theta ,$$

$$R_{00}=-3(H^2+\dot{H}),R_{0i}=R_{i0}=0,R_{ij}=a^2(3H^2+\dot{H}+2K/a^2){\gamma }_{ij},$$

$$R=6(2H^2+\dot{H}+K/a^2),$$

$$G_0^0=-3(H^2+K/a^2),G_i^0=G_0^i=0,G_j^i=-(3H^2+2\dot{H}+K/a^2){\delta }_j^i.$$

$$T_{\nu }^{\mu }=(\rho +P)u^{\mu }{u}_{\nu }+P{\delta }_{\nu }^{\mu },$$

$$H^2=\frac{8\pi G}{3}\rho -\frac{K}{a^2}$$

$$3H^2+2\dot{H}=-8\pi GP-\frac{K}{a^2}.$$

$$\frac{\ddot{a}}{a}=-\frac{4\pi G}{3}(\rho +3P).$$

$$\dot{\rho }+3H(\rho +P)=0,$$

$${\nabla }_{\mu }{G}_{\nu }^{\mu }\equiv \frac{\partial G_{\nu }^{\mu }}{\partial x^{\mu }}+{\Gamma }_{\alpha \mu }^{\mu }{G}_{\nu }^{\alpha }-{\Gamma }_{\nu \mu }^{\alpha }{G}_{\alpha }^{\mu }=0,$$

$${\Omega }_M+{\Omega }_K=1,\text{where}{\Omega }_M\equiv \frac{8\pi G\rho }{3H^2},{\Omega }_K\equiv -\frac{K}{(aH{)}^2}.$$

$${\Omega }_r^{(0)}=\frac{8\pi G{\rho }_r^{(0)}}{3H_0^2},{\Omega }_m^{(0)}=\frac{8\pi G{\rho }_m^{(0)}}{3H_0^2},{\Omega }_{DE}^{(0)}=\frac{8\pi G{\rho }_{DE}^{(0)}}{3H_0^2},{\Omega }_K^{(0)}=-\frac{K}{(a_0{H}_0{)}^2},$$

$$dt=\pm a(t)\frac{dr}{\sqrt{1-K{r}^2}}.$$

$${\int }_{t_1}^{t_0}\frac{dt}{a(t)}={\int }_0^{r_1}\frac{dr}{\sqrt{1-K{r}^2}}.$$

$$\frac{\delta t_1}{a(t_1)}=\frac{\delta t_0}{a(t_0)}.$$

$$\frac{{\nu }_0}{{\nu }_1}=\frac{a(t_1)}{a(t_0)}.$$

$$z\equiv \frac{{\lambda }_0}{{\lambda }_1}-1=\frac{a(t_0)}{a(t_1)}-1.$$

$$\dot{\mathbf{r}}=H\mathbf{r}+a\dot{\mathbf{x}}.$$

$$v\equiv \dot{\mathbf{r}}\cdot \mathbf{r}/r=Hr+{\mathbf{v}}_p\cdot \mathbf{r}/r.$$

$$v\simeq H_{0}r,$$

$$H_0=100hkmse{c}^{-1}Mp{c}^{-1},$$

$${\Omega }_M\equiv {\Omega }_m+{\Omega }_r+{\Omega }_{DE}=1-{\Omega }_K,$$

Full text

Universidade de São Paulo

Instituto de Física

Vínculos Observacionais sobre a Equação de Estado da Energia Escura

Williams Jonata Miranda Ribeiro

Orientador: Prof. Dr. Elcio Abdalla

Dissertação de mestrado apresentada ao Instituto de Física para a obtenção do título de Mestre em Ciências

Banca Examinadora:

Prof. Dr. Elcio Abdalla (IF-USP)

Prof. Dr. Gastão Cesar Bierrenbach Lima Neto (IAG-USP)

Prof. Dr. Carlos Alexandre Wuensche de Souza (INPE)

São Paulo

2018

See pages 1 of fichacatalografica.pdf

University of São Paulo

Institute of Physics

Observational Constraints on the Dark Energy Equation of State

Williams Jonata Miranda Ribeiro

Advisor: Prof. Dr. Elcio Abdalla

Dissertation presented to the Institute of Physics of the University of São Paulo in partial fulfillment of the requirements for the degree of Master of Science

Examining Committee:

Prof. Dr. Elcio Abdalla (IF-USP)

Prof. Dr. Gastão Cesar Bierrenbach Lima Neto (IAG-USP)

Prof. Dr. Carlos Alexandre Wuensche de Souza (INPE)

São Paulo

2018

Acknowledgments

First, I would like to thank God for renewing my strength in the most difficult moments I had been through in my master’s. Without His power, I would not have gone so far.

Also, I would like to thank my advisor Dr. Elcio Abdalla for accepting me to be part of his research group, even though in the beginning I had no prior knowledge about general relativity and cosmology at all.

I am very thankful to my whole family for all the support during these years of study and specially to my father Francisco Junior for his encouragement, which was essential for my personal growth.

I could not forget to mention the friends I made during these years in São Paulo (and the ones I already knew) who have assisted me in many ways in my journey, helping me either with physics or with conversations and laughs that kept me going: Bonifácio Lima, Christianne Moraes, Kariny Bonjorno, Naim Comar, Renan Boschetti, Riis Rhavia, Rodrigo Pinheiro and Thiago Silveira.

I would like to specially thank the postdoc André Costa for all the assistance, tips, discussions and even time invested on my work that made me to progress so much. A gratefully thanks to my friend Caroline Guandalin for the interesting discussions about cosmology and for lending me her computer to perform the simulations for my work (thank you very much!!!). Also, thanks to the colleagues Alessandro Marins and Leonardo Werneck for the programming tips that helped me to improve my code and my computer skills.

My sincere acknowledgement to the science agency CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico - Brasil. This dissertation would not be possible without its financial support.

Abstract

Late-time cosmic acceleration is one of the most interesting unsolved puzzles in modern cosmology. The explanation most accepted nowadays, dark energy, raises questions about its own nature, e.g. what exactly is dark energy, and implications to the observations, e.g. how to handle fine tuning problem and coincidence problem. Hence, dark energy evolution through cosmic history, together with its equation of state, are subjects of research in many current experiments. In this dissertation, using Markov Chain Monte Carlo sampling, we try to constrain the evolution of the dark energy equation of state in a nearly model-independent approach by combining different datasets coming from observations of baryon acoustic oscillations, cosmic chronometers, cosmic microwave background anisotropies and type Ia supernovae. We found no strong evidence that could indicate deviations from $\Lambda$CDM model, which is the standard model in cosmology accepted today.

Keywords: Cosmology. Dark Energy. Equation of State. Cosmological Parameters. Monte Carlo.

Resumo

A aceleração cósmica atual é um dos mais interessantes enigmas não resolvidos da cosmologia moderna. A explicação mais aceita hoje em dia, a energia escura, levanta questões acerca de sua própria natureza, como o que é exatamente a energia escura, e as implicações para as observações, por exemplo como lidar com o problema do ajuste fino e o problema da coincidência. Por isso, a evolução da energia escura durante a história cósmica, juntamente com sua equação de estado, são objetos de pesquisa em inúmeros experimentos atuais. Nesta dissertação, usando amostragem por cadeias de Markov de Monte Carlo, tentamos restringir a evolução da equação de estado da energia escura em uma abordagem quase independente de modelo ao combinar diferentes conjuntos de dados provenientes de observações de oscilações acústicas de bárions, cronômetros cósmicos, anisotropias da radiação cósmica de fundo e supernovas tipo Ia. Não encontramos evidências fortes que pudessem indicar desvios do modelo $\Lambda$CDM, o qual é o modelo padrão aceito hoje na cosmologia.

Palavras-Chave: Cosmologia. Energia Escura. Equação de Estado. Parâmetros Cosmológicos. Monte Carlo.

Contents

1. Introduction

2. 1 Overview

3. 1.1 Homogeneous and isotropic Universe

4. 1.2 Cosmological redshift, Hubble law, critical density and equation of state

5. 1.3 Cosmic distances

6. 1.3.1 Luminosity distance

7. 1.3.2 Angular diameter distance

8. 1.4 Cosmic components

9. 1.4.1 Relativistic species

10. 1.4.2 Non-relativistic matter

11. 1.4.3 Dark energy

12. 1.5 Radiation-matter equality, Universe dominated by dark energy and Hubble parameter evolution

13. 2 Datasets description

14. 2.1 Baryon acoustic oscillations

15. 2.1.1 Overview

16. 2.1.2 Correlation function and power spectrum

17. 2.1.3 Perturbation evolution in the primordial plasma

18. 2.1.4 BAO acoustic scale and relative BAO distance

19. 2.2 Cosmic microwave background radiation

20. 2.2.1 Overview

21. 2.2.2 Anisotropy power spectrum and Einstein-Boltzmann equations

22. 2.2.3 Optical depth

23. 2.2.4 Primordial fluctuations: scalar spectral index and amplitude perturbation

24. 2.3 Cosmic chronometers

25. 2.3.1 Overview

26. 2.3.2 Differential age method from galaxy surveys

27. 2.4 Supernovae

28. 2.4.1 Overview

29. 2.4.2 Magnitude system

30. 2.4.3 Luminosity distance, distance modulus and B-band generalized model

31. 3 Statistical analysis

32. 3.1 Likelihood function and Bayes’ theorem

33. 3.2 Parameter estimation

34. 3.3 Markov Chain Monte Carlo method: emcee

35. 3.4 Burn-in and convergence

36. 3.5 Principal components analysis

37. 4 Results

38. 4.1 BAO+Clocks+Supernovae

39. 4.2 CMB

40. 4.3 BAO+Clocks+CMB+Supernovae

41. 4.4 Comparison of datasets, results and models

42. Conclusions

Introduction

In 1998, a remarkable discovery set up a new era in modern cosmology: through observations of type Ia supernovae, it was realized that the Universe is going through a phase of accelerated expansion [1, 2]. As our common sense leads us to the idea that the expansion is decelerated due to the effects of gravitation (which was the general belief until 1997), many theories were (and are still being) developed in order to explain this acceleration. Among them, the currently accepted theory is that there exists some kind of “exotic fluid” with negative pressure, named dark energy, that contributes to approximately $70\%$ of our Universe and that is the responsible for this accelerated phase. Other explanations include models of dark energy with repulsive gravity [3], modifications to general relativity [4], models of vacuum decaying [5, 6], interacting dark matter-dark energy models [7, 8, 9, 10], and many others as explored in [11]. Nowadays, besides type Ia supernovae, many kinds of independent observations confirm the accelerated phase: cosmic microwave background anisotropies [12], large-scale structure [13], baryon acoustic oscillations [14], cosmic chronometers [15], etc. Thus, in this work we are going to assume dark energy as the correct explanation for the current Universe dynamics.

One of the ways of describing dark energy is considering it as a barotropic fluid (a fluid whose density depends only on its pressure) with equation of state $w$, which is defined as the ratio between the pressure and the density of the fluid. The model that better matches the observations is named $\Lambda$CDM model, in which $w=-1$ during all cosmic history. Despite its success in explaining many independent measurements, it is still not the final history for the fact that it does not solve problems like the cosmological constant problem [16] (also called fine tuning problem) and the coincidence problem [17]. Other well known models in literature are the $w$CDM model, which assigns a value different of $-1$ but still constant for $w$ and the CPL model (Chevallier-Polarski-Linder), which considers $w$ as a dynamical dark energy equation of state that varies with time [18, 19]. There are also parametrizations that try to capture the evolution of $w$ either by performing Taylor expansions or by using something called kink parametrization, which takes into account rapid transitions of dark energy equation of state (something not allowed by a conventional Taylor expansion) [20, 21]. Several constraints in these models can be found in [22, 23, 24, 25].

There is no fundamental theory that can predict whether dark energy equation of state has a static value or evolves with time. In addition to the fact that this quantity cannot be measured directly, strong degeneracy with cosmological parameters makes its evolution hard to predict, specially for high redshifts [26]. Also, although one can write the dark energy equation of state as a function of the background expansion (and consequently as a function of luminosity distance), accurate measurements of the latter coming from cosmic chronometers and type Ia supernovae cannot give good constrains on $w$ because in this approach there is an additional dependence on the derivative of the background expansion and because we need estimates of the present density of pressureless matter that only comes from large-scale structure methods, which actually measures the present density of clustered matter (both kinds of matter are not necessarily the same) [11, 27].

The importance of dark energy equation of state lies in the fact that its evolution determines the current state of the Universe, e.g. an accelerated expansion requires $w<-1/3$ [27], its eventual fate, e.g. a phantom dark energy defined by $w<-1$ leads to a finite-time singularity called big-rip singularity [28], and the evolution of dark energy itself via continuity equation. Nowadays, the main collaboration for studying the nature of dark energy is the Dark Energy Survey (DES) research project [29], but there are other projects being designed for studying the dark sector such as Euclid [30], BINGO [31, 32], J-PAS [33] and DESI [34].

Instead of testing specific models, a more robust way of extracting information from the evolution of $w$ is using different observational datasets in order to reconstruct the evolution of the dark energy equation of state in a model-independent approach (as much as possible) and check if the datasets lead to an agreement with $\Lambda$CDM model or not. To this end, we are going to use a method known as principal components analysis (PCA) [35, 36, 37] along with MCMC sampling applied to several datasets in order to obtain estimates for dark energy equation of state in a few redshifts.

Chapter 1 Overview

The initial chapter defines the general terms and concepts used in modern cosmology by using the framework of the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, which agrees with observations that display the Universe as homogeneous and isotropic at large scales. We also construct the tools in such a way that the need for dark matter and dark energy becomes clear for explaining many independent measurements.

1.1 Homogeneous and isotropic Universe

We consider the most general line element that describes a homogeneous and isotropic 4-dimensional spacetime, the FLRW metric, given by (using units of $c=1$) [27]

[TABLE]

where $g_{\mu\nu}$ is the metric tensor, $a(t)$ is the scale factor which depends on the cosmic time $t$ and $d\sigma^{2}$ represents the spatial part of the FLRW metric with constant curvature K, given by

[TABLE]

For curvature, $K=-1$, $K=0$ and $K=+1$ represent open, flat and closed geometries, respectively. We are also using spherical coordinates $(x^{1},x^{2},x^{3})=(r,\theta,\phi)$ and Einstein’s summation convention which states that terms with equal upper and lower indices are summed over.

Besides cosmic time $t$, we also define another useful “time” named conformal time $\eta$ and defined by

[TABLE]

which is the amount of time a photon would take to travel from where we are to the furthest observable distance provided the universe ceased expanding. It is interesting sometimes to work with the conformal time instead of cosmic time because it may simplify the evolution equations. From now on, a prime (′) indicates derivative with respect to the conformal time.

The Universe evolution dynamics is described by the Einstein equations. The first step to obtain them is to evaluate the connections from the metric tensor $g_{\mu\nu}$,

[TABLE]

We also define the Ricci tensor and the Ricci scalar (or scalar curvature), respectively, as

[TABLE]

and

[TABLE]

From (1.5) and (1.6) we can evaluate the Einstein tensor

[TABLE]

from where we obtain the cosmological dynamics by solving the Einstein equations

[TABLE]

where $T^{\mu}_{\hphantom{\mu}\nu}$ represents the energy-momentum tensor of matter components.

From the FLRW metric, the only non-vanishing connection coefficients are

[TABLE]

where $H\equiv\dot{a}/a$, a dot representing a derivative with respect to cosmic time $t$. The quantity $H$, called Hubble parameter, describes the rate of expansion of the Universe.

Using (1.5), (1.6) and (1.7), we can evaluate the components from the Ricci tensor, Ricci scalar and the components from the Einstein tensor, which are given by, respectively,

[TABLE]

[TABLE]

Considering a FLRW space-time imposes a constraint on the energy-momentum tensor, which should be represented by a perfect fluid in the form

[TABLE]

where $u^{\mu}=(-1,0,0,0)$ is the four-velocity of the fluid in comoving coordinates, with $\rho$ and $P$ representing the energy density and pressure, respectively. The $(00)$ and $(ii)$ components of the Einstein equations (1.8) can be obtained using (1.11) and (1.12), from which one obtains the Friedmann equations

[TABLE]

and

[TABLE]

The importance of these equations lies in the fact that they relate the evolution dynamics of the Universe to its own constitution and space geometry. Combining these equations by eliminating the curvature term yields

[TABLE]

Multiplying (1.13) by $a^{2}$, differentiating and using (1.15), one yields

[TABLE]

which is the continuity equation for the cosmological components. This equation may also be derived from relations satisfied by the Einstein tensor, the Bianchi identities

[TABLE]

$\nabla_{\mu}$ denoting the covariant derivative, which leads to the conservation of the energy-momentum tensor $\nabla_{\mu}T^{\mu}_{\hphantom{\mu}\nu}=0$ that implies in (1.16).

We may rewrite equation (1.13) in a way that makes it easier to study the evolution of cosmological parameters,

[TABLE]

Equation (1.18) defines a constraint among cosmological parameters, such that the amount of each component in the Universe must respect it at any moment of the cosmic evolution [38].

Referring to values today with the superscript $(0)$, let us define the density parameters $\Omega$. They are given by relativistic particles, non-relativistic matter, dark energy and curvature, respectively, as

[TABLE]

in which we can identify relativistic particles as electromagnetic radiation ($\gamma$) and neutrinos ($\nu$), non-relativistic matter as baryons ($b$) and cold dark matter ($c$), and dark energy (which we shall introduce later) as a cosmological constant ($\Lambda$) or a more general fluid ($X$). These cosmological parameters are of great interest, because if we can accurately measure or infer them from observations then it might be possible to understand the current state of the Universe, its initial moments and predict its future.

1.2 Cosmological redshift, Hubble law, critical density and equation of state

At first glance, Friedmann equations give us no information about the evolution of the scale factor in the FLRW metric: is it increasing, decreasing or is it a constant? This information comes from observations of distant galaxies that are seen to recede from us based on the displacement of their spectral lines to the red part of the spectrum [39]. Thereby, here we relate these shifts on the spectral lines to the scale factor.

In the FLRW metric, let us define our position at the center of coordinates, without loss of generality, and consider a light ray coming toward us in the radial direction. This light ray obeys the equation $ds^{2}=0$, such that (1.1) implies

[TABLE]

Considering a light ray coming towards us from a distant source at position $r_{1}$ at cosmic time $t_{1}$, as the distance decreases the cosmic time increases and therefore we must choose the minus sign in (1.20). Defining our position as $r=0$ at time $t_{0}$, one finds

[TABLE]

Taking (1.21) in a differential form and using the fact that the radial coordinate $r_{1}$ is time-independent, we see that the interval $\delta t_{1}$ between departure of subsequent light signals is related to the interval $\delta t_{0}$ between arrivals of these light signals by

[TABLE]

Considering these “signals” as subsequent wave fronts, the emitted frequency is $\nu_{1}=1/\delta t_{1}$ and the observed frequency is $\nu_{0}=1/\delta t_{0}$, therefore

[TABLE]

If the scale factor is increasing, then $\nu_{0}<\nu_{1}\Rightarrow\lambda_{0}>\lambda_{1}$, showing that the wavelength of photons stretches during its travel to us due to the Universe expansion. This is called redshift and is exactly what we observe today. Conventionally, the redshift $z$ is related to the scale factor by

[TABLE]

In an expanding Universe, the relation between the physical distance $\bm{r}$ and the comoving distance $\bm{x}$ to an object is given by $\bm{r}=a(t)\bm{x}$. Taking the derivative of the equation $\bm{r}=a(t)\bm{x}$ with respect to cosmic time $t$, yields

[TABLE]

The fist term in (1.25) represents the Hubble flow, while the second term represents movements with respect to the local Hubble flow and is called peculiar velocity $\bm{v}_{p}$. In the special case in which an object moves exclusively due to the Hubble flow (that is, the object moves only because of the expansion), the comoving distance does not change and the second term vanishes. In the general case, the projection of the object velocity in the radial direction is given by

[TABLE]

In most cases the peculiar velocity of galaxies does not exceed $10^{6}\,m/s$ [27]. Therefore, if the second term in (1.26) is negligible in comparison with the first one, then

[TABLE]

in which the replacement of $H$ by the present value $H_{0}$ is justified in the limit of low redshifts ($z\ll 1$). Equation (1.27) represents the Hubble law, discovered in 1929 by Edwin Hubble when he plotted the recessional velocity $v$ of galaxies versus distance $r$ to them [39]. Although his data were scarce and noisy, Hubble correctly concluded that the Universe was expanding. The Hubble parameter today $H_{0}$ is usually written in the form

[TABLE]

in which the parameter $h$ describes the uncertainty on the value of $H_{0}$ and we shall name it reduced Hubble parameter. The original Hubble measurements are displayed in Fig. 1.1.

An important quantity of interest in cosmology is the critical density, which is defined as the total energy density necessary for the Universe to be spatially flat [38]. From (1.18), one may rewrite it as

[TABLE]

which is related to the total equivalent mass density by $\Omega_{M}=8\pi G\rho/(3H^{2})$. Thus, for any spatially flat cosmological model, one requires $\Omega_{M}=1$ and is common to write the corresponding mass density as a critical density, which is given by

[TABLE]

The critical density today can be written as

[TABLE]

where we have used $G=6.67\times 10^{-8}cm^{3}g^{-1}sec^{-2}$. Therefore, if the total density of the Universe is greater than $\rho_{c}^{(0)}$ then the Universe is spatially closed and if the total density is smaller than $\rho_{c}^{(0)}$ then we live in a spatially open Universe.

Another useful definition is the equation of state of a fluid (in cosmology, modeled as a dilute gas), which is the ratio of its pressure $P$ to its density $\rho$ and is given by (in natural units)

[TABLE]

such that, as we shall see from thermodynamics, relativistic particles and non-relativistic matter have equation of state $w=1/3$ and $w=0$, respectively. This is an important quantity that governs the evolution of different species in the Universe and, when related to dark energy, it tells us if the Universe is going through an accelerated expansion or not and if it will expand forever or collapse in the near future [27].

Several parametrizations for dark energy equation of state have been proposed so far. The main ones are based on the expansion

[TABLE]

where $w_{n}$ are parameters to be adjusted using observational datasets and $x_{n}$ are phenomenological dependencies on redshift (or scale factor). The cases most considered in the literature are

[TABLE]

In case (i), if $w_{0}=-1$ we have a $\Lambda$CDM model and if $w_{0}\neq-1$ we have a $w$CDM model. In case (ii), the model in which $n\leq 1$ was studied in [40, 41]. Also at linear order, case (iii) is called CPL model [18, 19] and case (iv) was firstly introduced in [42]. The evolution of dark energy equation of state is the main objective of study in this work.

1.3 Cosmic distances

In order to relate observations to quantities of interest (e.g. matter density, curvature,…), it is important to define distances in cosmology. However, defining distance measures in the Universe can be confusing because light received by us from a distant galaxy was emitted when the Universe was younger and, therefore, the physical distance has increased since the time of emission because of the finite velocity $c$ of light. For that matter, defining cosmic distances in terms of the comoving distance may be a better approach.

As usual, we will work in a FLRW spacetime (1.1) by setting $r=\sin\chi$ ($K=+1$), $r=\chi$ ($K=0$) and $r=\sinh\chi$ ($K=-1$) in (1.2). Therefore, the 3-dimensional space line-element can be expressed as

[TABLE]

where

[TABLE]

The function (1.39) can be written in an unified way

[TABLE]

in which one might recover the flat case by taking the limit $K\rightarrow-0$.

Now, let us show that the comoving distance $d_{c}$ is a function of the Hubble parameter and, thereby, the cosmological parameters. The trajectory of a photon traveling along the $\chi$ direction follows the geodesic equation $ds^{2}=-c^{2}dt^{2}+a^{2}(t)d\chi^{2}=0$ (we have now restored $c$ in in the line element). Similarly to the derivation of the cosmological redshift, let us consider the case where light is emitted at time $t=t_{1}$ in position $\chi=\chi_{1}$ (redshift $z$), reaching an observer at time $t=t_{0}$ at position $\chi=0$ (redshift $z=0$). Integrating the geodesic equation, one recovers the comoving distance

[TABLE]

Combining the definition of the Hubble parameter with (1.24), one obtains the relation $dt=-dz/[H(1+z)]$ and therefore the comoving distance is given by

[TABLE]

1.3.1 Luminosity distance

We define now the cosmic distance most important in the study of type Ia supernovae, which was the first observation that determined the accelerated expansion rate of the Universe. This distance indicator is called luminosity distance ($d_{L}$) and is defined by

[TABLE]

where $L_{s}$ is the absolute luminosity of a source (in SI units, measured in $W=Js^{-1}$) and $\mathcal{F}$ is the observed flux (in SI units, measured in $Wm^{-2}$).

As one can see, the observed luminosity $L_{0}$ (detected at $\chi=0$ and $z=0$) is different from the absolute luminosity $L_{s}$ from the source (emitted at $\chi$ with the redshift $z$). By definition, the flux is $\mathcal{F}\equiv L_{0}/S$, where $S=4\pi(a_{0}f_{K}(\chi))^{2}$ is the area of a sphere at $z=0$ considering the space curvature. Then the luminosity distance (1.43) is given by

[TABLE]

In order to rewrite the ratio $L_{s}/L_{0}$ in terms of redshift, let us define the energy of a pulse emitted at the time-interval $\Delta t_{1}$ as $\Delta E_{1}$, which implies that the absolute luminosity is $L_{s}=\Delta E_{1}/\Delta t_{1}$. On the same way, the observed luminosity is $L_{0}=\Delta E_{0}/\Delta t_{0}$, where $\Delta E_{0}$ is the energy of light detected at the time-interval $\Delta t_{0}$. As the energy of a photon is directly proportional to its frequency $\nu$ (inversely proportional to $\lambda$) we have that $\Delta E_{1}/\Delta E_{0}=\lambda_{0}/\lambda_{1}=1+z$. Moreover, as the light signal propagates with the same velocity $c=\lambda/\Delta t$ through its path, then $\lambda_{1}/\Delta t_{1}=\lambda_{0}/\Delta t_{0}$, where $\lambda_{1}$ and $\lambda_{0}$ represent the wavelength of light at the emission and detection points, respectively. This implies that $\Delta t_{0}/\Delta t_{1}=\lambda_{0}/\lambda_{1}=1+z$. Hence we find that

[TABLE]

and luminosity distance reduces to

[TABLE]

Combining (1.40) and (1.42) with (1.46), then luminosity distance can finally be expressed as

[TABLE]

where $\Omega_{K}^{(0)}=-Kc^{2}/(a_{0}H_{0})^{2}$, as given by (1.19). As one can see, our expression for luminosity distance is strongly dependent on the cosmology adopted because there is a dependence with the cosmological parameters ($E(z)$ depends on the $\Omega_{i}^{(0)}$’s, as we shall see) that affects the evolution of $d_{L}$.

1.3.2 Angular diameter distance

Another useful definition of cosmic distance is the angular diameter distance, which is very popular in CMB and BAO data analysis and is defined by

[TABLE]

where $\Delta\theta$ is the angle subtended by an object of actual size $\Delta x$ orthogonal to the line of sight.

In order to express (1.48) in terms of the FLRW metric, suppose we have two radial null geodesics (light paths) meeting at the observer at time $t_{0}$ with angular separation $\Delta\theta$, which have been emitted from a source of size $\Delta x$ at time $t_{1}$ at a comoving coordinate $\chi$ (we also assume that $\phi=$ constant along the photon paths). From the angular part of the FLRW metric we have

[TABLE]

Therefore, we may identify the angular diameter distance as

[TABLE]

which also displays strong dependence on the cosmology adopted just like luminosity distance. The evolution of $d_{L}$ and $d_{A}$ for some cosmologies is shown in Fig. 1.2.

It is interesting to notice that combining (1.47) and (1.50) yields the relation

[TABLE]

which is called reciprocity or duality or Etherington relation [43], valid for all cosmological models based on Riemannian geometry as long as flux is conserved. This relation has been tested by many authors [44] and it seems to be in good agreement with data.

1.4 Cosmic components

In this section we shall study the species that are supposed to form the Universe today. These are generally classified into relativistic particles, non-relativistic matter and dark energy. There is an additional component, presumably a scalar field, that is supposed to have dominated the Universe during the inflationary era, but this component is not of our immediate interest here. Therefore, we will focus on the main components and analyze its thermodynamic properties of interest for the thermal history of the Universe.

During cosmic ages, many processes have happened so rapidly that the equilibrium was achieved for most of the time, with different kinds of particles sharing the same temperature. We wish to measure quantities like density and pressure in terms of this equilibrium temperature. For that matter, it is necessary to introduce the occupation number or distribution function or phase space occupancy of a species, which counts the number of particles in a given region in phase space around position $\bm{x}$ and momentum $\bm{p}$ [45]. Considering a particle with momentum $\bm{p}$ and mass $m$, from special relativity the energy of this particle is $E=\sqrt{p^{2}+m^{2}}$, where $p\equiv\left|\bm{p}\right|$. Thus, the distribution function of some species in equilibrium at temperature $T$ is given by

[TABLE]

where $\mu$ is the chemical potential for each of the species. The plus and minus signs are defined depending on which kind of particles we are dealing with: plus representing Fermi-Dirac statistics and minus the Bose-Einstein statistics.

Equation (1.52) defines the general case where the distribution may depend on the position of the species, but in an homogeneous Universe we can relax this assumption and write $f=f(|\bm{p}|)\equiv f(p)$. By Heisenberg’s principle, no particle can be localized in a region of phase space smaller than $(2\pi\hbar)^{3}$, therefore it is a fundamental size which defines the number of phase space elements in the volume $d^{3}x\,d^{3}p$ as $d^{3}x\,d^{3}p/(2\pi\hbar)^{3}$. By defining $g_{\ast}$ as the number of internal degrees of freedom (e.g. spin states), the energy density $\rho$ and pressure $P$ are given by

[TABLE]

In (1.53) and (1.4), the integrals are not evaluated over $d^{3}x$ because the energy density and pressure are defined as quantities per unit volume. In the first equality of (1.4) we have used the fact that the pressure per unit number density of particles is given by $pv/3$ ($v$ is the particle velocity), and in the second equality we used the relation $v=p/E$ (using units of $c=1$) from special relativity (which can be derived combining the equations $E=mc^{2}/\sqrt{1-v^{2}/c^{2}}$ and $p=mv/\sqrt{1-v^{2}/c^{2}}$). For the final expressions, we have adopted $\hbar=1$. In what follows, we particularize (1.53) and (1.4) for different particle species.

1.4.1 Relativistic species

The limit of relativistic species is equivalent to consider $T\gg m$ in (1.53) and (1.4), i.e. $m\rightarrow 0$. For non-degenerate particles ($T\gg\mu$) one obtains

[TABLE]

where we have used $\int_{0}^{\infty}dx\,x^{3}/(e^{x}-1)=\pi^{4}/15$ and $\int_{0}^{\infty}dx\,x^{3}/(e^{x}+1)=7\pi^{4}/120$. Equation (1.56) shows that for relativistic particles without degeneracies the equation of state is $w=1/3$.

The main relativistic particle that one might want to study is the photon, which is a boson. For this kind of particle, $g_{\ast}=2$ to account for the two spin states. Also, we could have solved (1.53) and (1.4) for photons simply by considering the chemical potential as zero, which is expected theoretically because in the early Universe photon number is not conserved (e.g. electrons and positrons can annihilate to produce photons) [45]. Moreover, one might safely ignore the chemical potential because precision measurements on the CMB spectrum (which shall be described in the next chapter) constrain the chemical potential to $\mu/T<9\times 10^{-5}$ [46]. This leads the CMB photon density to be

[TABLE]

The COBE satellite showed that the CMB spectrum is very close to the spectrum of a black-body radiation with a temperature of $T_{\gamma}=2.725\pm 0.002\,K$ [47], as one can see in Fig. 1.3, which has been improved with WMAP satellite to the value $T_{\gamma}=2.72548\pm 0.00057\,K$ [48]. Sticking to the value measured by COBE satellite and using the conversion factor $1\,K^{4}=1.279\times 10^{-35}g\,cm^{-3}$, the energy density of CMB photons today is $\rho_{\gamma}^{(0)}=4.641\times 10^{-34}g\,cm^{-3}$. Therefore, the photon density parameter is

[TABLE]

in which we have used (1.31). The reduced Hubble parameter has a current value close to $h=0.67$ [49], which implies $\Omega_{\gamma}^{(0)}\simeq 5.5\times 10^{-5}$. The amount of radiation in the present Universe is very low because, as we shall see later, radiation density scales as $a^{-4}$ so that it dilutes faster than other components. The scaling of radiation density leads to a scaling in the temperature of CMB photons of $T\propto a^{-1}\,(=1+z)$, which evolves this way even after the photons went out of equilibrium with matter [50].

Another relativistic particle to be added to the cosmic inventory is the neutrino, which has a very small mass. Neutrinos are fermionic particles with zero chemical potential and there are three types of species in the standard model (electron neutrino $\nu_{e}$, muon neutrino $\nu_{\mu}$ and tau neutrino $\nu_{\tau}$). As each species has one spin degree of freedom and we have to account for anti-neutrinos, from (1.55) the neutrino energy density is given by

[TABLE]

where $N_{eff}$ is the effective number of neutrino species and $T_{\nu}$ is the background neutrino temperature, predicted to have the value $T_{\nu}=(4/11)^{1/3}\,T_{\gamma}\simeq 1.945\,K$ (this relation comes from the conservation of entropy before and after the annihilation of electrons and positrons [27]). Unlike cosmic microwave background, which is measured with great accuracy, cosmic neutrino background has only indirect evidences due to the fact that neutrinos interact very weakly [51].

Although the effective number of neutrinos in the standard model is $N_{eff}=3$, the presence of relativistic degrees of freedom changes this value slightly to the value $N_{eff}=3.04$ [52]. As the temperature of neutrinos and photons are linked via the relation $T_{\nu}/T_{\gamma}=(4/11)^{1/3}$, by using (1.57) and (1.59) one shows that neutrino density and photon density are linked via $\rho_{\nu}=N_{eff}(7/8)(4/11)^{4/3}\rho_{\gamma}$. Hence the radiation density parameter today, which accounts for photons and relativistic neutrinos, is given by

[TABLE]

with $\Omega_{\gamma}^{(0)}$ given by (1.58). Considering again $h=0.67$ and $N_{eff}=3.04$, one obtains $\Omega_{r}^{(0)}\simeq 9.3\times 10^{-5}$, which shows that today radiation is very diluted.

1.4.2 Non-relativistic matter

Now we consider the case of non-relativistic particles ($T\ll m$). From (1.53) and (1.4), one gets

[TABLE]

which are valid for both bosons and fermions. Equation (1.62) is in the form of an equation of state, similar to (1.32), and thus one infers that $w\simeq 0$, as expected, because $T/m\ll 1$. The case $w=0$ is also called dust. The above expressions show that non-relativistic matter is not simply described only in terms of temperature (as relativistic particles do), and therefore we need to measure the density of non-relativistic particles (baryons and dark matter) directly from observations [27].

Baryons compose the visible matter we can observe in the Universe: mostly protons and neutrons (although we also include electrons in this classification even knowing they are leptons, because protons and neutrons are so much more massive than electrons that virtually all the mass in atoms is in the baryons). In literature there are four main techniques used to measure baryon density: observation of baryons from gas in group of galaxies [53]; baryon counting by spectra analysis of distant quasars [54]; sensitivity to baryon density on the amount of light elements produced during the Big Bang nucleosynthesis epoch [55]; effects of change in baryon density on the CMB anisotropy spectrum [56]. The tightest constraint on baryon density comes from recent measurements of Planck satellite [49], which constrains baryon density to the value

[TABLE]

at $68\%$ confidence level. Using $h=0.67$, $\Omega_{b}^{(0)}\simeq 0.05$ for the central value in (1.63).

All of these different techniques are in quite good agreement [57], ranging baryon density from $2-5\%$ with respect to the critical density. However, total matter density in the Universe is higher than this and, therefore, there must be some kind of matter which is not baryonic to account for this discrepancy.

In 1933, Zwicky [58] studied the Coma cluster by measuring radial velocities of galaxies (which by the time were called nebulae). He found a surprising result: galaxy velocities were highly dispersed to be explained by the “visible” stellar material, which indicated that the cluster density was much higher than density inferred from luminous matter only. His conclusion was that there was some kind of “missing matter” to account for the matter unobserved, which he dubbed dark matter. Another evidence comes from 1970’s when rotation curves of spiral galaxies showed that there was indeed some “missing matter” in these objects by the observation that rotation curves were flat to very large radii from the galactic center, when in fact the behavior expected was a declining one [59]. This fact is evidenced in Fig. 1.4.

These evidences and recent studies on the field (e.g. gravitational lensing [61], large-scale structure [62], Bullet Cluster [63]) also support the existence of dark matter, suggesting it as another kind of non-relativistic matter that differs from baryons because it does not interact electromagnetically but only gravitationally. In addition, the most probable kind of dark matter is cold dark matter (CDM), which considers dark matter as non-relativistic due to the fact that if dark matter was hot (relativistic), then structure formation would not have achieved the level of development observed today because relativistic particles hardly cluster, streaming out of overdense regions easily, and baryons would not have enough potential wells to fall into and form structures.

There are a few candidates to the origin of dark matter, which are divided into astrophysical candidates (e.g. white dwarfs, black holes and neutron stars) and particle candidates (e.g. axions and Weakly Interacting Massive Particles). We also cannot rule out new physics as non-standard gravitational effects unpredicted by general relativity, although this is unlikely due to the huge amount of independent observations that point to the existence of dark matter. The best constraint we have to date comes again from Planck satellite [49], which constrains dark matter density to the value

[TABLE]

at $68\%$ confidence level. By using $h=0.67$, $\Omega_{c}^{(0)}\simeq 0.26$ for the central value in (1.64), which displays dark matter as a dominating component over baryonic matter.

1.4.3 Dark energy

Observing the previous values calculated for the density of cosmic components, it is interesting to note that baryons and dark matter account for approximately $30\%$ of the total energy budget of the Universe, while photons and neutrinos are negligible. Moreover, CMB measurements [64] along with inflationary scenarios [65] put tight constraints on the contemporary Universe curvature, $|\Omega_{K}^{(0)}|\lesssim 0.01$, displaying our Universe as flat and showing that the total density is close to the critical density. As one can see, there is a lack of some cosmic component to account for the remaining $70\%$ and this is where dark energy comes in: it gives an explanation for the lacking matter in the Universe and for the late-time cosmic acceleration measured years ago by observing type Ia supernovae [1, 2]. Dark energy evidence due to type Ia supernovae can be clearly seen by analyzing Fig. 1.5 and Fig. 1.6.

Planck constraint on the dark energy density, for a cosmological constant (which shall be explained latter), is given by [49]

[TABLE]

at $68\%$ confidence level, which agrees with the missing $70\%$.

1.5 Radiation-matter equality, Universe dominated by dark energy and Hubble parameter evolution

Let us consider a Universe dominated, at some point of its evolution, by a single fluid (component) with equation of state $w$. If $w$ is a constant, one can analytically solve for the evolution of the density $\rho$ and scale factor $a$ for a flat Universe. Using (1.13) and (1.16), one obtains the solutions (with $K=0$)

[TABLE]

where $t_{i}$ is a constant. For instance, if one considers radiation as the dominant component then the total density scales as $a^{-4}$ and the scale factor as $(t-t_{i})^{1/2}$. On the other hand, in a matter-dominated era the pressure is negligible and thus the total density scales as $a^{-3}$ and the scale factor as $(t-t_{i})^{2/3}$.

Consider the primitive Universe, with dominance of radiation (density $\rho_{r}$ and pressure $P_{r}=\rho_{r}/3$) and non-relativistic matter (density $\rho_{m}$ and pressure $P_{m}=0$). As these components scale as $\rho_{r}\propto a^{-4}$ and $\rho_{m}\propto a^{-3}$, we have

[TABLE]

The redshift of the transition from radiation-dominance to matter-dominance, which corresponds to the radiation-matter equality ($\rho_{r}=\rho_{m}$), is

[TABLE]

where $\Omega_{r}^{(0)}$ is given by (1.58) and (1.60). If one considers the effective number of neutrino species $N_{eff}=3.04$, we get

[TABLE]

and, using Planck result $\Omega_{m}^{(0)}h^{2}=0.1415\pm 0.0019$ [49], one obtains $z_{eq}\simeq 3389$ for the central value, which displays the radiation-matter equality as a previous event to the decoupling of photons, as we shall see.

As mentioned in the introduction, many observations show that the Universe is going through an accelerated phase nowadays. In order to explain this framework from the point of view of a component dominance, which is supposed to drive this acceleration, we require $\ddot{a}>0$ in equation (1.15), which leads to

[TABLE]

with $\rho$ assumed to be positive. Therefore, a component that follows this requirement and dominates the Universe today could be a simple explanation for the accelerated expansion. Such an unusual component with negative pressure is what we call dark energy. This is our motivation here to consider dark energy as a good candidate to explain the late-time acceleration and thereby to study its equation of state evolution.

An interesting case is the so-called cosmological constant, defined when $w=-1$, which implies from (1.16) that $\rho$ is a constant during all epochs. From the first Friedmann equation, one can infer that $H$ is a constant in a flat Universe and therefore the scale factor evolves exponentially as $a\propto\exp(Ht)$. Actually, this model might give an explanation for dark energy because no one knows if the acceleration will last forever or end in a near future [27].

From the discussion above, let us consider dark energy as a fluid with equation of state $w_{DE}=P_{DE}/\rho_{DE}$, which satisfies the equation

[TABLE]

By integrating it using the relation $dt=-dz/[H(1+z)]$, one obtains the way dark energy density evolves, given by

[TABLE]

such that the integral form is kept because we are considering the general case in which $w$ may be time-dependent.

Expanding the first Friedmann equation in the components,

[TABLE]

as already mentioned, one can rewrite (1.74) as a constraint equation for the cosmological parameters today, given by

[TABLE]

If we combine (1.18) and (1.74), we get one of the main equations used in cosmology,

[TABLE]

which displays the Hubble parameter evolution as a function of redshift, the cosmological parameters today and dark energy equation of state.

Chapter 2 Datasets description

In order to constrain the evolution of the dark energy equation of state, in this section we shall describe the observational datasets we are going to use in Chapter 4 and the theoretical background required to model these observations. From now on, even whether not mentioned, we will stick to a flat geometry as testified by many observations [12, 56, 67].

2.1 Baryon acoustic oscillations

2.1.1 Overview

The observable Universe is not composed simply by randomly-positioned galaxies all over space, but there are many structures (e.g. clusters, super-clusters, voids) arranged in a very specific way (also called cosmic web). One of the main objectives of cosmology is to understand how these structures were developed and came to be what we observe today.

In the early Universe, temperature was so high and the Universe so dense that one can consider matter and radiation coupled as a single fluid, named photon-baryon fluid. Such a high temperature prevented formation of atoms because electrons were constantly being “casted out” of hydrogen atoms, and also prevented photons of propagating freely through space because the mean free path was too low due to scattering with electrons. Moreover, the photon-baryon fluid is also a plasma.

The standard vision of gravity is that regions in space more dense than its surroundings tend to attract more matter than underdense regions, becoming even more dense, but this is not exactly what happens in the early Universe. Instead, small overdensities of baryonic matter over the whole space, coming from the inflationary era, try to grow due to gravitational instability, but they are soon counterbalanced by the pressure imbalance of the plasma caused by the photons and heat trapped in it. This pressure is so high that overdensities just oscillate on its amplitude instead of growing by gravity, with spherical sound waves being emitted from each overdensity and propagating through the Universe, similar to sound waves propagating in a fluid.

This situation goes on until approximately 380.000 years after the Big Bang, at the point the Universe has expanded enough such that matter-radiation decoupling occurs, allowing electrons and nuclei to form neutral atoms. The plasma pressure is then released and the sound waves are frozen in, forming an overdense spherical shell of characteristic comoving radius of approximately $150$ Mpc, with center on the initial overdensity. Nevertheless, these oscillations in the primordial plasma, named baryon acoustic oscillations (BAO), end up but leave an imprint due to the spherical shell mainly in the galaxy distribution, causing two galaxies to be more likely separated by this characteristic scale[14].

A scheme on the evolution of perturbations in the primordial era for some redshifts is displayed in Fig. 2.1, where we plot the density of each component times the square of the characteristic comoving radius versus the comoving radius, such that the mass in the overdensity is the area under the curve. Initial point-like overdensities for each component (dark matter, gas containing nuclei and electrons, photons and neutrinos) evolve, having in mind that these perturbations occur in many places over the Universe and their effects sum linearly because fluctuations are very small. Due to the fact that initial perturbations are adiabatic, in the beginning all species are perturbed almost by the same amount (species at $z=82507$) because the energy perturbation of neutrinos and photons is $4/3$ bigger than dark matter and gas. The first component to decouple from the others and free-stream is the neutrino because it interacts very weakly and, as a relativistic particle, it does not cluster. As dark matter interacts only gravitationally, it stands still because has no intrinsic motion for the fact that it is cold dark matter. The gas and photons are coupled to each other as already described and the spherical sound wave starts to propagate with origin in the initial overdensity (species at $z=6824$). These component perturbations grow almost nothing during radiation era: gas and photon perturbations do not grow because of the coupling and dark matter perturbations due to Meszaros effect [68].

As the Universe evolves, the photon-gas perturbation continues to propagate, the neutrinos have completely streamed away and dark matter perturbation remains still but gets larger due to gravitational forces that tend to attract more material (species at $z=1440$). In the meantime, the Universe is cooling to the point that photons decouple from gas and start to stream away just as neutrinos did, while baryon perturbation stops propagating and starts to grow freely (species at $z=848$). What remains is a dark matter perturbation at the original center and a baryon perturbation in a shell of radius 150 Mpc from the center (species at $z=478$), which start to attract each other and grow quickly due to combined gravitational forces (species at $z=79$). Finally, we are left with similar profiles for dark matter and gas that enhance the acoustic peak (species at $z=10$). As galaxies form in regions that are overdense, now becomes intuitive the reason why galaxies are more probable to be found separated by this characteristic scale.

In order to have a good idea about the way one obtains the acoustic scale from galaxy distribution, we need to look at Fig. 2.2. Representing galaxies by dots, in the left hand panel the characteristic scale is clearly visible because there are many voids and the galaxy rings are dense. On the other way, the right hand panel displays a more realistic situation where galaxy rings are less dense and the voids are filled with other galaxy rings, which visually masks the characteristic scale. Therefore, the only way of recovering the characteristic BAO radius is using statistics, which requires mapping enormous volumes of the sky to detect the BAO signal that is very weak at large scales but yet measurable from galaxy redshift surveys[69].

2.1.2 Correlation function and power spectrum

The tool that allows one to measure the BAO peak is the correlation function, while cosmological information is better extracted from the power spectrum. In order to understand both, one needs to define the density contrast [71]

[TABLE]

which shows us deviations from the background density $\rho_{0}$. Now we define the correlation function as the spatial average

[TABLE]

made over the whole volume by taking many points ($a$ and $b$) in pairs to perform the evaluation. In this average, we are using the same separation distance $r_{ab}$ considered at various locations, which is called sample average. The correlation function tells us if there are any correlations between overdensities separated by the distance $r_{ab}$ that could indicate some mechanism creating dependence on the distance, quantifying the excess clustering on a given scale. Also, when the correlation function $\xi(\bm{r})$ depends only on the separation $\bm{r}$ and not on the locations $\bm{r}_{a}$ and $\bm{r}_{b}$, we have statistical homogeneity (the system has the same statistical properties everywhere). In this case, the correlation function as a sample average can be written as

[TABLE]

where the integration is made over all possible positions.

In practice, for obtaining the correlation function one has to compare the real catalog of galaxies with a mock galaxy distribution generated randomly. We can therefore obtain an estimate of the correlation function as

[TABLE]

where DD represents the number of galaxies at distance $r$ counted by an observer at a galaxy in the real catalog, while DR also represents the number of galaxies at distance $r$ but now counted in the random catalog. In the literature, there are other estimators for the correlation function which have been compared rigorously in [72]. The first BAO detection was made by Eisenstein et al. [14] while analyzing a spectroscopic sample of 46748 luminous red galaxies (LRG) measured from the Sloan Digital Sky Survey (SDSS), where a peak corresponding to the BAO scale was found at $100\,h^{-1}$ Mpc, as we can see in Fig. 2.3.

Sometimes it is useful to derive some quantities in the Fourier space. For that matter, we will apply the Fourier transform and its inverse, respectively, with the convention

[TABLE]

When studying perturbation theory in cosmology, one finds that it is interesting to work with the density contrast in Fourier space because this decouples the evolution equations, making the equations to evolve independently [27]. Since the average of a perturbed variable is zero, the first important quantity related to perturbations is a quadratic function. Therefore, we define the power spectrum as

[TABLE]

definition that applies to any perturbed variable in Fourier space (as we will see, temperature perturbations in the CMB follow the same pattern).

The Fourier transform of the density contrast of a density field $\delta(\bm{x})$ is given by

[TABLE]

and, by using the definition of the power spectrum in the form

[TABLE]

where the normalization $A$ is identified with the volume V considered, we have

[TABLE]

With the change $\bm{r}=\bm{x}-\bm{y}$, it follows that

[TABLE]

with $\xi(\bm{r})$ defined by (2.3), which tells us that the power spectrum is the Fourier transform of the correlation function (Fourier pair). Therefore, it follows conversely that

[TABLE]

If the correlation function does not depend on the direction $\bm{r}$ but only on the modulus $r=|\bm{r}|$ and, accordingly, the power spectrum depends only on $k=|\bm{k}|$, then the system has spatial isotropy and the power spectrum simplifies to

[TABLE]

Because the correlation function and the power spectrum are related by a Fourier transform via (2.11), features displayed in one have consequences on the other. As a simple example, if $\xi(\bm{r})$ has the form of a delta function centered at a scale $r_{*}$, the result is an oscillation pattern in the power spectrum. Although the correlation function of a real survey is not a delta but a small bump, this feature is still preserved and the acoustic peak in $\xi(\bm{r})$ induces oscillations in the power spectrum, which characterizes the baryon acoustic oscillations, as shown in Fig. 2.4.

2.1.3 Perturbation evolution in the primordial plasma

In order to get a quantitative picture of the baryon acoustic oscillations, one needs to consider the fact that baryons are tightly coupled to photons in the primordial plasma ($\delta_{b}\sim\delta_{\gamma}$). This implies that perturbations in both components are similar until recombination and analyzing photon features is the same as analyzing baryon features. Thus, let us define the Fourier transform of $\delta T/T$, the temperature perturbation of photons, as

[TABLE]

in which $\mu\equiv(\bm{k}\cdot\bm{p})/k$ and $\bm{p}$ is the photon direction. From now on we will work in the conformal time.

Instead of working with $\Theta$ directly, one might use the $l$th multipole moment of the temperature field, defined by

[TABLE]

where $\mathcal{P}_{l}$ is the $l$th Legendre polynomial. The multipole moments are an average over all photon directions and the first ones have specific names: $l=0$ is the monopole, $l=1$ is the dipole, $l=2$ is the quadrupole, and so on. Besides, photon perturbations can be characterized either by $\Theta(k,\mu,\eta)$ or by all the moments $\Theta_{l}(k,\eta)$.

The set of equations that describe the evolution of photon and baryon perturbations in the early Universe comes from the general Boltzmann equation

[TABLE]

which relates changes in the distribution function of species to possible collision terms that may also depend on the distribution function (e.g. Compton scattering). After long calculations and linearization of the perturbations, the final evolution equations relevant to BAO analysis are [45]

[TABLE]

[TABLE]

[TABLE]

where $v_{b}$ is the baryon velocity, $R_{s}\equiv(3\rho_{b})/(4\rho_{\gamma})$ is the ratio of baryon to photon density, $\tau$ is the optical depth defined by $\tau(\eta)\equiv\int_{\eta}^{\eta_{0}}d\eta^{\prime}\,n_{e}\,\sigma_{T}\,a$ ($n_{e}$ is the electron density and $\sigma_{T}$ is the Thomson cross section), and $\Phi$ and $\Psi$ are perturbations about the flat FLRW metric, corresponding to the Newtonian gauge

[TABLE]

in which $\Psi$ is the Newtonian potential and $\Phi$ is the perturbation in the spatial curvature.

In order to derive the BAO evolution equation, it is necessary to show that multipole terms higher than the dipole are negligible in the tightly coupled regime. This limit corresponds to the scattering rate of photons being much larger than the expansion rate ($\tau>>1$). The idea now is to turn the differential equation (2.17) into an infinite set of coupled equations for $\Theta_{l}(\eta)$ with the advantage that higher moments are very small. Thus, we multiply (2.17) by $\mathcal{P}_{l}(\mu)$ and integrate over $\mu$. Combining with (2.15), equation (2.17) for $l>2$ yields

[TABLE]

To solve the integral, we use the recurrence relation for Legendre polynomials

[TABLE]

to obtain

[TABLE]

Now we analyze the order of magnitude of each term. As $\tau>>1$, the first term on the left (of order $\Theta_{l}/\eta$) is much smaller than the term on the right (of order $\tau\Theta_{l}/\eta$). If we neglect the $\Theta_{l+1}$ term for a moment, it follows that in the tightly coupled limit

[TABLE]

and, for horizon size modes $k\eta\sim 1$, it implies that $\Theta_{l}<<\Theta_{l-1}$, which justifies neglecting the $\Theta_{l+1}$ term. The whole assumption is valid for $l>1$, such that the monopole and dipole terms dominate over all other multipoles in the primordial plasma.

The next step is to multiply (2.17) by $\mathcal{P}_{0}(\mu)$ and $\mathcal{P}_{1}(\mu)$ and integrate over $\mu$. One obtains

[TABLE]

and

[TABLE]

where we are using the fact that higher multipoles can be neglected. The above equations are supplemented by the equations governing baryon density perturbations, (2.18) and (2.19). One might rewrite the baryon velocity equation, (2.19), as

[TABLE]

The second term in (2.27) is much smaller than the first one because it is proportional to a factor of order $\tau^{-1}$. Thus, at lowest order, $v_{b}=-3i\Theta_{1}$. One might use this information to expand everywhere in the second term using this low order expression, obtaining

[TABLE]

In order to get rid of $v_{b}$, the above expression is inserted into (2.26). By rearranging the equation, one obtains

[TABLE]

With the pair (2.25) and (2.29), we have two first-order coupled equations characterizing the monopole and the dipole in the tightly-coupled limit. It is interesting to turn then into a second-order equation by differentiating (2.25) and eliminating $\Theta^{\prime}_{1}$ using (2.29), which leads to

[TABLE]

At last, to disappear with $\Theta_{1}$ we use (2.25) and we finally obtain

[TABLE]

where $\mathcal{H}\equiv a^{\prime}/a$ is the Hubble parameter in the conformal time, $F(k,\eta)$ is defined as a forcing function and $c_{s}$ is the propagation velocity (or sound speed) of the fluid given by

[TABLE]

The sound speed depends on the baryon density and the photon density, which implies that in the primordial Universe, when photon density was much larger that baryon density, sound waves propagated at the relativistic velocity $c_{s}=c/\sqrt{3}$. As the Universe continued the expansion, baryons became more important and this velocity was reduced (baryons make the fluid heavier, lowering the sound speed). Although (2.31) is written for the monopole, it is also valid for baryon density perturbations because $\Theta_{0}\sim\delta_{b}$ due to the photon-baryon coupling. Besides, equation (2.31) is a wave equation with friction and source terms, which is the “heart” of the baryon acoustic oscillations. As the terms for $\Theta_{0}$ and $\Phi$ have very similar forms, and noticing that $R^{\prime}_{s}=\mathcal{H}R_{s}$, one may rewrite (2.31) as

[TABLE]

Because (2.33) is an inhomogeneous second-order differential equation, it is necessary to use Green’s function method to solve it. First we find the solutions to the homogeneous part and then we use these to find the particular solution. In principle, there is a damping term in (2.33) and we should solve it considering this, with the right-hand side equal to zero. However, in the primordial Universe the pressure of the photon-baryon fluid induces oscillations with time scale much shorter than the expansion of the Universe (also $R_{s}$ is small around this epoch) [45], and therefore it is a good first approximation to discard the damping term and consider only the oscillating solutions. Thus, the homogeneous solution is

[TABLE]

with

[TABLE]

where we define the sound horizon as

[TABLE]

which, since $c_{s}$ is the sound wave propagation velocity, is the comoving distance traveled by the sound waves until the epoch $\eta$.

The general solution to a second-order equation is a linear combination of the homogeneous solutions and a particular solution. One might construct the particular solution by integrating the source term weighted by the Green’s function that is a combination of the homogeneous solutions. The general result is

[TABLE]

In this equation, we are considering $R_{s}$ very small except in the rapidly varying sines and cosines. In the $f_{1}$ term, for example, $kr_{s}$ is evaluated with $c_{s}$ in its complete form (2.32). As the constants $c_{1}$ and $c_{2}$ are fixed by the initial conditions ($\Theta_{0}$ and $\Phi$ constants, $\Theta^{\prime}_{0}=\Phi^{\prime}=0$ at $\eta=0$), $c_{1}$ must vanish and $c_{2}=\Theta_{0}(0)+\Phi(0)$. In the integrand, the denominator is $-kc_{s}(\tilde{\eta})\rightarrow-k/\sqrt{3}$ in the limit we are considering and the numerator reduces to $-\sin{[k(r_{s}-\tilde{r}_{s})]}$, leading to

[TABLE]

Equation (2.1.3) is an approximated solution that describes the oscillations in the baryon-photon fluid in the tightly coupled limit. This solution matches the exact solution with very good agreement, getting the peak locations correctly and the heights fairly well, as one can see in Fig. 2.5. Therefore, all we need are the external gravitational potentials from dark matter and then we can calculate the effect of these potentials on the anisotropies. Also, now we have an accurate way of finding the frequency of the oscillations and the location of the acoustic peaks. In the limit that the second term in (2.1.3) is negligible, the cosine term dominates and therefore the peaks should appear at the positions

[TABLE]

The dipole term is also non-negligible at this point of the cosmic history, and it will be useful in the next section because it contributes to the CMB spectrum. This one is obtained by differentiating (2.1.3) and inserting it in (2.25), which leads to

[TABLE]

2.1.4 BAO acoustic scale and relative BAO distance

The BAO acoustic scale, which corresponds to the sound horizon in a specific $\eta$, can be formally defined as the comoving distance traveled by sound waves in the photon-baryon fluid until baryons are released from Compton drag of photons [27]. This is called drag epoch (not to be confused with recombination epoch), defined as the epoch at which baryons “stop noticing” photons, occurring at redshift $z_{drag}$. The sound horizon at $z=z_{drag}$ is

[TABLE]

where

[TABLE]

is the square of the effective sound speed of sound waves in the primordial plasma, as mentioned before. There is a fitting formula for the redshift $z_{drag}$ due to Eisenstein and Hu [75], given by

[TABLE]

where

[TABLE]

$\omega_{b}\equiv\Omega_{b}^{(0)}h^{2}$ and $\omega_{m}\equiv\Omega_{m}^{(0)}h^{2}$.

Equation (2.41) can be solved analytically if one takes into account the fact that dark energy is negligible for $z>z_{drag}$. Thus, one obtains

[TABLE]

where $R_{s}(a)=(3\omega_{b}/4\omega_{\gamma})a$ and $a_{eq}=(1+z_{eq})^{-1}$ is the scale factor at the radiation-matter equality (1.70). Integrating the above equation, yields

[TABLE]

where $R_{s}^{(drag)}\equiv R_{s}(a_{drag})$ and $R_{s}^{(eq)}\equiv R(a_{eq})$. Recent measurements from Planck Satellite [49] constrain the BAO sound horizon value with great accuracy (approximately $1\%$) to be $r_{s}(z_{drag})=(147.60\pm 0.43)$ Mpc, with $z_{drag}=(1059.57\pm 0.47)$, which is the basic information for using BAO as a “cosmological ruler” in order to constrain cosmological parameters, especially the dark energy equation of state [76].

The power spectrum contains information about structures in the Universe and may be obtained analyzing the angular and redshift distribution of galaxies, which divides modes into perpendicular ($k_{\bot}$) and parallel ($k_{\parallel}$) to the line of sight. It can be showed [27] that the quantities

[TABLE]

and

[TABLE]

can, in principle, be measured independently and give good estimates for the evolution of the angular diameter distance $d_{A}$ and the Hubble parameter $H(z)$ [77]. Just as the modes, the angle $\theta_{s}(z)$ corresponds to observations perpendicular to the line of sight and the redshift difference $\delta z_{s}(z)$ to observations in the radial direction.

The BAO data so far is not sufficient for making good measurements on $\theta_{s}(z)$ and $\delta z_{s}(z)$ independently [78]. A way around this is to define a combined distance scale ratio that takes into account the combination of two spatial dimensions perpendicular to the line of sight and one dimension along the line of sight, defined by

[TABLE]

In literature, observational data is displayed in a slightly different manner from (2.49), as it follows. The main information that specifies a determined galaxy survey is the relative BAO distance, defined by

[TABLE]

where

[TABLE]

is the related effective distance. Combining (2.50), (2.46) and (1.50), the explicit form of the relative BAO distance in a flat spatial geometry is

[TABLE]

Another useful descriptor for BAO data is the acoustic parameter $A(z)$, which is independent of $h$ and is defined by [14]

[TABLE]

We shall be using BAO observational data displayed in Table 2.1. These are the main BAO data found in literature, which are used extensively for testing models and they also make part of parameter estimation codes for cosmology like CosmoMC [79] and Monte Python [80].

As WiggleZ data are correlated, we also show in Table 2.2 the inverse covariance matrix for the three WiggleZ measurements.

2.2 Cosmic microwave background radiation

2.2.1 Overview

In 1963, Penzias and Wilson [82] accidentally discovered a form of radiation that could be detected in all sky directions, which was invariant whether they changed the position of the detectors. This signal, dubbed Cosmic Microwave Background (CMB) radiation, was first predicted by George Gamow and Robert Dicke in the 1940’s, independently, and is one of the cosmic relics from which one can extract information about the primordial Universe. This discovery reinforced the Hot Big Bang model because if the Universe started in a hot, dense and opaque state, then one of the consequences would be a microwave background radiation [66].

As previously mentioned, the basic idea about CMB radiation is that in the primordial Universe there was a strong coupling between matter and radiation until the epoch where, after the Universe had expanded enough ($z\approx 1090$), photons decoupled from baryonic matter, propagating freely ever since and having almost null interaction with matter on its path. These photons are the ones we observe today in the form of microwave radiation.

There are three closely related epochs that need to be distinguished and that occur around this redshift: the recombination epoch, defined as the time in which baryons stop being ionized and become neutral by combining with electrons; the photon decoupling epoch, in which photon scattering rate with electrons becomes smaller than Hubble parameter (which is the expansion rate of the Universe), turning the Universe from opaque to transparent; and the epoch of last scattering surface (Fig. 2.6), defined as the time in which a CMB photon last scattered off an electron. Once the expansion rate becomes larger than the scattering rate, it is very unlikely that a photon will scatter again, which makes the last scattering epoch very close to the photon decoupling epoch. That is why in the last section these moments in the history of the Universe were taken to be the same, but here they are properly defined.

The CMB measurements were performed for the first time by COBE satellite in 1992 [47] and refined by other satellite measurements such as WMAP in 2003 [83] and Planck in 2013 [84]. The difference among these experiments can be seen in Fig. 2.7. Even though experiment improvements were huge, some of the important discoveries provided by COBE data are still valid and can be elucidated in the following terms:

To any direction in the sky (angular coordinates ($\theta,\phi$)), CMB spectrum is very close to the one of an ideal black-body (Fig. 1.3) with mean temperature

[TABLE] 2. 2.

There is a dipole distortion in the temperature maps due to Doppler effect caused by the motion of COBE satellite relative to the frame of reference where CMB is isotropic (Fig. 2.8), which means that different hemispheres from the sky are slightly blueshifted (or redshifted) to higher (or smaller) temperatures; 3. 3.

After dipole distortion removal, the remaining temperature fluctuations are very small in amplitude. This can be better understood if one defines the dimensionless temperature fluctuation across the sky as

[TABLE]

Without dipole distortion, the root mean square temperature fluctuation (which is an average over all points in the sky except the ones in the region contaminated by our own galaxy foreground emission) for COBE data is

[TABLE]

indicating an extraordinary closeness to isotropy.

Due to the huge complexity of the equations involved in generating the full CMB temperature spectrum (there are nine coupled equations that account for all components and for the metric, called Einstein-Boltzmann equations), numerical simulations are necessary to solve the problem. For this purpose, open-source codes such as CAMB [87] and CLASS [88] were created, which are written in Fortran and C, respectively. Even with these difficulties, we can elucidate in general lines the equations that lead to CMB anisotropies and also study important parameters that describe some characteristics of the spectrum.

2.2.2 Anisotropy power spectrum and Einstein-Boltzmann equations

Let’s consider temperature fluctuations $\delta T/T$ across the sky as measured from a particular experiment. Since each point in the sky has a unique temperature $T(\theta,\phi)$ defined on the surface of the celestial sphere, it is intuitive to expand the fluctuations in spherical harmonics as

[TABLE]

in which $Y_{\ell m}(\theta,\phi)$ represent spherical harmonic functions. Useful information can only be extracted by using statistics since we have a temperature distribution over the whole sky. Thereby, the most important statistical descriptor of $\delta T/T$ is the 2D correlation function $C(\theta)$ (and its corresponding Fourier pair, the 2D power spectrum $C_{\ell}$) likewise BAO information uses the 3D correlation function $\xi(r)$ (and 3D power spectrum $P(k)$).

For this purpose, consider two different directions in the sky at the last scattering surface. These points have directions $\hat{n}$ and $\hat{n}^{\prime}$ separated by the angle $\theta$ such that $\cos{\theta}=\hat{n}\cdot\hat{n}^{\prime}$. The correlation $C(\theta)$ between temperature fluctuations for pair of directions is defined as the product of the fluctuations averaged over all points separated by the same angle $\theta$, in the form

[TABLE]

If it was possible to obtain from measurements the precise form of $C(\theta)$ for angles from $\theta=0$ to $\theta=180\degree$, we would have a complete statistical description for the fluctuations over all scales. Unfortunately, there are experimental limitations that make it impossible because only a limited range of angular scales are accessible with precision [66]. Using (2.57) and (2.58), the correlation function can be written as

[TABLE]

where $P_{\ell}$ are Legendre polynomials.

The important thing in the expansion (2.59) are the $C_{\ell}$ coefficients, because now we can characterize the correlation function $C(\theta)$ by its $C_{\ell}$ multipole moments, which are the way CMB data are usually displayed. Generally speaking, a multipole moment is a measure of temperature fluctuations at angular scale $\theta\sim 180\degree/\ell$ [45]. Thus, it does not matter whether one prefers to use $\theta$ or $\ell$ to plot the anisotropy spectrum.

In order to obtain the $C_{\ell}$ coefficients, one needs to recall the definition of temperature perturbation $\Theta$ as mentioned in the last section, which comes from the equation

[TABLE]

As one can see, the temperature field is defined at every point in space ($\bm{x}$) and time ($\eta$), but our observations can only be performed here ($\bm{x}_{0}$) and now ($\eta_{0}$). The information we receive is based exclusively on the direction of the incoming photons ($\hat{p}$), which means we observe temperature fluctuations in all directions but we cannot, at least at first sight, infer that the temperature in a specific direction is due to a photon that came right from the last scattering surface because a photon might have its path altered due to sources in the way (there are other effects to account for potential wells in the path between last scattering surface and us, such as Sunyaev-Zel’dovich and Sachs-Wolfe [50]). Instead of working with angular coordinates ($\theta,\phi$), we will stick to the photon direction coordinate $\hat{p}$ and rewrite the expansion (2.57) as

[TABLE]

All the relevant information from the temperature field $T$ is encapsulated by the amplitudes $a_{\ell m}$ (which are related to $C_{\ell}$, as we shall see). In order to invert (2.61) and write $a_{\ell m}$ as a function of $\Theta$, we can use the orthogonality condition for the spherical harmonics, normalized via

[TABLE]

in which $\Omega$ is the solid angle covered by $\hat{p}$. Multiplying (2.61) by $Y^{*}_{\ell m}(\hat{p})$, integrating it and using (2.62), the result for $a_{\ell m}$ is

[TABLE]

where we have written $\Theta(\bm{x})$ in terms of its Fourier transform $\Theta(\bm{k})$, which is the one we have equations to work with.

The $a_{\ell m}$’s cannot be exactly predicted, and as usual it is their statistical properties we are interested in order to extract information. For this purpose, the distribution from which they are drawn is the main goal of CMB measurements [45]. As the coefficients $a_{\ell m}$ are assumed to be statistically independent [27], we define the mean and variance of the $a_{\ell m}$’s, respectively, as

[TABLE]

where the variance $C_{\ell}=\langle\left|a_{\ell m}\right|^{2}\rangle$ is called CMB temperature power spectrum.

We can now find an expression for $C_{\ell}$ by using (2.63) and (2.64). For that matter, we need first to evaluate $\langle\Theta(\bm{k},\hat{p})\Theta^{*}(\bm{k}^{\prime},\hat{p}^{\prime})\rangle$, where $\eta$ dependence is implicit. To solve the problem, we need to rewrite the temperature perturbation as $\delta\times(\Theta/\delta)$, in which dark matter perturbation $\delta$ does not depend on $\hat{p}$. As the ratio $\Theta/\delta$ does not depend on the initial amplitude of a mode, it can be taken out of the average. Thus,

[TABLE]

where in the second equality it was considered the fact that different modes are uncorrelated [27] and that the ratio $\Theta/\delta$ depends only on the magnitude of $\bm{k}$ and the product $\hat{k}\cdot\hat{p}$ [45]. Therefore, combining (2.63), (2.64) and (2.2.2) yields

[TABLE]

One can expand $\Theta(k,\hat{k}\cdot\hat{p})$ and $\Theta^{*}(k,\hat{k}\cdot\hat{p}^{\prime})$ using the inverse of (2.15), $\Theta(k,\hat{k}\cdot\hat{p})=\sum_{\ell}(-i)^{\ell}(2\ell+1)\,\mathcal{P}_{\ell}(\hat{k}\cdot\hat{p})\,\Theta_{\ell}(k)$, which leads to

[TABLE]

The two angular integrals are equal to $4\pi Y_{\ell m}(\hat{k})/(2\ell+1)$ (or the complex conjugate) in the case $\ell^{\prime}=\ell$ and $\ell^{\prime\prime}=\ell$, otherwise they are zero. The remaining $\left|Y_{\ell m}\right|^{2}$ combines with the angular part of the $d^{3}k$ integral (which is equal to one), leading to

[TABLE]

If one uses (2.68) to evaluate CMB power spectrum, it is necessary to know the form of the matter power spectrum, dark matter overdensity and the multipole moments as a function of scale, while performing an integral over all Fourier modes. Approximate solutions for these terms can be found in [45] and [50], but there is no general solution for the problem and, in practice, one needs a code for solving numerically the set of Einstein-Boltzmann equations. As an example, for $\Lambda$CDM model we have nine coupled differential equations describing interactions among components and the metric of the Universe (scalar perturbations are the main source of CMB anisotropies). These equations are summarized as [45]

[TABLE]

and, if dark energy is not $\Lambda$CDM, a new equation must be added to account for dark energy perturbations and possible interactions with other components.

From these ones, three of them were already used in BAO section in a compact form (eqs. (2.17), (2.18) and (2.19)). The term $\Pi=\Theta_{2}+\Theta_{P2}+\Theta_{P0}$ and $\Theta_{P}$ are terms due to polarization of the CMB, $\delta$ and $v$ represent CDM overdensity and velocity, respectively, and $\mathcal{N}$ is the neutrino perturbation (it is the analogue of the photon perturbation variable $\Theta$). It is now clear that either one uses approximations together with equation (2.68) or solves numerically the full set of Einstein-Boltzmann equations to find $\Theta$ (in order to get more accurate results) and combines with (2.63) and (2.64) to obtain finally the CMB power spectrum. In the literature, temperature power spectrum is usually plotted as

[TABLE]

in units of $\mu K^{2}$, as displayed in Fig 2.9.

Our approach for evaluating CMB power spectrum is numerical: we use CAMB code [79, 89], a code publicly available for over a decade, very well tested and improved by the community. In order to predict CMB power spectrum with high accuracy, we also take into account lensing in the code, which improves the smoothing effect on the acoustic peaks by $5\%$. This is the correct procedure as pointed out in [49].

2.2.3 Optical depth

An important cosmological parameter related to CMB spectrum is the Thomson scattering optical depth due to reionization $\tau$. Before giving a definition of it, it is useful to introduce some concepts. The Universe has gone through important phases in its history, some of them having hydrogen as the main factor related to these phase transitions. We can cite, for example, recombination (as already described earlier) and reionization epoch, the latter being the period in which the Universe goes from being neutral to becoming ionized due to the formation of the first generation of galaxies, which emitted ultraviolet radiation and started ionizing large intergalactic gas clouds [90].

Even though in this epoch there was a considerable amount of helium, the significant physics can be understood simply by considering just hydrogen in its neutral form (for which we shall use the letter $H$) and its ionized form (a single proton, letter $p$). Charge neutrality requires the number density of free protons and free electrons to be the same, $n_{e}=n_{p}$. In order to quantify the ionization level of the Universe, one defines free electron fraction (or fractional ionization) as

[TABLE]

in which the subscript $b$ stands for baryons. The quantity $X_{e}$ ranges from $X_{e}=1$ when baryons are fully ionized to $X_{e}=0$ when the Universe is entirely neutral.

In the reionization period, the interactions happened primarily between photons and electrons through Thomson scattering, which made these particle to exchange energy and momentum. Thus, the mean free path of a photon in this epoch (the mean distance it travels before scattering from a electron) is

[TABLE]

where $\sigma_{T}$ is Thomson scattering cross-section ($\sigma_{T}=6.65\times 10^{-29}\,m^{2}$). Due to the fact that photons travel at the speed of light $c$, their scattering rate is

[TABLE]

Photons we observe today came from the last scattering surface and, in its path toward us, may have suffered scattering in the reionization epoch. Therefore, at the time interval $t\rightarrow t+dt$ the probability that a photon is scattered off an electron is $dP=\Gamma(t)dt$. The definition of optical depth to reionization is the integral of the scattering rate over all cosmic history, which is equivalent to the expected number of scatterings it has suffered since earlier times. Mathematically,

[TABLE]

where $n_{e}$ is the number density of free electrons produced by reionization. Even though the integral is over all cosmic history, we do not need to worry about the fact that there is a small residual ionization fraction coming from recombination because this one is of order $\sim 10^{-3}$, being safely neglected [91].

Using the relation $dt=da/(aH)=-dz/[(1+z)H]$ and (2.79), one can rewrite $\tau$ as

[TABLE]

Considering the fact that reionization happens at a matter-dominated epoch, $H(z)\sim(1+z)^{3/2}$ and therefore

[TABLE]

Equation (2.84) can be used to evaluate $\tau$ either by using some differential equation describing the evolution of $X_{e}(z)$ [45] or by parametrizations of $X_{e}$ as a function of $y\equiv(1+z)^{3/2}$, as it is standard procedure in CAMB code [91]. Recent measurements for $\tau$ come from Planck satellite, yielding the value $\tau=0.078\pm 0.019$ [49].

2.2.4 Primordial fluctuations: scalar spectral index and amplitude perturbation

Inflation is the best scenario to explain why uncorrelated scales observed today have such similar temperatures, as we have seen, at the same time it solves other independent problems in cosmology and particle physics at once [92]. Besides, inflation is a mechanism for generating the primordial fluctuations that led to CMB anisotropies, which is why improvements on CMB measurements could confirm inflationary theories. One of the ways of checking whether CMB points to inflation is to work with the primordial power spectrum (scalar spectrum) that comes from the metric perturbations (2.20), which can be demonstrated to be [45]

[TABLE]

where $M_{Pl}$ is the Planck mass, $n_{s}$ is the scalar spectral index, $\delta_{H}$ is the scalar amplitude at horizon crossing $aH=k$, $D_{1}$ is the growth function which describes the growth of matter perturbations at late times [45], $\Omega_{m}^{(0)}$ is the matter density parameter today and $\epsilon$ is the slow-roll inflation parameter defined by

[TABLE]

which measures changes in the Hubble rate during inflationary era that are supposed to be small due to dominance of the potential energy.

If a spectrum has the property that $k^{3}P_{\Phi}(k)$ is a constant (i.e., has no dependence on the scale), then it is called a scale-invariant (also called scale-free or Harrison-Zel’dovich) spectrum, which is characterized by $n_{s}=1$. Therefore, the scalar index quantifies deviations from scale invariance and from perfect de Sitter limit [93]. For example, most models of inflation are based on the idea that a scalar field $\phi(\bm{x},t)$ (called inflaton) runs the inflationary era, “slowly rowling” to the minimum of the potential energy $V(\phi)$ during the short period that inflation takes place. In these models, it can be demonstrated [93] the relation

[TABLE]

where $\delta$ is another slow-roll inflation parameter, given by

[TABLE]

with $\phi$ representing now the zero-order term $\phi^{(0)}(t)$ from the space-time dependent scalar field.

From (2.87), it becomes obvious that an inflationary scenario predicts a deviation from scale invariance, since $\epsilon$ and $\delta$ are predicted to be small but non-zero. This is exactly what is observed: Planck measurements constrain the scalar spectral index to the value $n_{s}=0.9655\pm 0.0062$ [49], giving much credit to the inflationary scenario.

It is customary to model the almost scale-invariant spectrum as a power-law close to a reference scale $k_{0}$ (also called pivot scale or tilt) as [93]

[TABLE]

where $A_{s}$ is the amplitude of the scalar spectrum measured at the pivot scale $k_{0}=0.05\,\text{Mpc}^{-1}$, approximately in the middle of the logarithmic range of scales probed by Planck. The choice of $k_{0}$ value is made to reduce degeneracy between $n_{s}$ and $A_{s}$ [49].

Using (2.85), it is possible to identify the spectrum amplitude as $\delta_{H}\sim\mathcal{P}_{\mathcal{R}}(k_{0})^{1/2}=A_{s}^{1/2}$, which is constrained by Planck measurements to be $\ln{(10^{10}A_{s})}=3.089\pm 0.036$ [49], where $\ln$ is the natural logarithm (the logarithm of $A_{s}$ is also called log power of the primordial curvature perturbations). From (2.89), the scalar index can also be defined as

[TABLE]

We shall not work with tensor power spectrum that comes from tensor perturbations ($n_{T}=0$) neither with higher-order effects such as running of scalar ($dn_{s}/d\ln{k}=0$), since Planck measurements are weakly sensitive to these parameters [49].

While CMB theoretical power spectrum comes from CAMB code, CMB data comes from the most recent measurements of the Planck satellite. This information is encoded in the Planck likelihood code package (also called Plik), available at [94], which takes as input the theoretical power spectrum and a set of “nuisance” parameters describing unresolved foreground and instrument calibration, and outputs the log-likelihood (which shall be described in the next chapter). In this dissertation, we use two likelihood codes: Plik TT likelihood, which uses a high-$\ell$ ($30\leq\ell\leq 2508$) temperature-only dataset together with 16 “nuisance” parameters, and low TEB likelihood (lowP), which uses low-$\ell$ ($2\leq\ell\leq 29$) temperature plus polarization datasets with a single “nuisance” parameter. The latter likelihood code is important to break the large degeneracy between $\tau$ and $\ln{A_{s}}$ [95], while using polarization for high-$\ell$ multipoles is not necessary because it does not effectively improve parameter constraints [49].

2.3 Cosmic chronometers

2.3.1 Overview

Cosmic chronometers (or cosmic clocks) is a nearly model-independent tool developed for directly measuring the evolution of the Hubble parameter. This method is based on objects (in this case, galaxies) from which one may roughly infer the variation in the age of the Universe with the redshift [96]. Using these “cosmological clocks”, one may infer the age difference $\Delta t$ between samples of nearby early-type galaxies (ETG) from its spectral lines and also measure the redshift difference $\Delta z$ between them. The connection between these measurements and the Hubble parameter is made by identifying the derivative $dz/dt$ with the ratio $\Delta z/\Delta t$ and by using the straightforward equation

[TABLE]

to describe the Hubble parameter observationally. As this formula is very simple and has no integral nature (differently from methods based on integrated quantities, such as luminosity distance from supernovae), it has a good sensitivity for constraining a dynamical dark energy equation of state. Also, the way cosmic chronometers are applied, with a “differential age method”, is much more reliable than working with an “absolute age method”, which is much more susceptible to systematics and potential effects of galaxy evolution [97]. Therefore, the relevant physics is ranged between the redshifts where the differential is applied.

Although the idea of this method seems very simple, complications in estimating $H(z)$ are huge. For better exploration of the method, one needs to select samples of red elliptical galaxies in the form of passively evolving galaxies, which are defined as the ones with old stellar population, low-rate stellar formation and with stellar mass $\mathcal{M}$ above $10^{11}\mathcal{M}_{\odot}$ [98]. Besides, errors in the final estimates can be large (even more than $30\%$), requiring combination with other independent datasets to better constrain parameters.

One of the best ways to perform the analysis is to use the 4000 Å break (which we shall name D4000) present in ETG spectra, since it is linearly related to the age of the galaxy in the case where the galaxy has an old stellar population. This break models a discontinuity in the spectrum from these galaxies close to the wavelength $\lambda_{rest}=4000$Å, which occurs due to metal absortion lines with amplitude that scales linearly with metallicity Z and age of the stellar population [98]. Because star formation history (SFH) has a small influence only if one considers old passively evolving galaxies, the criteria adopted above must be respected. The 4000 Å break feature is displayed in Fig. 2.10.

2.3.2 Differential age method from galaxy surveys

The D4000 is defined as ratio between the continuum flux densities in a red band (4050-4250Å) and a blue band (3750-3950Å) around the 4000Å wavelength [99]. It is written as

[TABLE]

although reference [98] adopts a slightly different range for the bands in order to reduce errors due to dust reddening: 3850-3950Å and 4000-4100Å. From now on, this narrower band definition will be represented as $D4000_{n}$.

As mentioned, if the galaxy is old and passively evolving, there is a linear relation between $D4000_{n}$ and galaxy age (with constant metallicity) of the form

[TABLE]

where $A(Z,SFH)$ is a factor to be calculated by means of stellar population models and $B(Z,SFH)$ is a “nuisance” parameter that must vanish by applying the differential age method as

[TABLE]

Now becomes clear one of the advantages of the method: it directly traces the age evolution of ensembles of galaxies.

Thereby, Hubble evolution can be rewritten as a function of $D4000_{n}$ and $A(Z,SFH)$, which leads to

[TABLE]

As an example, measurements of $D4000_{n}$ for a sample of 11324 ETGs are displayed in Fig. 2.11 as gray dots (these data are the base for the calculation of 8 from the 31 cosmic chronometers measurements for $H(z)$ that will be presented later on). Orange and green dots represent the median value of $D4000_{n}$ in a determined redshift bin for different mass ranges: green for lower galaxy masses and orange for higher galaxy masses. The choice of the size of the bins depends on the amount of galaxies available to perform the statistic, and it has a crucial importance in the sampling. Also, for a redshift smaller than $0.3$, metallicity is a feature that can be estimated from the spectrum, while for a redshift larger than $0.3$ this information is no longer available and one needs to take a conservative choice for metallicity [98].

Another procedure that has to be made is to calibrate the relation between $D4000_{n}$ and galaxy age by using models of stellar population synthesis in order to evaluate $A(Z,SFH)$. It is important to confirm that there is no bias coming from the model, and for that matter reference [98] employs two independent samples of synthetic spectra (BC03 and MaStro), finding no significant difference between them. Other details about all the procedures to obtain estimates for $H(z)$ are fully described in reference [96] and [98].

The dataset for cosmic chronometers we are going to use is a compilation of 31 measurements for $H(z)$ in the redshift range $0<z<2$, as shown in Table 2.3.

2.4 Supernovae

2.4.1 Overview

A supernova is an event that occurs due to the death of certain types of stars, such that its brightness can be so intense as the brightness of a whole galaxy, as can be seen in Fig. 2.12. Historically, supernovae were divided in two classes, according to their spectra: Type I and Type II supernovae. Type I supernovae have no hydrogen absorption lines in their spectra, while Type II supernovae do [66]. For Type I supernovae there are three subdivisions: Type Ia (SN), characterized by an absorption line of singly ionized silicon; Type Ib, defined as supernovae with a line of helium in their spectra; and Type Ic, the ones that lack both lines. As Type Ib and Type Ic have similar properties to Type II supernovae (they are formed due to core collapse of the original star), we will focus on Type Ia supernovae because they have mechanisms of origin and evolution totally different of the previous ones and they are objects with “standardizable” properties.

It is believed that Type Ia supernovae occur when a white dwarf (a star whose gravity is counterbalanced with its electron degeneracy pressure) in a binary system accretes enough matter from its binary companion until it reaches Chandrasekhar limit [102], which renders the white dwarf unstable by increasing its temperature and density, ultimately leading to a thermonuclear explosion (this can be demonstrated by applying Fermi-Dirac statistics) [50]. Because the exploding star has always a mass close to this limit (approximately $1.44\,M_{\odot}$, see Fig. 2.13), supernovae absolute luminosity varies for only a small amount at peak brightness, making them one of the best distance indicators we have today in cosmology (we call objects with this property standardizable candles for fact that their luminosity can be standardized). Thus, the main information one needs to extract from a supernova is the light curve (absolute luminosity versus time since peak brightness).

The occurrence rate of supernovae in our galaxy is one per century, approximately. Even though it is a small rate, the fact that supernovae explosion is extremely luminous makes much easier to detect them even at very large distances (up to $z\sim 1$). For some days (and sometimes months), Type Ia supernovae light can even outshine light from its host galaxy.

Although variation in the absolute magnitude of different supernovae is small, this effect needs to be taken into account and supernovae need calibration in order to use their information as a cosmological probe. For example, the average luminosity of Type Ia supernovae is, at peak brightness, $L=4\times 10^{9}\,L_{\odot}$, while Cepheid methods for determining distances show that Type Ia supernovae peak luminosity lie in range $L\approx 3\rightarrow 5\times 10^{9}\,L_{\odot}$ [66]. In order to treat this variation in magnitude, M. Phillips [104] studied a sample of nine supernovae in the B, V and I bands of the spectrum, showing there is a linear relation between peak magnitude and the amount in magnitude that the light curve decays from its peak brightness at the time interval of 15 days ($\Delta m_{15}$). For this reason, it becomes evident there is a correlation between light curve width and peak luminosity: the higher the peak brightness, the broader the light curve (Fig. 2.14). One of the first calibrations using this method was performed at the end of the 1990s from a local sample ($z<<1$) with the Calán/Tololo survey [105].

2.4.2 Magnitude system

In order to relate luminosity distance (1.47) to observational data, one needs first to define the concept of magnitude used to express luminosity and flux. The magnitude system, created by the Greek astronomer Hipparchus in the second century BC, was divided into six classes depending on the apparent brightness of visible stars: the brightest ones were of first magnitude, the faintest ones were of sixth magnitude and the intermediate ones were of second magnitude, third magnitude and so on. Some time later, it was discovered that human eye responds to a large range of light intensity, actually better described in logarithmic scale, and this suggested to define the magnitude system in terms of a logarithmic function, as follows.

We define the bolometric apparent magnitude $m$ of a light source in terms of the source’s bolometric flux $f$ as

[TABLE]

where $f_{x}$ is the reference flux with value $f_{x}=2.53\times 10^{-8}\,W/m^{2}$. The negative sign ensures that an object with a large flux has small bolometric apparent magnitude (e.g. $f_{\odot}=1367\,\,W/m^{2}\rightarrow m_{\odot}=-26.8$). The reference flux was defined considering the star Vega, such that visible stars to the naked eye are typically in range $0<m<6$.

The second quantity to be defined is the bolometric absolute magnitude $M$ from a light source in terms of its luminosity $L$, given by

[TABLE]

where the reference luminosity is $L_{x}=78.7\,L_{\odot}$. The quantity $M$ is defined as the bolometric apparent magnitude that a light source would have if it were at a luminosity distance of $d_{L}=10\,$pc from the observer (e.g. for Type Ia supernovae $M\approx-19$ at the peak brightness and, for the Sun, $M_{\odot}=4.74$) [27, 66]. The reference $L_{x}$ is based on this definition, since it is the luminosity of a source at 10 parsecs of distance producing the reference flux $f_{x}$. The magnitude system seems awkward at first glance, but it is nothing more than the logarithmic measure of both flux and luminosity.

2.4.3 Luminosity distance, distance modulus and B-band generalized model

With these definitions, the relation between bolometric apparent magnitude, bolometric absolute magnitude and luminosity distance is given by, using (1.43),

[TABLE]

which may be rewritten as

[TABLE]

where $d_{L}$ must be given in Megaparsecs.

Since measured quantities such as flux and luminosity are usually quoted in terms of magnitudes, it is convenient to write luminosity distance in terms of distance modulus $m_{mod}$ to a light source, defined by

[TABLE]

Due to calibration of type Ia supernovae, as mentioned before, equation (2.100) must be generalized to the form ($d_{L}$ in Megaparsecs) [67]

[TABLE]

where $\alpha$, $\beta$ and $M_{B}$ are nuisance parameters ($M_{B}$ represents supernovae absolute magnitude), $m_{B}^{mod}$ is the apparent magnitude in the B-band (or observed peak magnitude in the rest-frame B band), $X_{1}$ describes the time stretching of the light curve and $\mathcal{C}$ represents the supernova color parameter at maximum brightness ($X_{1}$ and $\mathcal{C}$ are parameters that model the light curve). Besides, luminosity distance is defined in two different redshift measurements: $z_{hel}$ is the redshift in the Sun’s reference frame and $z_{CMB}$ represents the redshift in CMB reference frame. These redshift definitions must be accounted because luminosity distance is a combination of the real comoving distance $D(z)$ (measured when peculiar velocities and gravitational redshifts are subtracted) with a factor of $(1+z)$ coming from relativistic beaming and time dilation (each one contributing with $\sqrt{1+z}$) [108]. Therefore, luminosity distance must be written as $d_{L}=(1+z_{hel})D(z_{CMB})$.

It was found that supernovae absolute magnitude is dependent on host galaxy properties in a way still not fully understood [109], but that has important implications in cosmology such that a simple model must be assumed for $M_{B}$ in the form

[TABLE]

where $M_{stellar}$ represents the host stellar mass.

In our analysis, supernovae samples were taken from Joint Light-curve Analysis (JLA) dataset, which contains 740 supernovae distributed in the redshift range $0<z<1.3$ and that can be obtained in [110]. For JLA dataset, the covariance matrix is given by [111]

[TABLE]

where

[TABLE]

[TABLE]

The first three terms in (2.4.3) account for errors in redshift due to peculiar velocities, errors in the magnitude due to gravitational lensing and intrinsic errors that are not already accounted in the former terms, respectively. The fourth term accounts for errors in the observed magnitude, the fifth term takes into account errors in the stretching factor $X_{1}$, and the sixth accounts for colour $\mathcal{C}$ error, while the remaining terms represent their covariances. Terms in $\bm{D}_{stat}$ are all diagonal, while terms in $\bm{C}_{stat}+\bm{C}_{sys}$ have also non-diagonal elements accounting for the same thing as $\bm{D}_{stat}$ (except for the first three terms in (2.4.3)) and including systematics. All the $V$ terms in (2.105) are provided in [110].

Chapter 3 Statistical analysis

In this chapter, we review some statistical tools necessary to analyze our datasets: likelihood function, Bayes’ theorem, Markov chain Monte Carlo technique and principal component analysis. These tools are standard in modern cosmology, and their use is essential for a correct parameter estimation.

3.1 Likelihood function and Bayes’ theorem

Suppose we know that a random variable $x$ has a probability distribution function (PDF) $f(x|\theta)$ which is a function of an unknown parameter $\theta$. The “$|$” distinguishes the random variable $x$ from the parameter $\theta$. This PDF is defined as a conditional probability of having the data $x$ given the theoretical parameter $\theta$, where it is obvious that different $\theta$ will lead to different PDF’s. Besides, if one samples $x$ as $f(x|\theta)$ then it is more likely that $x$ has a value that maximizes $f(x|\theta)$, while if $x$ has a particular value, then it is better to choose $\theta$ such that it maximizes the occurrence of that $x$. Therefore, the best $\theta_{i}$ are the parameters that maximize $f(x_{1},x_{2},\dots,x_{n}|\theta_{1},\theta_{2},\dots,\theta_{m})$, where generally parameter space is large and one usually writes $f$ in the compact form $f(x_{i},\theta_{j})$. This leads us to the maximum likelihood method, which consists in finding parameters $\theta_{j}$ that maximize the likelihood function $f(x_{i},\theta_{j})$ by solving the equations

[TABLE]

In cosmology, sometimes data are pretty scarce and one has to use information from previous experiments in order to better constrain information. Besides, one is often worried about estimating the theoretical parameters $\theta_{j}$ from the data. The way to handle this problem is to use Bayesian statistics: instead of working with $f(x_{i}|\theta_{j})$, the probability of having the data given the model, it is more important to work with $L(\theta_{j}|x_{i})$, the probability of having the model given the data. These probabilities are related through the so-called Bayes’ theorem [112]:

[TABLE]

where data $x_{i}$ are denoted by $D$ and theoretical parameters $\theta_{j}$ are denoted by $T$. In this equation, $P(D|T)$ is the conditional probability of having the data given the theory, $P(T)$ and $P(D)$ are the probability of having the theory and the data, respectively, and $P(T|D)$ is the conditional probability of having the theory given the data. Thus,

[TABLE]

in which $p(\theta_{j})$ is called prior probability for the parameters $\theta_{j}$, $g(x_{i})$ is the PDF of the data $x_{i}$, $f(x_{i}|\theta_{j})$ is the likelihood function and $L(\theta_{j}|x_{i})$ (or simply $L(\theta_{j})$) is called posterior probability. The latter contains valuable information: it is the probability distribution for the theoretical parameters $\theta_{j}$ based on the observed data $x_{i}$ and knowledge from previous experiments.

Since $L(\theta_{j}|x_{i})$ is a probability distribution, one can normalize it to unity as

[TABLE]

which leads to

[TABLE]

where $n$ is the parameter dimension.

The term $g(x_{i})$ is called evidence and is just a normalization factor, having no influence on the estimation of $\theta_{j}$. On the other hand, the prior $p(\theta_{j})$ is an essential term and has large influence in the final estimation. For example, if a previous experiment convincingly excluded, say, $\Omega_{m}^{(0)}<0.1$, therefore one might safely use this information in the form $p(\Omega_{m}^{(0)}<0.1)=0$. As another example, if for some reason we believe that a parameter $y_{cal}$ has the value $y_{cal}=1.0000\pm 0.0025$ (this is actually the parameter describing absolute map calibration for Planck), then we could use for $p(y_{cal})$ a Gaussian prior with mean 1 and standard deviation 0.0025. In general, one usually excludes nonphysical values, e.g. $\Omega_{m}^{(0)}<0$ along with using flat priors, which is given by $p(\theta_{j})=1/\Delta\theta$ in a specified range $\Delta\theta$ where $\theta_{j}$ is defined and $p(\theta_{j})=0$ outside (such that $\int p(\theta_{j})\,d\theta_{j}=1$). If one uses only flat priors, then

[TABLE]

meaning that the posterior corresponds to the likelihood function $\mathcal{L}$ in this case.

3.2 Parameter estimation

After collecting data, it is important to translate them into a model, which is a theoretical background we assume to be true. The model has some parameters $\theta_{i}$ in it, which one wants to determine. Therefore, parameter estimation is the analysis performed to obtain estimates of the parameters, errors, or ideally the complete posterior probability distribution $L(\theta_{i})$, which takes into account some priors assumed. Thus, in order to find the maximum posterior estimators $\hat{\theta}_{i}$, equation (3.1) needs to be replaced by

[TABLE]

By discarding the evidence factor $g(x_{i})$, we need to normalize the posterior with an overall normalization $N$ through the integral

[TABLE]

where this one is extended to the whole parameter range. Using the normalized posterior (i.e. $L(\theta_{i})/N$, to which we shall continue to use $L(\theta_{i})$), we can evaluate the regions of confidence $R(\alpha)$ for the parameters, defined by

[TABLE]

where $R(\alpha)$ is the region for which the integral above yields $0<\alpha<1$. One usually chooses the values $\alpha=0.683,\,0.954,\,0.997$ (denoted as 1, 2 and 3$\sigma$, respectively).

It is common to have interest in some specific subset of parameters and consider the others as “nuisance” which we would be happy to get rid of. For example, supernovae theory predicts for distance modulus $m_{mod}$ the equation $m_{mod}=5\log_{10}d_{L}(z;\Omega_{m}^{(0)})+25-\alpha X_{1}+\beta\mathcal{C}+M_{B}$, from which the main parameter of interest is $\Omega_{m}^{(0)}$. The general posterior depends on the full set of parameters ($\Omega_{m}^{(0)}$, $\alpha$, $\beta$ and $M_{B}$) but we can turn it into a function only of $\Omega_{m}^{(0)}$ by integrating out all the other parameters,

[TABLE]

where the integration is over the domain of the “nuisance” parameters. This procedure is called marginalization and allows one to translate a multidimensional posterior into a one or two-dimensional probability distribution, which makes visualization of the problem much easier. For instance, if the maximum posterior estimator of $\Omega_{m}^{(0)}$ is 0.3 and

[TABLE]

in the region $\Omega_{m}^{(0)}=[0.27,0.33]$, the final result is quoted as $\Omega_{m}^{(0)}=0.30\pm 0.03$ at $68.3\%$ confidence level (or $1\sigma$).

Let us assume that the posterior $L(\bm{\theta}|\bm{x})$ ($\bm{\theta}$ and $\bm{x}$ denoting arrays of parameters $\theta_{i}$ and data $x_{i}$) is single-peaked. A common estimator for the parameters is the expected value (or mean) $\langle\bm{\theta}\rangle$ defined by

[TABLE]

and variance $\sigma^{2}$ defined by

[TABLE]

where we use Markov Chain Monte Carlo method (which shall be explained in the next section) to sample the posterior $L(\bm{\theta}|\bm{x})$, in which $\bm{\theta}^{(t)}$ are taken directly from the posterior by construction [112]. Also, it is common to work with standard deviation $\sigma$ instead of variance: $\sigma\equiv\sqrt{\sigma^{2}}$.

It is usual to approximate the likelihood function close to the peak with a multivariate Gaussian in parameter space, in the cases where the likelihood is single-peaked and well behaved (its peak is well located, not dispersed). This can be seen whether one performs a Taylor expansion close to the peak $\bm{\theta}_{0}$ of the log likelihood,

[TABLE]

which leads to

[TABLE]

where $H_{\alpha\beta}\equiv-\frac{\partial^{2}\ln{\mathcal{L}}}{\partial\bm{\theta}_{\alpha}\partial\bm{\theta}\beta}$ is called Hessian matrix and controls whether the estimates of $\theta_{\alpha}$ and $\theta_{\beta}$ are correlated or not [113]. Besides, in cosmology, for many problems of interest the observational data generally have a Gaussian distribution and therefore

[TABLE]

where $x_{i}$ is the observed data, $x_{i}^{th}(\bm{\theta})$ is the theoretical prediction for $x_{i}$ given a set of parameters $\bm{\theta}$ and $C_{ij}^{-1}$ is the inverse covariance matrix associated with $x_{i}$: $C_{ij}\equiv\langle(x_{i}-\langle x_{i}\rangle)(x_{j}-\langle x_{j}\rangle)^{T}\rangle$. The function inside square brackets (up to the constant $-1/2$) is usually denoted as the chi-square function

[TABLE]

and if there is no correlation between data, then (3.17) simplifies to

[TABLE]

where $\sigma_{i}^{2}$ stands for data variance. In this case, maximizing the likelihood function is equivalent to minimizing the $\chi^{2}$ function and, thus, one needs to find a set of parameters $\theta_{i}$ intrinsic to theory such that the latter agrees as best as it can with observations.

It is fundamental to combine different datasets for better constraining parameters because each observation has its own degenerate directions, i.e. directions in parameter space poorly constrained by the data. Thus, pairs of datasets have much stronger power constraints than each one separately because it is only their combination that mutually breaks parameter degeneracies [112]. For any datasets, one might consider them as coming from independent experiments such that their combination happens when one performs the product of the likelihoods. Equivalently, this is the same thing as adding up the chi-square functions from all the experiments. For the datasets described in the previous chapter,

[TABLE]

3.3 Markov Chain Monte Carlo method: emcee

Generally, problems in cosmology have large dimensionality in parameter space (up to 20, 30 parameters). This fact makes it unfeasible to evaluate the full likelihood function (and consequently the posterior) in a finite grid mainly because each point in the grid is in general expensive computationally, which would require a large amount of time to solve the problem completely. Besides, a considerable hypervolume inside the grid has a very small likelihood and it is thus of no interest, which turns this approach very inefficient.

In order to concentrate in the region of interest where the likelihood is high, we will work with Markov Chain Monte Carlo (MCMC) method. This one aims to generate a set of points that sample the target density, in this case the likelihood (which leads to the posterior) in an efficient manner. MCMC draws random points in parameter space using Markov process: the next point depends on the previous one but not in the entire sampling history. By sampling the posterior with MCMC, one might estimate the quantities of interest for the parameters such as mean, variance, etc.

The main idea behind MCMC method is to sample a new point $\bm{\theta}^{*}$ from the present point $\bm{\theta}$, accepting the new point in the chain with probability $p(\text{acceptance})$ that depends on the ratio between the new and old target density values ($p(\bm{\theta}^{*})$ and $p(\bm{\theta})$, respectively). The distribution which generates a new point based on the present one is called proposal distribution $q(\bm{\theta}^{*}|\bm{\theta})$, generally taken to be a Gaussian that matches the target distribution. On the same manner, one can also define a reversible proposal distribution $q(\bm{\theta}|\bm{\theta}^{*})$, which gives a distribution to reverse the chain (reversible distributions play an important role in MCMC algorithms, such as homogeneity of the chain [114]). A basic MCMC algorithm is the so-called Metropolis-Hastings (M-H) algorithm, where the probability of acceptance of a new point is

[TABLE]

and if the chain is reversible (symmetric proposal distribution) M-H simplifies to Metropolis algorithm

[TABLE]

The basic idea of Metropolis algorithm (which can be generalized for M-H) is that if the new point $\bm{\theta}^{*}$ has a target density $p(\bm{\theta}^{*})$ higher than the target density $p(\bm{\theta})$ of the present point $\bm{\theta}$, the new point is always accepted. However, if $\bm{\theta}^{*}$ has a target density smaller than the target density of $\bm{\theta}$, then the new point is accepted with probability given by (3.21), which means that in the latter case the point might be accepted or not and this makes MCMC exploration in parameter space much more efficient than grid approach. Detailed description of MCMC methods can be found in [115, 116, 114].

In cosmology, Monte Carlo methods are generally applied to data using known softwares such as CosmoMC [79] and Monte Python [80], which are based on Metropolis-Hastings. In this work, as an alternative approach, we will focus on a MCMC method that is mostly used in astrophysics and which has been extensively improved by the community: emcee, a pure-Python implementation which is an affine invariant MCMC ensemble sampler [117]. This code is designed for better exploration of parameter space, with improvements in efficiency and speed when compared to other standard MCMC methods. A faster sampler is preferable because reduces the number of likelihood evaluations, which generally have high computational cost, and thus accuracy is achieved in a minor number of steps.

Affine invariance means that by applying an affine transformation in the proposal distribution one is able to sample the target distribution efficiently even if this one is highly anisotropic, which would render sampling much more expensive computationally. For example, if one considers the skewed probability distribution in $\mathbb{R}^{2}$

[TABLE]

for small $\epsilon$, it follows that traditional MCMC methods have difficulty to sample this distribution. However, if one performs the affine transformation

[TABLE]

the initial challenging problem turns into the easier problem of sampling the distribution

[TABLE]

which is an isotropic density that even a simple Metropolis-Hastings would not have problems to solve. This example motivates an affine invariant sampler: an algorithm which samples target distributions equally well under all linear transformations, being insensitive to covariances among parameters [117]. Although target distributions in cosmology are not highly anisotropic, some parameters might be correlated and emcee is able to take this into account.

The base algorithm for emcee is called stretch move and converges faster than M-H methods, therefore with superior performance [117]. It is based on evolving an ensemble of $K$ walkers $S=\{X_{k}\}$ where the proposal distribution for one walker $k$ depends on the position of the $K-1$ walkers in the complementary ensemble $S_{[k]}=\{X_{j},\forall j\neq k\}$, in which position means the N-dimensional vector in parameter space characterizing each walker. Thus, in order to draw a new position $Y$ from the current position $X_{k}$, the proposal distribution is given by an affine transformation as

[TABLE]

where $X_{j}$ is a walker chosen randomly from $S_{[k]}$ and $Z$ is a random variable sampled from a distribution $g(Z=z)$. This distribution is important because if it satisfies the relation

[TABLE]

then the proposal (3.25) is symmetric, that is, the chain looks the same when run forward or backward in time [114], which can be stated as

[TABLE]

The appropriate acceptance probability for stretch move is

[TABLE]

in which $N$ is the parameter space dimension. Also, the choice made for $g(z)$ is [117]

[TABLE]

where $a$ is a scale parameter adjusted for better performance, set to 2. It is interesting to notice that the idea of using the position of other walkers to guide the next position of a specific walker is based on the fact that the complementary ensemble carries valuable information about the target density, which gives a way of adapting the trial move always aiming the target distribution (this is encoded in the affine invariant proposal), rendering the method very powerful. A more detailed description of affine invariant ensemble samplers can be found in [118].

3.4 Burn-in and convergence

It is expected that chains generated from MCMC methods fairly sample the target distribution once convergence to a stationary distribution has been achieved, that is, no matter where the parameters were initially sampled: the ensemble of walkers are led to the peak of the target distribution and their values oscillate close to the mean value for each parameter. Thus, the first steps in the chain should not be accounted for sampling in order to remove the dependence on the starting point [113]. This period in which walkers still have a “memory” from where they started is called burn-in period.

How long should we run different MCMC chains in order to get good estimates of the parameters? To answer this question, we will adopt a well-known convergence test called Gelman-Rubin convergence criterion [119], which basically compares variance between different chains in order to assess convergence. Suppose we are working with $m$ chains in our MCMC and that $n$ iterations were performed for each chain. Suppose also that burn-in period has already finished. Let $\theta_{i,j}$ represent a point in parameter space in position $i$ and chain $j$. We define the mean of each chain and the mean of all the chains, respectively, as

[TABLE]

where index $i$ represents points in a chain and $j$ runs over the chains. The chain-to-chain variance $B$ and mean within-chain variance $W$ are defined, respectively, as

[TABLE]

The total or true variance $\hat{\sigma}^{2}$ for the estimator $\bar{\theta}$ is defined as the weighted average $\hat{\sigma}^{2}\equiv W(n-1)/n+B$. One can also define the pooled variance

[TABLE]

which takes into account sampling variability of the estimator $\bar{\theta}$ [114]. In the initial sampling, when the distribution is over-dispersed, $B$ is large because chains are also dispersed and variance between them is high. At the same time, $V$ overestimates $W$ because of the $B$ term in (3.32). Also, $W$ is smaller than $\hat{\sigma}^{2}$ because exploration of parameter space is still insufficient (that is, $B$ is still high). However, as $n\rightarrow\infty$, chain-to-chain variance $B$ goes to zero and $W$ approaches true variance $\hat{\sigma}^{2}$, which renders $V/W$ close to 1. Thus, the ratio $\hat{R}\equiv\sqrt{V/W}$, named potential scale reduction factor, is the quantity of interest chosen to monitor convergence, in which a value close to one means that convergence has been achieved. The value of $\hat{R}$ used to stop sampling is chosen according to the level of precision desired, the number of parameters in the theory and the computational cost for every MCMC point, where usual choices are 1.1, 1.03 and 1.01 [120].

3.5 Principal components analysis

Before describing the next tool, it is useful to motivate its application. In order to constrain the dark energy equation of state, we are going to use an interpolation that makes $w(z)$ to be a smooth and continuous function because $w(z)$ is integrated over redshifts whenever $H(z)$ shows up. To this purpose, dark energy equation of state is sampled in six redshift bins, $z_{i}\in\{0.0,0.25,0.50,0.85,1.25,2.0\}$, with corresponding values for the equation of state $w_{i}\,(i=1,2,\dots,6)$. Including more than six bins does not effectively improve the constraints [121]. The last redshift ($z=2$) is chosen according to the fact that this is the epoch in which dark energy started to dominate the expansion of the Universe [96]. We construct $w(z)$ as

[TABLE]

where “spline” means we are using cubic spline interpolation between $z_{i}$’s. The parameters $w_{i}$ are MCMC parameters and must be constrained by data on the same way that other cosmological parameters do (e.g. $\Omega_{b}h^{2}$, $\Omega_{c}h^{2}$, …).

Once chains containing $w_{i}$’s are obtained, after marginalizing over all the remaining parameters, we must face the fact that estimates for $w_{i}$’s are correlated because their covariance matrix is not diagonal. Therefore, one needs to perform a rotation in parameter space, which rewrites $w_{i}$’s in a different basis, in order to obtain new estimates that have a diagonal covariance matrix. The procedure adopted to this purpose is called principal components analysis (PCA) [35, 36, 37] and it is described as follows.

Let us define $\mathcal{F}$ as the Fisher information matrix [27, 113] of a set of estimators $\hat{\theta}$ ($w_{i}$’s in our problem) of parameters $\theta$ measured from observations, assuming that $\mathcal{F}$ can be approximated by the inverse covariance matrix [36] of the estimators $\hat{\theta}$ as

[TABLE]

The decorrelation matrix $\mathcal{W}$ is the matrix that diagonalizes the Fisher matrix $\mathcal{F}$ (or $\mathcal{C}^{-1}$) and thus satisfies the relation

[TABLE]

in which $\Lambda$ is the diagonal matrix containing the eigenvalues of $\mathcal{F}$ and the rows of $\mathcal{W}$ are the eigenvectors of $\mathcal{F}$.

The quatity of interest is $\hat{\theta}^{\prime}=\mathcal{W}\hat{\theta}$ because this new set of $w_{i}$’s is uncorrelated due to the fact that its covariance matrix is, using (3.34) and (3.35),

[TABLE]

which is diagonal and represents the variances of the $w_{i}$’s (elements of $\Lambda^{-1}$ are necessarily positive). The set $\hat{\theta}^{\prime}=\mathcal{W}\hat{\theta}$ is called principal components and the rows of $\mathcal{W}$ are weights (or window functions) that relate the principal components to the initial set $\hat{\theta}$. It is interesting to notice that one might scale the decorrelated quantities $\mathcal{W}\hat{\theta}$ to unity variance by multiplying them by the square root of the diagonal elements of $\Lambda$. Thus, relation (3.35) yields

[TABLE]

As pointed out in [36], there are infinite choices of decorrelation matrices satisfying (3.37). For example, choosing an orthogonal matrix $\mathcal{O}$ and performing a rotation in $\mathcal{W}$ leads to another decorrelation matrix: $\mathcal{V}=\mathcal{O}\mathcal{W}\Rightarrow\mathcal{V}^{\text{T}}\mathcal{V}=\mathcal{W}^{\text{T}}\mathcal{O}^{\text{T}}\mathcal{O}\mathcal{W}=\mathcal{W}^{\text{T}}\mathcal{W}=\mathcal{F}$. There is an interesting choice for $\mathcal{W}$ which renders its weights (rows) positive almost everywhere, with very small negative contributions [36]. This choice is the square root of the Fisher matrix $\widetilde{\mathcal{W}}\equiv\mathcal{F}^{1/2}\equiv\mathcal{C}^{-1/2}$, defined as

[TABLE]

where we normalize $\widetilde{\mathcal{W}}$ such that its rows sum to unity. With this choice, the covariance matrix of the new parameters $\widetilde{\mathcal{W}}\hat{\theta}$ is diagonal [35]. We adopt this choice in our work.

Chapter 4 Results

This chapter displays analyses performed via MCMC sampling using emcee for different combinations of datasets. The first analysis accounts for BAO, cosmic clocks (CC) and supernovae (SN) datasets combined because these ones have lower constraining power than CMB data. The second analysis accounts only CMB data for constraining parameters and the third analysis considers all datasets combined in order to improve constraints. In the last section, a comparison with previous results is performed .

We developed a code in Python for sampling six cosmological parameters, six dark energy equation of state parameters, fifteen “nuisance” parameters from Planck 2015 data and four “nuisance” parameters from JLA dataset. In addition, as we are working with a time-dependent dark energy equation of state $w(z)$, when using CMB data we employ a modified version of CAMB which uses a Parameterized Post-Friedmann (PPF) [122, 123, 124] prescription for dark energy in order to cross the phantom divide [27], i.e. $w=-1$, multiple times. Parameter details are described in Table 4.1.

In results presented in the next sections, we have already removed burn-in and the Gelman-Rubin criterion adopted is $\hat{R}<1.01$ for the set BAO+CC+SN, $\hat{R}<1.6$ for CMB data and $\hat{R}<1.3$ for the set BAO+CC+CMB+SN. These values were chosen according to the size of parameter space and computational cost of each set (which is extremely large for the last two sets). The number of points generated in each MCMC chain is of order $\mathcal{O}(10^{6})$ for BAO+CC+SN and BAO+CC+CMB+SN sets, while for CMB set is $\mathcal{O}(3\times 10^{5})$. Also, we generate estimates of the posterior mean of each parameter of interest together with confidence intervals ($68\%$ and $95\%$) and best fit values using GetDist Python package [125].

4.1 BAO+Clocks+Supernovae

The first investigation refers to BAO, cosmic clocks and supernovae datasets combined. The sampling results are displayed in Fig. 4.1, Fig 4.2. and Table 4.2, which show that even by combining these datasets they do not have enough power to constrain parameters completely. This is due to the fact that datasets used in this sampling are very scarce for high redshifts (only cosmic clocks has information close to $z=2.0$, while the highest redshit for BAO is $z=0.73$ and for supernovae is $z=1.3$). In particular, we tested the last three $w_{i}$’s even by adopting a larger prior than [-10,8] (we even used [-50,48]), but they are also not fully constrained in these conditions, which is why we have not applied PCA to dark energy equation of state using this set. On the other hand, convergence has been achieved in a high level and our cosmological parameter estimates for $\Omega_{b}h^{2}$, $\Omega_{c}h^{2}$ and $H_{0}$ are in good agreement (at $1\sigma$) with latest observations from Planck [49].

4.2 CMB

We shall now analyze the impact of CMB data in the parameters. Results are showed in Fig. 4.3, Fig. 4.4, Fig. 4.5 and Table 4.3. Although this sampling has not fully converged, our cosmological parameter estimates are consistent with Planck data [49] at $1\sigma$ (Fig. 4.3), displaying no deviations when compared to a $\Lambda$CDM model (although CMB measurements are in strong tension with Hubble Space Telescope local measurement for the Hubble rate today, $H_{0}=73.48\pm 1.66\,\text{km}\,\text{s}^{-1}\text{Mpc}^{-1}$ [126], a still unsolved problem in cosmology). When we look at the $w_{i}$’s (Fig. 4.4), we see that Planck data alone seems to be able to constrain all dark energy equation of state parameters, differently from the previous dataset analyzed (BAO+CC+SN). Also, when we look at dark energy equation of state evolution (Fig. 4.5), we see that Planck data alone leads to an agreement with $\Lambda$CDM model at $1\sigma$. It is interesting to notice that, as opposed to the previous section where the first $w_{i}$’s are well constrained, in this case $w_{1}$ parameter has an error bar larger than other $w_{i}$’s (except for $w_{6}$), which is expected since most of CMB information comes from a high redshift ($z\sim 1000$). If sampling was larger and convergence was achieved, then $w_{1}$ would probably have the largest error bar.

4.3 BAO+Clocks+CMB+Supernovae

In the main analysis, we combine all datasets: BAO, cosmic clocks, CMB and supernovae. Our parameter estimates are displayed in Fig. 4.6, Fig. 4.7, Fig. 4.8 and Table 4.4. This sampling has converged more than the one from the previous section, although it still has not fully converged yet. As one can see, all cosmological parameters are constrained and consistent with results from Planck [49] at $1\sigma$ (Fig. 4.6 and Table 4.4). Dark energy equation of state parameters are better constrained than in the previous section because we are now using data (BAO, Clocks, Supernovae) that are also able to constrain parameters at lower redshifts, which can be seen in Fig. 4.7 and Fig. 4.8. On the other hand, this time not all $w_{i}$’s agree with $\Lambda$CDM model at $1\sigma$: there is a small deviation in $w_{3}$ and $w_{4}$, which agree with $\Lambda$CDM model only at $2\sigma$. We believe that this deviation happens because sampling has not fully converged and therefore using a longer chain would solve the problem. As expected, $w_{2}$ is the dark energy equation of state parameter with the smallest error bar, which is explained due to the fact that the greatest amount of information is around the corresponding redshit ($z=0.25$). Another interesting fact is that even Planck data has not enough power to constrain the last dark energy equation of state parameter, $w_{6}$, inside our prior (in the last section, apparently all $w_{i}$’s were constrained, but now we know that this was just because of the lack of convergence). This can be seen in Fig. 4.7 as a “cut” in the 2D plot of $w_{5}\,\times\,w_{6}$, which prevents it of forming a “smooth edge”.

4.4 Comparison of datasets, results and models

As pointed out in the last chapter, the importance of using different datasets is that only the combination of independent data breaks degeneracies and allows us to better constrain parameters. This is very clear when we look at Fig. 4.9, which shows that the combination of all datasets improves the precision of parameter estimation (for example, in the 2D plot $w_{1}\,\times\,w_{2}$ the area estimated using BAO+CC+CMB+SN is smaller than the area estimated using just CMB data alone). We would expect to see the same area in the 2D plot $w_{5}\,\times\,w_{6}$ for both BAO+CC+CMB+SN and CMB data alone because at higher redshits, as we have seen, BAO+CC+SN does not have enough power to constrain $w_{i}$’s. Unfortunately, as both chains are not fully converged, we do not see this happening. The important thing is that BAO+CC+SN combination constrains $w_{i}$’s better at lower redshifts, while CMB data constrains $w_{i}$’s better at higher redshifts.

To picture the implications of our results, in Fig. 4.10 we display our estimations of $w_{i}$’s constraints at $2\sigma$ along with results from Said et al. [121], also at $2\sigma$, and a few alternative models to $\Lambda$CDM. Said et al. paper uses the following datasets combination: WMAP(Wilkinson Microwave Anisotropy Probe [127])+SNLS(Supernova Legacy Survey [128])+BAO(6 degree Field Galaxy Redshift Survey at $z=0.1$ [129] and WiggleZ at $z=0.44,0.60,0.73$ [130])+$H(z)$(dataset from Moresco et al. [131]). The constraints from Said et al. seem better than our results, despite the fact that the datasets used in our work are more recent (and more precise) and should, at first sight, better constrain the $w_{i}$ parameters. The only explanation we could think of is that the parametrization for dark energy equation of state may be somehow affecting the results, because while we use a cubic spline interpolation for connecting the $w_{i}$’s Said et al. uses hyperbolic tangent functions for this purpose. Although we have this difficulty, the results from both estimations can be used to rule out or indicate that theoretical models for $w(z)$ present in the literature might be wrong, as stated in Fig. 4.10.

Conclusions

The main objective of this work was to analyze how current observations coming from baryon acoustic oscillations, cosmic chronometers, cosmic microwave background anisotropies and type Ia supernovae could better constrain dark energy equation of state and maybe indicate deviations from $\Lambda$CDM model, since although this one seems the best theory to explain different observations, it suffers from intrinsic problems that lack of a consistent explanation.

In the fist chapter, we made a review of the background cosmology necessary to understand the basis of the theory. In the second chapter, we described in details the datasets used in our work and the theory behind them in order to picture how data are obtained and what kind of information can be extracted from them. In the third chapter, we concerned about statistical tools to perform the analyses and in the final chapter we explored the results obtained and its consequences to cosmology.

The analysis of the datasets could show us the effect of certain kinds of experiments in parameter estimation. BAO, cosmic chronometers and supernovae data have the power to constrain parameters of low to intermediate redshift, while CMB data constrain parameters better at higher redshifts. However, we did not have very good results for high redshifts even using Planck data, probably from the fact that dark energy does not have much influence in the early Universe in which CMB photons dominate. Maybe combination with other datasets that possess more information at high redshifts (e.g. quasars, redshift space distortion data, CMB lensing) could improve the constrains.

The results obtained in our parameter estimations using MCMC sampling were strongly influenced by the level of convergence achieved in each combination of datasets. The most reliable analysis is the one performed using BAO+CC+SN, because the Gelman-Rubin criterion adopted is very good ($\hat{R}<1.01$) and the corner plot is visually acceptable. However, the other two analyses might be biased in a negative way because Gelman-Rubin criterion adopted is not good and corner plot is obviously not very pleasant, for both CMB data alone and BAO+CC+CMB+SN. Despite this problem, it is customary to say that new physics might be presented when we have at least three or four sigmas of discrepancy, and therefore even if our results are biased we do not see possibility for new physics, which tells us that concordance with $\Lambda$CDM still remains.

Our work showed the inherent difficulties related to studying dark energy equation of state. Besides the fact that the datasets used do not have enough power to constrain $w(z)$ completely, we also need to face the problem of the computational cost of using MCMC methods for sampling, which are generally very expensive in cosmology (in our case because CMB theory generated by CAMB code is extremely complex) and generally demand the use of parallel programming with the help of large computer clusters. As an example, in the sampling with BAO+CC+CMB+SN each MCMC point took 3.5 seconds on average to be generated in the fastest computer we used in this work (CPU: Intel Core i7, cores: 8, memory: SSD 16GB), while we would probably need at least $\mathcal{O}(3\times 10^{6})$ points to finally achieve convergence. The effect of these problems are seen as huge error bars in the dark energy equation of state parameters (the last $w_{i}$’s), turning the evolution of the equation of state very hard to predict.

Using emcee instead of standard parameter estimations codes for cosmology, such as CosmoMC and Monte Python, was made with the intuit of confirming that we could recover old results even with a different sampler. This was successfully achieved, since we tested simple models such as $\Lambda$CDM and our cosmological parameter estimation really matched results found in literature. This is important to give us much more confidence in theses codes, which are generally used as a “black box” in most applications to cosmology.

Due to the computational cost problem, it is worth mentioning some alternative tools to MCMC sampling as possible extensions to this work. One of them is Gaussian Processes [135], which provides a surprisingly computationally effective method to reconstruct functions (this has already been applied to $w(z)$ [136, 137]) from data at low computational cost. The other one is Variational Inference [138], a tool which has been tested in astrophysics problems [139] and converts the inference problem into an optimization problem, performing parameter estimation much faster than standard MCMC methods.

As a way of helping cosmology students interested in MCMC methods, all the codes developed in this dissertation can be found at Github in reference [140].