Stellar Cusps and Core Profiles

A stellar cusp is an intrinsic density that rises toward a massive black hole; a core is a resolved flattening or deficit relative to an explicitly chosen outer profile. Neither is read directly from a surface-brightness slope. Relaxation can produce a cusp, a black-hole binary can scour a core, collisions can remove luminous giants without removing total mass, and recent star formation can add a population-specific excess. The useful comparison therefore follows the distribution function through projection and selection, asks whether the system had time to relax, and keeps tracer light separate from the unseen mass that drives the dynamics.

Required background. Classical Field Theory supplies the phase-space distribution and projection integrals used below.

Helpful background. Black-Hole Geometry and Mechanics identifies the central spacetime and capture boundary; Relativistic Matter and Compact Objects supplies the stellar and remnant species whose slopes can differ; Black-Hole Astrophysics places the profile in its observable source context; and Relativistic Two-Body Dynamics clarifies the orbital regimes represented by energy and angular momentum.

Relaxed cusps in a Kepler potential

Inside the influence radius, take the relative potential to be $\Psi(r)=GM_\bullet/r$ and define positive binding energy per unit mass by $\mathcal E=\Psi-v^2/2$. For an isotropic number distribution in position–velocity phase space, take $f(\mathcal E)=f_0\mathcal E^p$ for bound stars and $f=0$ for $\mathcal E\le0$. With $p\gt-1$ so that the velocity integral converges at escape, the number density is

$$n(r)=4\pi\int_0^{\sqrt{2\Psi}} f\!\left(\Psi-\frac{v^2}{2}\right)v^2\,dv.$$

Set $v=\sqrt{2\Psi}\,x$. The remaining integral over $x$ is dimensionless, so

$$n(r)\propto \Psi^{p+3/2}\propto r^{-(p+3/2)}.$$

The single-mass, isotropic Bahcall–Wolf zero-flow solution has $p=1/4$ and hence

$$f(\mathcal E)\propto\mathcal E^{1/4}, \qquad n(r)\propto r^{-7/4}.$$

“Zero flow” refers to the net number current of stars through binding-energy space, not to a vanishing outward transport of energy. Diffusion and drift nearly balance; stars still exchange energy. A physical solution with an inner sink has a small nonzero stellar current and differs near its boundaries. Bahcall and Wolf 1976, § III(c), eqs. (59)–(60), and § III(d), pp. 220–223 derives the power law and distinguishes the stellar current from the outward energy flow. With several masses, heavier objects generally acquire steeper profiles, and the branch depends on their abundance as well as their mass; see Alexander 2017, § 2.1.4, author manuscript pp. 12–14.

A scaling argument explains why the slope $7/4$ is natural. For $n\propto r^{-\gamma}$ in a Kepler potential, $\sigma\propto r^{-1/2}$ and the equal-mass relaxation time scales as $t_r\propto r^{\gamma-3/2}$, ignoring the weak radial dependence of the Coulomb logarithm. A shell of logarithmic radial width contains $N_r\sim n r^3$ stars, each with characteristic binding energy proportional to $1/r$. The outward energy-transport scale is therefore

$$\dot E_{\rm transport}\sim\frac{N_r E_\star}{t_r} \propto r^{7/2-2\gamma}.$$

Making this independent of radius gives $\gamma=7/4$. This is an interpretation of the steady heat flow, not a replacement for solving the drift–diffusion equation with its inner sink and outer stellar reservoir.

For a spherical power law $n\propto r^{-\gamma}$ with $1<\gamma<3$, line-of-sight projection gives $\Sigma(R)\propto R^{1-\gamma}$ away from inner and outer boundaries. Thus the ideal $\gamma=7/4$ cusp projects to slope $-3/4$. The extra one in the exponent is easy to lose when comparing a theoretical three-dimensional cusp to star counts or integrated light.

Comparison axes for cusps and cores

The mechanisms differ in more than their fitted slope.

Mechanism	Intrinsic expectation	Discriminating axis	Principal limitation
Collisional steady cusp	Species-dependent rise inside a relaxed region	Stellar age compared with energy-relaxation time	A formal steady solution need not be reached
Binary black-hole scouring	Phase-space deficit and tangentially biased survivors	Mass deficit, anisotropy, and merger history	The unobserved pre-merger profile is model dependent
Giant destruction or envelope stripping	Deficit in a luminous tracer	Spectral type and collision radius	Light can flatten while remnant mass remains cusped
Recent nuclear star formation	Age-selected disk, ring, or central excess	Age, metallicity, rotation, and geometry	Projection can imitate a spherical cusp
Unresolved observation	Point-spread-smeared central profile	Injection and recovery at the actual resolution	A fitted inner slope below the resolution radius is not measured

The Galactic Centre makes these distinctions concrete: the observed late-type giant profile is much flatter than a textbook single-mass relaxed cusp, yet that fact alone does not determine the distribution of faint main-sequence stars or stellar remnants Genzel et al. 2010, §2.4. Across galaxies, size, mass, population, and resolution vary together, so a demographic trend is not automatically an evolutionary sequence Neumayer et al. 2020, §§5.3–5.6, pp. 15–21.

An honest dynamical calculation must specify $M_\bullet$, the mass spectrum, initial $f(\mathcal E,J)$, anisotropy, capture sink, outer source, and any merger or star-formation history. Orbit-averaged Fokker–Planck, Monte Carlo, and direct N-body methods represent different approximations. Agreement in a projected slope is meaningful only after their phase-space boundaries and projection kernels have been aligned.

A numerical projection benchmark

Before evolving a cluster, test the projection and fitting analysis against a known cusp and a cored control. The small program below integrates along the line of sight after the change of variable $r=R\cosh u$. It fits projected slopes over $0.01\le R\le0.1$ with outer radius $r_{\rm out}=100$, for $\gamma=7/4$. The cored control replaces $r^{-\gamma}$ by $(r^2+r_c^2)^{-\gamma/2}$ with $r_c=0.03$.

Use Python 3.12.5 and NumPy 2.1.0. Save the following program as cusp_projection_demo.py. The script emits the resolution, cusp slope, its difference from a high-resolution reference, and the cored-control slope.

import numpy as np

def projected(R, gamma, core_radius, n, r_out=100.0):
    u = np.linspace(0.0, np.arccosh(r_out / R), n + 1)
    radius = R * np.cosh(u)
    density = (radius**2 + core_radius**2) ** (-0.5 * gamma)
    return np.trapezoid(2.0 * density * R * np.cosh(u), u)

def fitted_slope(n, core_radius):
    radii = np.geomspace(0.01, 0.1, 24)
    sigma = np.array([projected(R, 1.75, core_radius, n) for R in radii])
    return np.polyfit(np.log(radii), np.log(sigma), 1)[0]

reference = fitted_slope(65536, 0.0)
for resolution in (32, 64, 128):
    cusp = fitted_slope(resolution, 0.0)
    core = fitted_slope(resolution, 0.03)
    print(resolution, cusp, abs(cusp - reference), core)

Create the pinned environment and rerun the recipe with:

python3 -m venv cusp-demo
source cusp-demo/bin/activate
python3 -m pip install numpy==2.1.0
python3 cusp_projection_demo.py

The expected teaching-run values are:

Resolution N	Cusp slope result	Error versus reference	Cored slope result
32	-0.7513382341	3.09e-6	-0.3916452702
64	-0.7513359163	7.73e-7	-0.3916442500
128	-0.7513353367	1.93e-7	-0.3916439951

The quadrature error falls by about four when $N$ doubles, as expected for the trapezoidal rule on this smooth transformed integrand. The limiting cusp slope is slightly steeper than $-0.75$ because the model has a finite outer boundary. That offset is a physical truncation effect in this benchmark, not a numerical failure.

Computational boundaries and failure modes

This executable recipe is a projection and regression test. It is not a claimed production simulation of a galactic nucleus, does not solve a Fokker–Planck equation, and contains no capture, source, anisotropy, mass segregation, star formation, or merger physics. A production calculation would need a documented solver version, initial distributions, energy and angular-momentum grids or particle number, stochastic seeds, physical time mapping, capture boundary, and checkpoints.

At minimum, such a calculation should first recover the equal-mass $r^{-7/4}$ benchmark over a radial interval separated from both boundaries. It should also preserve a deliberately depleted-core control when relaxation is turned off. Refinement must increase energy and angular-momentum resolution or particle number, repeat random realizations, and compare intrinsic slopes, projected slopes, anisotropy, and capture flux. A three-row table of fitted slopes alone cannot establish convergence if the resolved interval moves with resolution.

Important failure modes include artificial pile-up at a reflecting inner grid edge, capture flux set by the first cell, N-body relaxation accelerated by too few particles, unrelaxed initial transients mistaken for a steady cusp, and fitting inside the point-spread or force-softening scale. A central hole in bright giants is also a counterexample to the inference “flat light means flat total mass.” The unique dynamics of dense nuclei, including collision and relaxation qualifications, are reviewed in Neumayer et al. 2020, §10.1, pp. 36–37.

The scope remains the stellar environment dynamically tied to the central black hole. Generic galaxy cores produced without reference to a nuclear black hole require their own formation and selection analysis and are not classified here by surface-brightness shape alone.

Exercises

The converged teaching result is approximately $-0.751335$, not exactly $-0.75$. A student proposes increasing $N$ until the difference from $-0.75$ vanishes. Diagnose the proposal and name one calculation that would actually test the explanation.

Solution

The proposal confuses discretization error with the finite-$r_{\rm out}$ model. The table shows quadrature changes shrinking toward a stable value that remains about $0.001335$ from the infinite-power-law slope. Hold the numerical tolerance fixed and increase $r_{\rm out}$, or narrow the fitted interval farther from the boundary; the fitted slope should then approach $-0.75$. This is a boundary-domain test, distinct from resolution refinement, and both are needed before interpreting a physical slope difference.

References

• Alexander, Tal. 2017. “Stellar Dynamics and Stellar Phenomena Near a Massive Black Hole.” Annual Review of Astronomy and Astrophysics 55: 17–57. DOI.
• Bahcall, John N., and Richard A. Wolf. 1976. “Star Distribution around a Massive Black Hole in a Globular Cluster.” The Astrophysical Journal 209: 214–232. DOI; Open PDF.
• Genzel, Reinhard, Frank Eisenhauer, and Stefan Gillessen. 2010. “The Galactic Center Massive Black Hole and Nuclear Star Cluster.” Reviews of Modern Physics 82: 3121–3195. DOI.
• Neumayer, Nadine, Anil Seth, and Torsten Böker. 2020. “Nuclear Star Clusters.” Astronomy and Astrophysics Review 28: 4. DOI.