Schwarzschild and Made-to-Measure Modeling

A galaxy sits in dynamical equilibrium, its billions of stars streaming along orbits no telescope can resolve one by one. Yet we want its mass — central black hole, stellar mass-to-light ratio, invisible dark halo — and the internal shape of its orbits, whether stars swing on radial plunges or circulate on tangential loops. The only data are a smooth image and, in each pixel, a distribution of line-of-sight velocities. How do you invert that into a six-dimensional-phase-space equilibrium? Martin Schwarzschild's 1979 answer was disarmingly constructive. Pick a trial gravitational potential; integrate a large library of orbits inside it; then find the non-negative amounts of light to place on each orbit so their superposed density and velocity structure reproduces the observations. Each orbit is a legitimate building block of an equilibrium — Jeans' theorem made numerical — and the weights are physical light, so they cannot go negative. If no non-negative combination fits, the trial potential is simply wrong, and that failure is how mass is measured. Made-to-Measure modeling, introduced by Syer and Tremaine in 1996, reaches the same equilibrium by a different road: a swarm of live particles orbits while slowly breathing their masses up or down under a force-of-change law until the model relaxes onto the data. Together these two techniques are the workhorses of modern galactic dynamics — the tools behind most dynamical black-hole masses and the only general method able to model triaxial galaxies, where the Jeans equations have no closure.

Key takeaways

Schwarzschild's method integrates a large library of stellar orbits in a fixed trial potential, then solves for non-negative orbit weights that reproduce the observed surface brightness and line-of-sight kinematics.
Made-to-Measure (M2M) instead evolves a live particle system, slowly adjusting particle weights via a force-of-change equation until the model matches the data, avoiding a stored orbit library.
Both are numerical realizations of Jeans' theorem: any steady-state distribution function is a superposition of orbits, so a galaxy's phase-space structure is a weighted sum over orbit families with $w_k\ge 0$ .
The fit is repeated over a grid of trial potentials (varying $M/L$ , black-hole mass $M_\bullet$ , halo, and flattening); the best potential minimizes a combined $\chi^2$ over image and kinematic moments.
These are the state-of-the-art tools for measuring black-hole masses, mass-to-light ratios, dark-matter fractions, and anisotropy $\beta(r)$ , including triaxial systems where Jeans modeling fails.
The chief limitation is the mass–anisotropy degeneracy, broken only by fitting the full LOSVD shape ( $h_3,h_4$ ) with controlled regularization.
Key idea: by Jeans' theorem every steady-state galaxy is a superposition of orbits; Schwarzschild and M2M modeling find the non-negative orbit (or particle) weights that reproduce the observed surface brightness and stellar kinematics in a trial potential, and the best potential is the one that fits best.

The core insight

By Jeans' theorem every steady-state galaxy is a superposition of orbits; Schwarzschild and M2M modeling determine the non-negative orbit weights (or particle weights) that reproduce the observed surface brightness and stellar kinematics in a trial potential, and the best potential is the one that fits best.

Where this lives today

Orbit-based modeling is the working currency of galaxy dynamics. Schwarzschild's method, in the axisymmetric and triaxial codes of Cappellari, van den Bosch, Valluri, and Vasiliev (the open-source `forstand`), delivers most dynamical black-hole masses in nearby galaxies and underpins the $M_\bullet$ – $\sigma$ relation. The M87 mass measured this way, $\sim 6.5\times10^{9}\,M_\odot$ , was strikingly confirmed by the Event Horizon Telescope's 2019 shadow image — an independent method landing on the same number. Integral-field surveys (ATLAS3D, SAMI, MaNGA) feed these models thousands of galaxies' $V,\sigma,h_3,h_4$ maps, turning orbit modeling into population statistics on anisotropy and shape. Made-to-Measure, meanwhile, is the tool of choice for the Milky Way itself: M2M models fit to Gaia proper motions, spectroscopic velocities, and infrared star counts (Portail, Wegg, Gerhard) reconstruct the mass, pattern speed, and 3D structure of the Galactic bar. The frontier now is triaxiality with dark matter and adaptive-optics kinematics inside the sphere of influence with ELT-class telescopes.

The physical picture

Think of the galaxy's light as a pile of sand you must reassemble, but each handful may only be spread along one fixed track — an orbit — inside a chosen potential. Schwarzschild's method asks: what non-negative amounts of sand on each of thousands of tracks reproduce both the observed image and the velocities pointing toward us in every pixel? Each orbit, integrated for many periods, deposits a stationary density in the cells it visits and contributes a characteristic LOSVD to each aperture. Summing orbits with weights is exactly building a distribution function as a sum over integrals of motion. The first figure shows the three orbit families that populate a general potential — loop (tube) orbits, box orbits, and precessing rosettes — each carrying a weight.

Figure

Three orbit families in one trial potential — a loop or tube orbit, a box orbit, and a precessing rosette — each carrying a non-negative weight wₖ; their weighted time-averaged densities and velocity signatures sum to the model galaxy, ρ(xᵢ) = Σₖ wₖ ρᵢₖ with wₖ ≥ 0

The potential search is the real physics. Orbit weights are so flexible that they adapt to almost any potential, so the discriminating power comes not from any single orbit but from whether a good non-negative fit exists at all, and how small the residual $\chi^2$ is across a grid of potentials. Add a point mass and refit; add a heavier halo and refit. The potential whose best superposition leaves the smallest $\chi^2$ is the measured one. The second figure is what a modeller reads off: the $\chi^2$ landscape versus trial black-hole mass, whose minimum weighs the black hole.

Figure

The potential search made concrete — reduced χ² of the best non-negative orbit fit versus trial black-hole mass for an M87-like model; the minimum at M• ≈ 6.5×10⁹ M⊙ selects the potential, and the Δχ² = 9 threshold marks the 3σ range of roughly 5.5 to 7.5 ×10⁹ M⊙

The core equations

The whole method rests on one fact: in a steady state the phase-space density $f$ does not change as you drift along an orbit. That is the steady-state collisionless Boltzmann equation (CBE), which every equilibrium galaxy obeys and the orbit library is built to satisfy term by term.

\frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla f - \nabla\Phi\cdot\frac{\partial f}{\partial \mathbf{v}} = 0 \qquad (1)

where: $f(\mathbf{x},\mathbf{v})$ — phase-space distribution function; $\Phi(\mathbf{x})$ — trial gravitational potential; $\mathbf{v}$ — velocity. In steady state the first term vanishes, leaving $f$ constant along the phase-space flow — that is, along each orbit. Jeans' theorem then says $f$ depends on phase space only through isolating integrals of motion, so each orbit is an elementary equilibrium and any galaxy $f$ is a superposition of orbits. Schwarzschild discretizes this: replace the continuous $f$ by a finite library indexed by $k$ , each orbit carrying a weight.

f(\mathbf{x},\mathbf{v}) = \sum_{k} w_k\, \delta_{\text{orbit }k}(\mathbf{x},\mathbf{v}),\qquad w_k \ge 0 \qquad (2)

where: $w_k$ — non-negative weight (physical light or mass) on orbit $k$ ; $k$ — index over the library. The inequality $w_k\ge 0$ is not a convenience but a physical constraint — an orbit cannot hold negative light — and it is what makes a bad potential detectable. Each orbit projects onto the data twice. Integrated for many periods it deposits a fractional time (a density) in each cell $i$ , and the observed luminosity density must equal the weighted sum of those contributions.

\rho(\mathbf{x}_i) = \sum_k w_k\, \rho_{ik}, \qquad \rho_{ik} = \text{fraction of time orbit } k \text{ spends in cell } i \qquad (3)

where: $\rho_{ik}$ — contribution of orbit $k$ to spatial cell $i$ ; $\rho(\mathbf{x}_i)$ — observed luminosity density in cell $i$ . One linear constraint per cell, all linear in the weights. The kinematic constraints work the same way on velocities: each orbit contributes a line-of-sight velocity distribution (LOSVD) to every aperture $j$ , and the weighted sum must match the observed.

\mathcal{L}_j(v_{\parallel}) = \sum_k w_k\, \ell_{jk}(v_{\parallel}) \qquad (4)

where: $\mathcal{L}_j$ — observed LOSVD in aperture $j$ ; $\ell_{jk}$ — LOSVD contributed by orbit $k$ ; $v_{\parallel}$ — line-of-sight velocity. In practice the LOSVD is compressed into a few Gauss–Hermite moments. The Gauss–Hermite series expands the LOSVD about a best-fitting Gaussian; its moments $h_3$ (skewness) and $h_4$ (kurtosis) are the diagnostics that separate radial from tangential orbits — the key to breaking the mass–anisotropy degeneracy.

\mathcal{L}(v)\propto \frac{e^{-y^2/2}}{\sigma}\Big[1 + h_3 H_3(y) + h_4 H_4(y)\Big],\quad y=\frac{v-V}{\sigma} \qquad (5)

where: $V,\sigma$ — mean velocity and dispersion; $h_3,h_4$ — skewness and kurtosis; $H_n$ — Hermite polynomials. Radially biased orbits give positive $h_4$ , tangential bias negative $h_4$ , so the LOSVD shape carries anisotropy information the dispersion alone hides. With more orbits than constraints the system is underdetermined, so a regularization term picks a smooth solution among the many that fit — the objective combines a $\chi^2$ data-mismatch with an entropy-like smoothing term.

\chi^2_{\text{tot}} = \sum_j \left(\frac{\text{model}_j - \text{data}_j}{\Delta_j}\right)^2 - \lambda\, S(\{w_k\}) \qquad (6)

where: $\Delta_j$ — uncertainty on constraint $j$ ; $S$ — entropy/smoothing regularizer; $\lambda$ — regularization strength. A common $S$ is the Boltzmann entropy $S=-\sum_k w_k\ln(w_k/\tilde w_k)$ relative to a smooth prior $\tilde w_k$ , suppressing spiky, noise-driven weights. Made-to-Measure reaches the same optimum without storing $\rho_{ik}$ : it evolves live particles and nudges each weight along the gradient of a merit function, giving the force-of-change equation.

\frac{\mathrm{d}w_k}{\mathrm{d}t} = -\varepsilon\, w_k \left[ \sum_j \frac{K_{jk}}{\Delta_j^2}\,\Delta_j^{\text{model}} \; - \; \mu\,\frac{\partial S}{\partial w_k}\right] \qquad (7)

where: $\varepsilon$ — learning rate; $K_{jk}$ — particle $k$ 's contribution to observable $j$ ; $\Delta_j^{\text{model}}$ — model-minus-data residual; $\mu$ — entropy weighting. The leading $w_k$ keeps weights positive and updates fractional, so M2M relaxes toward the same equilibrium Schwarzschild's NNLS solve reaches directly.

Derivation AFrom the steady-state Boltzmann equation to the orbit superposition

The goal of this derivation is to justify the central ansatz, equation (2): to show that a non-negative sum of time-averaged single-orbit densities is not merely a plausible fitting basis but an exact solution of the steady-state Boltzmann equation in the given potential, so every trial superposition is a genuine equilibrium and the only question left is which weights match the data.

Step 1. Define a single orbit's phase-space density as its time-average over one traversal (or many periods for a non-closed orbit), $f_k(\mathbf{x},\mathbf{v}) = \frac{1}{T_k}\int_0^{T_k}\delta[\mathbf{x}-\mathbf{x}_k(t)]\,\delta[\mathbf{v}-\mathbf{v}_k(t)]\,\mathrm{d}t$ . Why: an orbit visited over a long time is, on average, a static smear of density — this turns a trajectory into a mass distribution the telescope could see.

Step 2. Show each $f_k$ solves the CBE. The trajectory obeys Hamilton's equations, so $f_k$ is constant along the flow; the streaming operator $\mathcal{L}\equiv \mathbf{v}\cdot\nabla - \nabla\Phi\cdot\partial_{\mathbf{v}}$ gives $\mathcal{L}f_k = 0$ . Why: Jeans' theorem in action — a time-averaged orbit is a stationary distribution function by itself, because motion along it changes nothing about the average.

Step 3. Superpose. Because $\mathcal{L}$ is linear, $\mathcal{L}\!\big(\sum_k w_k f_k\big) = \sum_k w_k\,\mathcal{L}f_k = 0$ for any constant $w_k$ . Why: linearity is the whole engine — we are free to search over weights knowing every combination is automatically a valid equilibrium, so fitting never breaks the dynamics.

Step 4. Impose physicality and project. Requiring $w_k\ge 0$ with each $f_k\ge 0$ makes the total $f$ a genuine non-negative distribution function; projecting $f$ onto the sky yields the model image, equation (3), and the LOSVDs, equation (4). Why: this closes the loop — a valid, non-negative equilibrium whose projections are exactly the linear constraints the solver matches. Restricting to energy-only spherical orbits recovers the isotropic Eddington $f(E)$ , confirming the construction against a known result.

Go deeper · Derivation B · 1 minThe Made-to-Measure force-of-change from a merit function

Derivation BThe Made-to-Measure force-of-change from a merit function

The goal here is to derive the M2M update law, equation (7), as plain gradient ascent on a single merit function, so the seemingly ad hoc "nudge the heavy particles where the model is wrong" prescription is revealed as the optimization dual of Schwarzschild's constrained least squares — the two methods maximize the same objective, one by linear algebra and one by relaxation.

Step 1. Combine data agreement and the entropy prior into one functional, $F = -\tfrac12\chi^2 + \mu S$ , with $\chi^2$ the mismatch of equation (6) and $S$ the entropy. Why: turning "fit the data" and "stay smooth" into a single scalar to maximize lets a blind gradient march find the answer without ever inverting a matrix.

Step 2. Ascend the gradient with a $w_k$ prefactor, $\dot w_k = \varepsilon\, w_k\,\partial F/\partial w_k$ ; writing out $\partial F/\partial w_k$ recovers the bracket in equation (7). Why: the prefactor $w_k$ keeps every weight positive (a zero stays zero, a positive one moves fractionally) and makes the step scale-free, so light and heavy orbits adjust proportionally.

Step 3. Temporally smooth the observables over many orbital periods before forming the residual $\Delta_j^{\text{model}}$ . Why: a single particle's position jitters wildly along its orbit; averaging over periods makes the force respond to the orbit's steady contribution, so the swarm relaxes onto a stationary equilibrium instead of chasing transient positions. The fixed point $\dot w_k=0$ is exactly the stationary point of $F$ that Schwarzschild's NNLS solve locates.

Go deeper · Derivation C · 1 minEnclosed mass and the mass-anisotropy degeneracy

Derivation CEnclosed mass and the mass-anisotropy degeneracy

The goal of this derivation is to expose, from the second velocity moment of the Boltzmann equation, exactly what line-of-sight kinematics can and cannot constrain: the enclosed mass and the anisotropy enter the projected dispersion in one combination, so dispersions alone cannot separate them — and to show why fitting the full LOSVD shape breaks the tie.

Step 1. Define the anisotropy parameter $\beta(r) = 1 - \overline{v_\theta^2}/\overline{v_r^2}$ , positive for radial bias, negative for tangential. Why: $\beta$ is the single number encoding the internal orbit mix the weights $w_k$ determine — exactly the structure we are trying to measure, not assume.

Step 2. The radial Jeans equation — the second moment of equation (1) — reads $\frac{\mathrm{d}(\nu\overline{v_r^2})}{\mathrm{d}r} + \frac{2\beta}{r}\,\nu\overline{v_r^2} = -\nu\,\frac{GM(<r)}{r^2}$ , which rearranges to $\frac{GM(<r)}{r} = -\overline{v_r^2}\left(\frac{\mathrm{d}\ln\nu}{\mathrm{d}\ln r} + \frac{\mathrm{d}\ln\overline{v_r^2}}{\mathrm{d}\ln r} + 2\beta\right)$ . Why: the punchline in one line — $M$ and $\beta$ sit inside the same bracket, so a projected dispersion fixes the combination, not either separately.

Step 3. Read off the degeneracy: increasing $\beta$ (more radial orbits) lowers the required $M$ at fixed dispersion, so many $(M,\beta)$ pairs fit equally well. Why: this is why a dispersion profile alone cannot weigh a black hole or halo — the trap orbit modeling exists to escape.

Step 4. Break it with the LOSVD shape: radial bias gives positive $h_4$ , tangential bias negative. Why: Schwarzschild and M2M never assume $\beta$ — they fit $h_3,h_4$ aperture by aperture, pinning down $M$ and $\beta$ simultaneously wherever data quality allows. In the isotropic, flat-dispersion limit $\beta=0$ this reduces to $M(<r)=v_c^2 r/G$ .

Limiting cases

Case	Result	Meaning
Spherical isotropic limit	$f = f(E)$	Orbits reduce to a one-parameter family in energy; the method collapses to the classical Eddington inversion and agrees with the isotropic Jeans equation.
Axisymmetric limit	$f=f(E,L_z)$ or $f(E,L_z,I_3)$	Two or three integrals; short-axis tube orbits dominate. Standard for disks and rotating bulges with integral-field data.
Triaxial general case	$f=f(E,I_2,I_3)$	Box, short-axis tube, and long-axis tube families all required; only Schwarzschild and M2M can handle this — the Jeans equations lack a closure.
Isotropic isothermal	$M(<r)=v_c^2 r/G$	With $\beta=0$ and flat dispersion the enclosed mass reduces to the elementary circular-speed estimate — the sanity check on the machinery.
Collisionless validity	$t_{\text{relax}}\gg t_{\text{age}}$	The static- $f$ assumption needs the two-body relaxation time to exceed the system age, true for galaxies ( $N\sim 10^{11}$ ) but not for globular-cluster cores.

Common trapA dispersion profile alone can weigh the black hole or the halo

It cannot. Line-of-sight kinematics constrain the enclosed mass and the anisotropy only in the combination that appears in the Jeans equation, $GM/r \sim -\overline{v_r^2}(\cdots + 2\beta)$ — so a more radial orbit distribution ( $\beta$ up) mimics extra mass, and different potentials fit the dispersion equally well. This mass–anisotropy degeneracy is the single most important caveat in dynamical modeling. It is broken only by fitting the full LOSVD shape — the $h_3,h_4$ moments — which register radial versus tangential orbits through their distinct kurtosis, and only with adequate data and a tested regularization. Never quote a black-hole or halo mass without stating the assumed anisotropy freedom and the regularization; a number from dispersions alone is not a measurement but a family of them.

Worked examples

1 · Weighing M87's black hole — the sphere of influence (fully worked)

A Schwarzschild model of M87 must resolve the region where the black hole dominates. For $M_\bullet = 6.5\times10^{9}\,M_\odot$ , central dispersion $\sigma = 300\ \mathrm{km\,s^{-1}}$ , and distance $D=16.7\ \mathrm{Mpc}$ , find the sphere-of-influence radius and the angular resolution required.

Step 1. The sphere of influence is where the black hole's Keplerian speed matches the dispersion, $r_\bullet = GM_\bullet/\sigma^2$ . Use $G = 4.301\times10^{-3}\ \mathrm{pc}\,(\mathrm{km\,s^{-1}})^2\,M_\odot^{-1}$ so units close in parsecs.

Step 2. Substitute: $r_\bullet = (4.301\times10^{-3})(6.5\times10^{9})/(300)^2\ \mathrm{pc} = (2.796\times10^{7})/(9.0\times10^{4})\ \mathrm{pc} = \mathbf{311\ pc}$ .

Step 3. Convert to an angle, $\theta_\bullet = r_\bullet/D = 311\ \mathrm{pc}/(16.7\times10^{6}\ \mathrm{pc}) = 1.86\times10^{-5}\ \mathrm{rad}$ . Multiply by $206265$ arcsec per radian: $\theta_\bullet = \mathbf{3.8\ arcsec}$ .

Step 4. Interpret. A 3.8 arcsec sphere of influence is comfortably resolved by HST and adaptive-optics spectrographs (sub-arcsec PSF), so the kinematics inside $r_\bullet$ are well sampled and the $\chi^2(M_\bullet)$ curve has a sharp minimum. The black hole dominates within ~311 pc ≈ 3.8 arcsec, which the data resolve. Had $\theta_\bullet$ fallen below the PSF, $M_\bullet$ would be unconstrained regardless of library size.

2 · How anisotropy fakes a factor of two in mass (fully worked)

At $r=5\ \mathrm{kpc}$ in a spherical elliptical, let the tracer density fall as $\nu\propto r^{-2}$ ( $\mathrm{d}\ln\nu/\mathrm{d}\ln r = -2$ ) with radial dispersion $\overline{v_r^2}^{1/2}=200\ \mathrm{km\,s^{-1}}$ , roughly flat ( $\mathrm{d}\ln\overline{v_r^2}/\mathrm{d}\ln r\approx 0$ ). Compare the enclosed mass for isotropy ( $\beta=0$ ) versus radial anisotropy ( $\beta=0.5$ ).

Step 1. From the rearranged Jeans relation, $\dfrac{GM(<r)}{r} = -\overline{v_r^2}\,(-2 + 0 + 2\beta) = 2\overline{v_r^2}(1-\beta)$ .

Step 2. Isotropic $\beta=0$ : $M = 2\overline{v_r^2}\,r/G$ . With $G = 4.301\times10^{-6}\ \mathrm{kpc}\,(\mathrm{km\,s^{-1}})^2\,M_\odot^{-1}$ , $M = 2\,(200^2)(5)/(4.301\times10^{-6}) = \mathbf{9.3\times10^{10}\,M_\odot}$ .

Step 3. Radial case $\beta=0.5$ : the bracket $2(1-\beta)$ halves, so $M = \overline{v_r^2}\,r/G = \mathbf{4.7\times10^{10}\,M_\odot}$ — a factor of exactly two less, from the same dispersion.

Step 4. The lesson. The projected dispersion cannot tell these two galaxies apart; only the LOSVD shape can. Radial orbits at $\beta=0.5$ produce a positive $h_4$ , and it is by fitting $h_4$ aperture by aperture that a Schwarzschild model refuses the wrong $(M,\beta)$ pair. A 0.5 shift in $\beta$ moves the inferred mass by a factor of 2 — which is why anisotropy must be fit, not assumed.

3 · Sizing the linear system

A modern axisymmetric model might integrate $N_{\text{orb}}\approx 10^{4}$ orbits against $N_{\text{con}}\approx 10^{3}$ constraints ( $V,\sigma,h_3,h_4$ over an integral-field map plus brightness cells). Over- or under-determined?

Step 1. With $N_{\text{orb}}\gg N_{\text{con}}$ the non-negative least-squares system has far more unknowns than equations, so it is under-determined — infinitely many weight sets fit within errors.

Step 2. That is why the regularizer $\lambda S$ of equation (6) is mandatory: it picks the smoothest weight set, and the recovered structure can depend on $\lambda$ . Orbit weights are never unique — regularization chooses among solutions, and its strength must be tested.

Physics of the portrait

The method carries one man's name because he saw, before the computers were really ready, that an equilibrium galaxy could be built rather than guessed — orbit by orbit, with nothing assumed about its distribution function.

Origins · 1 minMartin Schwarzschild and how this was discovered

Martin Schwarzschild

1912–1997 · numerical stellar and galactic dynamics

Martin Schwarzschild was born in Potsdam in 1912, the son of Karl Schwarzschild, who found the first exact solution of Einstein's field equations. Fleeing Nazi Germany, he settled at Princeton and became one of the twentieth century's foremost stellar astrophysicists — his 1958 Structure and Evolution of the Stars set the template for the field, and he was an early, fearless user of computers for physics. Late in his career he turned to a problem that had resisted every analytic attack: how to construct a self-consistent equilibrium model of a triaxial stellar system, the shape N-body simulations were then suggesting for elliptical galaxies. The Jeans equations were useless there — no closure exists for a triaxial potential. In his 1979 paper he cut the knot with a numerical idea of startling directness: choose the potential, integrate a library of orbits, and solve for the non-negative orbit occupation numbers that reproduce the model's own density. That a triaxial equilibrium could be assembled at all, from box and tube orbits alone, settled a real controversy about whether such galaxies could exist. The method was soon turned from a consistency demonstration into a measurement tool, and today "Schwarzschild modeling" is the standard route to dynamical black-hole masses and galaxy shapes. He completed what Jeans began: a galaxy's unseen distribution function, reconstructed from its orbits.

History

1915

— Jeans shows that a steady-state distribution function depends on phase space only through integrals of motion

1979

— Schwarzschild builds the first self-consistent orbit-superposition model of a triaxial galaxy

1993

— van der Marel and Franx parametrize the line-of-sight velocity distribution with Gauss-Hermite moments

1996

— Syer and Tremaine introduce made-to-measure N-body modeling

2011

— Gebhardt and collaborators weigh M87's black hole at about 6.6 billion solar masses by Schwarzschild modeling

Going further

Texts. Binney, J. & Tremaine, S., Galactic Dynamics, 2nd ed., §4.4–4.6 (Jeans' theorem and Schwarzschild's method) · Mo, van den Bosch & White, Galaxy Formation and Evolution, Ch. 5–6 · Carroll & Ostlie, An Introduction to Modern Astrophysics, 2nd ed., Ch. 24–25 · Schwarzschild, M. (1979), ApJ 232, 236 (founding paper) · Syer, D. & Tremaine, S. (1996), MNRAS 282, 223 (made-to-measure) · van der Marel & Franx (1993), ApJ 407, 525 (Gauss–Hermite LOSVD) · Cappellari, M. (2008), MNRAS 390, 71. Key consistency checks: (1) a non-negative orbit sum solves the CBE by linearity; (2) the second moment of the CBE reduces to the Jeans equation with the $M$ – $\beta$ combination; (3) the isotropic limit gives $M(<r)=v_c^2 r/G$ ; (4) the spherical isotropic case recovers the Eddington $f(E)$ . Worked numbers use $G=4.301\times10^{-3}\ \mathrm{pc}\,(\mathrm{km/s})^2 M_\odot^{-1}$ and the M87 anchors ( $M_\bullet\approx6.5\times10^{9}\,M_\odot$ , $\sigma\approx300\ \mathrm{km/s}$ , $D=16.7$ Mpc → $r_\bullet\approx311$ pc $\approx 3.8''$ ) and the factor-of-two anisotropy illustration. No galaxy-specific weights are invented.

Prerequisites. Collisionless Boltzmann Equation · Jeans' Theorem · Distribution Function and Phase Space · Orbits in Axisymmetric and Triaxial Potentials · Galactic Potential Models. Leads to. Jeans Equations · Action–Angle Variables · Gravitational N-Body Problem.

Ask the tutor: Curious how a wrong trial potential leaves no non-negative orbit combination that fits — ask me to walk through a toy two-orbit example where the NNLS residual refuses to vanish.

Why does a bad trial potential show up as a failed fit rather than just worse numbers?

Because the weights are constrained to be non-negative. In the right potential some non-negative combination of orbits reproduces both the image and the kinematics within the errors. In a wrong potential the orbit shapes are wrong, and forcing the image to match may leave no non-negative way to also match the velocities — the solver cannot drive the residual down. So a wrong potential announces itself as an irreducibly large $\chi^2$ , not a slightly worse but still-good fit. Scanning potentials and reading off where $\chi^2$ is minimized is exactly how $M_\bullet$ , $M/L$ , and halo parameters are measured, as the second figure shows.

What is the difference between Schwarzschild's method and Made-to-Measure?

Both find the non-negative orbit weights that reproduce the data, but by opposite routes. Schwarzschild integrates a fixed library, stores each orbit's contribution to every constraint in a matrix, and solves a one-shot non-negative least-squares problem. Made-to-Measure never stores that matrix: it runs a live particle simulation and adjusts each particle's mass by the force-of-change law, relaxing toward the fit. Schwarzschild is exact and transparent but memory-hungry; M2M scales to enormous particle numbers and self-consistent, evolving potentials, which is why it is favored for the Milky Way bar. At convergence they optimize the same merit function.

Do box and tube orbits actually matter, or are they jargon?

They are the physical alphabet the method spells galaxies with. In an axisymmetric potential, angular momentum about the symmetry axis is conserved, so stars circulate on short-axis tube orbits — loops that never reach the center. A triaxial potential has no such conservation law and admits box orbits that oscillate through the center, plus long-axis tubes. A triaxial galaxy's characteristic boxy shape requires box orbits to build it — you cannot assemble a triaxial equilibrium without them. That is why Jeans modeling, which needs a smooth closure, fails for triaxial systems, while Schwarzschild's explicit orbit library succeeds.

Why is regularization necessary, and what goes wrong if you get it wrong?

Because there are far more orbits than constraints, so infinitely many weight sets fit equally well — the problem is under-determined. Regularization (typically a maximum-entropy term) picks the smoothest solution, suppressing spiky, noise-driven weight patterns. The danger is real: too strong pulls the answer toward the prior (often isotropy), biasing the recovered $\beta(r)$ and even the best-fit potential; too weak overfits noise. Because recovered structure can depend on the strength $\lambda$ , a credible result reports how the answer moves as $\lambda$ is varied. Regularization is part of the measurement, and it must be disclosed.

How does adding a central point mass let these models weigh a black hole?

You include a Keplerian point mass $M_\bullet$ in the trial potential, integrate the orbit library, fit the kinematics, then repeat over a grid of $M_\bullet$ . Inside the sphere of influence, $r_\bullet=GM_\bullet/\sigma^2$ , stars move faster and the LOSVDs broaden in a way only the right $M_\bullet$ can reproduce, so $\chi^2(M_\bullet)$ has a minimum. The measurement only works if the data resolve $r_\bullet$ : for M87 that is about 3.8 arcsec, easily resolved, but for distant or lower-mass galaxies $r_\bullet$ can drop below the PSF, and then $M_\bullet$ is bounded, not measured. This is why adaptive optics and space telescopes are essential.

Why can't Jeans modeling just do everything Schwarzschild modeling does?

Jeans modeling takes velocity moments of the Boltzmann equation, but each moment introduces the next, so the chain never closes — you must assume a relation, typically a fixed anisotropy $\beta(r)$ . That assumption is exactly the quantity you most want to measure, and in a triaxial potential no consistent closure even exists. Schwarzschild modeling assumes nothing about $\beta$ : it reconstructs the full distribution function as an orbit superposition and reads $\beta(r)$ off the weights. The cost is enormous computation across a parameter grid, which is why fast Jeans-based methods (like JAM) stay popular for large surveys where the extra generality is not needed.