N cuerpos a gran escala y el árbol de Barnes–Hut

1. Motivación: de tres cuerpos a mil millones

En la lección anterior integramos tres cuerpos. Un cúmulo globular tiene $N \sim 10^{5}$–$10^{6}$ estrellas; una simulación cosmológica moderna, $N \gtrsim 10^{10}$ partículas. La aceleración de cada cuerpo,

exige $N - 1$ términos, luego cada evaluación de fuerzas cuesta $N(N-1)/2$ pares: $O(N^2)$. Con $N = 10^6$ son $5\times10^{11}$ interacciones por paso; a $10^{9}$ interacciones por segundo, ocho minutos por paso, y una simulación necesita $10^{4}$–$10^{6}$ pasos. Barnes y Hut (1986) mostraron que agrupando los cuerpos lejanos en «pseudo-partículas» el coste baja a $O(N\log N)$ con un error controlable: es el algoritmo que, con refinamientos (multipolos de orden alto, métodos multipolares rápidos, mallas), está detrás de casi toda la dinámica galáctica y cosmológica numérica.

En esta lección (y en src/lib/math/barneshut.ts) usamos unidades adimensionales con $G = 1$. Para el modelo de Plummer adoptamos además las unidades de Hénon: masa total $M = 1$, energía total $E = -\tfrac14$, con lo que el radio de escala es $a = 3\pi/16 \approx 0.589$.

2. La idea de Barnes y Hut: aproximación monopolar

2.1 Potencial de un grupo lejano

Sea un grupo de cuerpos de masas $m_j$ en posiciones $\mathbf r_j$, con masa total $M = \sum m_j$ y centro de masas $\mathbf R = \sum m_j\mathbf r_j/M$, contenido en una celda de lado $s$. Su potencial en un punto $\mathbf r$ a distancia $d = \lVert\mathbf r - \mathbf R\rVert$ del centro de masas, con $\mathbf x_j = \mathbf r_j - \mathbf R$ y $\mathbf d = \mathbf r - \mathbf R$, es

donde $P_\ell$ son los polinomios de Legendre y $\gamma_j$ el ángulo entre $\mathbf x_j$ y $\mathbf d$ (desarrollo válido para $x_j \lt d$). Los tres primeros términos son:

• $\ell = 0$ (monopolo): $-GM/d$.
• $\ell = 1$ (dipolo): $-\dfrac{G}{d^2}\sum_j m_j\,\mathbf x_j\cdot\hat{\mathbf d} = 0$, idénticamente nulo porque $\sum_j m_j\mathbf x_j = 0$ por definición del centro de masas.
• $\ell = 2$ (cuadrupolo): $-\dfrac{G}{d^3}\sum_j m_j\,x_j^2\,\tfrac12(3\cos^2\gamma_j - 1)$, de tamaño $\lesssim G M s^2/d^3$.

Por tanto, al sustituir el grupo por una sola masa $M$ en $\mathbf R$ (aproximación monopolar), el error relativo del potencial es $O\!\left((s/d)^2\right)$ y el de la fuerza, derivando, también $O\!\left((s/d)^2\right)$. Elegir el centro de masas como centro del desarrollo es lo que mata el dipolo y hace que un parámetro tan simple como $s/d$ controle el error cuadráticamente.

2.2 El árbol de octantes

Se encierra el sistema en un cubo y se subdivide recursivamente en 8 octantes hasta que cada celda contiene a lo sumo un cuerpo. Cada celda interna guarda su masa total y su centro de masas (se acumulan durante la inserción). La construcción cuesta $O(N\log N)$ en promedio y el árbol se reconstruye en cada paso (Barnes & Hut lo hacían así; es más barato que mantenerlo).

2.3 Criterio de apertura

Para calcular la fuerza sobre el cuerpo $i$ se recorre el árbol desde la raíz. En cada celda de lado $s$ cuyo centro de masas dista $d$ del cuerpo:

• si $s/d \lt \theta$ (la celda se ve «pequeña»), se usa su monopolo y no se desciende;
• si no, se abre la celda y se visitan sus hijos;
• una hoja con un solo cuerpo se trata exactamente.

$\theta$ es el parámetro de apertura. Con $\theta = 0$ ninguna celda se acepta y el algoritmo degenera en la suma directa (nuestro test lo comprueba a $10^{-9}$). Valores habituales: $\theta \approx 0.5$–$0.8$. El criterio original de Barnes y Hut usa el lado $s$ y la distancia al centro de masas; Salmon y Warren (1994) mostraron que puede fallar para celdas cuyo centro de masas está cerca del borde, y propusieron añadir a $d$ el desplazamiento del centro de masas respecto del centro geométrico; nuestra implementación usa el criterio original.

2.4 Coste $O(N\log N)$

Proposición. Para un cuerpo dado y un sistema de $N$ cuerpos distribuidos con densidad acotada, el número de celdas aceptadas con el criterio $s/d \lt \theta$ es $O(\theta^{-3}\log N)$.

Demostración (esbozo). Las celdas del nivel $k$ tienen lado $s_k = s_0/2^k$. Una celda de ese nivel se acepta si está a distancia $d \gt s_k/\theta$ y su padre fue abierto, es decir, $d \lt 2 s_k/\theta$. Las celdas aceptadas del nivel $k$ viven por tanto en una cáscara esférica de radios $s_k/\theta$ y $2 s_k/\theta$, de volumen $\tfrac{4\pi}{3}(8 - 1)(s_k/\theta)^3$, que contiene a lo sumo $\approx 29\,\theta^{-3}$ celdas de volumen $s_k^3$. Como hay $\approx\log_2 N^{1/3}\cdot$ const niveles no vacíos (la profundidad del árbol es $O(\log N)$ para densidad acotada), el total es $O(\theta^{-3}\log N)$. Sumando sobre los $N$ cuerpos: $O(N\theta^{-3}\log N)$. $\blacksquare$

La dependencia $\theta^{-3}$ es dura: pasar de $\theta = 1$ a $\theta = 0.5$ multiplica el coste por 8. En nuestra tabla (§5) el cociente medido entre $\theta = 1$ y $\theta = 0.5$ es $530/179 \approx 3$, menor que 8 porque a $N = 2000$ muchas interacciones son con hojas individuales, no con celdas.

3. Suavizado gravitatorio

3.1 Por qué hace falta

Con $N$ finito, dos cuerpos pueden acercarse arbitrariamente y la fuerza $\propto 1/r^2$ diverge: el paso de tiempo tendría que hacerse minúsculo y, peor, una «colisión» entre dos partículas de simulación no representa nada físico cuando cada partícula simboliza $10^{4}$ estrellas o una porción de materia oscura. Se sustituye el potencial de cada partícula por el de una esfera de Plummer de radio $\varepsilon$:

Es exactamente el potencial generado por la densidad de Plummer (1911), $\rho(r) = \dfrac{3m}{4\pi\varepsilon^3}\left(1 + r^2/\varepsilon^2\right)^{-5/2}$ (ejercicio 3). La fuerza es finita en todas partes: su máximo, $\lvert a\rvert_{\max} = \dfrac{2Gm}{3\sqrt3\,\varepsilon^2}$, ocurre en $r = \varepsilon/\sqrt2$ (ejercicio 4). Ese mismo $\varepsilon$ es el parámetro eps de accelerations en integrators.ts y de todo el módulo barneshut.ts.

3.2 Sesgo y ruido

El precio es un sesgo: para $r \gg \varepsilon$,

la fuerza queda subestimada en un $1.5\,(\varepsilon/r)^2$ ($1.5$ % a $r = 10\varepsilon$, $15$ % a $r = 3\varepsilon$: sesgoPlummer). Dehnen (2001) analiza el compromiso: $\varepsilon$ grande reduce el ruido de la fuerza (las fluctuaciones por el número finito de vecinos, que escalan como $N^{-1/2}$ para las fuerzas de corto alcance) pero aumenta el sesgo; existe un $\varepsilon$ óptimo que decrece con $N$ aproximadamente como $N^{-0.2}$–$N^{-0.3}$ para el perfil de Plummer, y los núcleos de soporte compacto (splines) son preferibles al de Plummer porque su sesgo se anula exactamente para $r \gt 2\varepsilon$. Una regla práctica: $\varepsilon$ del orden de la separación media entre vecinos en la región más densa; en nuestra esfera de Plummer con $N = 1000$ esa distancia media al vecino más próximo es $0.22$ (unidades de Hénon).

4. El teorema del virial

4.1 Enunciado y demostración

Sea $I = \sum_i m_i\lVert\mathbf r_i\rVert^2$ el momento de inercia escalar. Derivando dos veces,

donde $K$ es la energía cinética y $\mathcal V = \sum_i\mathbf r_i\cdot\mathbf F_i$ es el virial de Clausius. Para fuerzas por pares $\mathbf F_{ij} = -\mathbf F_{ji}$ dirigidas según $\mathbf r_{ij} = \mathbf r_i - \mathbf r_j$, agrupando pares,

Si $\mathbf F_{ij}$ deriva de un potencial homogéneo de grado $-1$ ($\Phi = -G m_i m_j/r_{ij}$), $\mathbf r_{ij}\cdot\mathbf F_{ij} = -r\,\partial_r\Phi = \Phi$, y $\mathcal V = W$, la energía potencial total. Promediando en el tiempo un sistema acotado ($\langle\ddot I\rangle = 0$):

$\blacksquare$ (Binney & Tremaine §4.8.)

4.2 El virial con suavizado

Con Plummer la fuerza no es homogénea de grado $-2$, y el virial de la fuerza difiere de la energía potencial suavizada:

Así que la relación correcta en equilibrio es $2K + \mathcal V_\varepsilon = 0$, no $2K + W_\varepsilon = 0$; la diferencia es positiva y de orden $(\varepsilon/r)^2\lvert W\rvert$. virial() devuelve las tres cantidades. En nuestra esfera de Plummer ($N = 1000$): con $\varepsilon = 0$, $K = 0.262$, $W = -0.521$, $2K/\lvert W\rvert = 1.006$ (la desviación del 0.6 % frente a 1 es ruido de muestreo; los valores teóricos son $K = \tfrac14$, $W = -\tfrac12$). Con $\varepsilon = 0.1$, $W_\varepsilon = -0.506$ y $\mathcal V_\varepsilon = -0.484$: la condición ingenua $2K/\lvert W_\varepsilon\rvert$ da $1.035$ y la correcta $2K/\lvert\mathcal V_\varepsilon\rvert$ da $1.08$: el sistema, generado en equilibrio newtoniano, está un 8 % «demasiado caliente» para la gravedad suavizada y se expandirá ligeramente al integrarlo. Con $\varepsilon = 0.3$ ($\approx a/2$) la discrepancia sube al 40 %: un suavizado así ya no simula el mismo sistema.

4.3 Uso como diagnóstico

En una simulación se monitoriza $Q = 2K/\lvert\mathcal V\rvert$ (además de $E$, $\mathbf P$ y $\mathbf L$):

• $Q$ oscila alrededor de 1 con periodo $\sim$ tiempo dinámico: sistema en equilibrio.
• $Q \gt 1$ sostenido: sistema caliente, se expande; $Q \lt 1$: colapsa. Un sistema en colapso frío ($Q \approx 0$) «relaja violentamente» hacia $Q \approx 1$ en pocos tiempos dinámicos (Lynden-Bell 1967).
• Deriva secular de $Q$ con $E$ constante: relajación por dos cuerpos (evaporación, colapso del núcleo), un efecto físico en cúmulos pero espurio en simulaciones de sistemas no colisionales, donde señala $\varepsilon$ demasiado pequeño o $N$ demasiado pequeño.

5. Resultados numéricos

Todos de barneshut.ts sobre esferas de Plummer generadas con modeloPlummer (semilla 7), $\varepsilon = 0.01$; el error relativo de cada cuerpo es $\lVert\mathbf a_{\text{BH}} - \mathbf a_{\text{directa}}\rVert/\lVert\mathbf a_{\text{directa}}\rVert$.

5.1 Coste y error en función de $N$ ($\theta = 0.5$)

Al multiplicar $N$ por 32 (250 → 8000) las interacciones por cuerpo se multiplican por 6.7, mientras que la suma directa se multiplica por 32. El crecimiento es más rápido que $\log N$ puro (que daría $\times 1.6$) porque a estos $N$ el perfil de Plummer, muy concentrado, hace que una fracción grande de interacciones sea con hojas; el régimen asintótico $O(\log N)$ solo se alcanza a $N$ mucho mayores, cuando la profundidad del árbol domina sobre las interacciones con hojas. El error máximo crece con $N$: siempre hay algún cuerpo para el que una celda queda justo en el límite del criterio; es el caso que motivó la corrección de Salmon y Warren.

5.2 Error en función de $\theta$ ($N = 2000$)

El error mediano escala aproximadamente como $\theta^{2.5}$ entre $0.2$ y $1$ (la predicción monopolar es $\theta^{2}$ por celda, pero el número de celdas aceptadas también crece con $\theta$ y los errores se suman con signos aleatorios). Con $\theta = 1.5$ hay cuerpos cuya fuerza es completamente errónea (error 150 %): el criterio ha aceptado celdas que contienen al propio vecino más cercano.

5.3 Qué precisión hace falta

Un error relativo de fuerza del $10^{-3}$ por partícula y por paso, con signos aleatorios, se acumula tras $n$ pasos como $\sqrt n\cdot 10^{-3}$: en $10^{4}$ pasos, un 10 %. Pero el ruido de discretización (usar $N$ partículas para representar un continuo) produce ya errores de fuerza del orden de $N^{-1/2}$ (Dehnen 2001): para $N = 10^{6}$, $10^{-3}$. No tiene sentido calcular la fuerza mucho mejor que ese ruido intrínseco: de ahí que $\theta \approx 0.5$–$0.7$ sea el estándar.

Qué observar. El simulador usa suma directa (tres cuerpos no justifican un árbol). Con el escenario «Pitagórico», activa el suavizado con $\varepsilon = 0.1$ y compara la evolución con la de $\varepsilon = 0$: los encuentros cercanos que expulsan a la masa 3 se amortiguan y el estado final cambia. Es la demostración más directa de que $\varepsilon$ es un parámetro físico del modelo, no un detalle numérico.

Problemas

1. (Cuadrupolo). Demuestra que para un grupo de masas con centro de masas en el origen el término cuadrupolar del potencial en $\mathbf d$ es $-\dfrac{G}{2d^3}\sum_j m_j\left(3(\mathbf x_j\cdot\hat{\mathbf d})^2 - x_j^2\right)$, y que para una celda cúbica de lado $s$ centrada en el centro de masas está acotado en valor absoluto por $\tfrac34\,G M s^2/d^3$. Concluye que el error relativo del potencial monopolar es $\le \tfrac34(s/d)^2$ y, con $\theta = 0.5$, inferior al 19 % en el peor caso.

2. (Dipolo nulo). Demuestra que si el desarrollo se hace alrededor de un punto $\mathbf C$ distinto del centro de masas (por ejemplo, el centro geométrico de la celda), aparece un término dipolar de error relativo $O(s/d)$, y calcula cuánto peor es con $\theta = 0.5$.

3. (Densidad de Plummer). Demuestra, con la ecuación de Poisson $\nabla^2\Phi = 4\pi G\rho$, que el potencial $\Phi_\varepsilon = -Gm/\sqrt{r^2+\varepsilon^2}$ corresponde a $\rho(r) = \dfrac{3m}{4\pi\varepsilon^3}\left(1 + r^2/\varepsilon^2\right)^{-5/2}$, y comprueba que $\int\rho\,dV = m$.

4. (Fuerza máxima). Demuestra que $\lvert\mathbf a_\varepsilon\rvert = Gm\,r(r^2+\varepsilon^2)^{-3/2}$ es máxima en $r = \varepsilon/\sqrt2$ con valor $\dfrac{2Gm}{3\sqrt3\,\varepsilon^2}$, y deduce una cota inferior del paso de tiempo de un leapfrog estable, $h \lesssim \varepsilon^{3/2}/\sqrt{Gm}$.

5. (Virial con suavizado). Demuestra la identidad $\mathcal V_\varepsilon = W_\varepsilon + \sum_{i<j}G m_i m_j\varepsilon^2(r_{ij}^2+\varepsilon^2)^{-3/2}$ y concluye que $\lvert\mathcal V_\varepsilon\rvert \lt \lvert W_\varepsilon\rvert$: un sistema en equilibrio newtoniano siempre parece «caliente» para la gravedad suavizada.

6. (Computacional: árbol en 2D). Implementa un árbol de cuadrantes para $N$ cuerpos en el plano, con criterio $s/d \lt \theta$, y mide (a) el número de interacciones por cuerpo para $N = 10^2, 10^3, 10^4$ con $\theta = 0.5$; (b) el error mediano frente a la suma directa para $\theta \in \{0.3, 0.5, 0.7, 1\}$. Compara con las tablas de §5.

import numpy as np
class Celda:
    def __init__(self, cx, cy, s):
        self.cx, self.cy, self.s = cx, cy, s
        self.m = 0.0; self.X = 0.0; self.Y = 0.0   # masa y centro de masas acumulados
        self.cuerpo = None; self.hijos = None
    def insertar(self, i, x, y, m):
        if self.m == 0 and self.hijos is None:          # hoja vacía
            self.cuerpo, self.m, self.X, self.Y = i, m, x, y; return
        if self.hijos is None:                          # hoja ocupada: subdividir
            self.hijos = [None]*4; j = self.cuerpo; self.cuerpo = None
            mj, xj, yj = self.m, self.X, self.Y; self.m = 0
            self._bajar(j, xj, yj, mj)
        self._bajar(i, x, y, m)
        self.X = (self.X*self.m + x*m)/(self.m+m); self.Y = (self.Y*self.m + y*m)/(self.m+m); self.m += m
    def _bajar(self, i, x, y, m):
        k = (x >= self.cx) + 2*(y >= self.cy); q = self.s/4
        if self.hijos[k] is None:
            self.hijos[k] = Celda(self.cx + (q if k & 1 else -q), self.cy + (q if k & 2 else -q), self.s/2)
        self.hijos[k].insertar(i, x, y, m)
def fuerza(c, i, x, y, theta, eps, cont):
    if c.m == 0 or c.cuerpo == i: return np.zeros(2)
    dx, dy = c.X - x, c.Y - y; d = np.hypot(dx, dy)
    if c.hijos is None or c.s < theta*d:
        cont[0] += 1; return c.m*np.array([dx, dy])/(d*d + eps*eps)**1.5
    return sum(fuerza(h, i, x, y, theta, eps, cont) for h in c.hijos if h is not None)

7. (Computacional: relajación violenta). Genera una esfera de Plummer con $N = 1000$ y reduce todas las velocidades a la mitad ($Q = 0.25$). Integra con leapfrog ($h = 0.01$, $\varepsilon = 0.05$, $\theta = 0.5$) durante $t = 20$ y grafica $Q(t) = 2K/\lvert\mathcal V_\varepsilon\rvert$, $E(t)$ y el radio que contiene la mitad de la masa. ¿En cuántos tiempos dinámicos se alcanza $Q \approx 1$?

Lecturas

• Barnes & Hut (1986), Nature 324, 446: el artículo original, de cuatro páginas, legible por completo.
• Aarseth (2003), Gravitational N-Body Simulations, cap. 6: códigos de árbol en el contexto de los métodos directos; cap. 1 para el planteamiento.
• Dehnen (2001), MNRAS 324, 273, y Dehnen & Read (2011), EPJ Plus 126, 55: la revisión moderna de suavizado, árboles, multipolos rápidos y pasos individuales.
• Binney & Tremaine (2008), Galactic Dynamics, §2.9 (resolución numérica de la ecuación de Poisson: árboles, mallas) y §4.8 (ecuaciones de Jeans y virial); §2.2.2 para la esfera de Plummer.
• Aarseth, Hénon & Wielen (1974), A&A 37, 183: el algoritmo de muestreo del modelo de Plummer que usa modeloPlummer.
• Salmon, J. K. & Warren, M. S. (1994). Skeletons from the treecode closet. J. Comput. Phys. 111, 136–155: los fallos del criterio original y los criterios con cota de error garantizada.