Reconstruction of primordial fluctuations from galaxy catalogs with higher order hamiltonian sampling

Abstract

Las diferentes estructuras que observamos en el Universo son el resultado de la evolución de las semillas primordiales, cuyo origen se cree que proviene de fluctuaciones cuánticas. En un período de inflación cósmica, estas fluctuaciones microscópicas fueron amplificadas en una fase de muy rápida expansión apenas 10−36 segundos después del Big Bang, dando lugar a un Universo prácticamente homogéneo y altamente gausiano. Sin embargo, en algunas regiones, la densidad era ligeramente superior a la media, produciéndose inestabilidades gravitacionales, haciendo que estos picos de densidad atrajesen más materia de sus alrededores. De esta manera, se originó la compleja red cósmica que observamos en surveys galácticos, la cual está compuesta por la agrupación de cúmulos de galaxias en nodos y filamentos, y grandes vacíos cósmicos. La evolución de las fluctuaciones primordiales se puede describir con las ecuaciones de un fluido, haciendo uso de la teoría de perturbaciones en un sistema de referencia Lagrangiano (teoría de perturbaciones Lagrangiana), que resulta ser muy útil para describir la formación de la red cósmica. Dentro de la estructura a gran escala podemos diferenciar dos regiones: aquella a grandes escalas, donde la evolución de las fluctuaciones primordiales sigue aproximadamente un crecimiento lineal, y a escalas intermedias y pequeñas, donde la gravedad es no lineal. En este régimen no lineal, la gravedad acopla distintos modos del espacio de Fourier por ser una interacción de largo alcance y altera la distribución gausiana primordial, introduciendo desviaciones importantes de una distribución simétrica normal. Esto hace muy complejo el análisis de la estructura a gran escala, requiriéndose, para ello, técnicas avanzadas que nos permitan caracterizar la red cósmica. En este contexto, escogemos una descripción Bayesiana, construyendo la función de probabilidad posterior que se quiere muestrear como el producto de un prior, describiendo la distribución continua de la materia oscura, y una likelihood, describiendo la distribución discreta de las galaxias. Esto nos permite, en un principio, obtener los campos de materia oscura que son compatibles con un catálogo de galaxias. Para ello, se requiere elegir unas funciones de probabilidad concretas, incluyendo un modelo de formación de estructuras que relacione el campo de densidad primordial con el evolucionado gravitacionalmente y una técnica concreta de muestreo estadístico. Comenzamos, en primer lugar, considerando un problema simplificado, en el cual asumimos un modelo de formación de estructuras donde el logaritmo del campo de densidad de la materia oscura sigue una distribución gausiana (modelo lognormal). La distribución de galaxias se asume que sigue una distribución poisoniana. En este proyecto se ha empleado el código ARGO, que utiliza una técnica de muestreo híbrida de cadenas de Markov basada en Hamiltonian sampling. Esta técnica tiene la ventaja de ser capaz de muestrear funciones de probabilidad arbitrarias (incluyendo no gausianas) y de forma eficiente, usando los gradientes de los potenciales determinados por la función posterior. Para ello, se define la energía potencial como el negativo del logaritmo de la función posterior que queremos muestrear y la energía cinética referida a unos momentos que son artificialmente introducidos en el método para posibilitar el muestreo aleatorio. Estos tienen la ventaja de ser descritos por una simple gausiana. El método consiste en muestrear la probabilidad conjunta de la señal que queremos reconstruir (los campos de materia oscura), que en la analogía hamiltoniana representan las posiciones; y los momentos. La marginalización del muestreo con respecto a los momentos, se logra despreciándolos de la iteración anterior en la cadena de Markov, introduciendo nuevos valores de los mismos en cada iteración de acuerdo a una gausiana con una matriz de covarianza determinada por la masa. El sistema de posiciones y momentos se evoluciona de acuerdo con las ecuaciones de Hamilton, que se resuelven de manera discreta y se van aceptando o rechazando en un paso de MetropolisHastings. Para la conservación del volumen del espacio de las fases se usa un esquema discreto de Leapfrog, hasta ahora considerado sólo hasta segundo orden. En este trabajo, hemos implementado por primera vez en el campo de la cosmología estructura a gran escala, esquemas de Leapfrog de mayor orden. En concreto, nos enfocamos en el cuarto orden, aunque el marco teórico que introducimos admite órdenes mayores. Esto nos permite evolucionar el sistema hamiltoniano con pasos mucho más amplios que a segundo orden y obtener una mayor aceptación de las iteraciones. Todo esto se traduce en un incremento de la eficiencia del método por un factor 18 en la convergencia de la cadena de Markov, en términos de tiempo computacional, requiriendo dos órdenes de magnitud menos iteraciones que en el caso anterior. Además, obtenemos una longitud de correlación entre los campos de materia muestreados de unas 10 iteraciones frente a las ∼ 300 requeridas en el método a segundo orden. Esto implica que podemos obtener muestras independientes cada 10 iteraciones una vez se alcanza la convergencia en unas 30 iteraciones. En un segundo estudio, usamos el código BIRTH, que incluye teoría de perturbación lagrangiana para la reconstrucción de las fluctuaciones primordiales, y lo aplicamos a un problema realista con un catálogo de galaxias emulando las galaxias rojas luminosas CMASS de BOSS final data release. Este catálogo incluye evolución cósmica en el cono de luz, velocidades peculiares de las galaxias, geometría del survey y función de selección radial. Para este código, obtenemos los mismos resultados en términos de eficiencia global que en el caso simplificado con el Leapfrog de cuarto orden. A partir de una reconstrucción del campo primordial de fluctuaciones obtenida con el código BIRTH, realizamos una simulación de N-cuerpos con el código PKDGRAV3 hasta z = 0.57, que corresponde con el redshift medio del catálogo de CMASS. Comprobamos que las estructuras finales de materia oscura guardan una relación espacial muy similar a las condiciones iniciales obtenidas a partir del catálogo. Esto nos muestra una de las aplicaciones interesantes del método desarrollado en este trabajo de fin de máster, dado que nos permitiría una comparación más directa entre los datos y los modelos. La eficiencia obtenida con el método Hamiltoniano desarrollado en este trabajo nos permitiría acometer muchos más estudios en un análisis Bayesiano global, tales como incluir un tratamiento totalmente no lineal de la gravedad en el proceso iterativo de reconstrucción.

The different structures we observe in the Universe result from the cosmic evolution starting from some primordial seeds. These are believed to have originated from some quantum fluctuations. During the inflationary epoch, they were stretched and frozen just 10−36 seconds after the Big Bang, giving rise to an almost homogeneous Universe closely Gaussian distributed. However, these fluctuations were slightly larger to the mean in some regions, causing gravitational instabilities, so that density peaks were attractors to their surrounding matter. In this way, a complex non-Gaussian cosmic web emerged, composed by the clustering of galaxies in knots, filaments and sheets; and cosmic voids, where the mean density is very low. This makes the analysis of the large scale structure very complex, so advanced techniques are required to characterize this cosmic web. We choose, in this context, a Bayesian framework, which allows us to approach the problem from a robust statistical perspective. We start considering a simplified problem, assuming a lognormal prior model for the density field and a Poisson likelihood for the galaxy distribution. Using the ARGO code we are able to sample the posterior resulting from the product of the prior and likelihood with a hybrid Markov Chain Monte Carlo method: the Hamiltonian sampling technique. In this way, we obtain a dark matter field in each Markov iteration, which is statistically compatible with the input galaxy distribution, given the model. In the Hamiltonian analogy, the dark matter fields represent the positions, and the momenta are artificially introduced to enable sampling arbitrary posterior distribution functions. They are sampled from a Gaussian distribution with a given mass. The Hamiltonian system is evolved through the discretized equations of motions, where the final positions and momenta are accepted or rejected according to a Metropolis-Hastings step. Until now, a second order Leapfrog scheme has been used to solve this problem. We explore in this thesis higher order discretizations, focusing, in particular, on fourth order Leapfrog. The result is a factor 18 faster convergence in terms of computational time, with respect to the original second order Leapfrog scheme. Moreover, we obtain a correlation length of about 10 iterations, as opposed to ∼ 300 with the old scheme. This implies that we can obtain independent samples each 10th iterations once we reach convergence, which is about 30 iterations, two orders of magnitude less than with the second order algorithm. In a second study, we used the BIRTH code, which includes Lagrangian perturbation theory for the reconstruction of the primordial fluctuations. We apply this code, which also relies on Hamiltonian sampling, to a realistic mock galaxy catalog resembling the BOSS-CMASS final data release. This catalog includes cosmic evolution in the light-cone, peculiar velocities, survey geometry, radial selection function. We obtain the same results in terms of global efficiency, as in the simplified case study, when we implement the fourth order Leapfrog scheme. Using one of the reconstructed primordial fluctuations with the BIRTH code, we make a constrained N-body simulation with PKDGRAV3 to z = 0.57, which corresponds to the mean redshift of CMASS galaxies. We observe a high spatial resemblance between the obtained dark matter structures and the initial conditions obtained from the catalog. This shows one of the interesting applications of the method developed in this thesis, as it permits a more direct comparison between data and models. The efficiency of the higher order Hamiltonian method developed here will enable us to tackle many more complex studies within a global Bayesian framework, such as including a full non-linear gravity solver within the iterative reconstruction process.