Closed-Form of the Local Galactic Potential and Stellar Distribution Function from Gaia DR3
arXiv:2609.09011v1 [astro-ph.GA] 8 Sep 2026
Indranil Das∗ Illinois Center for Advanced Studies of the Universe & Dept. of Physics University of Illinois Urbana-Champaign Urbana, IL 61801, USA [email protected]
Adam Kamoski∗ Department of Physics University of Massachusetts Boston NSF Institute for AI and Fundamental Interactions [email protected]
Dora Demiri European University of Tirana Tirana, Albania [email protected]
Brianna Isola University of New Hampshire Durham, NH 03824, USA [email protected]
Hanieh Karimi University of New Hampshire Durham, NH 03824, USA [email protected]@unh.edu
Dmitrii S. Zagorulia Lebedev Physical Institute Russian Academy of Sciences Moscow 119991, Russia [email protected]
Abstract The local dark matter density determines the strength of the signal expected in direct-detection experiments, yet published estimates from stellar motions disagree by more than their errors [1–6], and the most recent machine-learning analysis of Gaia data finds a local density consistent with zero [7]. According to Jeans’ theorem, a distribution function built from integrals of motion satisfies the collisionless Boltzmann equation (CBE) trivially for any choice of potential [8, 9], so a search that simultaneously fits the distribution function and the potential to the CBE identifies neither. Our pipeline instead estimates the distribution function in isolation, linearizing the equation in terms of accelerations and allowing for direct measurement of the local force field, and then fits closed forms to that field via symbolic regression. Throughout, we find that the usable information lies not in the CBE residual but in the stellar number counts, the observable most distorted by survey selection. Along the vertical profile, our recovered potential agrees with the classical self-gravitating isothermal disc.
1
Introduction
The Gaia DR3 catalog supplies six-dimensional phase-space coordinates for millions of stars [10], motivating machine-learning recovery of the Galactic potential Φ and stellar distribution function f from a single snapshot [7, 11, 12]. Symbolic regression [13, 14] provides an interpretable, differentiable output that can be checked against independent measurements. The local dark matter density inferred from such analyses is an input to direct-detection experiments, yet published values span 0.005–0.020 M⊙ pc−3 with error bars that do not overlap [1–6], and the most recent neural-network recovery reports a value consistent with zero [7]. The difficulty is identifiability: the equation relating ∗ Equal contribution.
observation to target is the stationary CBE, v · ∇x f − ∇x Φ · ∇v f = 0,
(1)
and when f and Φ are fit to it jointly, it fails to uniquely identify either; a model can match the data exactly without constraining Φ. Our pipeline avoids this by fixing f before solving for Φ, applied to 9.9 × 106 stars within 1 kpc. This work proceeds in three parts: In Sec. 2 we provide a quantified analysis of the degeneracy, which consistently puts more than 96% of the constraint on Φ in the spatial density. In Sec. 3 we directly measure the acceleration field, with no assumed potential, and we find closed forms for both unknowns under a Poisson-positivity constraint, with a controlled experiment showing the constraint is necessary. In Sec. 4 we recover a distribution function supervised so that the degenerate solution is no longer the optimum but the uninformative baseline.
2
Degeneracy of the potential function
Proposition 1 (Minimizing the residual cannot identify Φ). For any static and axisymmetric Φ, E = 12 |v|2 + Φ and Lz = Rvϕ are integrals of motion. Here, E is the specific orbital energy and Lz is the z-component of angular momentum. Then for any smooth h, the choice f = h(E, Lz ) satisfies Eq. (1) exactly. In particular f = const satisfies Eq. (1) for any Φ. The minimum of the residual is therefore zero on a set containing every potential. This follows from Jeans’ theorem [8, 9]. A residual-only symbolic search quickly converges to a simple but physically unreasonable zero-loss solution (complexity 2, both components constant). We therefore instead profile χ2 over the uniform volume density ρDM , which we scan on the interval ρDM ∈ [0, 20] × 10−3 M ⊙ pc−3 on mocks with known truth (Table A1). When the distribution function is allowed two free exponential components, the difference between exclusion and inclusion of the spatial density ν(z) is reflected in a χ2 span of 0.7 vs. one of 73. For the baseline mock the decomposition is CBE residual 0%, conditional velocity shape ∼1%, and spatial density ∼99%. ν(z) remains significant under every configuration tested, carrying more than 96% of the constraint in each: tracer temperature (12-60 km s−1 ), height range (0.5-2 kpc), and distribution-function family, including King-like truncated models (Table A1). This is because for an isothermal tracer ν depends 2 on Φ exponentially (ν ∝ e−Φ/σ ) [1, 15, 16], whereas the velocity shape enters only through moments that f can absorb. Since ν is the observable corrupted by survey completeness, the selection function must be modeled explicitly. When a completeness gradient S ∝ e−|z|/ℓ is present in a mock and omitted from the model, the recovered surface density shifts by σ 2 /2πGℓ, at a loss floor indistinguishable from the clean run with S = 1 (App. A). Proposition 2 (Two-integral models are meridionally isotropic). For static axisymmetric Φ, every f (E, Lz ) satisfies σR = σz exactly at every point. 2 Proof. E depends on vR and vz only through vR + vz2 , and Lz not at all, so both are invariant under 2 the exchange vR ↔ vz . Therefore, so is f , and ⟨vR ⟩ = ⟨vz2 ⟩ at fixed (R, z).
The two-integral family that minimizes the residual is the same family constrained to σR = σz . Hence, a measured anisotropy σR ̸= σz distinguishes our solution from it. This does not distinguish it from every stationary solution, since a three-integral f (E, Lz , I3 ) satisfies Eq. (1) identically and is generically anisotropic. An exactly separable third integral requires Φ of Stäckel form, which we neither impose nor test for (App. C). Therefore, we expect a non-zero residual.
3
Method
Selection. The 1 kpc sphere holds 9.94 × 106 stars with RUWE< 1.4, parallax S/N > 10, ≥ 4 radial-velocity transits, and velocity uncertainty < 10 km s−1 , split 70/15/15 before any fitting (in-sample selection picks the wrong front member on mocks; App. C). Gradients of ln f are evaluated on a 2.5 × 106 -star subsample; 1.73 × 106 of those fall in the 94 acceleration cells, and the symbolic search for ln f is fitted on 2.0 × 104 rows drawn from the 2.21 × 106 stars within 0.9 kpc. We use R0 = 8.122 kpc and z⊙ = 0, the latter by construction of the catalogue frame (App. B). The observed number density per unit volume falls to 5.1% of its value at 0.15 kpc 2
vc = 220 km/s vc = 240 km/s
0.10
−0.60
R = 7.55
R = 8.32
R = 7.74
R = 8.51
R = 7.93
R = 8.70
0.03
R = 8.12
0.02
−0.70
aϕ [cm s−1 yr−1]
az [cm s−1 yr−1]
aR [cm s−1 yr−1]
0.05 −0.65
0.00 −0.05
0.01
0.00
−0.10 −0.75 −0.15
7.6
7.8
8.0
8.2
8.4
8.6
−0.01
−0.6
−0.4
−0.2
R [kpc]
0.0
z [kpc]
0.2
0.4
0.6
−0.6
−0.4
−0.2
0.0
0.2
0.4
0.6
z [kpc]
Figure 1: The acceleration field a, measured cell by cell with no functional form assumed for Φ. Each point is one spatial cell solved from thousands of linear equations, Eq. (2); error bars on the rightmost plot are the cell-fit standard errors. Left: the radial component ar at the midplane, with constant-vc curves for scale. Center: the vertical component az , colored by radius, changing sign at the dynamical midplane. Right: the azimuthal component aϕ , which an axisymmetric Φ requires to vanish and which instead comes out systematically positive at the per-cent level of aR . across the sphere, with |d ln S/dd| peaking at 10.9 kpc−1 against a physical |∂z ln ν| of 3.0 kpc−1 (Fig. A2). We fit µ = V ·S(d)A(ℓ, b)ν(R, z) as a Poisson GLM. The three factors are identifiable up to two overall multiplicative constants because they are different functions of the same point: stars at equal distance in different directions lie at different heights. It returns a 2.88 kpc radial scale length and a two-scale-height vertical profile whose local scale height at 300 pc is 334 pc, neither constrained to a literature value. The complete pipeline and code are available at https: //github.com/briannaisola/GaiaSR. Empirical f . We factorize f = ν p(v|x) and fit the conditional with a 96-component full-covariance Gaussian mixture, whose log-gradient is closed-form and agrees withR central finite differences to 10−6 R 3 3 in all six coordinates. The conditional is selection-free: S(x)f / S(x)f d v = f / f d v, since S has no velocity dependence, so components are selected on held-out conditional log-likelihood. Refitting on a disjoint set of 5 × 105 stars moves ∂R ln f by 0.815 against its own root mean square (rms) of 1.673; refitting with Gaia’s errors injected a second time moves it by 0.736. The two shifts are comparable in size, so the limiting factor is the density estimator, and the implied residual floor is 48.7 km s−1 kpc−1 for any expression carrying its gradients. The azimuthal streaming term carries 93% of the total rms (222.2 against 67.4 and 41.7 km s−1 kpc−1 ), is measured over only a 14◦ baseline in azimuth, and cannot be balanced by any axisymmetric force; dropping it moves dvc /dR from −9.0 to −1.1 against a literature −1.7 ± 0.1 [17]. Acceleration Measurement. Eq. (1) is linear in a. In cylindrical coordinates with g = ln f , grouped by whether a term involves the potential, v2
v v
vR ∂R g + vz ∂z g + Rϕ ∂vR g − RR ϕ ∂vϕ g + Wi ·a = 0, (2) | {z } Ki with Wi = ∇v g and the azimuthal streaming term dropped as above. The centrifugal and Coriolis terms arise from the rotation of the cylindrical basis as a star moves, and are fixed by the star’s own coordinates. We solve 94 acceleration cells (7 radial bins × 14 vertical bins −4 star-count-cut cells) by Huber-reweighted least squares, assuming no functional form for Φ (Fig. 1). The field comes out at ∼ 0.6 cm s−1 yr−1 , the acceleration undetectable in any single star over the√mission lifetime but recoverable from the collective statistics of many. From it follow vc = −aR R, Σ(< |z|) = |az |/2πG, and the dynamical midplane where az changes sign. The whole volume as one cell gives vc = 231.4 km s−1 against 229 ± 3 [17]; the cells give Σ(< 0.5 kpc) = 44.0 M⊙ pc−2 , between the 41 and 65 ± 6 measured at 0.35 and 0.8 kpc [18], and place the midplane 18 pc from the Sun against a photometric 20.8 ± 0.3 pc [19]. The quoted errors are least-squares statistical errors on millions of stars, and they understate the real uncertainty: nothing below about a per cent here is resolved. Solving for aϕ rather than imposing aϕ = 0 bounds the departure from axisymmetry at 1.3% of |aR |, comparable to that floor, and we return to it in the limitations. Symbolic Φ. We search with PySR [14, 20]. Each cell enters the training set twice, tagged by an indicator, so a single Φ must reproduce both force components and the fitted field is curl-free by 3
0.15
R = 8.32
R = 7.74
R = 8.51
R = 7.93
R = 8.70
symbolic Φ
−0.62
vc = 229 km/s measured
R = 8.12
−0.64
0.05
aR [cm s−1 yr−1]
az [cm s−1 yr−1]
0.10
R = 7.55
0.00 −0.05
−0.66
−0.68
−0.70
−0.10
−0.72
−0.15
−0.6
−0.4
−0.2
0.0
0.2
0.4
0.6
7.6
7.8
8.0
z [kpc]
8.2
8.4
8.6
R [kpc]
Figure 2: Two plots representing the closed-form potential Φ(R, z) of Eq. (3), containing eleven nodes and two fitted constants, against the independently measured field of Fig. 1. Left: vertical acceleration az , points measured and lines the symbolic Φ, colored by radius R. Right: radial acceleration aR at the midplane, with a constant-vc curve at the literature value for comparison. The fit is to the measured field alone. construction. 5.24 ln R ± 8z 2 (∇2 Φ = ±16) generate equally valid acceleration fields a that the data cannot tell apart, and only the sign of the Laplacian distinguishes them. Hence, a Poisson positivity constraint becomes necessary, entering the objective as LΦ = σ −2 (pred − atgt )2 + λ[min(∇2 Φ, 0)]2 with λ = 25, one-sided so that it vanishes on the admissible set and cannot bias the choice among physical candidates. From the Pareto front we take the cheapest expression within a factor of two of the best admissible loss, subject to (1) ρ > −0.005 throughout the (R, z) domain where Poisson positivity is enforced and (2) a plausibility window ρ(R0 , 0) ∈ [0, 0.5] M⊙ pc−3 , vc ∈ [150, 320] km s−1 . This window is wide enough to exclude only nonsense: all 21 of the 26 front rows surviving Poisson positivity already lie inside it, and removing it returns the same expression. Selection is therefore by positivity and parsimony, and the window did not manufacture the vc and ρ(R0 , 0) reported below. Closed-form f . By Prop. 1, the residual cannot serve as the objective. We instead regress g onto P the empirical ln fˆ and its five gradients. The loss is Lf = ⟨((g − ln fˆ)/Sg )2 ⟩ + 15 i ⟨((∂i g − ∂i ln fˆ)/Si )2 ⟩, averaged over stars, where each S is the rms of the target it normalizes. The gradients are supervised because Eq. (2) absorbs them, and because an expression can track a function’s values closely while its derivatives wander. A constant g scores exactly 2.000, one from the value term and one from the five gradient terms, whose model gradients vanish. An expression carrying no information about the target can get a score as low as 2.000 under this loss. The same expression scores 0 under a residual objective, where it wins. Selection is on R2 against ln fˆ, taken over Pareto front rows that are finite on more than 99% of held-out stars. The residual and the anisotropy are computed afterwards and enter no decision. Two of the five supervised gradients are ∂vR ln fˆ and ∂vz ln fˆ, which carry σR and σz for a locally Gaussian conditional. The anisotropy reported below is therefore derived from supervised quantities. Front members of comparable loss return ratios between 0.91 and 9.90 (Table 1), so reproducing the supervised gradients pointwise does not by itself reproduce the ellipsoid they imply.
4
Results
The potential comes out at complexity 11 as Φ(R, z) = 5.3618 ln R + 0.39193 ln cosh z ,
(3)
with R, z in kpc and Φ in (100 km s−1 )2 . The recovered expression contains a ln cosh z vertical dependence characteristic of the potential of a self-gravitating isothermal sheet, whose density profile ν ∝ sech2 (z/2h) was derived by Spitzer [21]. Notably, this functional form emerged from the free operator set rather than being imposed a priori. 4
Table 1: Pareto front for ln f , held out; every fifth Table 2: Held-out results. σR /σz is measured member plus the two endpoints, from 37 (full from the same sample, not literature; the rest are front in App. C). Row two is g = const, the exact independent [1, 17–19, 22]. Two of the five exglobal optimum of a residual objective, present ternal comparisons fail, marked †. Uncertainties on our own front and rejected on R2 . ϵ∇ is the are discussed in the text and are not the statistical mean relative gradient error. errors. ∗: measured from the same sample. ‡: the supervising mixture on the same axisymmetric c loss R2 ϵ∇ CBE σR /σz target. CBE residual in km s−1 kpc−1 . 1 2.016 −8×10−3 1.00 6.2 2 2.000 −3×10−5 1.00 0.00 9 1.352 0.510 0.90 41.7 10 1.313 −2×109 3×104 1×108 20 0.879 0.745 0.76 51.9 28 0.691 0.876 0.74 41.7 39 0.603 0.917 0.70 36.1 44 0.589 0.923 0.69 38.2
1.14 0.91 8.11 9.90 1.37 1.87 2.13 1.88
Quantity
This work
Reference
σR /σz 1.83–2.13 1.893∗ −1 vc (R0 ) 231.6 km s 229 ± 3 Σ(< 0.5 kpc) 44.0 41–65 (0.35-0.8 kpc) midplane offset 18 pc 20.8 ± 0.3 dvc /dR† −1.1 −1.7 ± 0.1 ρ(R0 , 0)† 0.048 0.084–0.10 CBE resid. rel. to Φ=0 R2 (ln f )
38.2 0.204 0.923
58.6‡ 1.000 —
The distribution function comes out at complexity 44 as ln f = 2.407 − ln cosh(8.523 vz ) − ln cosh z ln cosh(vϕ ln R) − ln cosh (ln cosh vϕ − 1.625)(ln cosh(5.251 vz ) − R + 1.385) − ln cosh ln cosh(2.109 vR vϕ ) .
(4)
The selected expression is non-separable in the velocity components: vR enters through the product vR vϕ , while the large-|vz | behavior is approximately exponential. The recovered field reproduces three of the five external comparisons in Table 2. Table 1 shows that the complexity-2 constant solution achieves zero CBE residual but no meaningful fit to the empirical distribution (R2 = −3 × 10−5 ), demonstrating the failure of residual-only model selection. R2 then rises with complexity, though not monotonically along the whole front, and the meridional anisotropy settles into the range 1.83–2.13 over the top twelve members, mean 1.94, bracketing the measured 1.893. The maximum-R2 member reads 1.88, but we quote the spread, since no member reaches the observed ratio to better than the 5% scatter of the front. The recovered Φ and ln f predict this ellipsoid because two of the supervised gradients carry it (Sec. 3). That they come out at the right value shows the fit preserved the velocity shape it was given. Φ = 0 gives exactly 1.000, so 0.204 means the potential terms cancel most of the variance of the remaining, axisymmetric part of Eq. (2). This residual is evaluated with the azimuthal streaming term vϕ R ∂ϕ g removed, which no axisymmetric force can balance and which alone carries 93% of the signal (Sec. 3). With it retained the relative residual is 0.76. The value 38.2 lies below the 58.6 scored on the same axisymmetric target by the mixture that supervised it, because a smooth expression cannot reproduce estimator noise. The Φ residual is constant to 0.5% across the train, validation, and test splits, from two fitted constants against 1.7 × 106 evaluation stars, though the symbolic search itself saw only 2 × 104 . Ablations (App. C). Separable potentials Φ = ΦR (R) + Φz (z) occupy complexity 9–10 and are 43% worse than the coupled form, which predicts az falling 22% across the sample where a separable form predicts no change; since Ez is conserved only if ∂ 2 Φ/∂R ∂z = 0, the recovered coupling is also why Ez is not available as a√third integral here. Among vertical profiles fitted freely, ln cosh (χ2 /dof = 275.4) is matched by z 2 + h2 − h (275.8) and beaten by nothing, while any form with a midplane kink is 10× worse. Rescaling inputs by measured scales moves the residual to 23.69 and R2 to 0.937 at identical complexity. The observed anisotropy is stable to ±1% across apertures from |z| < 0.05 to |z| < 0.30 kpc. Where it departs. The local density is ρ = 0.048 against a literature 0.084–0.10 M⊙ pc−3 [1, 22]. Evaluating ρ = −(4πG)−1 [∂R aR +aR /R+∂z az ] on the measured field, with no symbolic expression involved, gives 0.0493: the fit reproduces the field to 3%, so the deficit is already present in the accelerations and enters upstream through the density estimator, exactly where the mocks place it, 5
a smoothed ν moving the recovered split from its clean value (Σ, ρ) = (50, 0.010) to (66, 0.003) because Σ integrates the vertical force while ρ differentiates it (App. A).
5
Conclusion
We present closed-form models of the local Galactic potential and stellar distribution function inferred from the Gaia DR3 catalog. Controlled mock tests show that stellar number counts provide the main constraint on the potential, whereas jointly fitting f and Φ to the stationary collisionless Boltzmann equation alone remains degenerate. From the completeness-corrected distribution function, we reconstruct the acceleration field to which the symbolic potential is fitted under a Poisson-positivity constraint. The recovered potential contains a ln cosh term characteristic of a self-gravitating isothermal sheet. Three of five external benchmarks are broadly reproduced, and the top Pareto-front models bracket the observed velocity anisotropy. Steady state and axisymmetry are assumed within 1 kpc, despite a non-zero measured azimuthal acceleration. The ln cosh term has a fixed vertical scale of 2h = 1 kpc. Although the surface density is approximately recovered, the local density is underestimated by about a factor of two, a deficit already present in the empirical acceleration field. The selection model has an 11% systematic in ∂z ln ν that dominates the uncertainty in Φ.
Acknowledgements We thank Akshay Ghalsasi for discussions, and Miles Cranmer and Jose M. Munoz for posing the problem. We thank IAIFI for inspiration and for facilitating the work through computing resources provided during during and after the Hackathon where the problem was posed. We thank NCSA and ACCESS/PSC for computing resources. This work has made use of data from the European Space Agency (ESA) mission Gaia (https: //www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
References [1] J. I. Read. The local dark matter density. J. Phys. G, 41:063101, 2014. [2] K. Schutz, T. Lin, B. R. Safdi, and C.-L. Wu. Constraining a thin dark matter disk with Gaia. PRL, 121: 081101, 2018. [3] A. Widmark and G. Monari. The dynamical matter density in the solar neighbourhood inferred from Gaia DR1. MNRAS, 482:262–277, 2019. [4] R. Guo, C. Liu, S. Mao, X.-X. Xue, R. J. Long, and L. Zhang. Measuring the local dark matter density with LAMOST DR5 and Gaia DR2. MNRAS, 495:4828–4844, 2020. [5] J.-B. Salomon, O. Bienaymé, C. Reylé, A. C. Robin, and B. Famaey. Kinematics and dynamics of Gaia red clump stars: revisiting north-south asymmetries and dark matter density at large heights. A&A, 643: A75, 2020. [6] A. Widmark, C. F. P. Laporte, P. F. de Salas, and G. Monari. Weighing the Galactic disk using phase-space spirals II. most stringent constraints on a thin dark disk using Gaia EDR3. A&A, 653:A86, 2021. [7] T. Kalda and G. M. Green. Deep Potential: Recovering the gravitational potential and local pattern speed in the solar neighborhood with GDR3 using normalizing flows. arXiv:2507.03742, 2025. [8] J. H. Jeans. On the theory of star-streaming and the structure of the universe. MNRAS, 76:70–84, 1915. [9] J. Binney and S. Tremaine. Galactic Dynamics. Princeton University Press, 2nd edition, 2008. [10] Gaia Collaboration, A. Vallenari, et al. Gaia Data Release 3: Summary of the content and survey properties. A&A, 674:A1, 2023. [11] G. M. Green, Y.-S. Ting, and H. Kamdar. Deep Potential: Recovering the gravitational potential from a snapshot of phase space. ApJ, 942:26, 2023. [12] J. An, A. P. Naik, N. W. Evans, and C. Burrage. Charting galactic accelerations: when and how to extract a unique potential from the distribution function. MNRAS, 506:5721–5730, 2021. [13] M. Cranmer, A. Sanchez-Gonzalez, P. Battaglia, R. Xu, K. Cranmer, D. Spergel, and S. Ho. Discovering symbolic models from deep learning with inductive biases. In NeurIPS, 2020. [14] M. Cranmer. Interpretable machine learning for science with PySR and SymbolicRegression.jl. arXiv:2305.01582, 2023. [15] K. Kuijken and G. Gilmore. The mass distribution in the galactic disc I. MNRAS, 239:571–603, 1989.
6
[16] J. Bovy and H.-W. Rix. A direct dynamical measurement of the Milky Way’s disk surface density profile. ApJ, 779:115, 2013. [17] A.-C. Eilers, D. W. Hogg, H.-W. Rix, and M. K. Ness. The circular velocity curve of the Milky Way from 5 to 25 kpc. ApJ, 871:120, 2019. [18] J. Holmberg and C. Flynn. The local surface density of disc matter mapped by Hipparcos. MNRAS, 352: 440–446, 2004. [19] M. Bennett and J. Bovy. Vertical waves in the solar neighbourhood in Gaia DR2. MNRAS, 482:1417–1425, 2019. [20] M. Cranmer. PySR: Template expressions and differential operators. https://ai.damtp.cam.ac.uk/ pysr/examples/, 2025. [21] L. Spitzer. The dynamics of the interstellar medium. III. Galactic distribution. ApJ, 95:329–344, 1942. [22] C. F. McKee, A. Parravano, and D. J. Hollenbach. Stars, gas, and dark matter in the solar neighborhood. ApJ, 814:13, 2015.
7
A
Mock anatomy of the degeneracy
Setup. Our mock inputs are Φ(z) = 2πGΣ(2h) ln cosh(z/2h) + 2πGρDM z 2 , Σ = 48 M⊙ pc−2 , h = 0.20 kpc, and ρDM = 0.010 M⊙ pc−3 . When run clean without the injected completeness gradient, we recover (Σ, ρDM ) = (50, 0.010) from this mock. The true distribution function is 2 2 F (E) = 0.7 e−E/σ1 + 0.3 e−E/σ2 , with (σ1 , σ2 ) = (18, 40) km s−1 and E = v 2 /2 + Φ. P 2 Profiling procedure. For each trial ρDM we fix Φ and optimize F (as i ai e−E/σi , ai ≥ 0, P ai = 1) to minimize χ2 against the conditional dispersion and kurtosis (σz , κ) in 10 height bins over 0.05–0.95 kpc, with assigned errors of 2% and 3%; when the spatial density is included the profile scores additionally against ν(z) normalized at z = 0 with 2% errors. Eight random restarts ensure convergence. These ∆χ2 spans are computed at fixed assumed errors and are heuristic; the percentages depend on the relative error scaling, though the ordering ν ≫ P (v|z) is insensitive to it over the tested 1–5% range. Table A1 gives the full sweep behind the decomposition of Sec. 2.
Table A1: The sensitivity hierarchy across mock configurations: ∆χ2 span over ρDM from the conditional velocity shape alone against the shape plus spatial density, with the distribution function free (three exponential components unless noted). Spans are calculated over the interval ρDM ∈ [0, 20] × 10−3 M ⊙ pc−3 . The ν(z) share is 1 − spanv /span+ν . Each block varies one axis about the baseline and the two sweeps were run independently, so the baseline row differs between them. The ν(z) share stays above 96% throughout and above 98% for every aperture |z|max ≥ 0.75 kpc. Span P (v|z)
Span + ν(z)
ν(z) share
Tracer temperature (σ1 , σ2 ) = (12, 30) km s−1 (σ1 , σ2 ) = (18, 40) [baseline] (σ1 , σ2 ) = (25, 50) (σ1 , σ2 ) = (35, 60) isothermal (σ = 15)
0.3 0.6 0.5 0.2 0.0
91 57 29 16 281
99.7% 98.9% 98.3% 98.8% 100%
Height range |z|max = 0.5 kpc |z|max = 0.75 kpc |z|max = 1.0 kpc [baseline] |z|max = 1.5 kpc |z|max = 2.0 kpc
0.1 0.3 0.7 0.9 0.7
3.0 20 73 209 438
96.7% 98.5% 99.0% 99.6% 99.8%
Distribution-function family F with 2 free components F with 5 free components King-like (truncated at 80 km s−1 )
0.7 0.4 0.0
73 73 303
99.0% 99.5% 100%
Configuration
An injected completeness gradient. Drawing the baseline mock through S ∝ e−|z|/ℓ with ℓ = 0.6 kpc and fitting with S ignored, the recovered potential converges on Φeff = Φ + σ 2 |z|/ℓ: an isothermal tracer cannot distinguish an unmodeled completeness gradient from a thin sheet of surface density ∆Σ = σ 2 /2πGℓ = 25 M⊙ pc−2 . The recovered surface density reads 74 against the predicted 75, the clean recovery plus ∆Σ; the scale height collapses to ∼120 pc; and ρDM scatters over 0.003–0.027 with complexity and fitting range where the clean recovery is stable at 0.010. The loss floor is indistinguishable from the clean run’s, so no goodness-of-fit or Pareto criterion detects the bias. This is why Sec. 3 models S explicitly. A smoothed density estimate. Estimating ν with a Gaussian-smoothed histogram (1.2 cells) deforms ln f by only 4% at the midplane, where its curvature peaks, yet pulls the recovered disc–halo split from (Σ, ρ) = (50, 0.010) to (66, 0.003); the smoothed surface is fit better than the truth fits it. Σ integrates the vertical force while ρ differentiates it, so smoothing ν moves density out of the halo term and into the sheet. That is the ρ(R0 , 0) deficit of Sec. 4, reproduced here where the truth is known. 8
F fixed at truth F free, P(vz|z) only F free, P(vz|z) + (z)
60 2( DM)
[%]
(a)
span = 107.3
40 20 0
span = 0.8
0
5
10
DM [10 3 M
15
pc 3]
20
25
share of constraint on
80
(b) 100
99%
80 60 40 20 0
0%
1%
CBE (z) P(vz|z) residual (vel. shape) (star counts)
Figure A1: Spans are calculated over the interval ρDM ∈ [0, 25] × 10−3 M ⊙ pc−3 , accounting for their slight deviation from corresponding spans reported in Table A1. Freeing the distribution function destroys the constraint on ρDM ; adding the star counts restores it. Left: ∆χ2 over ρDM for the baseline mock, with the distribution function fixed at truth (gray), free and scored on the conditional velocity shape alone (red), and free with the spatial density added (blue). Right: the same three profiles as fractions of the total. |b| < 15 ∘
volume complete
10
normalised number density
observed n(d) / n(0.15 kpc)
15 ∘ < |b| < 45 ∘ 0
10−1
0.0
0.2
0.4
0.6
0.8
|b| > 45 ∘ 100
10−1
1.0
0.0
heliocentric distance d [kpc]
0.2
0.4
0.6
0.8
1.0
heliocentric distance d [kpc]
Figure A2: Survey completeness dominates the observed density. Left: observed number density per unit volume against heliocentric distance, normalized at 0.15 kpc. It falls by more than an order of magnitude across the ball, driven by a selection gradient |d ln S/dd| of 10.9 kpc−1 against a physical |∂z ln ν| of 3.0 kpc−1 . Right: the same profile split by galactic latitude. The three lines separate because the distance falloff mixes the isotropic selection S(d) with the height-dependent stratification ν(R, z); the GLM separates them because stars at equal distance in different directions lie at different heights.
B
Data and the selection model
We adopt R0 = 8.122 kpc and z⊙ = 0: the catalogue is built in a frame centred on the Sun, so the dynamical midplane offset of Sec. 3 is measured against a fixed z = 0. The GLM µ = V ·S(d) A(ℓ, b) ν(R, z) fits the three factors jointly. The density factor returns a 334 pc tracer scale height and a 2.88 kpc scale length, neither tied to literature values. Dust is not separable in distance and direction, so we refit on the dust-poor |b| > 20◦ sight lines: ∂z ln ν(0.3 kpc) moves from −2.99 to −2.67 kpc−1 , an 11% systematic carried forward, and by the hierarchy of Sec. 2 that systematic is what limits Φ.
C
Ablations and model-selection protocol
Separable potentials. Constraining Φ = ΦR (R) + Φz (z) costs 43% in loss at complexity 9–10 against the coupled form of Eq. (3). The coupling accounts for the 22% decline of az across the sampled radial range, which a separable form cannot produce. The same cross term removes Ez as a candidate third integral, since Ez is conserved only when ∂ 2 Φ/∂R∂z = 0; this excludes the 9
−208
stars mixture
−210
150 55
28
50
26
100
−218
45
⟨vRvz⟩ [km2/s2]
−216
σz [km/s]
−214
σR [km/s]
⟨vϕ⟩ [km/s]
−212
24
22
−220 40
50 0 −50 −100 −150
−222 20
−224
−200
35 −0.50 −0.25
0.00
z [kpc]
0.25
0.50
−0.50 −0.25
0.00
0.25
0.50
−0.50 −0.25
z [kpc]
0.00
z [kpc]
0.25
0.50
−0.50 −0.25
0.00
0.25
0.50
z [kpc]
Figure A3: The 96-component mixture against conditional velocity moments it was never shown: held-out stars (points) versus the mixture (line), as functions of height. The mixture is fitted to individual stars, so the velocity moments serve as an independent check. The mixture reproduces them near the midplane and drifts from the dispersions beyond |z| ∼ 0.4 kpc, one source of the estimator floor discussed in Sec. 3. separable case but does not exclude the Stäckel family as a whole, so a non-zero CBE residual is expected here (Sec. 2). Vertical profile family. Fitting the vertical √ term freely among candidate profiles: ln cosh reaches χ2 /dof = 275.4, the softened modulus z 2 + h2 − h ties at 275.8, and every form with a midplane kink (e.g. |z|) is an order of magnitude worse. The data favor the isothermal-sheet family but do not distinguish its two smooth parameterizations, which agree to the width of the fitted region. Input scaling. Rescaling the inputs by the measured scale height and scale length moves the residual from 38.15 to 23.69 km s−1 kpc−1 and R2 from 0.923 to 0.937 at identical complexity: unit choices are worth as much as ∼10 complexity points, so all searches run in scaled variables. Aperture stability. The measured anisotropy σR /σz is stable to ±1% across apertures from |z| < 0.05 to |z| < 0.30 kpc, so its value is not an artifact of the aperture; the spread we report in Sec. 4 comes from the choice of front member. Why selection is held out. The 70/15/15 split of Sec. 3 is fixed before fitting because in-sample selection picks the wrong front member. On mocks with known truth, in-sample χ2 prefers a spuriously curved member over the straight one that is true. Held-out χ2 reverses the choice at every complexity (in-sample χ2 /N of 1.21 against held-out 1.53 for the line, 0.89 against 1.71 for the exponential). An exponential mimics a line over a bounded range, and only held-out scoring resists it. Every selection in this work follows the rule: mixture components on held-out conditional likelihood, ln f on held-out R2 , Φ on the screened front. Compute. The acceleration solve is linear least squares over 94 cells and runs in minutes on a single CPU node. Each PySR search (potential and distribution function) runs on one multi-core CPU node in a few hours; the full set of ablations reported here is ∼103 CPU-core-hours.
10