Source-linked AI summary

Closed-Form of the Local Galactic Potential and Stellar Distribution Function from Gaia DR3

Indranil Das, Adam Kamoski, Dora Demiri, Brianna Isola, Hanieh Karimi, Dmitrii S. Zagorulia

arXiv:2609.09011v1astro-ph.GAastro-ph.IMcs.LG

TL;DR

Published local dark matter density estimates disagree, despite the density’s importance for direct-detection experiments. The paper separates distribution-function estimation from potential fitting, reconstructs accelerations, and fits closed forms; it finds that spatial density supplies most of the potential constraint and recovers a vertical ln cosh form characteristic of an isothermal sheet.

  • Problem

    Published local dark matter density estimates disagree beyond their error bars, and the most recent Gaia neural-network analysis reports a value consistent with zero, despite the density’s importance for direct-detection signals.

    Source: §1 · p. 1
  • Method

    The pipeline estimates f before fitting Φ, reconstructs the acceleration field, and fits a symbolic potential under a Poisson-positivity constraint.

    Source: §1 · p. 2, §3 · p. 4
  • Results

    More than 96% of the constraint on Φ comes from spatial density, and the recovered potential contains a ln cosh z term characteristic of a self-gravitating isothermal sheet.

    Source: §1 · p. 2, §4 · p. 4
  • Takeaways & Limitations

    The paper concludes that stellar number counts provide the main constraint on the potential and that survey completeness must be modeled explicitly because it distorts those counts.

    Source: §5 · p. 6

    The azimuthal streaming term is measured over only a 14° azimuthal baseline and cannot be balanced by an axisymmetric force.

    Source: §3 · p. 3

Abstract

from arXiv · show

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, and the most recent machine-learning analysis of Gaia data finds a local density consistent with zero. According to Jeans' theorem, a distribution function built from integrals of motion satisfies the collisionless Boltzmann equation (CBE) trivially for any choice of potential, 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

Gaia DR3 enables recovery of the Galactic potential and stellar distribution function from stellar phase-space data, but local dark matter estimates disagree and joint fitting to the stationary CBE is non-identifying. The paper addresses this degeneracy by separating distribution-function estimation from potential fitting and examining where the potential’s constraint comes from.

1 Introduction

  • Published local dark matter densities span 0.005–0.020 M⊙pc−3 with non-overlapping error bars, while the most recent neural-network estimate is consistent with zero.The local density determines the expected signal strength in direct-detection experiments.Source: §1 · p. 1
  • Gaia DR3 provides six-dimensional phase-space coordinates for millions of stars, motivating machine-learning recovery of Φ and f from a single snapshot.The paper identifies symbolic regression as an interpretable, differentiable output that can be checked against independent measurements.Source: §1 · p. 1
  • Jointly fitting f and Φ to the stationary CBE does not uniquely identify either, allowing a model to match the data without constraining Φ.The pipeline instead fixes f before solving for Φ and applies it to 9.9 × 10^6 stars within 1 kpc.Source: §1 · p. 2
  • More than 96% of the constraint on Φ comes from spatial density, according to the paper’s quantified analysis of the degeneracy.The work also measures the acceleration field without assuming a potential, fits closed forms under a Poisson-positivity constraint, and supervises distribution-function recovery.Source: §1 · p. 2

2 Degeneracy of the potential function

For any static, axisymmetric potential, distribution functions built from energy and angular momentum can make the stationary CBE residual vanish, so residual minimization alone cannot identify the potential. The paper’s mock tests find that spatial density carries nearly all the constraint, while two-integral models have a specific anisotropy limitation.

2 Degeneracy of the potential function

  • For any static, axisymmetric Φ, f = h(E, Lz) satisfies the CBE exactly, so the residual minimum is zero for a set containing every potential.This includes the constant distribution function.Source: §2 · p. 2
  • Spatial density ν(z) carries more than 96% of the potential constraint across tested tracer temperatures, height ranges, and distribution-function families.For an isothermal tracer, ν depends exponentially on Φ; because survey completeness distorts ν, the selection function must be modeled explicitly.Source: §2 · p. 2
  • On the baseline mock, the constraint decomposes into 0% from the CBE residual, ∼1% from conditional velocity shape, and ∼99% from spatial density.With two free exponential components, including ν(z) yields a χ2 span of 73 versus 0.7 when it is excluded.Source: §2 · p. 2
  • Every two-integral f(E, Lz) model has σR = σz, but anisotropy does not distinguish the solution from all stationary models because three-integral distribution functions can be anisotropic.An exactly separable third integral requires a Stäckel-form potential, which the paper neither imposes nor tests.Source: §2 · p. 2

3 Method

The pipeline models the selection-corrected stellar distribution, measures accelerations in spatial cells without assuming a potential, then fits a closed-form potential under a Poisson-positivity constraint.

3 Method

  • The selection gradient peaks at 10.9 kpc−1 versus a physical vertical-density gradient of 3.0 kpc−1; a Poisson GLM estimates a 2.88 kpc radial scale length and 334 pc local scale height at 300 pc.The GLM factors the observed intensity into distance selection, angular selection, and spatial density, identifiable up to two overall multiplicative constants.Source: §3 · p. 3
  • A 96-component full-covariance Gaussian mixture models the velocity conditional, while refits show density estimation limits the gradients and imply a 48.7 km s−1 kpc−1 residual floor.The conditional is selection-free because S has no velocity dependence; refitting shifts ∂R ln f by 0.815 against an rms of 1.673.Source: §3 · p. 3
  • They solve 94 acceleration cells (7 radial bins × 14 vertical bins − 4 star-count-cut cells) by Huber-reweighted least squares without assuming a functional form for Φ, recovering a field of ∼0.6 cm s−1 yr−1, undetectable in any single star over the mission lifetime but recoverable from collective statistics of many.Source: §3 · p. 3
  • 231.4 km s−1 for whole-volume vc agrees with 229±3 [17]; cells yield Σ(<0.5 kpc)=44.0 M⊙pc−2 and place the midplane 18 pc from the Sun.The surface density lies between measurements of 41 and 65 ± 6 M⊙pc−2 at 0.35 and 0.8 kpc; the midplane is compared with a photometric 20.8 ± 0.3 pc.Source: §3 · p. 3
  • Least-squares errors understate uncertainty: structure below about 1% is unresolved, and the fitted azimuthal acceleration departs from zero by 1.3% of |aR|.The reported errors are statistical errors from millions of stars, and the departure is comparable to the stated uncertainty floor.Source: §3 · p. 3
  • Symbolic regression fits Φ to the independently measured field; a one-sided Poisson-positivity penalty and parsimony select the closed form.The fit uses both force components, and the positivity constraint addresses equally valid fields distinguished only by the Laplacian sign.Source: §3 · pp. 3–4, Figure 2 · p. 4
  • The distribution-function search supervises values and five gradients of ln f because the CBE residual cannot identify Φ; selection uses R2 against ln f.A constant expression can score 2.000 under the gradient-and-value loss, although it scores zero under a residual objective; the residual and anisotropy are computed afterward.Source: §3 · p. 4

4 Results

The recovered potential has a vertical ln cosh form characteristic of a self-gravitating isothermal sheet, while the closed-form distribution function is non-separable in velocity. The fitted field reproduces three of five external comparisons, but the local density remains underestimated.

4 Results

  • The complexity-11 potential contains a vertical ln cosh z dependence characteristic of a self-gravitating isothermal sheet, emerging from the free operator set.The supplied expression includes a 0.39193 ln cosh z term.Source: §4 · p. 4
  • The complexity-44 distribution function is non-separable in velocity: vR enters through vRvϕ, while large-|vz| behavior is approximately exponential.Source: §4 · p. 5
  • Three of five external comparisons are broadly reproduced, and top Pareto-front models bracket the measured σR/σz=1.893 with a 1.83–2.13 spread.Table 2 reports held-out results; σR/σz is measured from the same sample, while the other comparisons are independent.Source: §4 · p. 5, Table 2 · p. 5
  • Removing azimuthal streaming yields a relative residual of 0.204, versus 0.76 with it retained; no axisymmetric force can balance that term, which carries 93% of the signal.Source: §4 · p. 5
  • Separable potentials are 43% worse than the coupled form, which predicts az falling 22% across the sample where the separable form predicts no change.The recovered coupling means Ez is not available as a third integral here.Source: §4 · p. 5
  • The observed anisotropy remains stable to ±1% across apertures from |z| < 0.05 to |z| < 0.30 kpc.Source: §4 · p. 5
  • The local density is 0.048 versus a literature range of 0.084–0.10 M⊙pc−3, and smoothing ν shifts the recovered (Σ, ρ) split from (50, 0.010) to (66, 0.003).The density estimate from the measured field is 0.0493, close to the fitted value; the passage attributes the deficit upstream to the density estimator.Source: §4 · pp. 5–6

5 Conclusion

Using Gaia DR3, the paper reconstructs closed-form local Galactic potential and stellar distribution function, finding that stellar number counts provide the main potential constraint while fitting the stationary CBE alone remains degenerate. The recovered potential contains a characteristic ln cosh term, but its local density is underestimated and its scope is limited by assumptions and selection uncertainty.

5 Conclusion

  • Controlled mock tests identify stellar number counts as the main constraint on the potential, while joint fitting of f and Φ to the stationary CBE alone remains degenerate.Source: §5 · p. 6
  • The pipeline reconstructs the acceleration field from a completeness-corrected distribution function, then fits a symbolic potential under a Poisson-positivity constraint.Source: §5 · p. 6
  • The recovered potential contains a ln cosh term characteristic of a self-gravitating isothermal sheet, with a fixed vertical scale of 2h = 1 kpc.Source: §5 · p. 6
  • Three of five external benchmarks are broadly reproduced, and the top Pareto-front models bracket the observed velocity anisotropy.Source: §5 · p. 6
  • The analysis assumes steady state and axisymmetry within 1 kpc despite measuring a non-zero azimuthal acceleration.Source: §5 · p. 6
  • Although surface density is approximately recovered, local density is underestimated by about a factor of two; an 11% systematic in ∂z ln ν dominates the uncertainty in Φ.Source: §5 · p. 6

A Mock anatomy of the degeneracy

Mock tests show that stellar spatial density, rather than conditional velocity shape, carries nearly all the constraint on local dark matter density when the distribution function is free. Unmodeled completeness gradients and density smoothing can bias the recovered disc–halo split even when fit quality does not reveal the problem.

A Mock anatomy of the degeneracy

  • The clean mock recovers (Σ, ρDM) = (50, 0.010) from inputs (48 M⊙pc−2, 0.010 M⊙pc−3) before a completeness gradient is injected.The mock uses an isothermal-sheet plus dark-matter potential and a two-component exponential distribution function.Source: A Mock anatomy of the degeneracy · p. 8
  • The spatial-density share of the dark-matter constraint stays above 96% across mock configurations and above 98% for apertures |z|max ≥ 0.75 kpc.Table A1 compares conditional-velocity-shape-only and shape-plus-density Δχ² spans with the distribution function free.Source: Table A1 · p. 8
  • Ignoring a completeness gradient makes an isothermal tracer mistake it for a thin sheet, shifting recovered surface density to 74 M⊙pc−2 and ρDM across 0.003–0.027 M⊙pc−3.The clean mock recovery is stable at ρDM = 0.010 M⊙pc−3, but the biased run has an indistinguishable loss floor.Source: A Mock anatomy of the degeneracy · p. 8
  • Smoothing the density estimate by 1.2 cells shifts the mock disc–halo split from (Σ, ρ) = (50, 0.010) to (66, 0.003).The passage attributes the shift to Σ integrating vertical force while ρ differentiates it, moving density from the halo term into the sheet.Source: A Mock anatomy of the degeneracy · p. 8
  • Freeing the distribution function destroys the ρDM constraint from conditional velocity shape alone, while adding star counts restores it.Figure A1 compares the fixed-truth distribution function with a free one, with and without spatial density.Source: Figure A1 · p. 9

B Data and the selection model

The selection model jointly fits survey completeness, angular dependence, and tracer density, while carrying a dust-related systematic into the potential inference.

B Data and the selection model

  • The Poisson GLM jointly fits μ = V·S(d)A(ℓ, b)ν(R, z), separating factors because equal-distance stars in different directions lie at different heights.The adopted catalogue frame fixes z⊙ = 0, so the dynamical midplane offset is measured against that reference.Source: B Data and the selection model · p. 9
  • The fitted density factor gives a 334 pc tracer scale height and a 2.88 kpc radial scale length, neither fixed to literature values.Source: B Data and the selection model · p. 9
  • Refitting on dust-poor |b| > 20° sight lines shifts ∂z ln ν(0.3 kpc) from −2.99 to −2.67 kpc−1, an 11% systematic carried forward as a limit on Φ.Source: B Data and the selection model · p. 9

C Ablations and model-selection protocol

The ablations favor a coupled potential with a smooth, isothermal-sheet-like vertical profile, while input scaling and held-out selection materially affect model evaluation. The mixture’s velocity moments also drift from the dispersions beyond about 0.4 kpc.

C Ablations and model-selection protocol

  • Separable potentials incur 43% more loss at complexity 9–10 than the coupled form, which captures a 22% decline in az across the sampled radial range.The coupling makes Ez unavailable as a third integral here because conservation requires ∂2Φ/∂R∂z = 0.Source: C Ablations and model-selection protocol · pp. 9–10
  • The 96-component mixture reproduces held-out velocity moments near the midplane but drifts from dispersions beyond |z| ∼ 0.4 kpc.Because the mixture is fitted to individual stars, the moments provide an independent check and the drift contributes to the estimator floor.Source: Figure A3 · p. 10
  • The ln cosh vertical profile reaches χ2/dof = 275.4, statistically tied with the softened-modulus form at 275.8; profiles with a midplane kink are an order of magnitude worse.The data favor the isothermal-sheet family but do not distinguish its two smooth parameterizations.Source: C Ablations and model-selection protocol · p. 10
  • Rescaling inputs by measured scale height and length lowers the residual from 38.15 to 23.69 km s−1 kpc−1 and raises R2 from 0.923 to 0.937 at identical complexity.The authors report that unit choices are worth about 10 complexity points and use scaled variables in all searches.Source: C Ablations and model-selection protocol · p. 10
  • The measured anisotropy σR/σz remains stable to ±1% across apertures |z| < 0.05–0.30 kpc.The reported spread in Sec. 4 instead comes from the choice of front member.Source: C Ablations and model-selection protocol · p. 10
  • Held-out scoring reverses the in-sample preference for a spuriously curved mock member at every tested complexity, motivating fixed 70/15/15 splits before fitting.For the line, χ2/N is 1.21 in-sample versus 1.53 held-out; for the exponential it is 0.89 versus 1.71.Source: C Ablations and model-selection protocol · p. 10
  • The acceleration solve uses linear least squares over 94 cells and runs in minutes on a single CPU node.Each symbolic search takes a few hours on one multi-core CPU node; the reported ablations total about 103 CPU-core-hours.Source: C Ablations and model-selection protocol · p. 10
Loading 2609.09011v1…