Formación de estructura en distintos modelos cosmológicos

Trabajo de Fin de Grado Doble Grado en Física y Matemáticas Formación de estructura en distintos modelos cosmológicos Autora: Beatriz González González Tutor: José Oñorbe Bernis 2 de julio de 2021 A mi madre, a mi padre, a mi hermana. La astronomía obliga al alma a mirar hacia arriba y nos lleva de este mundo a otro. Platón (427-347 a. C.) The history of astronomy is a history of receding horizons. Edwin P. Hubble (1889-1953) Resumen En el paradigma actual de formación de estructuras cosmológicas a gran escala, la interacción gravitatoria hace que las partículas de materia oscura se agrupen y colapsen en objetos virializados, conocidos como halos. Estos crean los pozos de potencial donde se forman las galaxias. En este trabajo se estudia la función de masa de halo como propiedad fundamental del modelo cosmológico. Esta función da la densidad numérica de halos por intervalo de masa. Para ello, se ha desarrollado un código numérico para calcular las predicciones de la función de masa de halo de los modelos analíticos y semianalíticos más relevantes en el área. Además, se analizan catálogos de halos utilizando los resultados de un conjunto de simulaciones de N cuerpos a redshift 0. Se muestra que los modelos analíticos, que asumen una evolución lineal de las perturbaciones, tienen una precisión limitada. Se estudia la formación de halos en regiones del universo con distinta densidad, encontrando una dependencia entre este parámetro y la función de masa de halo, especialmente para los halos más masivos. Finalmente, en base a estos resultados, se discute la posible parametrización de la dependencia de la función de masa de halo respecto a la densidad como una simple dependencia temporal de una región con densidad media. Índice general I Introducción: Fundamento teórico 1 1. De la armonía de las esferas al universo en expansión 3 2. El modelo cosmológico 6 2.1. Materia y energía del universo . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Parámetros fundamentales del modelo . . . . . . . . . . . . . . . . . . . . . . . 9 2.3. Modelo ΛCDM.................................... 11 3. Teoría de formación de estructura 12 3.1. Fluctuacioneslineales................................ 12 3.2. La relación entre la función de correlación y el espectro de potencias . . . . . 15 3.3. Colapsoesférico................................... 19 4. La función de masa de halo 21 4.1. Formalismo de Press-Schecter . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.2. Más allá del formalismo de Press-Schechter: modelos semianalíticos . . . . . . 24 II Simulaciones 27 5. Simulaciones y resultados 29 5.1. Simulaciones cosmológicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2. SimulacionesLegacy ................................ 31 5.3. Resultados ...................................... 33 5.3.1. Filtrado de subhalos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.3.2. Errores: límites de resolución e incertidumbres . . . . . . . . . . . . . 34 5.3.3. Simulaciones estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.3.4. SimulacionesZoom............................. 40 vii 2 El modelo cosmológico La cosmología moderna se construye a partir del principio cosmológico, según el cual no existe ninguna localización privilegiada en el universo. A grandes escalas, del orden de cientos de megapársecs1, el universo es homogéneo e isótropo, esto es, presenta las mismas propiedades en todas las localizaciones y hacia todas las direcciones. 2.1. Materia y energía del universo La métrica más general que satisface la homogeneidad y la isotropía requeridas por el principio cosmológico es la métrica de Friedmann-Lemaître-Robertson-Walker (FLRW), caracterizada por el elemento de línea en unidades en las que c= 1 [6]: ds2=−dt2+a2(t)dr2 1−Kr2+r2(dθ2+ sin2θdφ2)[1] El factor de escala a(t)y la curvatura espacial Kson las propiedades fundamentales de la métrica FLRW. En concreto, la curvatura indica si el universo tiene una geometría abierta, plana o cerrada, correspondientes a valores de K=−1,0ó1respectivamente. Por su parte, a(t)describe la expansión o contracción del universo. El factor de escala indica la proporción de la distancia propia entre dos objetos en un tiempo trespecto a dicha distancia en un tiempo 1El pársec es una unidad de longitud utilizada en astronomía: 1pc = 3.2616 años luz = 3.0857 ·1016m 6 de referencia t0: a(t) = d(t) d(t0)[2] Este tiempo de referencia t0se fija en el presente, de manera que a(t0)=1. La distancia d(t0)recibe el nombre de distancia comóvil. Esta permanece constante con el tiempo si los dos objetos se mueven al ritmo de la expansión del universo, a diferencia de la distancia propia d(t). En cosmología se suele emplear lo que se conoce como corrimiento al rojo o, usando su término inglés, ‘redshift’. En un universo en expansión, la longitud de onda de radiación electromagnética recibida λoincrementa respecto a la emitida λ(t)en tiempo un t, esto es, se desplaza hacia el rojo en el espectro electromagnético. El redshift cósmico se mide con el parámetro z. Si a0se corresponde con el factor de escala en el presente, se tiene que: 1 + z=λ λ(t)=a0 a(t)[3] Llegado este punto, es importante hacer una aclaración en torno al redshift, el factor de escala y el tiempo. El redshift zse mide directamente de las galaxias y los cuásares lejanos. Asumida la teoría del Big Bang, el factor de escala representa una interpretación del redshift como expansión del universo, que es válida para cualquier modelo cosmológico, en el contexto de la relatividad general y de un universo homogéneo e isótropo. Por último, para relacionar el factor de escala con un tiempo ten años concretos, es necesario especificar el modelo cosmológico de manera exacta. Por ello, uno de los principales objetivos de la astrofísica teórica es determinar el modelo cosmológico que explique las observaciones del universo y permita fijar la escala de tiempo. La base teórica sobre la que se estudian los posibles modelos cosmológicos no es otra que la Teoría de la Relatividad General de Einstein. Esta conduce a que la geometría del espaciotiempo y, en consecuencia, la interacción gravitatoria de la materia dependen del contenido de materia-energía del universo. En concreto, las ecuaciones de campo de Einstein vienen dadas por [7], [8]: Rµν −1 2gµνR+ Λgµν = 8πGTµν [4] 7 Donde Rµν es el tensor de curvatura de Ricci; R, el escalar de curvatura de Ricci; Λ, la constante cosmológica; gµν, el tensor métrico; G, la constante de gravitación universal y Tµν el tensor momento-energía. De manera resumida, la ecuación anterior relaciona la geometría del espacio-tiempo, caracterizada por el tensor métrico gµν, con la densidad y flujo de energía y momento del espacio-tiempo dada por Tµν. En base al principio cosmológico, resulta adecuado describir el universo como un fluido ideal de densidad de energía ρ(t)y presión p(t), donde un observador que se mueva conjuntamente con dicho fluido (es decir, un fluido comóvil) verá el universo perfectamente isótropo. Con esto y a partir de la ecuación 4 de Einstein se llega a las ecuaciones de Friedmann [2]: ˙a a2 =8πG 3ρ+Λ 3−K a2[5] ¨a a=−4πG 3(ρ+ 3p) + Λ 3[6] Al igual que la ecuación 4, estas ecuaciones describen la expansión del universo y muestran que la evolución del universo está gobernada por su contenido de materia-energía [6]. En general, la presión y la densidad de un fluido están relacionadas por una ecuación de estado p=p(ρ). Para los fluidos estudiados en cosmología, se puede asumir que la presión es linealmente proporcional a la presión: p=αρ . En el caso de los tres fluidos fundamentales que se estudian en este trabajo (materia, radiación y energía de vacío), la constante αtoma los siguientes valores [6]. Para la materia, una buena aproximación a escalas grandes es considerarla acolisional, lo que significa que su presión es nula y α= 0. Para la radiación, que sigue un régimen relativista, se tiene que α= 1/3. Por último, para la energía de vacío, asociada a fluctuaciones cuánticas de vacío, la constante de proporcionalidad resulta α=−1. Una de las claves para entender la evolución del universo es que la densidad de cada fluido depende de su ecuación de estado. La conservación de la energía y la expansión adiabática del universo en la métrica FLRW conduce a la siguiente ecuación [9]: 8 ˙ρ+ 3 ˙a a(p+ρ) = 0 [7] De aquí, se puede deducir que la densidad de energía de materia decae como ρM∼a−3, mientras que la densidad de energía de radiación lo hace como ρR∼a−4. En el caso de la energía de vacío, su densidad de energía permanece constante: ρvac ∼const. La física cuántica predice una energía de vacío de unos 100 órdenes de magnitud menor que la observada [10]. Esto es uno de los grandes problemas de la física moderna y lleva a considerar un tipo de energía más abundante en el universo: la energía oscura, ligada a la ya mencionada constante cosmológica Λ. Esta conexión entre la energía de vacío y la constante cosmológica se remonta a mediados del siglo pasado [11]. Aún a día de hoy, la energía oscura entraña uno de los grandes misterios de la ciencia. 2.2. Parámetros fundamentales del modelo Uno de los parámetros más importantes de la cosmología es el parámetro de Hubble H(t), que describe el ritmo de expansión del universo para cierto tiempo y se define como: H(t) = ˙a(t) a(t)[8] Se denota H0al valor correspondiente al tiempo presente y recibe el nombre de constante de Hubble, constante de proporcionalidad en la ley de Hubble, que establece una relación entre la distancia entre galaxias y la velocidad a la que se alejan. Históricamente [12], el valor de H0 se suele expresar en función de otro parámetro h: H0= 100 hkm s−1Mpc−1[9] En los últimos años, se ha medido a partir de diversos métodos que el parámetro hse encuentra alrededor de un valor de 0.7. Concretamente, los resultados obtenidos a partir de los datos recogidos por las sondas WMAP9 [13] y Planck [14] determinaron un valor de h= 0.6950. 9 Por otro lado, se define la densidad crítica como la densidad media necesaria para que el universo frene su expansión (debido a la atracción gravitatoria) y permanezca estable: ρc=3H2 0 8πG [10] Usando el valor de la constante de Hubble que dan las observaciones (ver más arriba), esta densidad crítica es equivalente a unos 6 protones por metro cúbico, una cantidad a priori muy baja pero considerablemente superior a la densidad actual de materia que dan las observaciones, cercana 0.25 protones por metro cúbico [15]. ¿Cómo se explica esta diferencia? Primero, es conveniente definir el parámetro de densidad Ω0, igual al cociente entre la densidad actual y la densidad crítica ρc: Ω0=ρ(t0) ρc [11] Este parámetro sirve para determinar la geometría del universo: si Ω0>1, el universo será cerrado; si Ω0= 1, plano; y si Ω0<1, abierto. De particular interés para los desarrollos siguientes resultan los parámetros de densidad relativos a día de hoy de la materia ΩM, radiación ΩR, energía oscura ΩΛy curvatura ΩK, calculados como Ωi=ρi/ρ0. ΩM=8πGρM 3H2 0 ΩR=8πGρR 3H2 0 [12] ΩΛ=Λ 3H2 0 ΩK=−K a2 0H2 0 [13] A partir de estas densidades relativas y de la primera ecuación de Friedmann, se puede expresar el parámetro de Hubble en función su valor actual. Así, se llega a la siguiente expresión general en función del redshift z: H(z) = H0ΩR(1 + z)4+ ΩM(1 + z)3+ ΩK(1 + z)2+ ΩΛ1/2[14] 10 En definitiva, la idea principal a la hora de construir un modelo cosmológico es que el universo está caracterizado por su curvatura y formado principalmente por radiación, materia y energía oscura. Sus proporciones determinan el destino del universo: podría expandirse eternamente, colapsar o incluso tender asintóticamente a un estado límite entre estos dos casos. 2.3. Modelo ΛCDM En el marco general de modelo cosmológico, el conocido como Modelo Estándar o Modelo ΛCMD es el que presenta mejor concordancia con los datos observacionales [16]. Basado en tan solo seis parámetros fundamentales, se trata de un modelo sencillo, capaz de explicar de manera satisfactoria un gran número de observaciones a muchas escalas y distintos tiempo. Múltiples estudios e investigaciones, como por ejemplo los realizados a partir de los datos recogidos por las sondas WMAP9 [13] y Planck [14], indican que nuestro universo actual tiene una geometría plana (ΩK= 0) y está compuesto esencialmente de energía oscura (ΩΛ≈0.7) y materia (ΩM≈0.3), pues el término de radiación es despreciable. En el caso de la materia, cabe destacar que tan solo el 16 %, aproximadamente, es materia bariónica, esto es, materia formada por protones, neutrones y electrones. La gran mayoría es materia oscura, un tipo de materia acolisional (solo interactúa gravitacionalmente de forma significativa) que no emite radiación electromagnética y no se puede observar. Según el modelo ΛCDM, la materia oscura es fría, o dicho de otro modo, no relativista. Por todo ello, se suelen emplear los parámetros de densidad relativos a la materia bariónica y a la materia oscura, que toman valores de Ωb≈0.05 Ωm≈0.25, respectivamente. Esto explica la diferencia presentada más arriba entre la densidad crítica y la densidad de materia bariónica del universo. Esta última representa menos del 5 %del contenido del universo. En definitiva, la materia que se puede observar directamente constiye una parte muy pequeña del contenido de nuestro universo a día de hoy. Este está dominado por la expansión acelerada producida por una forma de energía desconocida: la energía oscura. 11 3 Teoría de formación de estructura La teoría de formación de estructura trata de explicar la transformación del universo homogéneo inicial, que muestra el fondo cósmico de microondas, en nuestro universo presente, donde la materia se reúne en galaxias y cúmulos de galaxias. Este proceso tiene su origen en pequeñas fluctuaciones primordiales de la densidad. Los principales modelos teóricos apunta a que estas fluctuaciones derivan de pequeñas inhomogeneidades en el plasma primigenio, originadas por fluctuaciones cuánticas ([17] y [6]), lo que parece concordar con las fluctuaciones observadas en el Fondo Cósmico de Microondas [18]. Con el tiempo, estas perturbaciones iniciales crecieron gracias a la inestabilidad gravitacional, dando lugar a objetos colapsados conocidos como halos de materia, que crearon los pozos de potencial donde se formaron las galaxias. 3.1. Fluctuaciones lineales Considerando que estas perturbaciones iniciales son pequeñas, se puede trabajar en el régimen lineal y obtener una visión general de la formación de estructura. Para ello, se define el contraste de densidad, una función adimensional que describe la sobredensidad en una región del espacio y en una época determinados. δ(x, t) = ρ(x, t)−¯ρ(t) ¯ρ(t)[15] 12 En esta ecuación, ρ(x, t)es el campo de densidad de materia en la posición xy el tiempo t y¯ρ(t)es la densidad media del universo de fondo, dependiente del modelo cosmológico. Gracias a la simplificación del régimen lineal, se puede emplear un tratamiento newtoniano para resolver la evolución de las perturbaciones iniciales en un universo en expansión bajo las siguiente condiciones. La materia del universo debe encontrarse en la forma de un fluido predominantemente autogravitante y acolisional, que, además, sea no relativista. Este es el caso de la materia oscura fría y la materia bariónica a grandes escalas, donde la hidrodinámica pierde importancia. En este contexto, las ecuaciones que gobiernan la densidad ρy la velocidad ude dicho fluido ideal son: Ecuación de continuidad: ∂δ ∂t +∇·[(1 + δ)u]=0 [16] Ecuación de Euler: ∂u ∂t + 2H(t)u+ (u·∇)u=−∇Φ a2[17] Ecuación de Poisson: ∇2Φ = 4πG¯ρa2δ[18] Cabe destacar que estas ecuaciones están expresadas en función de coordenadas comóviles x, esto es, coordenadas que se expanden al mismo ritmo que el universo. Así, si ryvson, respectivamente, la posición propia y la velocidad del fluido, se tiene que x=r/a yu= v−(˙a/a)r=a˙x, donde a=a(t)es el factor de escala. Por otra parte, Φ(x)es el potencial gravitatorio y H(t), el parámetro de Hubble, introducido en el Capítulo 2. Particularizando las ecuaciones anteriores para el caso de fluctuaciones de densidad pequeñas, δ << 1, se obtiene la aproximación lineal. Considerando únicamente los términos de primer orden de δyu, se obtiene la siguiente ecuación: ∂2δ ∂t2+ 2H(t)∂δ ∂t = 4πG¯ρδ [19] Se trata de una ecuación diferencial homogénea y lineal de segundo orden en el tiempo, lo que tiene dos consecuencias [19]. En primer lugar, la solución puede expresarse como la suma de dos funciones: δ(x, t) = δ+(x, t) + δ−(x, t). En segundo lugar, dado que la ecuación solo tiene derivadas respecto al tiempo, la evolución temporal será independiente de la posición x 13 y cada término se puede separar en una parte espacial y otra temporal: δ(x, t)=∆+(x)D+(t)+∆−(x)D−(t)[20] Aquí, D+(t)yD−(t)son los denominados modos de crecimiento y decaimiento y ∆+(x) y∆−(x)representan la distribución espacial del campo de materia inicial. La expresión para los modos de crecimiento y decaimiento dependerá de los parámetros del modelo cosmológico (Capítulo 2). Las estructuras crecerán de distinta manera en cada universo FLRW: el contenido de materia-energía del universo será clave en su evolución. Por ejemplo, para un universo Einstein-de-Sitter (EdS) donde ΩM= 1 y la curvatura es plana, se tiene que D+∝t2/3∝a(t)yD−∝t−1. En este caso, el contraste de densidad crece en armonía con la expansión del universo. Teniendo en cuenta que las fluctuaciones primordiales de densidad sucedieron hace millones de años, para estudiar la formación de estructuras en el momento presente solo se considera el modo de crecimiento, que pasa a denominarse simplemente D(t). En consecuencia, si ties el tiempo en el origen del Big Bang, la expresión para el contraste de densidad queda como: δ(x, t) = δ(x, ti)D(t) D(ti)[21] Para un universo plano con constante cosmológica, es decir, donde ΩΛ6= 0 yΩM<1, el modo de crecimiento puede expresarse mediante la siguiente integral [20]: D(z) = H(z) H0Z∞ z dz0(1 + z0) H3(z0)Z∞ 0 dz0(1 + z0) H3(z0)−1 [22] donde el factor entre corchetes representa la normalización para que D= 1 en z= 0. En definitiva, si se consiguen medir las fluctuaciones de densidad hoy en día y estas se comparan con las obtenidas a partir del Fondo Cósmico de Microondas, se deberían poder constreñir los parámetros cosmológicos del modelo cosmológico que rige nuestro universo. 14 3.2. La relación entre la función de correlación y el espectro de potencias Suponiendo que las fluctuaciones iniciales tienen un origen cuántico, el modelo actual más simple predice que, desde un punto de vista matemático, el contraste de densidad δ(x, ti)era un campo aleatorio gaussiano homogéneo e isótropo en el origen del Big Bang. En consecuencia, se define la función de correlación de dos puntos como ξ(r) = hδ(x1)δ(x2)i, que, debido a la homogeneidad y la isotropía impuestas por el principio cosmológico, solo depende de la distancia entre dichos puntos r=|x1−x2|. El significado físico de la función de correlación es el siguiente: si ξ(r)<0, el campo de densidad estará menos concentrado que la media; mientras que si ξ(r)>0, estará más concentrado. A menudo, resulta adecuado considerar δ(x, t)como una superposición lineal de ondas simples, lo que, en un universo plano como el nuestro, se consigue mediante el análisis de Fourier. Así, la transformada de Fourier del contraste de densidad viene dada por: ˜ δ(k, t) = 1 (2π)3/2Zdxδ(x, t)eik·x[23] A partir de aquí, se puede calcular la función de correlación del campo que representa la transformada de Fourier: h˜ δ(k)˜ δ∗(k’)i=h1 (2π)3Zδ(x)eik·xdxδ(y)e−ik’·ydyi[24] Intercambiando el orden entre las integrales y el promedio, se tiene que: h˜ δ(k)˜ δ∗(k’)i=1 (2π)3Zhδ(x)δ(y)ieik·xdxe−ik’·ydy=1 (2π)3Zξ(y−x)eik·xdxe−ik’·ydy[25] 15 En el capítulo anterior, se estudió que, en general, los modelos actuales conducen a que el contraste de densidades suavizado (ecuación 31) tiene una distribución gaussiana. A partir de las propiedades estadísticas estudiadas de las fluctuaciones de densidad, la relación entre la función de correlación y el espectro de potencias (ecuaciones 28 y 29) permite obtener su varianza σM, dada por la ecuación [32]. Así, sin más que usar la definición de la función de distribución gaussiana, se obtiene la siguiente expresión para la probabilidad: P(δM≥δc) = 1 σM√2πZδc ∞ dδ0exp −δ02 2σ2 M[37] Llegado este punto, el formalismo de Press-Schechter se encuentra con un problema: la probabilidad P(> M)es menor a la esperada. Al tomar el límite M→0, se obtiene que P(> M) = 0.5, lo que significa que solo la mitad de la materia del universo se encuentra en halos colapsados. Esto se debe al hecho de que, al ser el contraste de densidad una distribución gaussiana, la probabilidad de que la densidad sea menor a la de la media es P(δ < 0) = 0.5. Como solo las regiones inicialmente ‘sobredensas’ acaban en objetos colapsados, no se tienen en cuenta las zona ‘subdensas’. Sin embargo, la evolución muestra que las regiones vacías se hacen más vacías y las densas, más densas, lo que implica que la materia de las regiones ‘subdensas’ acaba cayendo en objetos masivos. Para salvar este obstáculo, se multiplica por un factor de 2 (conocido como factor compensatorio o, en inglés como ’fudge factor’ [27]) . A continuación, se muestra la deducción de la función de masa de halo en el formalismo de Press-Scheter. Para ello, se toman dos vías para definir la densidad numérica de halos cuyas masas se encuentran en un rango de [M, M +dM]. El objetivo es llegar a una igualdad que permita encontrar una expresión para la función de masa de halo. Por un lado, bajo la hipótesis de ergodicidad, la probabilidad P(> M)también puede verse como la fracción de materia del universo que se encuentra en un objeto colapsado de masa > M. Haciendo uso de esta definición, ∂P ∂M dM representa la fracción de masa encerrada en halos con una masa entre MyM+dM. Al multiplicar por la densidad media ¯ρ, se calcula la 22 masa total por unidad de volumen que se encuentra en dichos halos. Sin más que dividir por la masa, se obtiene el número de halos con masas en un rango de [M, M +dM]por unidad de volumen. Por otro lado, se denomina n(M, z)a la densidad numérica de halos con masa > M, es decir, la cantidad acumulada de halos con masa mayor que Mpor unidad de volumen. De este modo, la densidad numérica de halos con masa en un intervalo [M, M +dM]se puede expresar como dn dM dM. Recordando el primer razonamiento y añadiendo el factor 2, se tiene la siguiente igualdad. dn dM dM = 2 ¯ρ M ∂P ∂M dM [38] La función de masa de halo de Press-Scheter, que calcula el número de halos por unidad de volumen e intervalo de masa, viene dada por dn dM . Despejando y derivando en la ecuación 37, se deduce que, si ν=δc σM: dn dM =−¯ρ Mr2 πexp −ν2 2ν σM ∂σM ∂M [39] En este punto, es importante remarcar que esta fórmula permite relacionar la función de masa de halo con la varianza del campo de perturbaciones (dada por la ecuación 32) y, en consecuencia, con el espectro de potencias P(k). Tal y como se explica en el Capítulo 3, el espectro de potencias es una de las bases del modelo cosmológico actual. De aquí se deduce, que la función de masa de halo es también una característica esencial del modelo cosmológico. Más adelante, se verá en detalle cómo de precisa es esta fórmula. 23 4.2. Más allá del formalismo de Press-Schechter: modelos semianalíticos A partir de ciertas suposiciones sencillas, el formalismo de Press-Schechter consigue dar una expresión analítica para la función de masa de halo. Sin embargo, este modelo está limitado por las aproximaciones lineales y, a lo largo de los años, se ha comprobado que sus predicciones pierden precisión al compararlas con las obtenidas a partir de simulaciones cosmológicas. En vistas a mejorar dicha precisión de la función de masa de halo, se han ido desarrollando nuevos modelos semianalíticos donde se utilizan los resultados de simulaciones para calibrar la función de masa de halo. En el Capítulo 5 se hablará más en profundidad de las simulaciones cosmológicas y su uso para confrontar los resultados teóricos. A raíz de esto, se clarificará el porqué del desarrollo de modelos semianalíticos. Esta sección se limita a dar una pequeña introducción a dos de los primeros modelos semianalíticos y a las nuevas expresiones de la función de masa de halo. En este contexto, se suele utilizar la denominada función de masa o de multiplicidad f(σ)para expresar la función de masa de halo. Concretamente, esta función representa la fracción de masa que ha colapsado para formar halos por unidad de intervalo en ln σ−1[28]: f(σ) = M ¯ρ dn dln σ−1[40] Con esto, se puede reescribir la función de masa de halo de manera general, dependiente de la función de masas: dn dM =¯ρ M2f(σ)dln σ−1 dln M[41] En el caso de la función de masa de halo obtenida en el formalismo de Press-Schechter, la función de multiplicidad tiene la siguiente forma. Recordando que ν=δc σM: 24 fPS(ν) = r2 πνexp −ν2 2[42] Por otro lado, en las últimas décadas, se han dado múltiples definiciones de esta función de masas motivadas por resultados numéricos. En 1999, Sheth y Tormen presentaron la siguiente función de multiplicidad [29]: fST (σ) = Ar2a π1 + σ2 aδ2 cpδc σexp −aδ2 c 2σ2[43] Donde A= 0.3222,a= 0.707 yp= 0.3. Unos años más tardes, Jenkins [28], que trabajó con diferentes simulaciones cosmológicas con materia oscura fría a redshifts entre z= 5 yz= 0, encontró que la siguiente fórmula proporcionaba un buen ajuste, independientemente de la cosmología: fJ(σ) = 0.315 exp −|ln σ−1+ 0.61|3.8[44] 25 II Simulaciones 5 Simulaciones y resultados Desde el origen de la astrofísica, las observaciones astronómicas han sido, junto con las predicciones teóricas, la base fundamental de nuestro conocimiento sobre el cosmos. La comprensión de estas observaciones requiere teorías analíticas precisas. No obstante, las grandes escalas de tiempo y espacio que implica el estudio de la formación de estructuras suponen un problema para la comparación entre los modelos teóricos y los datos observacionales. En las últimas décadas, gracias al desarrollo de la tecnología computacional, las simulaciones cosmológicas se han convertido en una herramienta indispensable para comprender la evolución del universo en general y la formación de estructura en particular. Estas sirven para mejorar las predicciones analíticas de los modelos y, por tanto, extraer información más precisas de los datos observacionales. 5.1. Simulaciones cosmológicas Para simular la evolución de la materia [30], clave en la formación de estructura, uno de los métodos más utilizados es el de N cuerpos. Este método transforma en discreta la descripción continua del comportamiento de las partículas de materia, dado por la ecuación acolisional de Boltzmann, que describe la evolución en el espacio de fases de la función de densidad de la materia bajo la influencia de un potencial gravitatorio Φ. Para ello, se colocan N partículas de masa midistribuidas en determinadas posiciones rien un volumen prefijado (en general, un cubo con condiciones periódicas de contorno) y con velocidades iniciales vi. Esto puede interpretarse como un muestreo Monte Carlo del espacio de fases. A partir de aquí, el problema se 29 reduce a resolver la ecuación de Poisson para el potencial gravitatorio, que modela la interacción gravitatoria mutua entre las N partículas. Este punto es uno de los grandes desafíos a la hora de realizar simulaciones eficientemente. Para calcular las fuerzas gravitatorias se pueden emplear dos métodos. Por un lado, los métodos basados en resolver la forma integral de la ecuación de Poisson, Φ(r) = −GRdr0ρr0/|r0−r|, utilizan aproximaciones para acelerar la suma de interacciones entre partículas, como por ejemplo, el método de la aproximación de árbol [31]. Por otro lado, los métodos basados en resolver la forma diferencial de la ecuación de Poisson, ∇2Φ=4πGρ(r), utilizan la transformada de Fourier para trabajar en un retículo en el espacio de Fourier. Cabe destacar, que las simulaciones más modernas emplean métodos híbridos con el fin de alcanzar una eficiencia alta. Las simulaciones cosmológicas se realizan en el marco del modelo cosmológico, explicado en el Capítulo 2, a partir de ciertas condiciones iniciales. Para que haya concordancia con las observaciones, se impone un universo de fondo con una curvatura plana [32] y con una determinada composición de materia oscura, materia bariónica, energía oscura y radiación. Además, se deben especificar las perturbaciones primordiales sobre este universo de fondo. Para ello, dado que las predicciones de la teoría de inflación conducen a una distribución gaussiana para las fluctuaciones de densidad, basta establecer el espectro de potencias P(k), tal y como se detalla en el Capítulo 3, donde se muestra que este toma la forma de P(k) = AknT2(k). El último ingrediente para completar las simulaciones cosmológicas son las condiciones iniciales de las posiciones y velocidades de las partículas. Habitualmente, se pueden emplear dos tipos distintos. Por un lado, se estudian simulaciones donde se han impuesto grandes volúmenes periódicos con partículas distribuidas de manera uniformemente aleatoria. Por otro lado, se examinan diversas simulaciones zoom, donde la región de interés muestra una alta resolución frente al resto del volumen simulado. Una vez ejecutada la simulación, el siguiente paso consiste en identificar los halos de materia, objetos virializados de alta densidad. Formalmente, habría que recorrer cada punto comprobando si se cumple el Teorema del Virial, pero esto supondría un alto coste computacional. En lugar de esto, se emplean algoritmos que buscan los picos de densidad. A partir de estos, se definen los centros de los halos. Suponiendo que tienen simetría esférica, se va aumen- 30 tando su radio hasta que la densidad es unas 200 veces la densidad crítica ρc(ecuación 10), correspondiente al valor que predice el modelo de colapso esférico para que se forme un halo (Sección 3.3). Con ello, se asigna a cada halo un radio y una masa, denominados comúnmente como radio virial Rvir y masa virial Mvir del halo. El algoritmo generalmente utilizado se llama ‘Friends-Of-Friends’ (amigos de amigos en inglés), abreviado como FoF. Este define una ‘longitud de amistad’ la partir de la que impone dos reglas. La primera de ellas consiste en que las partículas a una distancia ≤lson partículas amigas y pertenecen al mismo halo. La segunda exige que las partículas también sean amigas de sus amigas. Así, de cada simulación se obtiene catálogo a distintos redshifts que registra las diferentes propiedades de los halos, tales como su masa, su radio, sus coordenadas o si se trata de un halo o un subhalo. 5.2. Simulaciones Legacy Las simulaciones analizadas en este trabajo, denominadas a partir de ahora simulaciones Legacy, han sido ejecutadas con el código Gadget-4 [33]. El modelo cosmológico utilizado en las simulaciones Legacy es consistente con los últimos resultados obtenidos a partir de los datos recogidos por las sondas WMAP9 [13] y Planck [14]. En la Tabla 1, se muestran los valores utilizados para los parámetros cosmológicos fundamentales, presentados en los Capítulos 2 y 3. ΩΛΩmΩbh σ8n 0.7150 0.2850 0.0450 0.6950 0.828 0.9632 Tabla 1: Parámetros cosmológicos usados en las simulaciones La función de transferencia empleada en las simulaciones Legacy fue generada usando CAMB [34], mientras que las condiciones iniciales se han creado mediante el código MUSIC [35]. El conjunto de simulaciones Legacy consta de dos tipos distintos de simulaciones. Por un lado, se han realizado simulaciones de 20483partículas en dos volúmenes diferentes, consis- 31 las incertidumbres de los valores calculados de la función de masa de halo de las simulaciones L1600 y L100 resultan despreciables para los análisis que se van a realizar a continuación. En la siguiente sección, se compararán los resultados de las simulaciones estándar L1600 y L100 con las predicciones analíticas. 5.3.3. Simulaciones estándar Las dos simulaciones estándar a analizar, L1600 y L100 tienen una diferencia clara: el volumen de la primera es muy superior al de la segunda. Esto va a permitir estudiar la función de masa de halo en un amplio rango de masa. En esta sección, se compararán los resultados de ambas simulaciones entre sí y con la predicción analítica de Press-Schechter. Además, también se incluirán las funciones semianalíticas de masa de halos descritas en el Capítulo 4: las funciones de Sheth-Tormen y Jenkins. Calculada la densidad numérica de halos por intervalo de masa en cada simulación, se observa en el panel superior de la Figura 4, que los resultados de las simulaciones L1600 (círculos naranjas) y L100 (cuadrados azules) parecen concordar con las predicciones analíticas. Se puede apreciar que las tres funciones de masa de halo representadas, esto es, las líneas de Press-Schechter (negra continua), Sheth-Tormen (verde discontinua) y Jenkins (roja punteada), guardan cierta similitud entre sí. Para estudiar con detalle las diferencias entre estas funciones y las desviaciones de los datos experimentales respecto a ellas, se calcula el error relativo ε respecto a la función de Press-Schechter, tanto de los resultados de las simulaciones como de las funciones de Sheth-Tormen y Jenkins. Este error relativo se calcula de la siguiente manera. Si Fjes un valor de la función de masa de halo a una masa determinada, bien de una de las simulaciones o bien de las funciones propuestas por Sheth-Tormen o Jenkins, y FP S es el valor que toma la función de Press-Scheter evaluada en la misma masa, el error relativo εviene dado por: ε=Fj−FPS FPS [46] Así, en el panel inferior de la Figura 4, se observa este error relativo respecto al 0, que se corresponde con la función de Press-Scheter. Esto evidencia las diferencias existentes entre 38 los resultados de las simulaciones y las funciones analíticas, cosa difícil de apreciar en el panel superior pues utiliza la escala logarítmica en ambos ejes. Se puede apreciar que los resultados experimentales se ajustan a la función analítica de Press-Schechter con porcentajes de error menores al 20 %hasta masas de 1014 M. Más allá, se observa que los resultados de ambas simulaciones tienden alejarse de la predicción de Press-Schechter y toman valores más altos, en especial los resultados de L1600, cuyo último valor tiene un porcentaje de error del 480 %. Figura 4: Comparación entre la función de masa de halo teóricas y los resultados de las simulaciones. En el panel de arriba, se representa el log10 de la función de masa de halo frente al log10 de la masa. En el panel de bajo se representa la desviación respecto a la función analítica de Press-Schechter frente al log10 de la masa. La línea negra representa la función de Press-Schechter; la línea verde discontinua, la función de Sheth-Tormen y la línea roja punteada, la función de Jenkins. Los puntos naranjas se corresponden con los resultados de la simulación L1600. Los cuadrados azules se corresponden con los resultados de la simulación L100. Una posible explicación para este comportamiento reside en la estadística de las simulaciones. Al estar formados por más partículas, el número de halos de grandes masas (Mvir >1014 M) resultantes de las simulaciones es pequeño, pues están formados por muchas partículas. Esto hace que pequeñas variaciones en el número de estos halos puedan provocar grandes fluctua- 39 ciones y, a grandes masas, se observen desviaciones respecto a las funciones analíticas. Por esta misma razón, la región de grandes masas es la más difícil de modelar. Para ello, hace falta un volumen muy grande. No obstante, diversos estudios en los que también se analizan resultados obtenidos de simulaciones de N cuerpos, indican que la función de masa de halo del formalismo de Press- Schechter subestima el número de halos más masivos [36]. A raíz de esto, su principal objetivo fue obtener una función de masa de halo semianalítica, esto es, utilizando los resultados de simulaciones cosmológicas, tal y como se introdujo en la Sección 4.2. En este trabajo, se han introducido las funciones de masa semianalíticas de Sheth-Tormen (ecuación 43) y de Jenkins (ecuación 44), las cuales también se representan en la Figura 4. Aquí se aprecia como ambas predicen a alta masa una densidad numérica de halos creciente respecto a la función de Press-Schechter, de modo que los resultados de las simulaciones se aproximan con menor error relativo a estas predicciones semianalíticas, aunque siguen dando un número de halos mayor. Con este análisis, se pretende ilustrar que los modelos semianalíticos no son infalibles y solo son válidos hasta una precisión concreta. A pesar de los intentos por obtener una función de masa de halo exacta y universal, a día de hoy todavía existen obstáculos para ello. Por un lado, es necesario dar una definición precisa de la masa de halo [37]. Por otro lado, las simulaciones cosmológicas de N cuerpos trabajan con un volumen finito y un número de partículas restringido debido a los límites de la tecnología computacional. Esto hace que el rango de masas donde las funciones semianalíticas son válidas esté limitado [36]. Es más, tal y como se ha explicado anteriormente, las simulaciones presentan errores. Estos deben ser minimizados en vista a poder igualar la precisión de posibles observaciones [38]. 5.3.4. Simulaciones Zoom La evolución del universo muestra que las zonas densas tienden a ser más densas y las vacías, más vacías. Bajo esta premisa, cabe preguntarse si el número de halos existentes y sus propiedades se ve afectado por la densidad de la zona donde se formen. En esta sección se pretende estudiar la función de masa de halo en regiones con distinta densidad, desde zonas muy poco densas hasta zonas muy densas en comparación con la densidad media. 40 Con este fin, se emplean distintas simulaciones zoom, cuyas principales características se presentan en la Tabla 2. Todas ellas han sido ejecutadas en un volumen de alta resolución formado por una caja de lado 83 Mpc/h inmersa y centrada en una caja mayor de 1600 Mpc/h. Fuera de la caja de 83 Mpc/h hay, gradualmente, menos partículas y de menor resolución, es decir, de mayor masa. Esto significa que, al alejarse de la caja de 83 Mpc/h, hay cada vez menos partículas y estas son más pesadas. Así, cerca de los bordes, las partículas de alta resolución de la caja de 83 Mpc/h interactuarán con otras de menor resolución, lo que puede influir en el número de halos que se formen. Para evitar estos posibles efectos de borde por la diferencia de resolución, se toma como volumen total una caja de 80 Mpc/h centrada en la caja de 83 Mpc/h. El catálogo de halos de cada simulación proporciona la posición de cada halo, dando las coordenadas de su centro (x,y,z). De este modo, basta seleccionar los halos que se encuentren en una caja de 80 Mpc/h centrada en el volumen total. De nuevo, esto se realiza mediante un programa propio en lenguaje Python. A partir de aquí, el procedimiento para calcular el número de halos por unidad de volumen e intervalo de masa de cada una de las simulaciones zoom es similar al realizado para el caso anterior de las simulaciones estándar. Se realiza el filtrado de subhalos, se calcula el número de halos por unidad de volumen e intervalo de masa y se suprimen los datos por debajo del límite de resolución, es decir, por debajo de la masa correspondiente a 70 partículas. En el panel superior de la Figura 5, se representan los resultados obtenidos para las distintas simulaciones zoom junto con la función analítica de Press-Schechter y las funciones semianalíticas de Sheth Tormen y Jenkins. Además, se añaden tanto los resultados de la simulación L100, que tiene un volumen y una resolución de un orden similar a las simulaciones zoom, como los de la simulación L1600, donde se forman halos más masivos. Ambas servirán para comparar los resultados de las simulaciones zoom y estándar. De nuevo, se agrega un segundo panel en el que se muestran los errores relativos εrespecto a la función de Press-Schechter. 41 Figura 5: Comparación entre las funciones de masa de halo analíticas y semianalíticas, la simulacies L100 y L1600, y las simulaciones zoom. En el panel de arriba, se representa el log10 de la función de masa de halo frente al log10 de la masa. En el panel de bajo se representa la desviación respecto a la función analítica de Press-Schechter frente al log10 de la masa. La línea negra representa la función de Press-Schechter; la línea verde discontinua, la función de Sheth-Tormen y la línea roja punteada, la función de Jenkins. Los puntos naranjas se corresponden con los resultados de la simulación L1600. Los cuadrados azules se corresponden con los resultados de la simulación L100. El resto de marcadores se corresponden con los resultados de las simulaciones zoom: S+7σ(círculos verdes), S+2σ(rombos marrones), S+σ(triángulos hacia abajo amarillos), Smedia (triángulos rojos), S-σ(equis grises), S- 2σ(cruces lilas) y S-6σ(estrellas azules). Las simulaciones zoom de zonas de distintas densidades presentan una mayor diferencia a alta masa Si se comparan los resultados de la región de densidad media, Smedia, con la función de Press-Schechter, se observa que, al igual que en las simulaciones estándar, la densidad numérica por intervalo de masa de halos difiere de la función analítica entre 20-80 %en un gran rango de masas (Mvir <1014.5M). De nuevo, en el caso de los halos más masivos, se tiene un porcentaje de error alto, concretamente, del 233 %. Es más, también existe una diferencia de un orden similar respecto a las funciones semianalíticas de Sheth-Tormen y Jenkins. 42 Cabe esperar que esta región de densidad media tenga un comportamiento similar a las simulaciones estándar analizadas anteriormente. Al compararla con los resultados de la simulación L100, en un amplio rango de masas (Mvir <1014 M) difieren menos de un 20 %. Por encima de esta masa, las diferencias se hacen más notables, hasta el punto de que, en el caso de L100, no se forman halos de Mvir ≈1015 M, a diferencia de Smedia. Si, por consiguiente, se comparan los resultados de Smedia con los de la simulación L1600 para masas Mvir ≈1015 Mse observa que Smedia da un número de halos por unidad de volumen e intervalo de masa mayor. Esto puede deberse a las ya mencionadas fluctuaciones estadísticas de varianza cósmica, mayores a grandes masas. Llegado este punto, resulta de gran interés estudiar las fluctuaciones del número de halos por unidad de volumen e intervalo de masa debidas a la densidad de la zona donde se formen. Para ello, se toma la ya analizada simulación Smedia como referencia. De nuevo, la región de alta masa (Mvir >1014 M) presenta las diferencias más importantes. En general, la tendencia observada respecto a Smedia es la siguiente. En las regiones menos densas que la media se forma un número menor de halos por unidad de volumen e intervalo de masa; mientras que, en las zonas más densas que la media este número es mayor. En el caso de S+σ, se puede apreciar que al contrario de lo que se podría intuir, se forman menos halos muy masivos (Mvir ≈1015 M) que en Smedia. Nuevamente, esto puede ser ocasionado por las fluctuaciones estadísticas. Estas fluctuaciones también pueden ser la causa por la que, para las simulaciones S-2σy S-6σ, no se forman halos para masas en torno a Mvir ≈1015 M. En cuanto a las simulaciones de regiones más densas, S+2σy S+7σ, el número de halos más masivos por unidad de volumen e intervalo de masa que se forman tiene una diferencia entre 300-450 %respecto a los resultados de Smedia. En esta línea, resulta interesante estudiar la dependencia de la función de masa de halo con el redshift z. Para ello, se muestra la Figura 6, donde se representa la función de masa de halo de Press-Schechter a distinto zpara ilustrar esta dependencia. 43 Figura 6: Función de masa de halo de Press-Schechter a distinto redshift. Se representa el log10 de la función de masa de halo frente al log10 de la masa. Cada línea se corresponde con: z=0 (azul), z=2 (naranja), z=4 (verde), z=6 (rojo), z=8 (lila), z=10 (marrón). En vista de cómo evoluciona la función de masa de halo respecto al redshift, se puede dar la siguiente explicación a los comportamientos observados en las regiones de diferente densidad. Se ha estudiado que, a masas grandes, los resultados de las simulaciones de regiones más densas (S+7σy S+2σ) dan más halos que media. Esto puede deberse a que estas zonas tienden a hacerse más densas y evolucionan más rápidamente, como si, aparentemente, estuvieran a un redshift menor. En la misma línea, en cuanto a las simulaciones de regiones menos densas, el hecho de que no presenten halos muy masivos encaja con que las regiones menos densas tienden a vaciarse. 44 III Conclusiones 6 Discusión y conclusiones Dentro del paradigma cosmológico del modelo cosmológico en general y del modelo ΛCDM en particular, se ha estudiado el número de halos de materia oscura que se forman por unidad de volumen e intervalo de masa, lo que se conoce como la función de masa de halo. El estudio de la función de masa de halo resulta crucial a la hora de constreñir los parámetros cosmológicos característicos del modelo cosmológico que rige nuestro universo. Para ello, se ha introducido tanto el formalismo analítico de Press-Schechter, construido a partir de las propiedades estadísticas del campo de densidades, como los modelos semianalíticos de Sheth-Tormen y Jenkins, que utilizan los resultados de simulaciones de N cuerpos para ajustar la función de masa de halo. Mediante un código propio en lenguaje Python, se ha hecho un programa capaz de calcular estas tres funciones de masa de halo para distintos intervalos de masa y modelos cosmológicos. Además, se han analizado los resultados de distintas simulaciones cosmológicas de N cuerpos. De nuevo, mediante un código propio en lenguaje Python, se ha podido calcular el número de halos por unidad de volumen e intervalo de masa que se forma en cada una de las simulaciones. En primer lugar, se han analizado dos simulaciones estándar ejecutadas en volúmenes distintos con el fin de estudiar el número de halos por unidad de volumen e intervalo de masa en un rango amplio de masas. La comparación entre los resultados obtenidos y la predicción analítica de Press-Schechter refleja la falta de precisión de dicho modelo. En concreto, esto se hace muy evidente a alta masa, donde los resultados divergen notablemente de la función de Press-Schechter dando un número mayor de halos por unidad de volumen e intervalo de masa. En principio, esto puede deberse a posibles fluctuaciones debidas a la estadística de las 47 [34] A. Lewis and A. Challinor, “CAMB: Code for anisotropies in the microwave background,” Feb. 2011. [35] O. Hahn and T. Abel, “Multi-scale initial conditions for cosmological simulations,” Monthly Notices of the Royal Astronomical Society, vol. 415, p. 2101–2121, Jul 2011. [36] Z. Lukić, K. Heitmann, S. Habib, S. Bashinsky, and P. M. Ricker, “The Halo Mass Function: High-Redshift Evolution and Universality,” The Astrophysical Journal, vol. 671, p. 1160–1181, Dec 2007. [37] M. White, “The mass of a halo,” Astronomy & Astrophysics, vol. 367, p. 27–32, Feb 2001. [38] K. Heitmann, P. M. Ricker, M. S. Warren, and S. Habib, “Robustness of Cosmological Simulations. I. Large-Scale Structure,” The Astrophysical Journal Supplement Series, vol. 160, p. 28–58, Sep 2005. [39] D. Reed, J. Gardner, T. Quinn, J. Stadel, M. Fardal, G. Lake, and F. Governato, “Evolution of the Mass Function of Dark Matter Haloes,” Monthly Notices of the Royal Astronomical Society, vol. 346, p. 565–572, Nov 2003. [40] M. S. Warren, K. Abazajian, D. E. Holz, and L. Teodoro, “Precision Determination of the Mass Function of Dark Matter Halos,” The Astrophysical Journal, vol. 646, p. 881–885, Aug 2006. [41] J. Tinker, A. V. Kravtsov, A. Klypin, K. Abazajian, M. Warren, G. Yepes, S. Gottlöber, and D. E. Holz, “Toward a halo mass function for precision cosmology: The limits of universality,” The Astrophysical Journal, vol. 688, p. 709–728, Dec 2008. [42] M. Crocce, P. Fosalba, F. J. Castander, and E. Gaztañaga, “Simulating the Universe with MICE: The abundance of massive clusters,” Monthly Notices of the Royal Astronomical Society, vol. 403, p. 1353–1367, Apr 2010. [43] S. Bhattacharya, K. Heitmann, M. White, Z. Lukić, C. Wagner, and S. Habib, “Mass Function Predictions Beyond ΛCDM,” The Astrophysical Journal, vol. 732, p. 122, Apr 2011. [44] S. Bocquet, K. Heitmann, S. Habib, E. Lawrence, T. Uram, N. Frontiere, A. Pope, and H. Finkel, “The Mira-Titan Universe. III. Emulation of the Halo Mass Function,” 2020. 54 [45] J. S. Bullock and M. Boylan-Kolchin, “Small-Scale Challenges to the ΛCDM paradigm,” Annual Review of Astronomy and Astrophysics, vol. 55, p. 343–387, Aug 2017. [46] B. P. Moster, R. S. Somerville, C. Maulbetsch, F. C. van den Bosch, A. V. Macciò, T. Naab, and L. Oser, “Constraints on the relationship between stellar mass and halo mass at low and high redshift,” The Astrophysical Journal, vol. 710, p. 903–923, Jan 2010. [47] P. S. Behroozi, R. H. Wechsler, and C. Conroy, “The Average Star Formation Histories of Galaxies in Dark Matter Halos from z= 0-8,” The Astrophysical Journal, vol. 770, p. 57, May 2013. 55 Agradecimientos Primero de todo, quiero agradecer a mi tutor, José Oñorbe Bernis, por la dedicación y la constante ayuda que me ha brindado para la elaboración de este Trabajo de Fin de Grado. Durante meses me ha guiado por el emocionante mundo de la astrofísica y la cosmología. También, me gustaría agradecer a los que siempre han estado conmigo. A mis padres, por el incansable apoyo en tiempos difíciles. A mi hermana, por escucharme, por entenderme, por todo. Por último, quiero agradecer a Puccini, Debussy y Ravel por ser mi particular banda sonora en este viaje de la armonía de las esferas al universo en expansión. 57