A general spectral solver for the axisymmetric Jeans equations: fast galaxy modelling with arbitrary anisotropy

TL;DR

This paper addresses the computational bottleneck in axisymmetric Jeans modelling by introducing a general spectral solver that can handle arbitrary anisotropy $\beta(r,\theta)$. The method solves directly for the radial velocity dispersion $\bar v_r^2$ on a logarithmic grid using Chebyshev collocation, and enforces asymptotic behavior via a Robin boundary condition, enabling exponential convergence with modest grid sizes. It recovers intrinsic moments with sub-percent accuracy and is orders of magnitude faster than traditional high-accuracy quadratures, with implementation designed for GPUs and integration into the JamPy package. Validations against analytic isotropic solutions and comparisons with JAM_sph show the solver's accuracy across regimes, including high central cusps and varying anisotropy, making it a robust, scalable tool for exploring large parameter spaces in galaxy dynamics.

Abstract

Dynamical modelling is a fundamental tool for measuring galaxy masses and density profiles in the era of large integral-field spectroscopic surveys and Bayesian inference. Solutions based on the Jeans equations are popular due to their robustness and computational efficiency. However, traditional semi-analytic Jeans solvers often require restrictive assumptions about the velocity anisotropy to remain computationally tractable. This paper presents a new spectral solver for the axisymmetric Jeans equations designed to overcome these limitations. I first illustrate, using orbit integrations in realistic potentials, that spherical alignment of the velocity ellipsoid is a physically well-motivated approximation for galaxy modelling. The new method employs a spectral technique to solve the Jeans partial differential equations directly. Two design choices are critical for accuracy and speed: (i) solving for the slowly-varying velocity dispersion rather than the rapidly varying pressure, and (ii) imposing a Robin boundary condition to enforce the asymptotic decay on a finite domain. This formulation supports arbitrary anisotropy distributions beta(r, theta) while simultaneously increasing computational speed by orders of magnitude compared to standard high-accuracy quadratures. Validated against exact analytic benchmarks, the solver recovers intrinsic moments with sub-percent accuracy. The implementation will be included in the public JamPy package and is structured to be optimally suited for massive parallelization on specialized hardware such as GPUs, enabling the rigorous exploration of complex parameter spaces.

Figures (7)

• Figure 1: Comparison of orbital envelopes in a separable potential versus a realistic galaxy potential. Left: Orbits in a Perfect Ellipsoid potential deZeeuw1985. The orbital boundaries rigidly follow prolate spheroidal coordinates (dashed lines), forcing the velocity ellipsoid to become cylindrically aligned in the core. Right: Orbits in a realistic multi-component potential (Sersic1968$n=4$ bulge + NFW halo Navarro1996nfw + Supermassive Black Hole). While the orbits do not strictly follow a single coordinate system, their envelopes are remarkably well approximated by spherical coordinates (dashed lines) at all radii, from the BH-dominated centre ($0.1 R_{\rm e}$) to the halo-dominated outskirts ($3 R_e$). This justifies the use of spherically aligned velocity ellipsoids for modelling real galaxies.
• Figure 2: Performance of the spectral JAM solver against the exact analytic isotropic Satoh model ($\beta=0$). Top panels: Maps of the intrinsic radial and azimuthal velocity dispersions, $(\overline{v_r^2})^{1/2}$ and $(\overline{v_\phi^2})^{1/2}$, computed on a $32 \times 24$ grid. Bottom panels: Relative percentage error between the spectral and analytic solutions. The residuals are essentially structureless and remain negligible ($<0.2$ per cent) throughout the meridional plane.
• Figure 3: Numerical convergence of the spectral solver as a function of grid resolution. Imposing a Dirichlet boundary condition ($\overline{v_r^2}=0$) at a finite radius $r_\mathrm{max}=3\max\{\sigma_k\}$ causes the error to saturate (red solid line with squares). The adoption of a Keplerian Robin boundary condition (blue dashed line with circles) restores exponential convergence, allowing the error to decrease toward machine precision. An exponential fit (green dotted line) to the data for $N_\theta \ge 7$ demonstrates the rapid convergence rate. Remarkably, the method achieves $\sim 1\%$ accuracy with a grid as coarse as $14 \times 7$ points.
• Figure 4: Validation of the spectral solver using a realistic composite galaxy model (Sérsic $n=4$ stellar bulge and a cuspy gNFW dark matter halo with inner slope $\gamma=-2.2$). Top panels: Intrinsic radial and azimuthal velocity dispersion maps $(\overline{v_r^2})^{1/2}$ and $(\overline{v_\phi^2})^{1/2}$. Bottom panels: Relative percentage error between the spectral solution and the exact analytic MGE solution derived in \ref{['sec:analytic_solution']}. Despite the high dynamic range and central cusp, the maximum error remains remarkably low ($<0.2$ per cent).
• Figure 5: Numerical convergence analysis for the Sérsic + gNFW model. As observed in the Satoh test, the Dirichlet boundary condition (red squares) suffers from truncation errors that cause the accuracy to plateau. The Keplerian Robin boundary condition (blue circles) restores exponential convergence, reaching a precision of $10^{-4}$. The slightly shallower slope relative to the Satoh case (\ref{['fig:satoh_convergence']}) is expected given the steep density gradients inherent to the $n=4$ profile and the gNFW cusp.