Modelo de expansión del universo en Python

Introducción

La idea de que el universo se expande nació casi a la fuerza de las ecuaciones. En 1922, Alexander Friedmann derivó la primera solución cosmológica no estacionaria de la relatividad general aplicando el principio cosmológico (homogeneidad e isotropía a gran escala) a las ecuaciones de Einstein (Friedmann, 1922). Cinco años después, Lemaître (1927) llegó a la misma conclusión por una vía independiente y, sobre todo, la conectó con las observaciones: si el universo se expande, las galaxias deben alejarse de nosotros con velocidades proporcionales a su distancia. Hubble (1929) lo confirmó experimentalmente, y Penzias y Wilson (1965) cerraron el círculo al detectar el fondo cósmico de microondas, fósil térmico de aquella expansión primitiva.

Hoy, un siglo después, los parámetros que rigen esa expansión se conocen con una precisión inimaginable para aquellos pioneros: la edad del universo se mide con un error inferior al 1 % (Planck Collaboration, 2020). Esa precisión permite algo que entonces no era posible: construir un modelo simplificado del universo y ver hasta dónde puede llevarnos antes de quebrarse.

Este trabajo aborda esa construcción en tres pasos. Primero, se deriva la ecuación de aceleración de Friedmann a partir del principio cosmológico, separando el contenido del universo en sus tres componentes (materia, radiación y energía oscura) y obteniendo la dependencia explícita de cada densidad con el factor de escala. Segundo, se resuelven las ecuaciones en Python con las constantes y sus errores 1σ tomados de Planck Collaboration (2020). Tercero, se integra el modelo hacia el pasado y hacia el futuro para distintos valores del parámetro ω de la ecuación de estado de la energía oscura, lo que permite leer en una sola figura el papel de cada componente y, en particular, el destino del universo en función de ω. Se verá cómo ω, un número adimensional, conduce a escenarios en los que el universo se expande, se acelera o termina en un Big Rip.

1. Marco teórico

1.1. El principio cosmológico y la métrica FLRW

Toda la cosmología estándar descansa sobre una hipótesis aparentemente inocente: a escalas suficientemente grandes (del orden de cientos de Mpc), el universo es homogéneo (sus propiedades no dependen del lugar) e isótropo (no dependen de la dirección). Esta hipótesis no es gratuita: los mapas de galaxias del SDSS y el espectro casi perfecto del fondo cósmico de microondas la confirman observacionalmente con una precisión de una parte en 10⁵ (Planck Collaboration, 2020). La única métrica espaciotemporal compatible con esta simetría es la de Friedmann–Lemaître–Robertson–Walker (FLRW), cuyo elemento de línea adopta la forma:

ds² = −c² dt² + a²(t) · [dr² / (1 − k r²) + r² (dθ² + sen²θ dφ²)]

Aquí t es el tiempo cósmico que mide un observador comóvil con el flujo de la expansión; (r, θ, φ) son coordenadas espaciales comóviles, que no cambian para objetos que solo se mueven con la expansión; y a(t) es el factor de escala, la magnitud central de toda la cosmología: cuantifica cuánto se han estirado las distancias respecto a un momento de referencia. Por convención se elige a(t₀) = 1 en el presente, de modo que a < 1 corresponde al pasado y a > 1 al futuro. El parámetro k ∈ {−1, 0, +1} describe la curvatura espacial: hiperbólica, plana o esférica, respectivamente. Las medidas de Planck son compatibles con k = 0 dentro de errores muy pequeños (|Ω_k| < 0,005), así que en adelante se adopta un universo espacialmente plano.

1.2. Derivación de la ecuación de aceleración

La relatividad general relaciona la geometría del espacio-tiempo con su contenido material a través de las ecuaciones de campo de Einstein:

G_μν = (8πG / c⁴) · T_μν

Sustituyendo la métrica FLRW en el tensor de Einstein G_μν y suponiendo que el contenido del universo se comporta como un fluido perfecto (caracterizado únicamente por una densidad de energía ρ y una presión isótropa p), las diez ecuaciones de campo se reducen a solo dos ecuaciones independientes. La componente G₀₀ da la primera ecuación de Friedmann:

H²(t) ≡ (ȧ / a)² = (8πG / 3) · ρ − k c² / a²

Y la combinación con las componentes espaciales (equivalente a usar la conservación covariante ∇_μT^μν = 0) da la segunda ecuación de Friedmann o ecuación de aceleración:

ä / a = −(4πG / 3) · (ρ + 3p / c²)

Esta segunda ecuación es la que se utiliza en el cálculo en Python y merece una lectura cuidadosa. El miembro izquierdo es la aceleración relativa del factor de escala: si es positivo, el universo se expande cada vez más rápido; si es negativo, se desacelera. El miembro derecho tiene siempre el signo opuesto a (ρ + 3p/c²). Por tanto, la única forma de que el universo se acelere es que su contenido cumpla (ρ + 3p/c²) < 0, lo que exige una presión negativa muy grande. Ningún fluido conocido (gas, polvo, fotones) cumple esa condición: la aceleración observada implica obligatoriamente una componente exótica, la energía oscura.

Multiplicando la ecuación de aceleración por a(t) se obtiene la forma:

d²a / dt² = −(4πG / 3) · [ρ_M(a) + 2 ρ_R(a) + (1 + 3ω) ρ_E(a)] · a(t)

donde se han separado las tres componentes según la ecuación de estado p = ωρc². El paso siguiente es averiguar cómo escala cada densidad con a(t).

1.3. Las tres componentes del universo y su evolución con a(t)

La conservación covariante del tensor energía-momento, aplicada a la métrica FLRW, da una sola ecuación escalar para cada componente, la ecuación de continuidad cosmológica:

dρ / dt + 3H (1 + ω) ρ = 0

Su lectura es directa: la densidad cambia por dos motivos, la dilución por el aumento de volumen (factor 3H) y el trabajo que la presión realiza sobre ese volumen al expandirse (factor adicional 3Hω). Separando variables e integrando se obtiene la dependencia genérica:

ρ(a) = ρ₀ · a^−3(1+ω)

Con esta fórmula y el valor adecuado de ω para cada componente se derivan las tres relaciones que rigen la evolución del universo.

I. Materia (ω_M = 0). Engloba todo lo que se comporta de forma no relativista: estrellas, gas, polvo, planetas y, sobre todo, materia oscura, que aporta más del 80 % de esta categoría. Su presión a escalas cosmológicas es completamente despreciable (la velocidad media de las partículas es mucho menor que c), por lo que ω_M = 0:

ρ_M(a) = ρ_M,0 / a³

La interpretación es la conservación del número de partículas: en un volumen comóvil siempre hay el mismo número de galaxias, así que la densidad cae como el inverso del volumen físico, que crece como a³.

II. Radiación (ω_R = 1/3). Incluye los fotones del fondo cósmico de microondas y los neutrinos mientras son relativistas. Su ecuación de estado es p_R = ρ_Rc²/3 (un resultado clásico de la relatividad especial para un gas de partículas ultrarrelativistas), de modo que ω_R = 1/3:

ρ_R(a) = ρ_R,0 / a⁴

Aparece una potencia adicional respecto a la materia porque, además de la dilución volumétrica (a⁻³), cada fotón pierde energía al desplazarse al rojo en un factor a⁻¹. El factor 2 que multiplica a ρ_R en la ecuación de aceleración se origina precisamente al sustituir 3p_R/c² = ρ_R: (ρ_R + 3p_R/c²) = 2ρ_R.

III. Energía oscura (ω_E = ω). Es la componente conceptualmente más delicada. No se sabe qué es, pero sí cómo se comporta a primer orden: su presión es muy negativa, p_E = ω ρ_E c², con ω ≈ −1. La interpretación canónica es la constante cosmológica Λ de Einstein, entendida como la densidad de energía del vacío cuántico, que daría exactamente ω = −1. Las restricciones observacionales más estrictas, que combinan los datos del CMB de Planck con las supernovas de tipo Ia (Planck Collaboration, 2020), limitan ω = −1,03 ± 0,03, compatible con Λ dentro del 3 %. Sustituyendo en la fórmula general:

ρ_E(a) = ρ_E,0 / a^3(1+ω) = ρ_E,0 / a^3+3ω

Para ω = −1 el exponente es cero y la densidad es estrictamente constante: cuando el universo se estira, se genera más vacío con la misma densidad (de ahí su nombre alternativo de energía del vacío). Para ω < −1 (la llamada “energía fantasma”), la densidad crece con la expansión, lo que conduce al desbocamiento conocido como Big Rip (Caldwell, Kamionkowski y Weinberg, 2003).

IV. La frontera ω = −1/3. Sustituyendo p_i = ω_i ρ_i c² en la ecuación de aceleración y dividiendo entre la densidad crítica ρ_c = 3H²/(8πG), se obtiene la condición de aceleración en forma adimensional:

Σ_i (1 + 3ω_i) Ω_i(a) < 0

• Para la materia, 1 + 3 · 0 = 1 (siempre positivo).
• Para la radiación, 1 + 3 · (1/3) = 2 (siempre positivo).
• Solo una componente con ω < −1/3 puede invertir el signo de la suma y producir aceleración.

El valor ω = −1/3 marca, por tanto, la frontera teórica entre dos escenarios cualitativamente distintos: por encima, el universo se desacelera como una pelota lanzada hacia arriba; por debajo, se acelera como si “algo invisible” empujara cada galaxia hacia fuera. Esta frontera reaparece explícitamente al comparar varios valores de ω en los resultados.

2. Metodología del cálculo

2.1. Integración numérica de la ecuación de Friedmann

La ecuación d²a/dt² = f(a, ω) es una ecuación diferencial ordinaria de segundo orden. Para integrarla se reescribe como un sistema de primer orden introduciendo la velocidad v ≡ da/dt:

v̇ = f(a, ω); ȧ = v

y se aplica un esquema de Euler explícito:

v_n+1 = v_n + Δt · f(a_n, ω) a_n+1 = a_n + Δt · v_n+1

Para la integración hacia el pasado se cambia el signo del paso. El paso temporal Δt = 10⁷ años es 1/1450 del tiempo de Hubble, lo bastante pequeño para que el error global de Euler, de primer orden O(Δt), quede por debajo del 0,1 % en el rango de integración. Como verificación, repetir el cálculo con Δt = 5 × 10⁶ años no altera la edad obtenida en más de 0,01 Gyr.

La integración parte de la condición actual (t = 0, a = 1) y termina, hacia el pasado, cuando el factor de escala alcanza un valor de referencia a_min = 0,01 (cuando el universo era una centésima del actual y la dinámica estaba dominada por la radiación); hacia el futuro, se prolonga hasta t = +10¹⁰ años en los análisis estándar y hasta t = +5 × 10¹⁰ años para estudiar el comportamiento asintótico.

2.2. Parámetros cosmológicos adoptados y conversión de unidades

Las constantes se expresan en un sistema mixto:

• La constante de Hubble, H₀, en años⁻¹.
• La constante de gravitación universal, G, en m³ kg⁻¹ año⁻².
• Las densidades, ρ_i, en kg m⁻³.

La elección de unidades no afecta al resultado físico, pero exige convertir los valores publicados. El valor de Planck H₀ = (67,4 ± 0,5) km/s/Mpc se transforma mediante:

H₀ [año⁻¹] = (H₀ [km/s/Mpc] / 1 Mpc [km]) · 1 año [s]

con 1 Mpc = 3,0857 × 10¹⁹ km y 1 año juliano = 3,15576 × 10⁷ s. El resultado es H₀ = 6,8930 × 10⁻¹¹ año⁻¹. Análogamente, la constante de gravitación universal en su valor CODATA 2018 (Tiesinga et al., 2021), 6,6743 × 10⁻¹¹ m³ kg⁻¹ s⁻², se reescala multiplicando por año², lo que da G = 6,6468 × 10⁴ m³ kg⁻¹ año⁻². Las densidades de cada componente se calculan a partir de los parámetros de densidad adimensionales Ω_i, definidos como ρ_i,0 = Ω_i · ρ_c, donde la densidad crítica es ρ_c = 3H₀²/(8πG) = 8,53 × 10⁻²⁷ kg/m³. Los valores Ω_M = 0,3153 ± 0,0073 y Ω_E = 0,6847 ± 0,0073 proceden directamente de Planck Collaboration (2020), mientras que Ω_R ≈ 9,24 × 10⁻⁵ se calcula a partir de la temperatura del fondo cósmico medida por COBE/FIRAS, T₀ = 2,725 K (Fixsen, 2009), añadiendo la contribución de los tres tipos de neutrinos relativistas. La Tabla 1 resume todos los valores con su error 1σ.

Parámetro	Valor adoptado	Unidades	Origen y error 1σ
Constante de Hubble H₀	6,8930 × 10⁻¹¹	año⁻¹	67,4 ± 0,5 km/s/Mpc (Planck)
Constante de gravitación G	6,6468 × 10⁴	m³ kg⁻¹ año⁻²	CODATA 2018 (error 10⁻¹⁰)
Densidad de materia ρM,0	(2,69 ± 0,07) × 10⁻²⁷	kg m⁻³	ΩM = 0,3153 ± 0,0073
Densidad de radiación ρR,0	(7,88 ± 0,12) × 10⁻³¹	kg m⁻³	ΩR ≈ 9,24 × 10⁻⁵ (Fixsen, 2009)
Densidad de energía oscura ρE,0	(5,84 ± 0,11) × 10⁻²⁷	kg m⁻³	ΩE = 0,6847 ± 0,0073
Parámetro ω (ecuación de estado)	−1,03 ± 0,03	—	Planck + SNe Ia

Tabla 1. Constantes y densidades adoptadas, con las unidades empleadas y las incertidumbres 1σ propagadas desde Planck Collaboration (2020) y CODATA 2018 (Tiesinga et al., 2021).

3. Resultados

3.1. Curvas a(t) para distintos valores de ω

La forma natural de hacer hablar al modelo es integrarlo con varios valores de ω y representar todas las soluciones a(t) en la misma gráfica. Para cubrir el espectro físico relevante se eligen cinco valores.

Antes de presentarlos conviene introducir la categoría a la que pertenecen dos de ellos. El término quintessence procede de la “quinta esencia” aristotélica (el quinto elemento, el éter celeste que se añadía a la tierra, el agua, el aire y el fuego). En cosmología moderna designa a la quinta componente del universo, después de los bariones, la materia oscura, los fotones y los neutrinos: una alternativa dinámica a la constante cosmológica como motor de la expansión acelerada. Mientras que Λ es un valor estrictamente constante (ω = −1 exacto), la quintessence se modela como un campo escalar φ(t) con energía potencial V(φ) que evoluciona en el tiempo. Su ecuación de estado cumple −1 < ω < −1/3 (siempre por encima de −1, a diferencia de la energía fantasma) y puede en principio variar con el tiempo, ω(t). Aquí se trata como constante a primer orden, pero variando su valor entre simulaciones se exploran modelos teóricos que predicen distintos ω efectivos.

• Energía fantasma marginal: ω = −1,10. Fuera del intervalo observacional de Planck (ω = −1,03 ± 0,03), pero útil para visualizar el régimen ω < −1, en el que la densidad de energía oscura crece con la expansión: ρ_E(a) = ρ_E,0 · a^+0,3.
• Constante cosmológica Λ canónica: ω = −1,00. Caso ΛCDM puro, dentro del intervalo de Planck. La densidad de energía oscura es estrictamente constante: ρ_E(a) = ρ_E,0.
• Quintessence próxima a −1: ω = −0,90. Modelo dinámico cercano a la constante cosmológica, aunque fuera del intervalo de Planck (a unos 4σ), cuya densidad se diluye lentamente con la expansión: ρ_E(a) = ρ_E,0 · a^−0,3.
• Quintessence fuerte: ω = −0,50. Claramente fuera de los límites de Planck; se incluye para visualizar el contraste. La dilución de la energía oscura es más rápida que en los casos anteriores: ρ_E(a) = ρ_E,0 · a^−1,5.
• Frontera teórica: ω = −1/3. Valor crítico entre aceleración y desaceleración. La densidad se diluye como la curvatura espacial efectiva: ρ_E(a) = ρ_E,0 · a⁻².

La primera lectura de la Figura 1 es que, para los tres valores más próximos a −1, las curvas son prácticamente indistinguibles en todo el rango pasado y solo empiezan a separarse a partir de t ≈ +3 Gyr. Esto refleja la elevadísima precisión con la que Planck restringe ω = −1,03 ± 0,03: cerca de ese valor, la historia pasada del universo es indistinguible de la del modelo ΛCDM puro. Por el contrario, las curvas con ω = −0,5 y ω = −1/3 se separan ya en el pasado (alcanzan a → 0 unos 2 Gyr más tarde, lo que da edades menores) y, sobre todo, en el futuro: con ω = −1/3 la expansión es asintóticamente lineal y el universo no se acelera nunca.

3.2. Futuro lejano y destino del universo

Para que el papel de ω se vea con claridad conviene prolongar la integración y representar a(t) en escala logarítmica (Figura 2). En esa escala, una expansión exponencial aparece como una recta, lo que facilita la comparación.

La Figura 2 transmite tres ideas:

1. Para ω = −1 (ΛCDM puro), el factor de escala crece exponencialmente, a(t) ∝ exp(H_∞ · t), con H_∞ = H₀ √Ω_E ≈ 5,7 × 10⁻¹¹ año⁻¹. Esa expansión lleva, en un tiempo del orden de varios tiempos de Hubble, a un universo de De Sitter en el que la mayor parte de las galaxias salen del horizonte y la temperatura del CMB tiende a un valor finito muy pequeño (Gibbons y Hawking, 1977).
2. Para ω = −0,9 (próximo a −1, aunque fuera del intervalo de Planck), la expansión sigue acelerándose, pero ya no es exponencial: la densidad de energía oscura se diluye lentamente y, a largo plazo, a(t) ∝ t^2/(3(1+ω)) = t^6,7, una ley de potencias muy rápida que en escala logarítmica se curva poco a poco por debajo de la recta de De Sitter.
3. Para ω < −1 (la curva ω = −1,10), la densidad de energía oscura crece con a, lo que retroalimenta la aceleración y produce un desbocamiento: el factor de escala diverge en un tiempo finito, el fenómeno conocido como Big Rip (Caldwell, Kamionkowski y Weinberg, 2003). Aplicando a los parámetros de Planck Collaboration (2020) la expresión analítica de dichos autores se obtiene un horizonte hasta la singularidad de ≈ 117 Gyr, en concordancia con el límite numérico en el que la integración diverge (t ≈ 116 Gyr; véase la Figura 3). Antes de esa singularidad, las estructuras gravitacionalmente ligadas se desligan progresivamente en orden inverso a su densidad media.

Los tres regímenes se resumen cuantitativamente en la Tabla 2, junto con la aceleración actual y el parámetro de desaceleración q₀.

ω adoptado	Régimen	q₀	ä hoy (año⁻²)	Destino
−1,10	Energía fantasma	−0,630 ± 0,012	+(2,99 ± 0,07) × 10⁻²¹	Big Rip
−1,00	ΛCDM (Λ pura)	−0,527 ± 0,011	+(2,50 ± 0,06) × 10⁻²¹	De Sitter asintótico
−0,90	Quintessence	−0,424 ± 0,010	+(2,02 ± 0,06) × 10⁻²¹	Expansión potencial
−0,50	Quintessence fuerte	−0,013 ± 0,006	+(6,4 ± 2,6) × 10⁻²³	Casi en la frontera
−1/3	Curvatura efectiva	+0,158 ± 0,004	−(7,5 ± 0,2) × 10⁻²²	Desaceleración perpetua

Tabla 2. Resumen de los cinco escenarios analizados. El parámetro de desaceleración se calculó como q₀ = ½Ω_M + Ω_R + ½(1 + 3ω)Ω_E; la aceleración actual ä es la salida directa de la ecuación de Friedmann en a = 1.

Un caso interesante: ω = −0,50. En la tabla se observa que, para ω = −0,50, el error relativo de ä es del 40 %, enorme comparado con el ~2–3 % de los demás casos. Se debe a que, cerca de ω = −1/3, la aceleración pasa por cero (es la frontera entre acelerar y desacelerar). En ese régimen, ä es la diferencia entre dos cantidades comparables, y los errores de cada una se vuelven comparables a la propia diferencia. No es un defecto del cálculo, es física: se está midiendo algo que es casi cero, así que la precisión relativa se degrada.

Como exploración complementaria del comportamiento asintótico, se prolongó la integración hasta t = +1,1598 × 10¹¹ años (≈ 116 Gyr), el límite a partir del cual el cálculo numérico deja de ser estable. El resultado (Figura 3) revela con claridad la asimetría dramática que el parámetro ω impone sobre el destino del universo.

Para ω = −0,5 y ω = −1/3, el factor de escala apenas alcanza a ≈ 4 tras más de cien mil millones de años de expansión, lo que evidencia un régimen subexponencial cada vez más lento. Para ω = −1, a(t) crece exponencialmente, en consonancia con la solución de De Sitter asintótica a(t) ∝ exp(H_∞ t); para ω = −0,9 el crecimiento es casi tan rápido (ley de potencias a ∝ t^6,7), y ambas superan un factor 10² al final del intervalo.

El caso ω = −1,1, en cambio, muestra la firma inequívoca de una singularidad finita: el factor de escala se dispara hasta a ~ 10⁹¹ en los últimos diez gigaaños del rango simulado y diverge de forma superexponencial. La interrupción del integrador con el error OverflowError: Numerical result out of range en t ≈ 1,16 × 10¹¹ años no es un artefacto del esquema de Euler, sino la traducción computacional de la singularidad real del modelo: el Big Rip predicho analíticamente por Caldwell, Kamionkowski y Weinberg (2003), cuyo tiempo característico para ω = −1,1 se estima en ≈ 117 Gyr.

Que la divergencia numérica se adelante ligeramente a la singularidad analítica es esperable: cuando las cantidades involucradas (a y su derivada ȧ) crecen tan rápido que algún producto intermedio supera el límite de los números de coma flotante de doble precisión (~10³⁰⁸), el integrador se detiene. Esto refuerza la lectura física: pequeñas desviaciones de ω por debajo de −1, aunque indistinguibles en las escalas observacionales accesibles hoy, condenarían al universo a una desintegración estructural en un horizonte cosmológico finito.

3.3. Evolución de las densidades: las dos transiciones cósmicas

Las simulaciones anteriores muestran cómo evoluciona el factor de escala a(t) bajo distintos escenarios de energía oscura, pero la dinámica subyacente (qué componente domina en cada época) queda oculta tras la suma global. Para hacer explícita esta composición, la Figura 4 representa la densidad de cada componente en función del factor de escala.

Las pendientes en escala log-log son las potencias derivadas en la sección 1.3: −3 para la materia, −4 para la radiación, 0 para la energía oscura con ω = −1 y −3(1 + ω) = −1,5 para ω = −0,5, incluida como contraste ilustrativo. Esta figura no procede de la integración numérica de las simulaciones, sino de la evaluación directa de las expresiones analíticas ρ_i(a) = ρ_i,0 · a^−3(1+ω_i) sobre una rejilla logarítmica del factor de escala. Esto permite identificar visualmente dos transiciones cósmicas con significado físico crítico:

• Igualdad materia–radiación: z_eq,MR = 3402 ± 26 (Planck Collaboration, 2020, Tabla 2), equivalente a un factor de escala a_eq,MR ≈ 2,94 × 10⁻⁴ y a un tiempo t ≈ 51.000 años después del Big Bang. Antes de ese punto, la radiación dominaba la dinámica (línea roja por encima de la naranja en la Figura 4) y la presión de los fotones impedía el colapso gravitacional de las inhomogeneidades primordiales. Después, la materia toma el control y se inaugura el régimen en el que las perturbaciones primordiales pueden crecer y dar lugar a las estructuras a gran escala que se observan hoy (Dodelson, 2003).
• Igualdad materia–energía oscura: a partir de Ω_M = 0,3153 ± 0,0073 y Ω_Λ = 0,6847 ± 0,0073 (Planck Collaboration, 2020), imponiendo ρ_M = ρ_Λ se obtiene a_eq,MΛ = (Ω_M/Ω_Λ)^1/3 = 0,772 ± 0,007, equivalente a un redshift z_eq,MΛ = 0,295 ± 0,012 y a un tiempo en el pasado de aproximadamente 3,5 Gyr. Para poner la escala temporal en contexto: en ese momento la Tierra ya estaba formada y la vida microbiana llevaba miles de millones de años evolucionando. Antes de ese cruce dominaba la materia (línea naranja por encima de la azul en la Figura 4) y el universo se desaceleraba bajo su propio peso gravitacional; después, la energía oscura toma el control y el universo entra en régimen de aceleración. Esta es, por tanto, la transición que explica físicamente por qué hoy se observa un universo en expansión acelerada (Riess et al., 1998; Perlmutter et al., 1999).

El problema de la coincidencia. El cruce materia–Λ ocurrió hace apenas el equivalente al 25,4 % de la edad del universo, lo que significa que vivimos en una época cosmológicamente excepcional: el momento mismo en que el universo cambia de naturaleza dinámica. Es el llamado problema de la coincidencia cósmica (Steinhardt, 1997; Zlatev, Wang y Steinhardt, 1999): ¿por qué toca observar el universo justo cuando sus dos componentes dominantes, materia y energía oscura, tienen densidades comparables, si esa coincidencia solo se da durante una fracción muy reducida de la historia cósmica? No hay consenso sobre la respuesta. Una interpretación es el principio antrópico: solo en esta época pueden existir observadores como nosotros, ya que antes las estructuras aún no se habían formado y más tarde la aceleración habrá disgregado las galaxias y la materia se habrá diluido en exceso. Otra, más sustantiva físicamente, es que la energía oscura no sea una constante cosmológica Λ inerte, sino un campo dinámico de tipo quintessence cuyo potencial V(φ) converja de forma natural hacia el régimen actual mediante una dinámica de tracking, lo que evitaría el ajuste fino que requiere Λ (Zlatev, Wang y Steinhardt, 1999). Las observaciones de las próximas dos décadas, especialmente las de Euclid y Roman, deberán arbitrar entre ambas hipótesis.

4. Discusión e interpretación física

Los resultados confirman, dentro del marco simplificado de este análisis, los tres rasgos esenciales del modelo cosmológico estándar.

I. Hubo un Big Bang. Para los cinco valores de ω considerados, la integración hacia el pasado alcanza a → 0 en un tiempo finito. La edad obtenida varía moderadamente con ω, entre aproximadamente 11,6 Gyr (ω = −1/3) y 14,0 Gyr (ω = −1,1), una dispersión cercana al 10 % en torno al valor central. Esta variación, aunque visible en la Figura 1, es modesta porque la dinámica del universo temprano está dominada por la radiación y la materia, no por la energía oscura, que solo controla la expansión durante los últimos ≈ 6 Gyr de la historia cósmica. Para el caso de referencia ω = −1 se obtiene 13,78 ± 0,11 Gyr, en concordancia, dentro de menos de 1σ, con el valor publicado por Planck Collaboration (2020) de 13,797 ± 0,023 Gyr.

II. La expansión actual es acelerada. En el modelo ΛCDM, la transición de un régimen de desaceleración a uno de aceleración (q = 0) se produce en a ≈ 0,61, hace aproximadamente 6 Gyr, cuando la densidad de materia cae por debajo del doble de la de energía oscura (ρ_M = 2ρ_Λ). Algo después, en a ≈ 0,77 (hace unos 3,5 Gyr), se alcanza la igualdad estricta de densidades, ρ_M = ρ_Λ, y la energía oscura pasa a ser la componente dominante de forma definitiva. Este resultado no es ajustable libremente en el modelo: lo impone el cociente Ω_M/Ω_Λ medido por Planck Collaboration (2020) y constituye una de las predicciones más sólidas del modelo ΛCDM. En el caso límite ω = −1/3, la aceleración nunca llega a producirse (q₀ = +0,158 según la Tabla 2): el universo se asentaría en una expansión asintóticamente lineal, en contradicción frontal con las observaciones de supernovas de tipo Ia que descubrieron la aceleración (Riess et al., 1998; Perlmutter et al., 1999) y que el modelo ΛCDM reproduce con naturalidad.

III. El destino futuro depende críticamente de ω. Para ω = −1 (constante cosmológica pura), el universo se expande exponencialmente para siempre, las galaxias lejanas salen del horizonte y la temperatura del CMB tiende a un valor finito y muy pequeño (Gibbons y Hawking, 1977). Para −1 < ω < −1/3 (quintessence), el universo también se expande para siempre y de forma acelerada, pero según una ley de potencias en lugar de una exponencial. Para ω < −1 (energía fantasma), la aceleración se desboca y el universo termina en un Big Rip en un tiempo finito; para ω = −1,1, este horizonte es de ≈ 117 Gyr (sección 3.2), de acuerdo con la solución analítica de Caldwell, Kamionkowski y Weinberg (2003). Medir ω con la precisión necesaria para distinguir entre estos destinos es el objetivo central de las misiones Euclid (Euclid Collaboration et al., 2025) y Roman (Spergel et al., 2015), y los resultados de DESI apuntan en esa dirección: tras las primeras tensiones de 2024 (DESI Collaboration, 2024), la segunda publicación de datos (DR2, 2025) prefiere una energía oscura dinámica (w₀ > −1, w_a < 0) con una significancia de 2,8–4,2σ según la muestra de supernovas que se combine (DESI Collaboration, 2025).

Conclusiones

Se ha derivado e integrado numéricamente la ecuación de aceleración de Friedmann para un universo plano dominado por materia, radiación y energía oscura. Con los parámetros de Planck Collaboration (2020) en sus unidades naturales y propagando sus errores 1σ, el modelo reproduce de forma satisfactoria los observables fundamentales del modelo ΛCDM: la edad del universo (13,78 ± 0,11 Gyr), el tiempo de Hubble (14,51 ± 0,11 Gyr) y el parámetro de desaceleración actual (q₀ = −0,527 ± 0,011), todos compatibles con los valores publicados dentro de menos de 1σ.

La exploración sistemática del parámetro ω confirma que ω = −1/3 actúa como frontera teórica entre dos universos cualitativamente distintos (q₀ > 0 frente a q₀ < 0) y que el destino último queda controlado por el signo de (ω + 1): expansión eterna para ω ≥ −1 (exponencial para ω = −1 y potencial para −1 < ω < −1/3) y Big Rip para ω < −1. La concordancia entre el horizonte numérico de inestabilidad (t ≈ 117 Gyr para ω = −1,1) y el tiempo de la singularidad predicho analíticamente por Caldwell, Kamionkowski y Weinberg (2003) constituye una verificación cruzada del esquema numérico.

La identificación gráfica de las dos transiciones cósmicas (igualdad materia–radiación en a ≈ 2,94 × 10⁻⁴ e igualdad materia–energía oscura en a ≈ 0,77) permite contextualizar la época actual: la última de estas transiciones ocurrió hace apenas 3,5 Gyr, el 25 % de la edad total del universo. El régimen de aceleración propiamente dicho (q < 0) comenzó algo antes, en a ≈ 0,61 (hace ≈ 6 Gyr). Esta circunstancia de vivir tan cerca del momento en que el universo cambia de naturaleza dinámica es el problema de la coincidencia cósmica (Zlatev, Wang y Steinhardt, 1999), uno de los argumentos observacionales que motivan la exploración de alternativas a la constante cosmológica, como la quintessence dinámica, cuya resolución corresponderá a las próximas generaciones de observatorios cosmológicos.

Anexo. Herramientas utilizadas

Este anexo describe las herramientas computacionales empleadas para producir los resultados numéricos y las figuras. Todo el desarrollo está contenido en un cuaderno Jupyter.

A.1. Entorno y dependencias. El análisis se realizó íntegramente en Python 3 dentro de un cuaderno Jupyter. Las únicas dependencias son tres bibliotecas del ecosistema científico estándar: math (módulo estándar), para constantes matemáticas como π; NumPy (Harris et al., 2020), para el manejo de arrays y operaciones vectorizadas; y Matplotlib (Hunter, 2007), para todas las figuras. No se utilizaron bibliotecas cosmológicas especializadas, de forma deliberada, por el objetivo del trabajo: construir el modelo desde primeros principios.

A.2. Esquema de integración numérica. La ecuación de aceleración de Friedmann d²a/dt² = f(a, ω) se integra con un esquema de Euler explícito de primer orden, reescribiendo la ecuación diferencial ordinaria de segundo orden como un sistema de dos ecuaciones de primer orden con la variable auxiliar v ≡ da/dt. En cada iteración del bucle while se actualiza primero la velocidad con la aceleración instantánea y, a continuación, el factor de escala con la velocidad recién calculada. El paso temporal Δt = 10⁷ años es 1/1450 del tiempo de Hubble, lo bastante pequeño para que el error global de Euler quede por debajo del 0,1 % en todo el rango. La integración hacia el pasado se obtiene cambiando el signo del paso. Las condiciones de parada son a = 0,01 hacia el pasado y t = 10¹⁰ años hacia el futuro en las simulaciones estándar; t = 5 × 10¹⁰ años y t = 1,1598 × 10¹¹ años en los análisis de futuro extendido.

A.3. Organización del cuaderno y flujo de trabajo. El cuaderno se organiza en tres bloques funcionales:

1. Bloque 1: constantes y bucles de integración originales. Celdas en las que se asignan los valores numéricos de H₀, G, ρ_M,0, ρ_R,0 y ρ_E,0, y los bucles while originales del cuaderno base (uno para el pasado y otro para el futuro), que producen la simulación individual para un único valor de ω y su visualización básica.
2. Bloque 2: función reutilizable y barrido sistemático sobre ω. Para no repetir el código cinco veces, los bucles del bloque 1 se encapsulan en la función simular(omega, t_futuro_max, a_pasado_min), que devuelve un array (N, 2) con el tiempo en Gyr y el factor de escala. A partir de ella se itera sobre la lista omegas = [−1.10, −1.00, −0.90, −0.50, −1/3] para producir las cinco curvas de la Figura 1 y, ampliando t_futuro_max, la Figura 2.
3. Bloque 3: verificaciones y figuras complementarias. Incluye el cálculo de magnitudes derivadas (edad del universo, q₀ y aceleración actual para cada ω); la verificación de convergencia del esquema de Euler, calculando la edad con cinco pasos temporales entre 10⁶ y 2 × 10⁷ años y comprobando que las cifras concuerdan dentro de 0,01 Gyr; la exploración del Big Rip mediante una integración prolongada hasta t = 1,1598 × 10¹¹ años para ω = −1,1; y la evaluación de las densidades analíticas ρ_i(a) sobre una rejilla logarítmica que genera la Figura 4.

A.4. Reproducibilidad. Todo el análisis es reproducible ejecutando secuencialmente las celdas del cuaderno. Los resultados intermedios de cada simulación se guardan automáticamente en ficheros de texto con el patrón expansionUniverso_<omega>.txt. No se requiere semilla aleatoria, porque la integración es completamente determinista, y el cálculo completo tarda menos de un minuto en un ordenador estándar.

Referencias

• Caldwell, R. R., Kamionkowski, M., & Weinberg, N. N. (2003). Phantom energy: Dark energy with w < −1 causes a cosmic doomsday. Physical Review Letters, 91(7), 071301.
• DESI Collaboration (2024). DESI 2024 VI: Cosmological constraints from the measurements of baryon acoustic oscillations. arXiv:2404.03002.
• DESI Collaboration (2025). DESI DR2 results II: Measurements of baryon acoustic oscillations and cosmological constraints. arXiv:2503.14738.
• Dodelson, S. (2003). Modern Cosmology. Academic Press.
• Euclid Collaboration, Mellier, Y., et al. (2025). Euclid. I. Overview of the Euclid mission. Astronomy & Astrophysics, 697, A1.
• Friedmann, A. (1922). Über die Krümmung des Raumes. Zeitschrift für Physik, 10(1), 377–386.
• Gibbons, G. W., & Hawking, S. W. (1977). Cosmological event horizons, thermodynamics, and particle creation. Physical Review D, 15(10), 2738–2751.
• Harris, C. R., Millman, K. J., van der Walt, S. J., et al. (2020). Array programming with NumPy. Nature, 585, 357–362.
• Hogg, D. W. (1999). Distance measures in cosmology. arXiv:astro-ph/9905116.
• Hubble, E. (1929). A relation between distance and radial velocity among extra-galactic nebulae. Proceedings of the National Academy of Sciences, 15(3), 168–173.
• Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. Computing in Science & Engineering, 9(3), 90–95.
• Lemaître, G. (1927). Un univers homogène de masse constante et de rayon croissant rendant compte de la vitesse radiale des nébuleuses extra-galactiques. Annales de la Société Scientifique de Bruxelles, A47, 49–59.
• Penzias, A. A., & Wilson, R. W. (1965). A measurement of excess antenna temperature at 4080 Mc/s. The Astrophysical Journal, 142, 419–421.
• Perlmutter, S., et al. (1999). Measurements of Ω and Λ from 42 high-redshift supernovae. The Astrophysical Journal, 517(2), 565–586.
• Planck Collaboration (2020). Planck 2018 results. VI. Cosmological parameters. Astronomy & Astrophysics, 641, A6.
• Riess, A. G., et al. (1998). Observational evidence from supernovae for an accelerating universe and a cosmological constant. The Astronomical Journal, 116(3), 1009–1038.
• Spergel, D., et al. (2015). Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 report. arXiv:1503.03757.
• Steinhardt, P. J. (1997). Cosmological challenges for the 21st century. En V. L. Fitch, D. R. Marlow y M. A. E. Dementi (Eds.), Critical Problems in Physics (pp. 123–147). Princeton University Press.
• Tiesinga, E., Mohr, P. J., Newell, D. B., & Taylor, B. N. (2021). CODATA recommended values of the fundamental physical constants: 2018. Reviews of Modern Physics, 93(2), 025010.
• Zlatev, I., Wang, L., & Steinhardt, P. J. (1999). Quintessence, cosmic coincidence, and the cosmological constant. Physical Review Letters, 82(5), 896–899.