Curvas de rotación y materia oscura

1. Motivación: la masa que no brilla

Si casi toda la masa de una galaxia espiral estuviera donde está su luz —en el bulbo y en el disco, que se extinguen a unos pocos kpc—, la velocidad circular fuera de esa región debería caer como la de los planetas: $v_c\propto r^{-1/2}$. No lo hace. Las curvas de rotación de galaxias espirales se mantienen planas, o incluso suben, hasta el último punto medido. Esa discrepancia, cuantificada por Vera Rubin y colaboradores en los años setenta y ochenta y por los mapas de hidrógeno neutro de Albert Bosma, es la evidencia dinámica más directa de la materia oscura. En esta lección aprendemos a medir, modelar y descomponer una curva de rotación, con la Vía Láctea como caso de estudio.

2. Qué mide una curva de rotación

2.1 De la velocidad radial a la velocidad circular

Un disco en rotación circular con velocidad $v_c(R)$, inclinado un ángulo $i$ respecto al plano del cielo ($i = 90°$ es de canto), proyecta sobre la visual, en un punto del disco a azimut $\phi$ medido desde el eje mayor, la velocidad

donde $v_{\rm sys}$ es la velocidad sistémica de la galaxia. A lo largo del eje mayor ($\phi = 0, \pi$) se obtiene $v_c(R) = |v_{\rm los}-v_{\rm sys}|/\sin i$. La incertidumbre en $i$ se propaga como $\delta v_c/v_c = \cot i\,\delta i$, por lo que Rubin, Ford y Thonnard escogieron galaxias de alta inclinación. Las velocidades se miden por efecto Doppler de líneas de emisión: H$\alpha$ y [N II] en el óptico (regiones H II), o la línea de 21 cm del hidrógeno neutro en radio, que se extiende mucho más allá del disco estelar.

2.2 La expectativa kepleriana

Si la masa $M$ está contenida dentro de un radio $R_M$, para $R$ mayor que $R_M$ se tiene $v_c = \sqrt{GM/R}$. Con la masa dinámica interior al radio solar de la Vía Láctea, $M(R_\odot) = v_c^2R_\odot/G = 9.9\times10^{10}\,M_\odot$ (con $v_c = 229$ km/s, $R_\odot = 8.122$ kpc), la velocidad en 25 kpc debería ser $229\sqrt{8.122/25} = 131$ km/s. Se mide $\approx 200$ km/s.

3. Las observaciones

3.1 Rubin, Ford y Thonnard (1980)

En ApJ 238, 471, Rubin, Ford y Thonnard publicaron espectros a lo largo del eje mayor de 21 galaxias Sc con radios de 4 a 122 kpc (para $H_0 = 50$ km/s/Mpc; con $H_0 = 70$ los radios se reducen en un factor 0.71), de NGC 4605 a UGC 2885. Sus conclusiones: todas las curvas suben rápidamente hasta $\sim 125$ km/s en $R\sim 5$ kpc y después más despacio; en general las curvas llegan al 83 % del radio isofotal $R_{25}$ y la mayoría siguen subiendo en el último punto medido. UGC 2885, la mayor de la muestra, ha completado menos de diez rotaciones en su parte externa desde el origen del universo y aun así muestra un patrón espiral regular. Ningún modelo con la masa distribuida como la luz reproduce estas curvas.

3.2 Bosma (1981) y el hidrógeno neutro

Los mapas de 21 cm de Bosma (tesis de 1978; AJ 86, 1825, 1981) extendieron las curvas más allá del disco óptico en 25 galaxias de distintos tipos: la mayoría permanecen planas hasta el radio máximo observado. Como el gas llega donde ya no hay estrellas, el argumento de que «la luz no sigue a la masa» se hizo inapelable. El índice $\alpha$ de la curva, $v_c\propto R^{\alpha}$, es $\approx 0$ en la región externa, mientras que la ley de Kepler exige $\alpha = -1/2$.

3.3 La Vía Láctea desde dentro: Eilers et al. (2019)

Medir la curva de rotación de nuestra propia galaxia es más difícil: estamos dentro del disco y no vemos el movimiento circular proyectado. Eilers, Hogg, Rix y Ness (2019) usaron las posiciones y velocidades tridimensionales de más de 23 000 gigantes rojas (Gaia DR2 + APOGEE + WISE + 2MASS) y la ecuación de Jeans (lección 3) con un potencial axisimétrico para obtener $v_c(R)$ entre 5 y 25 kpc (excluyen la región interior, dominada por la barra). Su resultado, con $R_\odot = 8.122\pm0.031$ kpc (GRAVITY Collaboration 2018):

con incertidumbres sistemáticas del 2–5 % en $v_c(R_\odot)$ y de 0.46 km/s/kpc en la pendiente. La curva declina suavemente: en 25 kpc vale $\approx 200$ km/s, muy por encima de los 131 km/s keplerianos. La masa dinámica interior a 25 kpc es $v_c^2R/G = 2.3\times10^{11}\,M_\odot$, unas cuatro veces la masa estelar total de la Galaxia ($5.4\times10^{10}\,M_\odot$ según McMillan 2017).

4. Descomposición: bulbo + disco + halo

4.1 Suma en cuadratura

Si el potencial es la suma de componentes, $\Phi = \sum_i\Phi_i$, también lo es la fuerza, y como $v_c^2 = R\,\partial\Phi/\partial R$,

Las velocidades circulares se suman en cuadratura, no linealmente.

4.2 El disco exponencial

El brillo superficial de los discos sigue $\Sigma(R) = \Sigma_0 e^{-R/R_d}$ (Freeman 1970), con escala $R_d$ y masa $M_d = 2\pi\Sigma_0R_d^2$. Un disco no es esférico y los teoremas de Newton no se aplican: la masa exterior sí ejerce fuerza. Para un disco infinitamente delgado, la velocidad circular en el plano es (Freeman 1970; Binney & Tremaine §2.6)

con $I_n$, $K_n$ las funciones de Bessel modificadas. La curva sube como cuerpo sólido, alcanza el máximo en $R\approx 2.2R_d$ y luego cae; en $R\gg R_d$ tiende a la kepleriana de masa $M_d$. (Las funciones de Bessel están en scipy.special; en el navegador se pueden usar las aproximaciones polinómicas de Abramowitz & Stegun §9.8, con error relativo menor que $10^{-7}$.)

4.3 Un modelo de tres componentes para la Vía Láctea

Tomemos un modelo de juguete con parámetros del orden de los ajustados en la literatura (McMillan 2017; Eilers et al. 2019): bulbo de Hernquist con $M_b = 10^{10}\,M_\odot$ y $a = 0.7$ kpc; disco exponencial con $M_d = 5\times10^{10}\,M_\odot$ y $R_d = 2.5$ kpc; halo NFW de Eilers et al. ($M_{\rm vir} = 7.25\times10^{11}\,M_\odot$, $c = 12.8$, $R_{\rm vir} = 189.3$ kpc). Las contribuciones, en km/s:

Sin halo, la curva cae de 185 a 104 km/s entre 8 y 25 kpc; con el halo NFW se mantiene entre 230 y 193, cerca de la medida. El ajuste no es perfecto (un modelo real incluye disco grueso, gas, y ajusta todos los parámetros a la vez), pero muestra el punto esencial: más allá de $\sim 14$ kpc la dinámica está dominada por la materia oscura (Eilers et al. 2019, §5.3), mientras que en el radio solar las estrellas aún aportan dos tercios de $v_c^2$.

El simulador reproduce este modelo de tres componentes. Compara el preset «Solo materia visible» con «Vía Láctea con halo» y la recta de Eilers et al. (2019); luego sube la masa del disco y baja la del halo buscando otra combinación que pase igual de cerca de los datos: es la degeneración de la §4.4.

4.4 La degeneración disco–halo

El ajuste a una curva de rotación no determina unívocamente el reparto de masa: un disco más masivo con un halo menos concentrado produce casi la misma $v_c(R)$ (la llamada degeneración disco–halo, o hipótesis del «disco máximo»). Para romperla se recurre a la dispersión vertical de velocidades del disco (que fija su densidad superficial local), a la relación masa–luminosidad de las poblaciones estelares, o —en la Vía Láctea— a la densidad local de materia oscura. Eilers et al. (2019) fijan las componentes estelares a un modelo previo y ajustan solo el halo; obtienen $\rho_{\rm dm}(R_\odot) = 0.30\pm0.03$ GeV cm$^{-3}$ (McMillan 2017: $0.40\pm0.04$; la diferencia refleja justamente esa degeneración).

5. La velocidad circular del Sol

5.1 Tres ingredientes

La velocidad del Sol respecto al centro galáctico tiene una componente circular $\Theta_0 = v_c(R_\odot)$ más el movimiento peculiar $(U,V,W)_\odot$ respecto al estándar local de reposo (LSR, la órbita circular en $R_\odot$). Schönrich, Binney y Dehnen (2010) obtuvieron a partir de Hipparcos $(U,V,W)_\odot = (11.1^{+0.69}_{-0.75},\,12.24^{+0.47}_{-0.47},\,7.25^{+0.37}_{-0.36})$ km/s, con sistemáticos adicionales de (1, 2, 0.5) km/s. $V_\odot = 12.24$ km/s es el exceso del Sol sobre la rotación circular.

5.2 El movimiento propio de Sgr A*

Reid y Brunthaler (2020) midieron con el VLBA, durante 18 años, el movimiento propio aparente de la radiofuente Sgr A* respecto a cuásares de fondo: $\mu_l = -6.411\pm0.008$ mas/año a lo largo del plano galáctico y $\mu_b = -0.219\pm0.007$ mas/año hacia el polo norte galáctico. Sgr A* está en reposo en el centro (su velocidad propia es menor que $\sim 1$ km/s); el movimiento aparente es el reflejo de la órbita del Sol. La velocidad tangencial total del Sol es

(el factor $4.74$ convierte mas/año × kpc en km/s: $1\ \text{kpc}\times1\ \text{mas} = 1\ \text{UA}$, y $1$ UA/año $= 4.74$ km/s). Con $R_\odot = 8.122$ kpc: $v_{\rm t} = 246.8$ km/s, y $\Theta_0 = v_{\rm t} - V_\odot = 234.6$ km/s. Con $R_\odot = 8.275\pm0.034$ kpc (GRAVITY Collaboration 2022): $v_{\rm t} = 251.5$ km/s y $\Theta_0 = 239.2$ km/s.

5.3 Propagación de incertidumbres y comparación

Como $v_{\rm t}\propto\mu_lR_\odot$, $\delta v_{\rm t}/v_{\rm t} = \sqrt{(\delta\mu_l/\mu_l)^2 + (\delta R_\odot/R_\odot)^2}$: con $\delta\mu_l/\mu_l = 0.12\,\%$ y $\delta R_\odot/R_\odot = 0.4\,\%$, $\delta v_{\rm t} = 1.0$ km/s. La incertidumbre dominante es la de $V_\odot$ ($\pm2$ km/s sistemáticos) y, sobre todo, la de $R_\odot$ entre autores: 8.122 (GRAVITY 2018) frente a 8.275 kpc (GRAVITY Collaboration 2022) cambia $\Theta_0$ en 4.6 km/s. Las distintas determinaciones se agrupan así:

Un valor de $\Theta_0\approx 230$–$240$ km/s, con $R_\odot\approx 8.1$–$8.3$ kpc, da un periodo orbital del Sol $2\pi R_\odot/v_{\rm t}\approx 200$ millones de años.

6. Análisis de error en una curva de rotación

• Inclinación. $\delta v_c/v_c = \cot i\,\delta i$. Para $i = 60°$ y $\delta i = 3°$ ($0.052$ rad): 3 %. Para $i = 30°$: 9 %.
• Distancia. Las velocidades no dependen de la distancia, pero los radios sí ($R\propto D$) y las masas también ($M = v_c^2R/G\propto D$).
• Movimientos no circulares. Barras, brazos espirales y alabeos producen desviaciones de 10–20 km/s; de ahí que Eilers et al. excluyan $R$ menor que 5 kpc.
• Asimetría de deriva. Las estrellas no rotan exactamente a $v_c$: su velocidad media azimutal $\bar v_\phi$ es menor en una cantidad que depende de la dispersión de velocidades. La ecuación de Jeans (lección 3) corrige este sesgo; ignorarlo subestima $v_c$ en $\sim 5$–10 km/s para las gigantes.

7. Problemas

1. Kepler frente a observación. Con $v_c(R_\odot) = 229$ km/s y $R_\odot = 8.122$ kpc, calcula la masa interior al Sol y la velocidad kepleriana esperada en 15, 20 y 25 kpc. Compara con la recta de Eilers et al.

2. Suma en cuadratura (demostración). Demuestra que si $\Phi = \Phi_1+\Phi_2$ entonces $v_c^2 = v_{c,1}^2+v_{c,2}^2$ y que, en cambio, las velocidades circulares no se suman linealmente. Muestra con el modelo de §4.3 en $R_\odot$ cuánto se equivocaría quien sumara linealmente.

3. Velocidad del Sol. A partir de $\mu_l = 6.411\pm0.008$ mas/año, $R_\odot = 8.275\pm0.034$ kpc y $V_\odot = 12.24\pm2$ km/s, calcula $\Theta_0$ con su incertidumbre y el periodo orbital del Sol.

4. Inclinación. Una galaxia tiene $v_{\rm los}-v_{\rm sys} = 180\pm3$ km/s en el eje mayor e $i = 55°\pm4°$. Obtén $v_c$ con su incertidumbre. ¿Qué pasa si $i$ fuera $20°$?

5. Disco máximo. Demuestra que, para el disco exponencial, $v_{c,d}(R)\to\sqrt{GM_d/R}$ cuando $R\gg R_d$ usando las expansiones asintóticas $I_n(y)K_n(y)\simeq\frac{1}{2y}\left[1 + \frac{(4n^2-1)\cdot(-1)}{8y^2}\cdot\ldots\right]$; basta con el término dominante de la diferencia $I_0K_0-I_1K_1\simeq\frac{1}{4y^3}$.

6. Computacional: ajuste disco + halo. Toma la recta de Eilers et al. muestreada cada kpc entre 5 y 25 kpc con error de 2 km/s. Fija el bulbo ($10^{10}\,M_\odot$, 0.7 kpc) y el disco ($R_d = 2.5$ kpc) y ajusta por mínimos cuadrados $M_d$, $M_{\rm vir}$ y $c$ del halo NFW (con $R_{\rm vir}$ ligado a $M_{\rm vir}$ por $M_{\rm vir} = \frac{4}{3}\pi R_{\rm vir}^3\cdot200\rho_{\rm crit}$, $\rho_{\rm crit} = 136\ M_\odot\,\text{kpc}^{-3}$). Dibuja las curvas y el mapa de $\chi^2$ en el plano $(M_d, M_{\rm vir})$ para visualizar la degeneración.

import numpy as np
from scipy.special import i0, i1, k0, k1
from scipy.optimize import least_squares
G = 4.30091e-6  # kpc (km/s)^2 / Msun
R = np.arange(5, 26.0); v_obs = 229.0 - 1.7*(R - 8.122); err = 2.0
def v_disk(R, Md, Rd):
    S0 = Md/(2*np.pi*Rd**2); y = R/(2*Rd)
    return np.sqrt(4*np.pi*G*S0*Rd*y**2*(i0(y)*k0(y) - i1(y)*k1(y)))
def v_hern(R, M, a): return np.sqrt(G*M*R/(R + a)**2)
def v_nfw(R, Mvir, c):
    Rvir = (3*Mvir/(4*np.pi*200*136.0))**(1/3); rs = Rvir/c
    g = lambda x: np.log1p(x) - x/(1 + x)
    return np.sqrt(G*Mvir*g(R/rs)/g(c)/R)
def model(p):
    Md, Mvir, c = p
    return np.sqrt(v_hern(R, 1e10, 0.7)**2 + v_disk(R, Md, 2.5)**2 + v_nfw(R, Mvir, c)**2)
res = least_squares(lambda p: (model(p) - v_obs)/err, x0=[5e10, 7e11, 12])
print(res.x, (res.fun**2).sum())

8. Autoevaluación

9. Lecturas

• Rubin, Ford & Thonnard (1980), ApJ 238, 471, doi:10.1086/158003. El artículo clásico; lee el resumen y la figura de las 21 curvas.
• Bosma (1981), AJ 86, 1825. Las curvas de 21 cm.
• Eilers, Hogg, Rix & Ness (2019), ApJ 871, 120, arXiv:1810.09466, §3 (método de Jeans), §5 (descomposición).
• Reid & Brunthaler (2020), ApJ 892, 39, arXiv:2001.04386; Schönrich, Binney & Dehnen (2010), MNRAS 403, 1829; GRAVITY Collaboration (2022), A&A 657, L12.
• Binney & Tremaine (2008), §2.6 (potenciales de discos, fórmula de Freeman) y §1.1.2 (la Galaxia).
• Bertone & Hooper (2018), Rev. Mod. Phys. 90, 045002: historia de la materia oscura, §IV para las curvas de rotación.