Las ecuaciones de estructura estelar

1. Motivación: una estrella como problema de contorno

Del Sol vemos la fotosfera: su luminosidad $L_\odot = 3.828\times10^{26}$ W, su radio $R_\odot = 6.957\times10^8$ m y, por la dinámica planetaria, su masa $M_\odot = 1.989\times10^{30}$ kg (valores nominales IAU 2015). Nada de eso nos dice directamente qué temperatura hay en el centro, qué densidad, ni cómo se transporta la energía. Sin embargo, desde Eddington (1926) sabemos que cuatro ecuaciones diferenciales ordinarias, con una ecuación de estado, una opacidad y una tasa de generación de energía, bastan para construir un modelo del interior completo. El modelo solar estándar resultante predice flujos de neutrinos y velocidades del sonido que la heliosismología y los detectores de neutrinos han verificado con precisión del orden del 1 %.

Esta lección deriva las cuatro ecuaciones, las condiciones de contorno y, antes de resolver nada, muestra lo lejos que llega el análisis dimensional: la temperatura central del Sol se estima con error del 12 % sin integrar una sola ecuación.

2. Hipótesis y variables

Supondremos una estrella esférica (sin rotación ni campo magnético apreciables), aislada y en estado cuasiestático: las magnitudes cambian en escalas de tiempo muy largas comparadas con la escala dinámica (§7). Las variables son funciones de la distancia al centro $r$:

3. Las cuatro ecuaciones

3.1 Continuidad de masa

La masa de una capa esférica de espesor $dr$ es $dm = 4\pi r^2\rho\,dr$:

3.2 Equilibrio hidrostático

Consideremos un cilindro de fluido de base $dA$ y altura $dr$ a distancia $r$ del centro. Sobre él actúan la presión en la base, $P(r)\,dA$, hacia fuera; la presión en la tapa, $P(r+dr)\,dA$, hacia dentro; y la gravedad. Por el teorema de Newton para distribuciones esféricas, solo la masa interior $m(r)$ contribuye, con aceleración $g = Gm/r^2$. La segunda ley de Newton para el elemento de masa $\rho\,dA\,dr$ es

En equilibrio ($\ddot r = 0$) obtenemos la ecuación de equilibrio hidrostático:

La presión decrece hacia fuera; el gradiente de presión sostiene cada capa contra el peso de la materia que tiene encima. Si (2) deja de cumplirse, la estrella responde en la escala de tiempo dinámica (§7), que para el Sol es de media hora: por eso una estrella «normal» está siempre en equilibrio hidrostático con enorme precisión.

3.3 Conservación de la energía

Sea $\varepsilon$ la energía liberada por reacciones nucleares por unidad de masa y tiempo, y $\varepsilon_\nu$ la que se llevan los neutrinos (que escapan sin interactuar). La capa $dm$ recibe $L(r)$ por su cara interna y emite $L(r+dr)$ por la externa; si además su entropía específica $s$ cambia, absorbe o cede calor $T\,ds/dt$ por unidad de masa. El primer principio da

En equilibrio térmico ($\partial s/\partial t = 0$) la luminosidad crece exactamente con la energía nuclear generada: $dL/dr = 4\pi r^2\rho\,\varepsilon$. El término de entropía es esencial en fases de contracción o expansión (presecuencia principal, gigantes) y se llama a veces «luminosidad gravitatoria».

3.4 Transporte de energía

La energía fluye de dentro hacia fuera porque hay un gradiente de temperatura. Hay tres mecanismos:

• Radiación. Los fotones difunden a través de un medio opaco. En la lección 2 se deriva que el flujo es $F = -\dfrac{4acT^3}{3\kappa\rho}\dfrac{dT}{dr}$, con $a = 4\sigma/c = 7.5657\times10^{-16}$ J m⁻³ K⁻⁴ la constante de radiación y $\sigma$ la de Stefan–Boltzmann. Como $L = 4\pi r^2 F$, $$\frac{dT}{dr} = -\frac{3\,\kappa\rho}{16\pi a c}\,\frac{L}{r^2T^3}. \tag{4a}$$
• Convección. Cuando el gradiente radiativo necesario es demasiado empinado, el gas se vuelve inestable (criterio de Schwarzschild, lección 2) y burbujas macroscópicas transportan el calor. En el interior profundo el gradiente real queda muy cerca del adiabático: $$\frac{dT}{dr} = \left(1-\frac1\gamma\right)\frac{T}{P}\frac{dP}{dr}, \tag{4b}$$ con $\gamma$ el exponente adiabático ($5/3$ para gas monoatómico ideal).
• Conducción. Despreciable en estrellas normales; domina en materia degenerada (enanas blancas, lección 4), donde se incorpora a (4a) como una opacidad conductiva efectiva.

Es cómodo escribir ambos casos con el gradiente logarítmico $\nabla \equiv d\ln T/d\ln P$:

La expresión de $\nabla_{\rm rad}$ se obtiene dividiendo (4a) por (2) (ejercicio 1).

3.5 Relaciones de cierre

Las ecuaciones (1)–(4) contienen cinco funciones incógnita ($m, P, T, L, \rho$) y tres «funciones materiales» que dependen de la composición química $X_i$ (fracciones en masa de hidrógeno $X$, helio $Y$ y metales $Z$):

Para gas ideal más radiación,

donde $k_B$ es la constante de Boltzmann, $m_u$ la unidad de masa atómica y $\mu$ el peso molecular medio (masa media por partícula en unidades de $m_u$). Para la mezcla solar superficial $X\approx0.71$, $Y\approx0.27$, $Z\approx0.02$, $\mu\approx0.61$.

4. Condiciones de contorno

El sistema (1)–(4) es de cuarto orden, así que hacen falta cuatro condiciones. No están todas en el mismo punto: es un problema de contorno en dos puntos.

En el centro ($r = 0$):

Ambas son consecuencia de la regularidad: la masa y la luminosidad encerradas en una esfera de radio nulo son nulas. De (1) y (3) se sigue que cerca del centro $m \approx \tfrac{4\pi}{3}\rho_c r^3$ y $L\approx\tfrac{4\pi}{3}\rho_c\varepsilon_c r^3$, y de (2) que $dP/dr \to 0$: la presión y la temperatura tienen un máximo plano en el centro, $P \approx P_c - \tfrac{2\pi}{3}G\rho_c^2 r^2$.

En la superficie ($r = R$): la aproximación más sencilla es la «superficie cero»,

que es excelente para el interior porque la presión y la temperatura superficiales son despreciables frente a las centrales ($P_{\rm fot}/P_c \sim 10^{-12}$). Un modelo más fino impone condiciones fotosféricas: en la capa de profundidad óptica $\tau = 2/3$ la temperatura es la efectiva, $T = T_{\rm ef}$ con $L = 4\pi R^2\sigma T_{\rm ef}^4$, y la presión se obtiene integrando (2) con $\tau = \int\kappa\rho\,dr$ hacia abajo desde $\tau = 0$: $P_{\rm fot}\approx \dfrac{2}{3}\dfrac{g}{\bar\kappa}$ con $g = GM/R^2$.

Teorema de Vogt–Russell. Fijadas la masa $M$ y la composición $X_i(m)$, el sistema (1)–(4) con (6)–(7) determina (en general de manera única) $R$, $L$ y toda la estructura interna. No es un teorema matemático riguroso (existen contraejemplos con soluciones múltiples), pero explica por qué la secuencia principal es una familia uniparamétrica: cada masa fija un punto del diagrama H–R.

5. Forma lagrangiana: la masa como variable independiente

Como $m(r)$ es monótona creciente, se puede usar $m$ en vez de $r$ como variable. Dividiendo (2)–(4) por (1):

Tres ventajas: (i) el dominio $0\le m\le M$ está fijado de antemano (el radio $R$ es una incógnita del problema, el dominio en $r$ no se conoce hasta resolverlo); (ii) la composición química se transporta con la masa, no con el radio, así que $X_i(m)$ es la descripción natural en evolución estelar; (iii) la ecuación de la presión ya no contiene $\rho$, lo que suaviza su comportamiento.

6. Estimaciones por análisis dimensional frente al modelo solar estándar

Sin resolver nada, (2) dice que $P$ cae de $P_c$ a $0$ en una distancia $R$, con $m\sim M$, $\rho\sim\bar\rho$:

Y de (5), despreciando la radiación, $T_c \sim \dfrac{\mu m_u}{k_B}\dfrac{P_c}{\rho_c}\sim \dfrac{\mu m_u}{k_B}\dfrac{GM}{R}$. Para el Sol, con los valores nominales IAU 2015 y $G = 6.6743\times10^{-11}$ m³ kg⁻¹ s⁻² (CODATA 2018):

Comparemos con el modelo solar estándar BP2000 (Bahcall, Pinsonneault & Basu 2001, tablas 5–6), cuyas predicciones heliosísmicas y de neutrinos han sido contrastadas con la observación:

La temperatura sale sorprendentemente bien (12 % de error) mientras que la presión y la densidad fallan por uno o dos órdenes de magnitud. La razón es instructiva: $T_c\propto P_c/\rho_c$, y la concentración central ($\rho_c/\bar\rho\approx108$) afecta por igual al numerador y al denominador. La temperatura central de una estrella queda fijada, en orden de magnitud, por $GM/R$ y nada más: es un teorema del virial disfrazado (§7).

Un detalle de consistencia: con la composición central del modelo ($X_c = 0.340$, $Y_c = 0.6405$, $Z_c = 0.0198$) la ecuación (5) da $\mu_c = 0.854$ y la temperatura de gas ideal $T = \mu_c m_u P_c/(k_B\rho_c) = 1.576\times10^7$ K, que coincide con el valor tabulado al 0.4 %. La presión de radiación central es $aT_c^4/3 = 1.5\times10^{13}$ Pa, apenas el 0.07 % del total: en el Sol la radiación no cuenta para el equilibrio, aunque sí para el transporte.

El politropo como puente. Si se cierra el sistema con $P = K\rho^{1+1/n}$ (lección 1), la ecuación de Lane–Emden da para $n = 3$: $\rho_c = 7.64\times10^4$ kg m⁻³ y $P_c = 1.25\times10^{16}$ Pa, ya dentro de un factor 2 del modelo completo. Explóralo en el simulador: ¿qué índice $n$ reproduce mejor la concentración central $\rho_c/\bar\rho = 108$ del Sol?

7. Dos resultados rigurosos: cota de la presión central y teorema del virial

Cota inferior de $P_c$. De (2) y (1), $\dfrac{dP}{dm} = -\dfrac{Gm}{4\pi r^4}$. Integrando desde el centro a la superficie y usando $r(m)\le R$:

Para el Sol, $P_c \ge 4.5\times10^{13}$ Pa. Es una cota débil (el valor real es 500 veces mayor) pero rigurosa: no supone nada sobre la ecuación de estado ni sobre el transporte. Con la hipótesis adicional de que $\bar\rho(r)$, la densidad media interior a $r$, no crece hacia fuera, se obtiene una cota mejor (ejercicio 3).

Teorema del virial. Multiplicando (2) por $4\pi r^3$ e integrando por partes entre $0$ y $R$ con $P(R) = 0$:

donde $\Omega\lt0$ es la energía potencial gravitatoria. Luego

Para gas ideal monoatómico $P/\rho = \tfrac23 u$, con $u$ la energía interna por unidad de masa, así que $2U = -\Omega$ y la energía total es $E = U + \Omega = -U = \Omega/2$. Dos consecuencias: (i) una estrella que pierde energía se calienta (capacidad calorífica negativa); (ii) la temperatura media satisface $\dfrac{3k_B}{\mu m_u}\langle T\rangle M = -\Omega \sim \dfrac{GM^2}{R}$, de donde $\langle T\rangle \sim \dfrac{\mu m_u GM}{3k_B R}\approx 4.6\times10^6$ K para el Sol. Esta es la justificación de la estimación dimensional de §6.

Escalas de tiempo (todas del Sol, con $M_\odot$, $R_\odot$, $L_\odot$ nominales):

• Dinámica, de (2) sin presión: $t_{\rm din}\sim\sqrt{R^3/GM} = 1.6\times10^3$ s ≈ 27 min.
• Térmica o de Kelvin–Helmholtz, de (10): $t_{\rm KH}\sim\dfrac{GM^2}{RL} = 9.9\times10^{14}$ s ≈ $3.1\times10^7$ años. Es la edad máxima que el Sol tendría si brillara solo por contracción: demasiado corta frente a la edad geológica, lo que obligó a buscar la fuente nuclear.
• Nuclear, suponiendo que el 10 % de la masa quema hidrógeno con eficiencia $0.007\,c^2$ (lección 3): $t_{\rm nuc}\sim\dfrac{0.1\times0.007\,Mc^2}{L} = 3.3\times10^{17}$ s ≈ $1.0\times10^{10}$ años.

La jerarquía $t_{\rm din}\ll t_{\rm KH}\ll t_{\rm nuc}$ es la que legitima las hipótesis de §2: en la escala nuclear, la estrella está en equilibrio hidrostático y (casi siempre) térmico.

8. Análisis numérico: por qué el sistema es difícil

1. Dos puntos de contorno. Un disparo (shooting) desde el centro con $P_c$, $T_c$ de prueba y desde la superficie con $R$, $L$ de prueba, ajustando en un punto de encuentro, es inestable: cerca de la superficie $P\to0$ y $T\to0$, y las soluciones vecinas divergen exponencialmente. Por eso los códigos modernos usan el método de Henyey (relajación de Newton–Raphson sobre toda la malla a la vez), que trata las cuatro condiciones simultáneamente.
2. Rigidez. $\varepsilon\propto T^{\nu}$ con $\nu\approx4$ (pp) o $\approx20$ (CNO), y $\kappa\propto\rho T^{-3.5}$: las escalas de variación de $L$ y $T$ difieren en órdenes de magnitud. Métodos explícitos (RK4) necesitan pasos diminutos; los implícitos (Henyey es implícito por construcción) no.
3. Singularidades en los extremos. En $r = 0$ hay términos $m/r^2$ de tipo $0/0$; se arranca con las series de §4. En la superficie $\rho\to0$ y $1/\rho$ diverge: se integra en $m$ (ecuaciones (8)) o se usa la atmósfera gris como condición de contorno.
4. Propagación de errores. Una incertidumbre relativa $\delta\kappa/\kappa$ en la opacidad cambia $\nabla_{\rm rad}$ en la misma proporción (ecuación 4): Bahcall, Serenelli & Basu (2005) mostraron que un aumento del $\sim10$ % de la opacidad entre $2\times10^6$ y $4.5\times10^6$ K bastaría para reconciliar los modelos de baja metalicidad con la heliosismología. La estructura es más sensible a la opacidad que a cualquier otro ingrediente microfísico.

Problemas

1. Gradiente radiativo logarítmico. Deduce $\nabla_{\rm rad} = \dfrac{3\kappa LP}{16\pi acGmT^4}$ a partir de (4a) y (2).

2. Comportamiento central. Demuestra que cerca del centro $P(r) = P_c - \dfrac{2\pi}{3}G\rho_c^2r^2 + O(r^4)$ y que $T(r) = T_c - \dfrac{2\pi}{3}G\rho_c^2\nabla_c\dfrac{T_c}{P_c}r^2 + O(r^4)$, con $\nabla_c$ el gradiente logarítmico en el centro.

3. Cota de Milne. Suponiendo que la densidad media interior $\bar\rho(r) = 3m/(4\pi r^3)$ no crece con $r$ (la estrella es más densa hacia el centro), demuestra que $P_c \ge \dfrac{3}{8\pi}\dfrac{GM^2}{R^4}$, el triple de (9), y que la igualdad se da para densidad uniforme.

4. Virial con presión de radiación. Sea $\beta = P_{\rm gas}/P$ constante en toda la estrella (modelo estándar de Eddington). Demuestra que $2U_{\rm gas} + U_{\rm rad} = -\Omega$, calcula la energía total $E$ en función de $\Omega$ y $\beta$, y discute el límite $\beta\to0$.

5. Temperatura media del Sol. Con el teorema del virial y $\Omega = -\tfrac{3}{5-n}\tfrac{GM^2}{R}$ para un politropo (lección 1, ejercicio 5), calcula $\langle T\rangle$ del Sol para $n = 3$ con $\mu = 0.61$, y compara con $T_c$ del BP2000.

6. (Computacional) El Sol como politropo n = 3 frente al modelo BP2000. Integra las ecuaciones (1)–(2) cerradas con $P = K\rho^{4/3}$ (sin usar Lane–Emden directamente) por disparo sobre $\rho_c$: fija $K$ de modo que $M = M_\odot$ y $R = R_\odot$ y compara $\rho_c$, $P_c$ y $T_c$ (con $\mu = 0.61$) con la tabla de §6. Estima el error de truncamiento del integrador duplicando el paso.

import numpy as np
G, Msun, Rsun, kB, mu_u = 6.6743e-11, 1.989e30, 6.957e8, 1.380649e-23, 1.66053907e-27

def integrar(rho_c, K, h=2e5):
    # y = [m, P];  rho = (P/K)^(3/4);  RK4 desde r0 con las series del §4
    r = h; P_c = K*rho_c**(4/3)
    y = np.array([4*np.pi/3*rho_c*r**3, P_c - 2*np.pi/3*G*rho_c**2*r**2])
    def f(r, y):
        rho = (max(y[1], 0)/K)**0.75
        return np.array([4*np.pi*r*r*rho, -G*y[0]*rho/r**2])
    while y[1] > 0:
        k1 = f(r, y); k2 = f(r+h/2, y+h/2*k1); k3 = f(r+h/2, y+h/2*k2); k4 = f(r+h, y+h*k3)
        y = y + h/6*(k1+2*k2+2*k3+k4); r += h
    return r, y[0]          # radio y masa (la última capa se sobrepasa en < h)

# Para n = 3, M no depende de rho_c: M(K) = 4 pi omega_3 (K/pi G)^1.5. Fijamos K por M = Msun...
K = np.pi*G*(Msun/(4*np.pi*2.01824))**(2/3)
# ... y rho_c por R = Rsun (biseccion: R decrece con rho_c)
lo, hi = 1e3, 1e7
for _ in range(40):
    rc = np.sqrt(lo*hi); R, M = integrar(rc, K)
    lo, hi = (rc, hi) if R > Rsun else (lo, rc)
print(rc, K*rc**(4/3), K*rc**(1/3)*0.61*mu_u/kB)   # rho_c, P_c, T_c

7. Fotosfera. A partir de (2) y $d\tau = -\kappa\rho\,dr$, demuestra que si $\kappa$ es constante en la atmósfera $P(\tau) = g\tau/\kappa$, y justifica la condición $P_{\rm fot} = 2g/(3\kappa)$ en $\tau = 2/3$. Evalúa $P_{\rm fot}$ para el Sol con $\kappa = 0.1$ m² kg⁻¹ y compárala con $P_c$.

Lecturas

• Kippenhahn, Weigert & Weiss (2012), Stellar Structure and Evolution, cap. 1–2 (ecuaciones mecánicas y virial), 4–5 (energía y transporte), 9–10 (condiciones de contorno y teorema de Vogt–Russell), 11 (método de Henyey).
• Carroll & Ostlie (2017), An Introduction to Modern Astrophysics, §10.1–10.5.
• Hansen, Kawaler & Trimble (2004), Stellar Interiors, cap. 1 (análisis dimensional y escalas de tiempo) y cap. 7 (métodos numéricos).
• Eddington (1926), The Internal Constitution of the Stars, cap. IV–VI: el origen del modelo estándar.
• Bahcall, Pinsonneault & Basu (2001), ApJ 555, 990, DOI 10.1086/321493 — modelo BP2000, tablas 5–6 (valores centrales usados aquí).
• Bahcall, Serenelli & Basu (2005), ApJ 621, L85, DOI 10.1086/428929; (2006), ApJS 165, 400, DOI 10.1086/504043 — sensibilidad del modelo a opacidades y abundancias.