Politropos y la ecuación de Lane–Emden

1. Planteamiento: dos ecuaciones y una relación de cierre

En la lección 0 vimos que la estructura de una estrella exige cuatro ecuaciones acopladas más tres funciones materiales. Aquí estudiamos el modelo más sencillo que conserva la mecánica completa: cerramos las dos ecuaciones mecánicas con una relación directa entre presión y densidad, sin temperatura ni transporte. A cambio de esa simplificación obtenemos una sola EDO adimensional cuya solución describe todas las estrellas de una familia, y que anticipa el resultado más célebre de la astrofísica estelar: la masa de Chandrasekhar.

Consideremos una estrella esférica, estática y no rotante. Sea $r$ la distancia al centro, $\rho(r)$ la densidad de masa, $P(r)$ la presión y $m(r)$ la masa encerrada en una esfera de radio $r$. El equilibrio hidrostático de una capa de espesor $dr$ y la conservación de la masa dan

donde $G$ es la constante de la gravitación. Despejando $m = -\dfrac{r^2}{G\rho}\dfrac{dP}{dr}$ y derivando respecto de $r$ obtenemos una sola ecuación de segundo orden, que no es otra cosa que la ecuación de Poisson $\nabla^2\Phi = 4\pi G\rho$ en simetría esférica (con $d\Phi/dr = Gm/r^2$ y $dP = -\rho\,d\Phi$):

La ecuación (1) tiene dos funciones incógnita, $P$ y $\rho$. En general se cierra con la ecuación de estado, el transporte de energía y la producción nuclear, que introducen la temperatura. Un politropo es el modelo en el que se postula una relación directa

con $K$ constante y $n$ el índice politrópico. Casos físicos importantes:

• $n = 3/2$: gas ideal en una estrella totalmente convectiva (adiabática con $\gamma = 5/3$), y gas de electrones degenerados no relativistas ($P \propto \rho^{5/3}$): enanas blancas poco masivas, enanas marrones.
• $n = 3$: gas de electrones degenerados ultrarrelativistas ($P\propto\rho^{4/3}$) y el modelo estándar de Eddington (presión de radiación en proporción constante a la total).
• $n = 0$: densidad uniforme (fluido incompresible). $n = 1$: aproximación tosca a las estrellas de neutrones.

2. Derivación de la ecuación de Lane–Emden

Escribimos la densidad como $\rho = \rho_c\,\theta^n$, donde $\rho_c$ es la densidad central y $\theta(r)$ es una función adimensional con $\theta(0) = 1$. Por (2),

Sustituyendo en (1):

El coeficiente tiene dimensiones de longitud al cuadrado. Definimos la escala de longitud

y llegamos a la ecuación de Lane–Emden (Lane 1870, Emden 1907):

con las condiciones de contorno $\theta(0) = 1$ (definición de $\rho_c$) y $\theta'(0) = 0$. Esta última se sigue de la regularidad: $dP/dr = -Gm\rho/r^2 \approx -\tfrac{4\pi}{3}G\rho_c^2\,r \to 0$ cuando $r\to0$. La superficie de la estrella es la primera raíz $\xi_1$ de $\theta$: allí $\rho = P = 0$, y el radio es

Observa la virtud de la adimensionalización: una sola función $\theta_n(\xi)$ describe todas las estrellas politrópicas de índice $n$, cualesquiera que sean $K$ y $\rho_c$. Las estrellas de igual $n$ son homólogas entre sí.

3. Masa, concentración y presión central

La masa encerrada es, usando (4) para integrar,

En la superficie, $M = 4\pi\alpha^3\rho_c\,\omega_n$ con $\omega_n \equiv -\xi_1^2\theta'(\xi_1)$. Dividiendo por $\tfrac{4}{3}\pi R^3 = \tfrac43\pi\alpha^3\xi_1^3$ se obtiene la concentración central

De (3), $K\rho_c^{1/n} = 4\pi G\alpha^2\rho_c/(n+1)$, de modo que $P_c = K\rho_c^{1+1/n} = 4\pi G\alpha^2\rho_c^2/(n+1)$. Eliminando $\alpha = R/\xi_1$ y $\rho_c$ mediante (6):

Valores calculados con src/lib/math/polytrope.ts y contrastados con Chandrasekhar (1939, tabla 4) y Horedt (2004). La concentración crece rapidísimo con $n$: los politropos de índice alto tienen un núcleo denso y una envoltura extensa y tenue.

El Sol como politropo. Con $M_\odot$ y $R_\odot$ nominales (IAU 2015), $\bar\rho_\odot = 1.41\times10^3\ \mathrm{kg/m^3}$. Para $n = 3$, (6) y (7) dan $\rho_c = 7.64\times10^4\ \mathrm{kg/m^3}$ y $P_c = 1.25\times10^{16}\ \mathrm{Pa}$, frente a $\rho_c = 1.527\times10^5\ \mathrm{kg/m^3}$ y $P_c = 2.342\times10^{16}\ \mathrm{Pa}$ del modelo solar estándar BP2000 (Bahcall, Pinsonneault & Basu 2001): un factor 2, nada mal para una sola ecuación sin temperatura. Con $n = 3/2$ el resultado es mucho peor ($\rho_c = 8.4\times10^3$, $P_c = 8.7\times10^{14}$ Pa): el Sol no es una estrella convectiva. El desacuerdo residual del $n = 3$ refleja que el núcleo solar tiene composición alterada por la combustión de hidrógeno ($\mu$ crece hacia el centro, lo que concentra la masa más de lo que permite un politropo) y que la envoltura convectiva ($n\approx1.5$) tiene otra estructura. La concentración real, $\rho_c/\bar\rho = 108$, corresponde a un $n$ efectivo de $\approx3.3$.

4. Comportamiento cerca del centro: la serie

En $\xi = 0$ el coeficiente $2/\xi$ es singular: $\xi=0$ es un punto singular regular de (4). Buscamos una serie par (la ecuación es invariante bajo $\xi\to-\xi$) $\theta = 1 + a_2\xi^2 + a_4\xi^4 + a_6\xi^6 + \dots$ Entonces

Igualando coeficientes con signo cambiado: $6a_2 = -1$, $20a_4 = -n a_2$, $42a_6 = -n a_4 - \tfrac{n(n-1)}{2}a_2^2$, es decir

El primer término dice que cerca del centro toda estrella politrópica es parabólica, $\rho \approx \rho_c(1 - n\xi^2/6)$, y que $\theta''(0) = -1/3$: el término $2\theta'/\xi$ tiene límite finito $-2/3$, aunque numéricamente sea $0/0$. Compárese con la lección 0: $P\approx P_c - \tfrac{2\pi}{3}G\rho_c^2r^2$ es exactamente $\theta^{n+1}\approx1-(n+1)\xi^2/6$ con (3).

5. Soluciones analíticas

Solo para $n = 0, 1, 5$ existen soluciones en forma cerrada que cumplan $\theta(0)=1,\ \theta'(0)=0$:

Para $n=1$, el cambio $u = \xi\theta$ transforma (4) en $u'' + u = 0$ (problema 1). Para $n = 5$, $\theta\gt0$ para todo $\xi$: el radio es infinito, pero la masa es finita, porque $-\xi^2\theta' = \tfrac{\xi^3}{3}(1+\xi^2/3)^{-3/2}\to\sqrt3$. Se puede demostrar que para $n\ge5$ no hay raíz finita y que para $0\le n\lt 5$ sí la hay; $n=5$ es la frontera.

6. Relación masa–radio y masa de Chandrasekhar

De $M = 4\pi\alpha^3\rho_c\omega_n$ y $R = \alpha\xi_1$ se tiene $\rho_c \propto M R^{-3}$. Llevando esto a (3), $R^2\propto K\rho_c^{(1-n)/n}\propto K\,(MR^{-3})^{(1-n)/n}$, de donde

• $n = 3/2$ (enanas blancas no relativistas): $R\propto M^{-1/3}$. ¡Más masa, menor radio!
• $n = 1$: $R$ independiente de $M$ (de hecho $R = \pi\alpha$ con $\alpha$ fijado por $K$).
• $n = 3$: el exponente de $R$ desaparece y (9) degenera: la masa queda fijada por $K$,

Para un gas de electrones completamente degenerado y ultrarrelativista (lección 4), $K = \dfrac{hc}{8}\left(\dfrac{3}{\pi}\right)^{1/3}\dfrac{1}{(\mu_e m_u)^{4/3}}$, donde $h$ es la constante de Planck, $c$ la velocidad de la luz, $m_u$ la unidad de masa atómica y $\mu_e$ el número de nucleones por electrón ($\mu_e = 2$ para He, C u O). Sustituyendo en (10) con $\omega_3 = 2.01824$ y constantes CODATA 2018:

(Con la masa del hidrógeno $m_H$ en lugar de $m_u$ se obtiene el «1.44» de muchos textos; Chandrasekhar 1931 obtuvo $0.91\,M_\odot$ con el $\mu_e = 2.5$ entonces en uso.) Interpretación: al aumentar la masa de una enana blanca, el radio disminuye (rama $n=3/2$), la densidad crece, los electrones se vuelven relativistas y el índice efectivo tiende a 3. En el límite, el radio tiende a cero para una masa finita: no hay enanas blancas en equilibrio con $M > M_{\mathrm{Ch}}$ (Chandrasekhar 1931, 1939 cap. XI). La lección 4 integra el modelo completo.

7. Resolución numérica

Para $n$ general, (4) se integra numéricamente. Tres decisiones merecen justificación:

1. Arranque con serie. No se puede evaluar el lado derecho $f(\xi,\theta,\theta') = -\theta^n - 2\theta'/\xi$ en $\xi = 0$. Arrancamos en $\xi_0 = 10^{-3}$ con (8): el error de truncamiento es $O(\xi_0^8)\sim10^{-24}$, muy por debajo del redondeo. Arrancar en $\xi_0$ con $\theta = 1$, $\theta'=0$ (sin la serie) introduce en cambio un error $O(\xi_0^2)$ que se propaga a $\xi_1$; con RK4 de orden 4 ese error domina para cualquier paso razonable.
2. Paso adaptativo. Para $n\to5$, $\xi_1$ crece sin cota ($\xi_1\approx 347$ para $n=4.95$), mientras que cerca del centro la solución varía en escalas $\sim1$. Usamos RK4 con duplicación de paso: la diferencia entre un paso $h$ y dos pasos $h/2$ estima el error local, $\varepsilon \approx |y_{h/2}-y_h|/15$, y $h$ se reescala por $(\mathrm{tol}/\varepsilon)^{1/5}$.
3. La raíz y $\theta\lt0$. Para $n$ no entero, $\theta^n$ no es real si $\theta\lt 0$. Detenemos la integración en el paso en que $\theta$ cambia de signo y refinamos $\xi_1$ con Newton sobre la longitud $s$ del último paso: $s_{k+1} = s_k - \theta(s_k)/\theta'(s_k)$, donde ambos valores salen del propio paso RK4. Como $\theta'(\xi_1)\ne0$ la raíz es simple y la convergencia es cuadrática.

La implementación (resolverLaneEmden en src/lib/math/polytrope.ts) reproduce $\xi_1$ para $n=0,1$ con error menor que $10^{-8}$ y las tablas de Chandrasekhar dentro de su última cifra (los tests de src/lib/math/__tests__/polytrope.test.ts lo comprueban).

Análisis de error. El error global de RK4 con paso adaptativo de tolerancia $\mathrm{tol}$ por paso es $O(\mathrm{tol}\cdot N)$ con $N$ el número de pasos; el test de convergencia para $n = 4.95$ (donde $N$ es mayor) muestra que endurecer la tolerancia de $10^{-11}$ a $10^{-14}$ cambia $\xi_1$ en menos de $10^{-7}$ relativo. La magnitud derivada más sensible es $\rho_c/\bar\rho = -\xi_1/(3\theta'(\xi_1))$: cerca de $n = 5$, $\theta'(\xi_1)\to0$ y un error absoluto $\delta$ en $\theta'$ se amplifica como $\delta/\theta'(\xi_1)^2$. Por eso la tabla de §3 da $\rho_c/\bar\rho$ con menos cifras para $n = 4.5$.

Explora en el simulador: (a) cómo $\xi_1$ y $\rho_c/\bar\rho$ se disparan al acercarse $n$ a 5; (b) que la curva de densidad $\theta^n$ (trazo discontinuo) está mucho más concentrada que $\theta$; (c) con el eje normalizado $\xi/\xi_1$, qué $n$ reproduce mejor la densidad central real del Sol ($\rho_c/\bar\rho = 108$).

Problemas

1. Solución $n=1$. Demuestra que con $u = \xi\theta$ la ecuación (4) con $n=1$ se convierte en $u''+u=0$ y deduce $\theta = \sin\xi/\xi$.

2. Esfera incompresible. Para $n = 0$ prueba que $\rho_c = \bar\rho$ y que $P_c = \dfrac{3}{8\pi}\dfrac{GM^2}{R^4}$.

3. Masa finita para $n=5$. Comprueba que $\theta = (1+\xi^2/3)^{-1/2}$ satisface (4) y que $M = 4\sqrt3\,\pi\alpha^3\rho_c$.

4. Transformación de homología. Prueba que si $\theta(\xi)$ resuelve (4), también lo hace $\tilde\theta(\xi) = A^{2/(n-1)}\theta(A\xi)$ para todo $A\gt0$ ($n\ne1$), aunque con $\tilde\theta(0) = A^{2/(n-1)}$. ¿Por qué entonces basta resolver (4) una vez para cada $n$?

5. Energía gravitatoria (Betti–Ritter). Demuestra que la energía potencial de un politropo es $\Omega = -\dfrac{3}{5-n}\dfrac{GM^2}{R}$. Comprueba el caso $n=0$ ($-\tfrac35 GM^2/R$) y explica por qué diverge en $n = 5$ aunque la masa sea finita. Indicación: usa el teorema del virial $\Omega = -3\int P\,dV$ y $d(P/\rho)= \tfrac{n+1}{n}\,K\rho^{1/n-1}d\rho$ (Chandrasekhar 1939, cap. IV, §9; Kippenhahn et al. 2012, §19.9).

6. Error del arranque (numérico). Usando el simulador o la librería, compara $\xi_1(n=3)$ obtenido (a) arrancando con (8) en $\xi_0 = 10^{-2}$ y (b) arrancando con $\theta=1$, $\theta'=0$ en el mismo $\xi_0$. Demuestra que la diferencia escala como $\xi_0^2$ y estima la constante. Indicación: el estado inicial (b) difiere del verdadero en $\delta\theta\approx\xi_0^2/6$, $\delta\theta'\approx\xi_0/3$; propaga esa perturbación con la ecuación variacional linealizada.

7. Relación masa–radio. A partir de (9), prueba que en la rama no relativista $M R^3 = \text{const}$ y calcula el cociente de radios $R(0.6\,M_\odot)/R(1.0\,M_\odot)$ de dos enanas blancas con la misma $K$. ¿Qué hipótesis de (9) falla al acercarse a $M_{\mathrm{Ch}}$?

8. (Computacional) Convergencia de la raíz. Implementa RK4 de paso fijo $h$ para (4) con $n = 3$, arrancando con (8) en $\xi_0 = 10^{-3}$, y detecta $\xi_1$ (a) por el último punto con $\theta\gt0$ y (b) interpolando linealmente entre los dos puntos que encierran la raíz. Mide numéricamente el orden de convergencia de (a) y de (b) y compáralo con el orden 4 del integrador. Explica por qué (b) no alcanza el orden 2 que cabría esperar de una interpolación lineal.

import numpy as np
n, xi1_ref = 3, 6.896848619
def f(x, y): return np.array([y[1], -max(y[0], 0)**n - 2*y[1]/x])
def lane_emden(h):
    x = 1e-3; y = np.array([1 - x*x/6 + n*x**4/120, -x/3 + n*x**3/30])
    while True:
        k1 = f(x, y); k2 = f(x+h/2, y+h/2*k1); k3 = f(x+h/2, y+h/2*k2); k4 = f(x+h, y+h*k3)
        yn = y + h/6*(k1+2*k2+2*k3+k4)
        if yn[0] <= 0:
            return x, x + h*y[0]/(y[0]-yn[0])     # (a) último punto positivo, (b) interpolación lineal
        x += h; y = yn
for h in [0.1, 0.05, 0.025, 0.0125]:
    a, b = lane_emden(h); print(f'h={h:<7} err_a={abs(a-xi1_ref):.2e}  err_b={abs(b-xi1_ref):.2e}')

Lecturas

• Chandrasekhar (1939), An Introduction to the Study of Stellar Structure, cap. IV (politropos: derivación, soluciones, tablas) y cap. XI (configuraciones degeneradas).
• Kippenhahn, Weigert & Weiss (2012), Stellar Structure and Evolution, cap. 19 (politropos, §19.1–19.9) y cap. 37 (enanas blancas).
• Horedt (2004), Polytropes: Applications in Astrophysics and Related Fields, cap. 2 (tablas de alta precisión de ξ₁ y θ′(ξ₁)).
• Hansen, Kawaler & Trimble (2004), Stellar Interiors, cap. 7 (politropos y métodos numéricos).
• Carroll & Ostlie (2017), §10.5.
• Chandrasekhar (1931), ApJ 74, 81 — el artículo original de dos páginas sobre la masa máxima.
• Bahcall, Pinsonneault & Basu (2001), ApJ 555, 990, DOI 10.1086/321493 — valores centrales del Sol usados en la comparación.