Jeans Anisotropic Models (JAM) for Stellar Dynamics

• Jeans Anisotropic Models (JAM) are dynamical techniques that solve the Jeans equations under anisotropic assumptions to recover gravitational potentials.
• JAM utilizes Multi-Gaussian Expansion (MGE) to predict projected kinematic observables, including line-of-sight velocities and proper motions.
• The framework constrains dark matter profiles, stellar mass-to-light ratios, and orbital anisotropy, achieving accuracy within 5–15% in simulations.

Jeans Anisotropic Models (JAM) constitute a class of dynamical modeling techniques for stellar systems and galaxies that solve the axisymmetric or spherical Jeans equations under explicit assumptions about the anisotropy of the velocity dispersion tensor. JAM models leverage Multi-Gaussian Expansion (MGE) parameterizations of observed (or simulated) tracer and mass distributions, enabling computationally efficient predictions of projected kinematic observables—including not only the line-of-sight velocity moments but also proper motions—under axisymmetry or spherical symmetry. This framework has become a foundational tool for reconstructing the underlying gravitational potential of galaxies and stellar systems, constraining dark matter profiles, supermassive black holes, and orbital anisotropy.

1. Mathematical Formulation and Foundations

The JAM approach is based on solving the steady-state Jeans equations, which follow from the second moments of the collisionless Boltzmann equation. When axisymmetry is assumed and cylindrical alignment of the velocity ellipsoid is imposed, the key equations in cylindrical coordinates $(R, \phi, z)$ are: $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$

$$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$

with $\nu(R, z)$ the intrinsic tracer density and $\Phi(R, z)$ the gravitational potential.

The traditional JAM procedure aligns the velocity ellipsoid with the cylindrical coordinate grid, assumes constant anisotropy $\beta_z = 1 - (\sigma_z^2/\sigma_R^2)$ (or allows this parameter to vary per Gaussian component in advanced implementations), and adopts parameterized profiles for the total mass distribution via MGE. For spherically aligned velocities, the equivalent Jeans equation is

$$\frac{\partial(\nu\overline{v_r^2})}{\partial r} + \frac{2\beta\,\nu\,\overline{v_r^2}}{r} = -\nu\frac{\partial\Phi}{\partial r}$$

with the anisotropy $$\beta=1-\langle v_\theta^2\rangle/\langle v_r^2\rangle$$.

A major advance is the ability to compute all projected second moments (including proper motions), as given in JAM implementations: $$\Sigma \langle v_\alpha v_\beta \rangle(x', y') = 4\pi^{3/2} G \int_0^1 \sum_k \sum_j \nu_{0k} q_j \rho_{0j} u^2 \mathcal{F}_{\alpha\beta} \frac{\exp\left(-\mathcal{A}[x'^2 + y'^2 (\mathcal{A}+\mathcal{B})/\mathcal{E}]\right)}{(1-\mathcal{C}u^2)\sqrt{\mathcal{E}[1-(1-q_j^2)u^2]}} \, du$$ where mathematics of the terms $$\nu_{0k}, q_j, \rho_{0j}, \mathcal{F}_{\alpha\beta}, \mathcal{A}, \mathcal{B}, \mathcal{C}, \mathcal{E}$$ are all determined by the MGE decompositions and anisotropy structure (Cappellari, 2012).

2. Implementation, Software, and Validation

The JAM methodology is implemented as an efficient, publicly available software package (Cappellari, 2012, Cappellari, 2015, Cappellari, 2019). Key computational features include:

• Direct, analytic computation of line-of-sight, radial, and tangential proper motion moments using quadrature, without recourse to special functions except in the generalized spherical case.
• Integration with surface-brightness-to-MGE routines, enabling seamless transformation from observed photometry to model input (Cappellari, 2012).
• Adaptability for multiple kinematic components, arbitrary inclination, axisymmetric or spherically aligned solutions, and NFW/gNFW dark matter profiles.
• Native support for modern programming environments (IDL, Python, C), facilitating large-scale parameter inference studies (e.g., Markov Chain Monte Carlo sampling).

Extensive validation against N-body and hydrodynamical simulations (e.g., Illustris for galaxies (Li et al., 2015), Auriga for dwarfs (2206.12121), TNG300 for clusters (Shi et al., 2024), and star particle sets for the Milky Way (Zhang et al., 14 May 2025)) consistently shows that JAM robustly recovers the total mass profile within critical radii (5–15% accuracy, depending on the system), although there is unavoidable degeneracy between stellar and dark matter contributions.

3. Anisotropy, Proper Motions, and Degeneracies

A key strength of JAM is explicit modeling of velocity anisotropy, which is essential for breaking the classical mass–anisotropy degeneracy. Modern versions provide for arbitrary, possibly radially varying, anisotropy profiles ($$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$0) or assign different $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$1 values per Gaussian in the MGE decomposition (Cappellari, 2015). Including proper motion second moments (radial and tangential) in the modeling—feasible with high-precision astrometry from Gaia or HST—dramatically constrains $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$2, as proper motions sample orthogonal combinations of the underlying distribution function (Cappellari, 2012, Cappellari, 2015, 2206.12121).

In contrast, if only line-of-sight information is used, the total mass profile is still robustly recovered within the tracer half-mass radius, but the breakdown between dark matter, stars, and the detailed shape of the distribution function suffers from a strong degeneracy (Li et al., 2015, 2206.12121). This is particularly severe in triaxial systems or low-resolution imaging, where JAM defaults to an oblate deprojection and can overestimate $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$3 by 18% for strongly prolate galaxies.

4. Physical Applications: Galaxies, Clusters, and the Milky Way

JAM models see wide-ranging application across galaxy types and scales:

• Early-type galaxies: JAM constrains mass-to-light ratios, total density slopes (e.g., $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$4), dark matter fractions, and black hole masses using integral field kinematics (Cappellari, 2012, Cappellari, 2019, Lu et al., 2023).
• Spiral and late-type galaxies: JAM is essential for capturing thick disc and bulge contributions, especially in dynamically hot systems where thin-disc (ADC) methods dramatically underestimate $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$5 (Kalinova et al., 2016).
• Dwarf galaxies: Mass within the half-light radius and, critically, inner dark matter content, can be constrained unbiasedly in ensemble (scatter $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$60.17 dex), highlighting the power and the caution required for systems out of equilibrium (subject to contraction or outflows) (2206.12121).
• The Milky Way: JAM, particularly the spherically aligned variant (JAM$$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$7), is applied with full phase-space, multi-tracer (KG/BHB) samples and Gaia kinematics to infer both the mass and shape (flattening $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$8) of the Galactic halo, with evidence for a radially varying dark halo shape—more oblate at small radii, rounder by $$\frac{\partial(\nu\langle v_R^2\rangle)}{\partial R} + \frac{\partial(\nu\langle v_R v_z\rangle)}{\partial z} + \nu\left(\frac{\langle v_R^2\rangle - \langle v_\phi^2\rangle}{R}\right) = -\nu \frac{\partial \Phi}{\partial R}$$9 kpc (Zhang et al., 14 May 2025).
• Galaxy clusters: Satellite kinematics modeled with JAM yield unbiased total masses within the tracer half-mass radius/boundary, even under redshift-space selection and realistic survey limitations (biases $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$00.06 dex, scatters $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$10.18 dex) (Shi et al., 2024).

JAM retains accuracy for the total density even in the presence of observed complexities such as edge-on dust, bars, and boxy/peanut bulges. Kinematic residuals from axisymmetric JAM fits can be exploited to identify and interpret signatures of such non-axisymmetric structures (Rutherford et al., 10 Sep 2025).

5. Extensions: Spherical, Cylindrical, and External Field Effects

JAM is available in multiple geometrical configurations, each appropriate for different physical regimes:

• Cylindrically-aligned (JAM$$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$2): Best suited for disk-dominated, fast-rotating systems. Velocity ellipsoid is aligned with $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$3.
• Spherically-aligned (JAM$$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$4): Natural for nearly spherical halos and systems where the orbits are not disk-supported (Cappellari, 2019). Solutions are provided by integrating along characteristics in $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$5, with explicit dependence on the anisotropy parameter $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$6 (see equations above).
• Spherical formalism: For full spherical symmetry, all three projected second moments (including two proper motion components) can be calculated efficiently (Cappellari, 2015).

In cosmological and galactic environments, the instability criteria underpinning the JAM framework can be modified by additional physics. For example, external tidal fields cause direction-dependent collapse thresholds and promote filamentary structure formation, requiring anisotropic Jeans masses based on the eigenvalues of the tidal tensor (Li, 2024). In modified gravity scenarios (e.g., $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$7 models), the linear stability analysis leads to a scale-dependent and model-dependent Jeans mass, directly testable with dense molecular cloud observations (Vainio et al., 2015).

6. Limitations and Interpretational Cautions

JAM models assume steady-state equilibrium and axisymmetry or spherical symmetry. When these assumptions are violated (e.g., strong outflows, recent mergers, or persistent contraction/expansion flows in simulated dwarfs), systematic biases arise: the total mass within the half-mass radius is still generally well constrained, but inner mass in quiescent, contracting systems is typically underestimated, and star-forming galaxies with strong feedback may show overestimates (2206.12121).

Further, the breakdown of degeneracy between luminous and dark matter components is limited if only line-of-sight information is used; precise measurements of proper motions or external constraints (e.g., stellar population synthesis M/L, spatially resolved IMF gradients) are required for unique decomposition (Lu et al., 2023). In triaxial or strongly barred/misaligned systems, axisymmetric JAM solutions may yield biased inferences regarding $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$8 by up to 18% (Li et al., 2015).

Careful masking of dust and exclusion of contaminated data is necessary, but as shown by (Rutherford et al., 10 Sep 2025), robust dynamical parameters can be recovered if affected regions are excluded, with dust masking affecting results at the $$\frac{\partial(\nu\langle v_R v_z\rangle)}{\partial R} + \frac{\partial(\nu\langle v_z^2\rangle)}{\partial z} + \frac{\nu\langle v_R v_z\rangle}{R} = -\nu \frac{\partial \Phi}{\partial z}$$910% level.

7. Contemporary and Future Prospects

The JAM framework continues to underpin the majority of dynamical modeling for nearby galaxies and the Milky Way, especially as massive integral-field and spectrographic surveys (e.g., MaNGA, DESI, GECKOS) deliver high-quality resolved kinematics and proper motion data. The methodology is a benchmark for testing more flexible approaches (e.g., Schwarzschild orbit superposition) and plays a pivotal role in advancing constraints on black hole masses, dark matter halo shapes, and the distribution function of stars. The inclusion of multi-tracer and multi-moment data (proper motions, multiple populations) is establishing JAM as the practical standard for a new era of dynamical inference at galactic and extragalactic scales.