BLAST: Beyond Limber Angular power Spectra Toolkit. A fast and efficient algorithm for 3x2 pt analysis

Resumen Técnico: Blast.jl – Un algoritmo rápido y eficiente para el análisis de 3 × 2 pt

Planteamiento del problema
La llegada de los sondeos fotométricos y espectroscópicos de próxima generación (por ejemplo, Euclid, LSST, Roman, DESI) promete mediciones de alta precisión de la Estructura a Gran Escala (LSS). Estos sondeos dependen en gran medida de las estadísticas "3 × 2 pt", que combinan el agrupamiento de galaxias (galaxy clustering), las correlaciones de cizalladura de materia oscura (shear-shear correlations) y las correlaciones cruzadas de lente gravitacional galaxia-galaxia. El modelado teórico de estas estadísticas requiere calcular los espectros de potencia angulares ($C_\ell$) integrando el producto del espectro de potencia de materia 3D con funciones de Bessel esféricas altamente oscilantes sobre tres dimensiones.

El enfoque estándar para reducir el coste computacional es la aproximación de Limber, que simplifica la integral 3D a una integral 1D al asumir que los modos en diferentes épocas no contribuyen a la integración a lo largo de la línea de visión. Sin embargo, esta aproximación introduce errores significativos, particularmente en escalas angulares grandes (bajo $\ell$). Para explotar plenamente el poder de constreñimiento de los futuros sondeos y controlar los sistemas, son esenciales los cálculos precisos no-Limber. No obstante, realizar estas integrales 3D completas es computacionalmente prohibitivo para los métodos de inferencia de parámetros como el Monte Carlo por Cadenas de Markov (MCMC), que requieren millones de evaluaciones. Por ejemplo, análisis recientes utilizando la aproximación de Limber requirieron 0.5 millones de horas de CPU para converger.

Metodología
Los autores presentan Blast.jl, un nuevo algoritmo diseñado para calcular espectros de potencia angulares no-Limber completos de manera eficiente sin asumir una tasa de crecimiento dependiente de la escala. La metodología central se basa en la descomposición del espectro de potencia de materia utilizando polinomios de Chebyshev.

1. Descomposición de Chebyshev: El algoritmo aproxima el espectro de potencia de materia suave $P(k, \chi_1, \chi_2)$ en el intervalo $[k_{min}, k_{max}]$ mediante una expansión de serie de Chebyshev. Esto permite la separación de la integral en componentes geométricos y componentes cosmológicos.
2. Precomputación de términos geométricos: La integral se reestructura de modo que la parte más costosa computacionalmente —la integración del producto de dos funciones de Bessel esféricas y polinomios de Chebyshev— sea independiente de los parámetros cosmológicos. Estos términos ($\tilde{T}^{AB}_{n; \ell}$) se computan una sola vez y se almacenan.
3. Evaluación eficiente: Para cualquier nuevo modelo cosmológico, el algoritmo solo necesita recalcular los coeficientes de Chebyshev del espectro de potencia, un proceso que es computacionalmente económico. Las integrales exteriores restantes sobre la distancia comóvil ($\chi$) y la razón $R = \chi_2/\chi_1$ se realizan utilizando reglas de cuadratura estándar (regla de Simpson para $\chi$ y cuadratura de Clenshaw-Curtis para $R$).
4. Transformación de coordenadas: Se emplea un cambio de variables de $(\chi_1, \chi_2)$ a $(\chi, R)$ para muestrear mejor la región donde $\chi_1 \approx \chi_2$, la cual contribuye de forma más significativa a la integral.
5. Corrección no lineal: El método trata los espectros de potencia lineales y no lineales por separado. La contribución no lineal, que es significativa solo en escalas pequeñas donde la aproximación de Limber es válida, se corrige utilizando un enfoque basado en Limber, mientras que la integral completa no-Limber se aplica al componente lineal.

Contribuciones clave

• Innovación algorítmica: Blast.jl introduce un enfoque de descomposición utilizando polinomios de Chebyshev que evita las restricciones de espaciado de rejilla logarítmica de los métodos basados en FFTLog (como FKEM y matter) y la resolución de ecuaciones diferenciales iterativas del método de Levin.
• Estrategia de precomputación: Al aislar las integrales geométricas independientes de la cosmología, el algoritmo reduce drásticamente el tiempo requerido para las evaluaciones repetidas durante la inferencia de parámetros.
• Flexibilidad: A diferencia de algunos métodos competidores, Blast.jl no asume una tasa de crecimiento independiente de la escala, lo que lo hace aplicable a cosmologías con neutrinos masivos o gravedad modificada.
• Código abierto: El código está disponible públicamente en GitHub, acompañado de un repositorio que detalla su uso.

Resultos
El algoritmo fue evaluado contra el desafío "N5K" (Non-local No-Nonsense Non-Limber Numerical Knock-out), organizado por la Colaboración de Ciencia de Energía Oscura de LSST para evaluar los métodos de integración no-Limber. Las pruebas se realizaron en un nodo de computación de alto rendimiento (Narval) con 64 hilos, comparando Blast.jl contra los ganadores del desafío: FKEM (CosmoLike), matter y Levin.

• Velocidad: Con el requisito de precisión fiducial ($\Delta\chi^2 < 0.2$ para $\ell < 200$), Blast.jl logró un tiempo de ejecución de 0.00512 segundos usando 64 núcleos. Esto es aproximadamente 10 veces más rápido que FKEM (el método más rápido anterior, 0.05871 s) y aproximadamente 15 veces más rápido cuando se optimiza a 3-2 núcleos (0.00360 s).
• Escalabilidad: Blast.jl demuestra una excelente escalabilidad con respecto al número de cubetas de redshift tomográfico. Mientras que otros métodos muestran un aumento en el tiempo de ejecución con el número de cubetas, el tiempo de ejecución de Blast.jl permanece efectivamente independiente del número de cubetas (oscilando entre 2 y 15).
• Precisión: El método cumple con los requisitos de precisión del desafío. Los autores señalan que la precisión puede ajustarse modificando el número de polinomios de Chebyshev y los puntos de muestreo sin impactar significativamente el tiempo de ejecución.
• Compromiso entre precisión y velocidad: Reducir la precisión de punto flotante de Float64 a Float32 mejoró el tiempo de ejecución por un factor de 2 con un impacto insignificante en la precisión ($\Delta\chi^2$ aumentó de 0.1054 a 0.1064).

Significado
El artículo afirma que Blast.jl representa el algoritmo de integración no-Limber más rápido conocido actualmente disponible, superando a los métodos del estado del arte por más de un orden de magnitud en el escenario fiducial. Su principal significado radica en hacer computacionalmente factible el análisis cosmológico utilizando estadísticas 3 × 2 pt para los sondeos de próxima generación sin depender de la errónea aproximación de Limber.

Los autores enfatizan que esta velocidad es crítica para permitir la inferencia de parámetros basada en MCMC con los masivos conjuntos de datos esperados de los futuros sondeos. Además, el diseño del algoritmo permite futuras extensiones, como la inclusión de distorsiones de espacio de velocidades (redshift space distortions), sesgo de magnificación y correcciones relativistas, así como la integración de la lente de la CMB. Los autores también destacan planes para incorporar diferenciación automática, lo que permitiría el uso de métodos de inferencia basados en gradientes, acelerando aún más los procesos de análisis cosmológico.