Accurately simulating core-collapse self-interacting dark matter halos

This paper investigates the numerical challenges and convergence properties of self-interacting dark matter (SIDM) N-body simulations during the gravothermal collapse phase, demonstrating that high-resolution simulations with strict energy conservation can accurately model the formation of compact halos capable of explaining observations like the GD-1 stellar stream.

The Gist

The Big Picture: The "Ghostly" Crowd

Imagine the universe is filled with invisible, ghostly particles called Dark Matter. In the standard story (Cold Dark Matter), these ghosts never bump into each other; they just float past one another like ghosts in a haunted house, only feeling gravity.

But what if these ghosts do bump into each other? This is the theory of Self-Interacting Dark Matter (SIDM). If they collide, they can swap energy, kind of like a crowded dance floor where people bumping into each other eventually spread out or bunch up.

This paper is about a specific, dramatic event in the life of a dark matter "cloud" (a halo): The Gravothermal Collapse.

The Story of the Halo: From a Crowd to a Black Hole

Think of a dark matter halo as a giant, fluffy cloud of gas.

1. The Dance Floor: Initially, the particles are spread out. They bump into each other, and heat (energy) flows from the hot center to the cooler edges.
2. The Core Expansion: At first, this makes the center puff up and get less dense.
3. The Reversal (The Collapse): Eventually, the center gets so crowded that the flow reverses. Now, the center starts sucking energy out of itself to keep the outer edges warm. The center gets hotter, moves faster, and shrinks tighter and tighter. It's like a crowd of people in a room suddenly deciding to huddle so tightly in the center that the room becomes a vacuum around them.
4. The Result: The center becomes incredibly dense—potentially dense enough to explain strange observations in the sky, like gaps in star streams (GD-1) or weird gravitational lensing effects.

The Problem: The "Video Game" Glitch

The authors are trying to simulate this collapse on a supercomputer. But simulating physics is like playing a video game with a physics engine. If the settings aren't perfect, the game glitches.

The paper is essentially a "How-To Guide for Not Breaking the Simulation." They found that if you aren't careful, your computer simulation will lie to you.

Here are the main "glitches" they fixed, explained with analogies:

1. The "Pixel Size" Problem (Softening Length)

In a simulation, you can't have particles touching perfectly; you need a tiny "cushion" around them so they don't crash the computer. This is called the softening length.

• The Glitch: If this cushion is too small, the particles get too close, and the computer's math gets shaky. It's like trying to calculate the distance between two people standing on the same square inch of floor; the math gets messy.
• The Fix: The authors found that if the cushion is too small, the simulation loses energy (like a leaky battery), causing the collapse to happen too fast. They learned you need a "Goldilocks" cushion: not too big, not too small.

2. The "Stopwatch" Problem (Time Steps)

The simulation moves forward in tiny slices of time.

• The Glitch: As the halo collapses, things happen faster and faster. The computer needs to take smaller and smaller "steps" in time to keep up. If you force the computer to take big steps (a "minimum time step") just to save time, it starts making mistakes.
• The Analogy: Imagine trying to film a hummingbird's wings. If you take a photo every second, you just see a blur. If you take a photo every millisecond, you see the wings clearly. If you force the camera to take photos only once a second, you miss the action entirely.
• The Fix: They found that forcing a "minimum time step" makes the simulation think the center is denser than it really is. It's a lie. To get the truth, you have to let the computer slow down and take tiny steps, even if it takes a long time to run.

3. The "Kernel" Problem (How Many Neighbors?)

To calculate how particles interact, the computer looks at a particle's neighbors.

• The Glitch: If the computer looks at too few neighbors, it misses the big picture. If it looks at too many, it blurs the details.
• The Fix: They found that the "neighbor count" needs to be high enough so that the particles aren't "smearing" each other out. If the simulation area is too fuzzy, the heat moves too fast, and the collapse happens artificially quickly.

The Satellite vs. The Isolated Halo

The authors ran two types of experiments:

1. The Hermit: A dark matter halo floating alone in space.
2. The Commuter: A dark matter halo orbiting a giant galaxy (like our Milky Way).

The Surprise: The "Commuter" halo collapses much faster.

• The Analogy: Imagine a person walking alone in a park (Hermit). They walk at their own pace. Now imagine that same person walking through a crowded, shaking subway station (Commuter). The crowd (tidal forces) bumps them around, heats them up, and pushes them toward the center. The "Commuter" halo gets stripped of its outer layers and collapses into a dense core much faster than the Hermit.

The "King Model" Solution

When the halo finally collapses into a tiny, dense ball, the authors found that a specific mathematical shape called the King Model (usually used for star clusters) describes the density perfectly.

• Why it matters: This gives astronomers a simple formula to use when they look at real data. Instead of needing a supercomputer to guess what a collapsed halo looks like, they can just plug the King Model into their equations and say, "Ah, that matches!"

The Bottom Line

This paper is a manual for astronomers who want to simulate the universe. They say:

• Don't rush: If you force the simulation to run fast (big time steps), you get fake results.
• Watch your energy: If your simulation loses energy, the collapse happens too fast.
• Be precise: To explain weird observations like the "GD-1" star stream, you need to simulate the very center of the halo with extreme precision.

They have even released their best, most detailed simulation data (with 50 million particles) for others to use as a "benchmark"—a gold standard to test their own simulations against.

In short: They figured out how to build a perfect digital microscope to watch dark matter collapse, ensuring that what we see on the screen is actually what happens in the real universe.

Technical Summary

1. Problem Statement

Self-interacting dark matter (SIDM) is a leading candidate to explain small-scale structure anomalies in the universe. A key prediction of SIDM is gravothermal collapse: as dark matter particles scatter, they transfer heat from the hot outer regions to the cold core, eventually causing the core to contract and reach extremely high densities.

While this mechanism offers a potential explanation for compact dark matter substructures observed in gravitational lensing and stellar streams (specifically the GD-1 stream in the Milky Way), accurately simulating the late-stage collapse phase is numerically challenging. Previous studies have struggled with:

• Convergence issues: Determining if simulation results are physical or artifacts of numerical resolution.
• Energy conservation: Small errors in energy conservation can artificially accelerate or decelerate the collapse.
• Parameter sensitivity: The impact of softening lengths, time steps, and kernel sizes on the evolution of the halo is not fully understood.

The authors aim to establish a robust numerical framework for simulating SIDM halos in the deep collapse phase, providing benchmarks for interpreting observational data.

2. Methodology

The authors utilized the cosmological N-body code OpenGadget3, which features a specialized module for SIDM that explicitly conserves energy and linear momentum.

• Simulation Setup:

  • Initial Conditions (ICs): Spherically symmetric halos following an NFW profile ($\rho \propto r^{-1}(1+r/r_s)^{-2}$) with a mass of $\approx 2.8 \times 10^9 M_\odot$ and high concentration ($c \approx 22$).
  • Resolutions: Simulations were run with varying particle counts, up to a high-resolution benchmark of $N = 5 \times 10^7$ particles.
  • Scenarios:
    1. Isolated Halos: To study intrinsic gravothermal evolution.
    2. Satellite Halos: Placed in an external potential mimicking the Milky Way (NFW halo + Hernquist bulge + 4 Miyamoto-Nagai disks) to study tidal effects, motivated by the GD-1 stellar stream perturber.
  • Cross-sections: Tested isotropic vs. forward-dominated scattering, and velocity-independent vs. velocity-dependent cross-sections (using a specific functional form $\sigma \propto (1 + (v/w)^2)^{-2}$).

• Numerical Parameter Investigation:
  The study systematically varied key numerical parameters to assess their impact on accuracy:

  • Gravitational softening length ($\epsilon$).
  • Adaptive vs. fixed time steps.
  • Minimum time step limits ($\Delta t_{min}$).
  • SIDM kernel size (number of neighbors, $N_{ngb}$).
  • Truncation radius of initial conditions ($r_{cut}$).

3. Key Contributions and Results

A. Numerical Convergence and Accuracy

• Energy Conservation is Critical: The evolution of the halo is extremely sensitive to energy conservation errors. A loss of total energy artificially speeds up the collapse, while a gain slows it down. The authors recommend keeping energy conservation errors well below 1%.
• Softening Length & Time Steps: A very small softening length combined with adaptive time steps can lead to energy loss due to non-time-symmetric time-step changes. However, a constant time step can mitigate this. The optimal softening length for the collapse phase may be smaller than standard CDM recommendations to resolve the shrinking core.
• Kernel Size ($h$) vs. Mean Free Path ($l$): The ratio $h/l$ is crucial. While a ratio slightly above 1 is acceptable, ratios exceeding 10 lead to artificially strong effective heat conduction, causing premature collapse. The authors recommend $h/l < 10$.
• Initial Condition Truncation: Sampling the initial halo only up to $10 r_s$ (vs. $15 r_s$) can artificially accelerate the collapse in isolated simulations. A larger truncation radius is preferred for controlled studies.
• Minimum Time Steps: Imposing a minimum time step ($\Delta t_{min}$) allows simulations to run longer but introduces significant errors in the late collapse phase. It can lead to unphysical increases in total energy and incorrect predictions of enclosed mass, particularly for lower-resolution runs.

B. Physical Evolution of SIDM Halos

• Isolated vs. Satellite Evolution:
  • Isolated: Follows standard gravothermal evolution (core expansion $\to$ collapse).
  • Satellite: Tidal forces significantly accelerate the collapse. Tidal heating injects energy, increasing velocity dispersion, followed by mass stripping.
  • Velocity Dependence: For velocity-dependent cross-sections, tidal stripping reduces the internal velocity dispersion, which increases the effective cross-section (since $\sigma \propto v^{-4}$ in their model). This creates a feedback loop that accelerates collapse much more dramatically in satellites than in isolated halos.
• Density Profiles: In the late collapse phase, the inner density profile of both isolated and satellite halos is well-described by a King model ($\rho \propto (1 + (r/r_c)^2)^{-3/2}$) for radii $r < 0.5$ kpc.
• Projected Mass: The projected enclosed mass initially decreases during core expansion but increases during collapse. However, in the very late stages, mass inflow/outflow equilibrium shifts, causing the enclosed mass at fixed radii to eventually decrease.

C. Benchmark Data

The authors have made the data from their highest-resolution simulation ($N = 5 \times 10^7$) publicly available. This dataset serves as a benchmark for testing other SIDM codes and numerical schemes.

4. Significance and Implications

• Observational Constraints: The study validates that current numerical methods, when tuned correctly (high resolution, strict energy conservation, appropriate kernel sizes), can accurately simulate SIDM halos to the densities required to explain observations of GD-1 stellar stream gaps and strong gravitational lensing by compact subhalos.
• Guidance for Cosmological Simulations: The paper provides specific recommendations for numerical parameters (e.g., $\alpha = 5 \times 10^{-4}$ for tree opening, $\eta = 0.025$ for gravity time steps) to avoid artificial acceleration of collapse in large-scale cosmological simulations.
• Limitations of Current Methods: The authors note that simulating the extreme late stages (e.g., collapse into a black hole) using pure N-body methods becomes computationally prohibitive due to shrinking time steps and length scales. They suggest that subgrid models (e.g., sink particles) or fluid-based approaches may be necessary for the deepest collapse phases.
• Physics Insights: The work highlights that velocity-dependent cross-sections combined with tidal forces can produce substructures that are significantly more compact than those predicted by velocity-independent models, offering a new pathway to reconcile SIDM theory with observational constraints on substructure abundance.

In conclusion, this paper establishes a rigorous numerical foundation for studying gravothermally collapsing SIDM halos, resolving previous ambiguities regarding numerical artifacts, and providing the tools necessary to test SIDM against high-precision astrophysical observations.