Potenciales galácticos: Poisson, Plummer, Hernquist, isotermo y NFW

1. Motivación: la gravedad de 10¹¹ estrellas

Una galaxia contiene del orden de $10^{11}$ estrellas, pero ninguna de ellas «siente» a sus vecinas una a una: como veremos en la lección 3, el tiempo que tardarían los encuentros binarios en desviar una órbita supera en muchos órdenes de magnitud la edad del universo. Cada estrella se mueve, con excelente aproximación, en el campo medio $\Phi(\mathbf x)$ creado por la distribución suave de masa $\rho(\mathbf x)$. Toda la dinámica galáctica descansa, por tanto, en una pregunta: dada $\rho$, ¿cuál es $\Phi$? Y en su recíproca, la que responde la materia oscura: dada la cinemática (y por tanto $\Phi$), ¿cuál es $\rho$?

En esta lección construimos la herramienta matemática —la ecuación de Poisson— y la aplicamos a los cuatro modelos esféricos que aparecen una y otra vez en la literatura: el de Plummer (cúmulos globulares), el de Hernquist (bulbos y elípticas), la esfera isoterma (curvas de rotación planas) y el perfil de Navarro–Frenk–White (halos de materia oscura en simulaciones cosmológicas).

2. La ecuación de Poisson

2.1 Del campo de una masa puntual al de una distribución continua

Sea $\rho(\mathbf x')$ la densidad de masa. La fuerza por unidad de masa en $\mathbf x$ es la suma de las contribuciones newtonianas de cada elemento $\rho(\mathbf x')\,d^3x'$:

donde $G$ es la constante de gravitación y $\Phi$ el potencial gravitatorio (energía potencial por unidad de masa). Para obtener una ecuación local tomamos la divergencia:

El integrando se anula para $\mathbf x\neq\mathbf x'$ (compruébalo: es la divergencia de un campo $\propto 1/r^2$ radial) y es singular en $\mathbf x = \mathbf x'$. Integrando sobre una bola $B_\epsilon$ de radio $\epsilon$ alrededor de $\mathbf x$ y aplicando el teorema de la divergencia,

de modo que la divergencia vale $-4\pi\,\delta^3(\mathbf x-\mathbf x')$ y

Esta es la ecuación de Poisson. En las regiones vacías se reduce a la ecuación de Laplace $\nabla^2\Phi = 0$. Integrada sobre un volumen $V$ da el teorema de Gauss gravitatorio:

2.2 Simetría esférica: los teoremas de Newton

Si $\rho = \rho(r)$, el laplaciano se reduce a su parte radial y la ecuación de Poisson queda

Aplicando Gauss a una esfera de radio $r$: $4\pi r^2\,\Phi'(r) = 4\pi G M(r)$, con $M(r) = 4\pi\int_0^r\rho(r')r'^2dr'$ la masa encerrada. Así:

donde $v_c$ es la velocidad circular: la de una órbita circular de radio $r$ (fuerza centrípeta $v_c^2/r$ igual a la atracción $GM(r)/r^2$). Este es el primer teorema de Newton (el campo en $r$ solo depende de la masa interior) y el segundo (la masa exterior no ejerce fuerza neta). Integrando desde el infinito con $\Phi(\infty) = 0$ y partiendo la integral en capas interiores y exteriores,

El segundo término es la contribución de las capas exteriores: no ejercen fuerza, pero sí contribuyen al potencial (una constante en cada capa).

2.3 Energía potencial total

Con $W = \tfrac12\int\rho\,\Phi\,d^3x$ y una integración por partes se llega a la forma útil

$W$ es la energía necesaria para dispersar el sistema al infinito, y aparecerá en el teorema del virial (lección 3).

3. Cuatro pares potencial–densidad

A cada modelo le exigimos lo mismo: dado $\Phi$, calcular $\rho$ por Poisson (o al revés), obtener $M(r)$ y $v_c(r)$, y entender el comportamiento asintótico. Todas las fórmulas están implementadas en src/lib/math/galaxias.ts, cuyo test comprueba numéricamente $\nabla^2\Phi = 4\pi G\rho$ en cada caso.

3.1 Plummer (1911)

Plummer propuso para los cúmulos globulares el potencial de una masa puntual «suavizada» a una escala $b$:

Derivando, $\Phi' = GMr/(r^2+b^2)^{3/2}$, luego $r^2\Phi' = GMr^3(r^2+b^2)^{-3/2}$ y

Dividiendo por $4\pi Gr^2$:

La densidad es finita en el centro ($\rho_0 = 3M/4\pi b^3$) y cae como $r^{-5}$: la masa total converge a $M$. El radio de media masa satisface $r_h^3/(r_h^2+b^2)^{3/2} = \tfrac12$, es decir $r_h = b/\sqrt{2^{2/3}-1} \approx 1.305\,b$. La velocidad circular crece linealmente cerca del centro (cuerpo sólido), alcanza su máximo en $r = \sqrt2\,b$ y decae como kepleriana ($r^{-1/2}$) a gran distancia. La energía potencial es $W = -3\pi GM^2/(32\,b)$ (problema 2).

Ejemplo. Un cúmulo globular con $M = 5\times10^5\,M_\odot$ y $b = 3\ \text{pc}$ tiene $\rho_0 = 4.4\times10^3\ M_\odot\,\text{pc}^{-3}$ y $v_{c,\max} = 16.6\ \text{km/s}$ en $r = 4.2\ \text{pc}$ (cálculo con densidadPlummer y velocidadCircularPlummer).

3.2 Hernquist (1990)

Para elípticas y bulbos, cuya luz sigue la ley $R^{1/4}$ de de Vaucouleurs, Hernquist encontró un modelo con toda la dinámica en forma cerrada:

Derivación. $\Phi' = GM/(r+a)^2$, luego $r^2\Phi' = GMr^2/(r+a)^2$ y $\frac{d}{dr}[r^2\Phi'] = GM\left[\frac{2r}{(r+a)^2}-\frac{2r^2}{(r+a)^3}\right] = \frac{2GMra}{(r+a)^3}$; dividiendo por $4\pi Gr^2$ se obtiene $\rho_H$. La masa encerrada sale de $M(r) = r^2\Phi'/G$. $\square$

La densidad diverge como $r^{-1}$ en el centro (una cúspide, integrable: $M(r)\propto r^2$) y cae como $r^{-4}$ fuera. El radio de media masa es $r_h = (1+\sqrt2)\,a$; la velocidad circular máxima, $v_c^2(a) = GM/4a$, se alcanza en $r = a$; y $W = -GM^2/6a$.

Ejemplo. Un bulbo con $M = 10^{10}\,M_\odot$ y $a = 0.7\ \text{kpc}$ tiene $v_{c,\max} = 124\ \text{km/s}$ en $r = 0.7\ \text{kpc}$ y aporta solo $67\ \text{km/s}$ en $r = 8\ \text{kpc}$: el bulbo no sostiene la curva de rotación del Sol.

3.3 La esfera isoterma singular

Si la distribución de velocidades es maxwelliana con dispersión unidimensional $\sigma$ (la misma a todo radio, de ahí «isoterma»), la ecuación de Jeans (lección 3) da $\rho \propto \exp(-\Phi/\sigma^2)$ y Poisson se convierte en una ecuación no lineal (de Emden) sin solución elemental. Sí tiene una solución singular de ley de potencia, que se comprueba por sustitución:

Verificación. $M(r) = 4\pi\int_0^r\frac{\sigma^2}{2\pi G}dr' = 2\sigma^2 r/G$; $v_c^2 = GM/r = 2\sigma^2$; $\Phi' = v_c^2/r = 2\sigma^2/r$. Y $\frac{1}{r^2}\frac{d}{dr}(r^2\cdot 2\sigma^2/r) = 2\sigma^2/r^2 = 4\pi G\rho$. $\square$

Es el modelo más simple de curva de rotación plana, pero la masa diverge linealmente: cualquier halo real debe truncarse. Lo que la esfera isoterma enseña es la correspondencia $v_c = \text{cte} \iff \rho\propto r^{-2}$.

3.4 Navarro, Frenk y White (1996, 1997)

En simulaciones de N cuerpos de materia oscura fría, Navarro, Frenk y White hallaron que los halos en equilibrio, desde escalas de galaxias enanas hasta cúmulos, se ajustan con el mismo perfil de dos pendientes:

con $r_s$ el radio de escala y $\rho_0$ una densidad característica. Como Hernquist, tiene cúspide $r^{-1}$, pero decae solo como $r^{-3}$: la masa diverge logarítmicamente y hay que truncar en el radio virial $R_{\rm vir}$. La masa encerrada se integra con el cambio $u = 1+x$:

El potencial sale de la fórmula general de §2.2; con $\int_r^\infty\rho\,r'\,dr' = \rho_0 r_s^2/(1+x)$:

donde la simplificación usa $\frac{x}{1+x}\cdot\frac{1}{x} = \frac{1}{1+x}$. La velocidad circular, $v_c^2 = 4\pi G\rho_0r_s^2\,g(x)/x$, alcanza su máximo donde $d[g(x)/x]/dx = 0$, es decir en $x_{\max} = 2.16258$ (raíz de $\ln(1+x) = x(1+2x)/(1+x)^2$, problema 4). Los halos se parametrizan por $M_{\rm vir}$ y la concentración $c \equiv R_{\rm vir}/r_s$, con

Ejemplo con datos reales. Eilers et al. (2019, ApJ 871, 120) ajustan a la curva de rotación de la Vía Láctea un halo NFW con $M_{\rm vir} = (7.25\pm0.26)\times10^{11}\,M_\odot$, $c = 12.8\pm0.3$ y $R_{\rm vir} = 189.3\pm2.2\ \text{kpc}$. Entonces $r_s = 14.8\ \text{kpc}$, $\rho_0 = 1.05\times10^7\ M_\odot\,\text{kpc}^{-3}$ (el artículo da $1.06\pm0.09$) y, en el radio solar $R_\odot = 8.122\ \text{kpc}$, la densidad de materia oscura es $\rho_{\rm dm}(R_\odot) = 0.30\ \text{GeV}\,\text{cm}^{-3}$, exactamente el valor que publican ($0.30\pm0.03$). El halo aporta $v_{c,\rm halo}(R_\odot) = 137\ \text{km/s}$, que sumado en cuadratura a las componentes estelares reproduce los $229\ \text{km/s}$ medidos (lección 2). El test galaxias.test.ts verifica estos tres números.

3.5 Resumen comparativo

La forma de $v_c(r)$ es la huella del perfil: $v_c\propto r$ si $\rho$ es constante, $v_c =$ cte si $\rho\propto r^{-2}$ y $v_c\propto r^{-1/2}$ fuera de toda la masa. Es el puente con la lección siguiente.

4. Análisis numérico: verificar Poisson con diferencias finitas

La función laplacianoRadial(Phi, r, h) de galaxias.ts evalúa $\frac{1}{r^2}\frac{d}{dr}(r^2\Phi')$ con diferencias centradas anidadas. Cada derivada centrada tiene error de truncamiento $O(h^2)$, pero el error de redondeo crece como $\epsilon_{\rm m}|\Phi|/h^2$ al anidar dos derivadas. El error total, $E(h)\approx C_1h^2 + C_2\epsilon_{\rm m}|\Phi|/h^2$, se minimiza en $h_{\rm opt}\propto\epsilon_{\rm m}^{1/4}r$: con doble precisión ($\epsilon_{\rm m}\approx 2\times10^{-16}$), $h/r\sim10^{-4}$, que es el valor por defecto. El test exige concordancia con $4\pi G\rho$ a cinco cifras; con $h/r = 10^{-2}$ domina el truncamiento y con $h/r = 10^{-7}$ domina el redondeo (problema 6).

Para la masa encerrada el test usa la regla de Simpson con $n = 2000$ subintervalos; su error $O(h^4)$ es despreciable salvo en los perfiles con cúspide, donde $\rho r^2\propto r$ cerca del origen no es suave en $r=0$ y conviene desplazar el límite inferior o usar la variable $u = \ln r$.

5. Problemas

1. Potencial de una capa. Demuestra, a partir de la fórmula general de §2.2, que una capa esférica delgada de masa $m$ y radio $R$ produce $\Phi = -Gm/R$ constante en su interior y $-Gm/r$ fuera.

2. Energía potencial de Plummer (demostración). Demuestra que $W = -\dfrac{3\pi GM^2}{32\,b}$.

3. Dispersión de velocidades virial. Usando $2K + W = 0$ con $K = \tfrac32 M\langle\sigma^2\rangle$ (isótropo), calcula $\langle\sigma^2\rangle$ para el cúmulo globular del ejemplo de §3.1 ($M = 5\times10^5\,M_\odot$, $b = 3$ pc).

4. Máximo de la curva NFW. Demuestra que $v_c$ de NFW es máxima en la raíz de $\ln(1+x) = \dfrac{x(1+2x)}{(1+x)^2}$ y resuélvela numéricamente por Newton–Raphson.

5. Masa de un halo NFW. Para el halo de Eilers et al. (2019) calcula $M(r)$ en 8.122, 20 y 50 kpc y compara con la masa total $v_c^2R/G$ que implica $v_c(R_\odot) = 229$ km/s.

6. Computacional: error del laplaciano discreto. Para el potencial de Hernquist con $M = 10^{11}\,M_\odot$, $a = 2$ kpc y $r = 2$ kpc, calcula el error relativo de laplacianoRadial frente a $4\pi G\rho$ para $h/r = 10^{-1}, 10^{-2}, \dots, 10^{-8}$ y representa $\log E$ frente a $\log h$. Identifica las pendientes $+2$ (truncamiento) y $-2$ (redondeo) y el óptimo.

import numpy as np
G=6.6743e-11; Ms=1.98841e30; kpc=3.0857e19
M=1e11*Ms; a=2*kpc; r=2*kpc
Phi=lambda x: -G*M/(x+a)
rho=lambda x: M/(2*np.pi)*a/(x*(x+a)**3)
def lap(r,h):
    d=lambda x:(Phi(x+h)-Phi(x-h))/(2*h)
    f=lambda x:x*x*d(x)
    return (f(r+h)-f(r-h))/(2*h)/r**2
for k in range(1,9):
    h=r*10.0**-k
    print(k, abs(lap(r,h)/(4*np.pi*G*rho(r))-1))

7. Modelo de dos componentes. Construye $v_c(r)$ para bulbo de Hernquist ($10^{10}\,M_\odot$, $a = 0.7$ kpc) más halo NFW de Eilers et al. y dibuja ambas contribuciones y su suma en cuadratura entre 0.1 y 100 kpc. ¿Qué falta para reproducir 229 km/s en 8.1 kpc?

6. Autoevaluación

7. Lecturas

• Binney & Tremaine (2008), Galactic Dynamics, §2.1 (teoría del potencial, teoremas de Newton) y §2.2 (modelos esféricos: Plummer, isócrono, Hernquist, NFW y la familia de dos potencias). Las tablas del capítulo 2 recogen las fórmulas de esta lección.
• Sparke & Gallagher (2007), Galaxies in the Universe, cap. 3, §3.1–3.2.
• Plummer (1911), MNRAS 71, 460; Hernquist (1990), ApJ 356, 359, doi:10.1086/168845; Navarro, Frenk & White (1996), ApJ 462, 563; (1997), ApJ 490, 493, doi:10.1086/304888.
• Eilers, Hogg, Rix & Ness (2019), ApJ 871, 120, doi:10.3847/1538-4357/aaf648 (arXiv:1810.09466), §5.3 para el ajuste NFW usado en los ejemplos.