The Rotation Curve of the Milky Way: State of the Art, the Keplerian Decline Debate, and Implications for Dark Matter

This review paper addresses a fundamental question in Galactic astronomy: Whether its outer rotation curve is flat, gently declining, or Keplerian is the question this review addresses. The paper states that "For four decades, the roughly flat rotation curve (RC) of the Milky Way (MW) stood as local evidence for an extended dark matter (DM) halo. Gaia DR3 has changed the picture: several analyses now find a velocity decline beyond R ≈ 15 kpc, and the most radical interpretation (Jiao et al., 2023) claims a nearly Keplerian fall-off (vc ∝ R−1/2) that excludes a flat RC at 3σ and implies a total dynamical mass of only ∼ 2 × 10 11 M⊙ —three to five times below pre-Gaia estimates."

The paper emphasizes the stakes: "Taken at face value, this would make the MW exceptional among comparable spirals and challenge both ΛCDM and MOND. It would also alter predictions for direct-detection experiments, although the local DM density is not fixed by virial mass alone."

The paper traces the history: "The revolutionary observational discovery of the 1970s and 1980s was that galaxies do not behave this way. Rubin et al. (1978) showed optically, and Bosma (1981) confirmed with 21 cm HI measurements extending to much larger radii, that the RCs of spiral galaxies remain flat or even slightly rising at radii well beyond the stellar disc."

Regarding the Milky Way's unique position, the paper notes: "It is the only galaxy whose individual stars can be observed in full six-dimensional phase space, and at the same time the one for which the rotation curve is hardest to extract, precisely because we observe it from within the disc we are trying to measure."

The Gaia mission's impact is described: "The Gaia mission (Gaia Collaboration et al., 2016), launched in December 2013, measures parallaxes and proper motions for ∼ 1.5 billion stars, with parallax precision σϖ ≈ 0.01–0.1 mas, and radial velocities for ∼ 33 million objects in DR3 (Katz et al., 2023). For the first time, the full six-dimensional phase-space distribution of millions of stars in the disc and halo is accessible, enabling RC determinations whose statistical errors are smaller than the systematic uncertainties."

Eilers et al. (2019) established the pre-DR3 benchmark using 23 000 RGB stars with APOGEE spectroscopy. Their results: vc(R0) = 229.0 ± 0.2 km s−1 (formal), dvc/dR = −1.7 ± 0.1 (formal) ± 0.46 (syst.) km s−1 kpc−1, Mvir = (7.25 ± 0.26) × 10 11 M⊙, c = 12.9. These results apply over R = 5–25 kpc.

Wang et al. (2023) applied the Lucy iterative deconvolution method to the full Gaia DR3 source catalogue... Their RC declines from vc ≈ 235 km s−1 at 8 kpc to vc ≈ 195 km s−1 at 28 kpc. Importantly, they found a significant north-south kinematic asymmetry at R ∼ 20–25 kpc—stars above and below the plane show mean radial velocities differing by ∼ 10–20 km s−1—directly indicating non-axisymmetric or non-equilibrium kinematics.

"Jiao et al. (2023) re-derived the RC of Wang et al. (2023) with a revised binning and a full systematic error budget... Both RCs decline with a very similar linear slope, β = −(2.18 ± 0.23) km s−1 kpc−1 for their own measurement and β = −(2.22 ± 0.20) km s−1 kpc−1 for that of Ou et al. (2024). They then fit a power-law model vc ∝ R γ to the data at R > 19 kpc, finding γ = −0.47 ± 0.15 for their RC and γ = −0.56+0.23−0.22 for the Ou et al. (2024) RC, both consistent with Keplerian (γ = −0.5) at 1σ. A flat RC (γ = 0) is excluded at the 3σ level."

The resulting mass estimate: Jiao et al. (2023) quote 2.06+0.24−0.13 × 10 11 M⊙ for their preferred total mass and a strict upper limit of 5.4 × 10 11 M⊙.

"Ou et al. (2024) used an axisymmetric Jeans analysis but neglected the vertical cross-term in equation (22)... Their RC extends to R ≃ 27.3 kpc, declining smoothly from vc ≈ 230 km s−1 at 5 kpc to ≈ 173 km s−1 at 27.3 kpc. Their preferred Einasto fit has M200 = 1.81+0.06−0.05 × 10 11 M⊙ and ρDM,⊙ = 0.447 ± 0.004 GeV cm−3. The much poorer gNFW fit gives M200 = 6.94 ± 0.12 × 10 11 M⊙ and ρDM,⊙ = 0.405 ± 0.004 GeV cm−3."

The paper derives the full formalism from first principles. The radial Jeans equation is:

v2c(R) = ⟨vϕ⟩2 + σ2ϕϕ − σ2RR − (R/ν) ∂(νσ2RR)/∂R − (R/ν) ∂(νσRz)/∂z

The asymmetric drift formula is derived as:

v2a = −σ2RR (∂ ln ν/∂ ln R + ∂ ln σ2RR/∂ ln R + 1 − κ2) − (R/ν) ∂(νσ2Rz)/∂z

The paper notes: "For plausible outer-disc tracer profiles the quadratic correction is va ∼ 30–55 km s−1... even a 30% error in va shifts vc by only ∼ 1–4 km s−1, below the ∼ 30 km s−1 outer-disc decline claimed by Jiao et al. (2023). The asymmetric drift correction by itself is therefore not a dominant systematic in the outer-disc RC inference."

The paper derives the NFW profile properties: MNFW(r) = 4πρs r3s [ln(1 + r/rs) − (r/rs)/(1 + r/rs)] and notes The peak of vc,NFW occurs at rmax ≈ 2.16 rs. For the Milky Way: "At R ≈ 20 kpc ≈ 1.35 rs the Galaxy is still inside the radius rmax ≈ 2.16 rs ≈ 32 kpc at which the halo-only vc,NFW curve peaks—the total RC (baryons + halo) has already begun declining, but at a rate much slower than Keplerian."

The paper presents its own illustrative analysis using a compilation of 35 measurements spanning 6.3 ≤ R/kpc ≤ 49.0. The phenomenological model is:

V(R) = A + B/[1 + (R/C) n]

Results: "The 68% credible intervals for the best-fit parameters are: A = 193.6+3.5−4.2 km s−1, B = 28.9+4.8−3.9 km s−1, C = 18.9+0.8−0.8 kpc, and n = 12+6−4. The derived inner-plateau velocity... is inferred as Vin = 222.5+1.5−1.5 km s−1, consistent with the canonical value vc(R0) ≈ 229 km s−1; the outer asymptote is Vout = A = 193.6+3.5−4.2 km s−1."

The mass deficit: The inferred Msph,30 = 2.62+0.09−0.10 × 10 11 M⊙ is ∼ 25% lower than the mass expected from a flat RC at the same radius. The discrepancy reaches ∆M ≈ 8.5 × 10 10 M⊙ at R = 30 kpc.

The Gaussian Process reconstruction confirms: a mass deficit of ∆M ∼ 8–9 × 10 10 M⊙ builds up for R ≳ 20 kpc and peaks at R ≈ 30 kpc.

The paper cautions: "The low goodness-of-fit statistic requires care. A reduced chi-square substantially below unity (χ2/ν ≈ 0.47 here) should not be over-interpreted as evidence of an exceptionally good model... the credible intervals quoted in this section should be regarded as approximate."

The paper estimates: At R = 20 kpc with σRR ≈ 35 km s−1, A ≈ 2.4... δ(v2a) ≈ 2060 (km s−1)2, δva ≈ 19 km s−1. The corresponding shift in the inferred vc is δvc ≈ δ(v2a)/(2vc) ≈ 4 km s−1.

The paper addresses a key criticism: "Klačka & Šturc (2025) showed that in the best-fit mass models of Jiao et al. (2023) and Ou et al. (2024), the spherically distributed DM halo accounts for ≳ 88% of the total radial gravitational acceleration at R = 27 kpc, and argued on this basis that the cylindrical Jeans analysis becomes inappropriate there. The paper clarifies: A spherical potential is a special case of an axisymmetric one, so the axisymmetric Jeans equations in cylindrical coordinates... remain mathematically valid regardless of whether the acceleration is dominated by the disc or by a spherical halo. However, the closure assumptions... that are benign where the disc dominates are much less well tested in the halo-dominated regime."

"The Sagittarius dwarf galaxy's multiple pericentre passages excite bending waves in the disc with amplitudes ∼ 1–2 kpc and associated non-circular motions of ∼ 10–20 km s−1 at R ∼ 20 kpc (Laporte et al., 2018). The LMC (MLMC ≈ 1–2 × 10 11 M⊙) creates a dynamical-friction wake and induces a reflex motion of the MW of order tens of km s−1."

FIRE-2 simulations (Ou et al., 2025) find combined biases of ∼ 10–40 km s−1 at R ∼ 20–25 kpc when the full simulation kinematics are passed through a Jeans-equation pipeline mimicking observational analyses. Also: "Koop et al. (2024)... finds that the difference between the true circular-speed curve and the Jeans-inferred one can reach ∼ 10–15% in MW-like systems, and cautions specifically that a radially truncated tracer density can produce a spuriously steep apparent decline at large R even when the underlying RC is essentially flat."

The paper derives the MOND prediction: In the deep-MOND regime (g ≪ a0), µ(g/a0) ≈ g/a0, so equation (62) becomes g2/a0 = gbar. For circular orbits (g = v2c/R): v4c = GMbar a0 — the baryonic Tully-Fisher relation.

The tension: At R ∼ 20 kpc the MW is in the transition toward the low-acceleration regime... gbar(R = 20 kpc) ≈ (0.17–0.29) a0. The paper notes: "Coquery & Blanchard (2025) found that a standard baryonic model with the standard MOND acceleration scale fails to reproduce the declining RC. When both the baryonic model and a0 were allowed to vary, the fit preferred a heavy stellar disc of order 10 11 M⊙ and an a0 consistent with zero, with a reported 95% upper bound a0 < 0.53 × 10−10 m s−2."

Caveats: the flat-RC prediction of equation (63) is strictly an asymptotic result... at R ∼ 20 kpc the MW is in the transition toward that regime, not asymptotically deep within it and the MW is not an isolated system: the external-field effect (EFE)... can produce a declining outer RC for galaxies embedded in the gravitational field of a larger structure.

Moffat et al. (2025) fit the Jiao et al. data with αMOG = 13 and λ = 25 kpc, obtaining Mtot ≈ 1.3 × 10 11 M⊙, attributing the residual mass to the hot CGM plasma.

"The Orphan-Chenab stream (Koposov et al., 2023) extends to r ∼ 60 kpc, constraining: MMW(< 50 kpc) = 3.77+0.26−0.19 × 10 11 M⊙, already exceeding the total mass from the Keplerian interpretation. The Sagittarius stream gives M(< 100 kpc) = 5.6 ± 0.4 × 10 11 M⊙ and Mvir = 9.0 ± 1.3 × 10 11 M⊙."

"Eadie & Jurić (2019) obtained a median M200 = 0.70 × 10 12 M⊙ with a 50% interval 0.62–0.81 × 10 12 M⊙, and M(< 100 kpc) = 0.53+0.21−0.12 × 10 12 M⊙. Posti & Helmi (2019): Mvir = (1.3 ± 0.3) × 10 12 M⊙."

Given (r0 = 770 kpc, vr = −110 km s−1, t0 = 13.8 Gyr), numerical solution gives MLG ≈ (4.7±0.6)×10 12 M⊙ (Benisty, 2024).

"Hayati et al. (2025)... obtain log10(Mvir/M⊙) = 12.20+0.16−0.14, i.e. Mvir ≈ 1.6 × 10 12 M⊙, in good agreement with the stream, GC, and other satellite-based estimates and in strong tension with the Jiao+23 Keplerian value."

The paper discusses a statistical test: "Bariego-Quintana & Llanes-Estrada (2026), who applied a uniform slope-fitting procedure to the DM-subtracted rotation curves of all 175 SPARC galaxies... they find that the SPARC sample is fully consistent with a null terminal slope — the binomial, Wilcoxon, and Z-test p-values all cluster around 0.9–1.0, providing no statistical evidence for a Keplerian decline in external spirals. Applying the same pipeline to the modern MW data... yields, by contrast, a clearly negative terminal slope of the DM contribution of −1.57 ± 0.11 km s−1 kpc−1 in the combined sample for R > 15 kpc."

The paper outlines: "In Scenario 1, the Keplerian decline is real and the outer-halo mass tracers are systematically biased upward... This would require multiple independent methods (streams, GCs, satellites, timing argument) to fail simultaneously in the same direction—unlikely but not logically excluded."

"Scenario 2 is a plausible intermediate interpretation: there is a genuine decline at R > 15 kpc, but its magnitude is amplified by systematic errors in the Jeans methodology. The true virial mass lies in the range (5–8) × 10 11 M⊙."

In Scenario 3, the decline is entirely an artefact of the cylindrical Jeans methodology, and the true RC is nearly flat at vc ≈ 220–230 km s−1 out to R ∼ 30 kpc, consistent with Mvir ≈ 10 12 M⊙.

The paper notes: "The local density cannot be inferred from a virial mass alone. Indeed, the two Ou et al. (2024) halo fits have very different M200 values but give ρDM,⊙ = 0.447 ± 0.004 GeV cm−3 (Einasto) and 0.405 ± 0.004 GeV cm−3 (gNFW)."

Regarding the escape speed: "For the canonical NFW halo (rs = 14.8 kpc, ρs = 1.06 × 10 7 M⊙ kpc−3), equation (71) gives a halo-only escape speed of vesc(R0) ≈ 447 km s−1; adding the baryonic contribution... raises the total to vesc(R0) ≈ 530–550 km s−1. Gaia measurements of high-velocity stars (Deason et al., 2019) give vesc(R0) = 528+24−25 km s−1."

For microlensing: "García-Bellido & Hawkins (2024) found that adopting the Jiao et al. density profile allows a MACHO fraction up to 20% for masses 10−3–1M⊙—relaxing the earlier tight constraints."

The paper's central conclusions are:

1. "Several Gaia-era analyses find a decline beyond R ≳ 15 kpc, but its magnitude and even the robustness of the outermost tail remain contested because the analyses share distance calibrations, tracer assumptions, and equilibrium modelling."

2. "The Keplerian interpretation is plausible but systematically challenged... FIRE-2 simulations find combined biases of ∼ 10–40 km s−1 at R > 20 kpc."

3. RC-based virial masses are strongly profile-dependent: for example, the two Ou et al. (2024) fits span M200 ≃ 1.8–6.9 × 10 11 M⊙ despite using the same RC.

4. An isolated MOND model predicts an asymptotically flat RC, so a confirmed Keplerian tail would be difficult to accommodate.

5. The local DM density is not determined by the virial mass alone.

The paper calls for: "Resolution requires better data and improved modelling: Gaia DR4, scheduled for 2 December 2026, WEAVE, 4MOST, and DESI for kinematics; SKA-MID for maser parallaxes; and a consistent three-dimensional Jeans analysis validated against cosmological simulations, using the spherical (or oblate) Jeans equation at large radii."

Improvements to AI Systems Based on This Paper

---

• Improvement: Build an AI system that explicitly models systematic uncertainties (e.g., asymmetric drift corrections, Jeans equation closure assumptions, tracer density truncation) as latent variables rather than fixed inputs. Use Bayesian hierarchical models that propagate these uncertainties through rotation curve inference.

• What it can do: Automatically flag when inferred rotation curves (and derived dark matter masses) are dominated by methodological choices rather than data. It can output probability distributions over virial masses that span the full range of plausible systematics (e.g., 1.8–6.9 × 1011 M⊙ for the same data), preventing overconfident claims like the Keplerian interpretation.

• Improvement: Train a neural network on cosmological zoom-in simulations (e.g., FIRE-2) to learn the mapping between true circular speed curves and Jeans-inferred ones, including biases from non-equilibrium features (bending waves, LMC wake, Sagittarius perturbations). Use this as a correction layer in observational pipelines.

• What it can do: Automatically correct observed rotation curves for known biases (10–40 km/s at R > 20 kpc) before fitting dark matter profiles. This would reconcile conflicting results (e.g., Jiao et al. vs. Eilers et al.) by providing a de-biased rotation curve that is more robust to non-axisymmetric kinematics.

• Improvement: Develop an AI system that jointly fits rotation curve data, stellar stream dynamics (Orphan-Chenab, Sagittarius), globular cluster kinematics, satellite galaxy orbits, and the timing argument within a single self-consistent Bayesian framework. Use normalizing flows to model the full phase-space distribution of tracers.

• What it can do: Automatically detect inconsistencies between rotation-curve-derived masses and independent dynamical tracers. For example, it would flag that a Keplerian decline (M tot ≈ 2×1011 M⊙) is incompatible with Orphan-Chenab stream constraints (M(3σ, forcing the system to either reject the Keplerian claim or identify which tracer assumptions are flawed.

• Improvement: Replace fixed halo profile assumptions (NFW, Einasto, gNFW) with a flexible, non-parametric density field (e.g., using Gaussian processes or deep generative models) that is constrained by rotation curves, escape velocity measurements, and local density estimates simultaneously.

• What it can do: Provide mass estimates that are not artificially narrowed by profile choice. The system would output a full posterior over the radial mass distribution, showing that M 200 is poorly constrained (spanning an order of magnitude) even with perfect rotation curve data, and explicitly quantify how much of the uncertainty comes from profile degeneracy versus data quality.

• Improvement: Train an anomaly-detection AI on simulated galaxies with known perturbations (satellite impacts, bar instabilities, spiral structure) to identify non-axisymmetric signatures in stellar kinematics (e.g., north-south asymmetries, radial velocity gradients).

• What it can do: Automatically flag regions of the rotation curve where equilibrium assumptions break down (e.g., the 10–20 km/s asymmetry at R 20–25 kpc found by Wang et al.). This would prevent over-interpretation of apparent declines that are actually due to bending waves or tidal interactions, rather than dark matter physics.

• Improvement: Build a model-selection AI that simultaneously fits ΛCDM, MOND (with external field effect), and MOG to all available Milky Way data (rotation curve, streams, escape speed, satellite orbits). Use Bayesian evidence ratios to quantify relative support.

• What it can do: Automatically compute the probability that a Keplerian decline is better explained by modified gravity (e.g., MOND with EFE) versus a low-mass dark matter halo. It would also identify which additional observations (e.g., outer halo stream velocities) would most efficiently break the degeneracy between theories.

• Improvement: Create an AI system that predicts local dark matter density (ρ DM,⊙) from global halo parameters using a learned mapping from cosmological simulations, rather than assuming a single profile. Include escape velocity and microlensing constraints as additional inputs.

• What it can do: Provide robust predictions for direct detection experiments that are insensitive to the virial mass controversy. For example, it would show that ρ DM,⊙ ≈ 0.4–0.45 GeV/cm3 regardless of whether M 200 = 2×1011 or 7×1011 M⊙, decoupling local experimental predictions from the outer rotation curve debate.

• Improvement: Train a system to predict how future datasets (Gaia DR4, WEAVE, 4MOST, DESI, SKA-MID masers) will reduce uncertainties on the rotation curve and dark matter mass. Use mock catalogs with realistic systematics to optimize observing strategies.

• What it can do: Automatically prioritize which measurements (e.g., maser parallaxes at R > 20 kpc vs. stellar stream velocities at R > 50 kpc) will most efficiently resolve the Keplerian decline debate. It can also forecast whether Gaia DR4 alone will be sufficient or whether complementary datasets are essential.

• Improvement: Build an NLP-based system that reads all rotation curve papers, extracts reported values and error bars, and automatically identifies shared assumptions (e.g., same distance calibration, same Jeans closure) that could induce correlated biases across seemingly independent measurements.

• What it can do: Produce a meta-analysis that correctly accounts for non-independent measurements. It would show that the apparent consensus on a declining rotation curve may be inflated by shared systematics, and quantify the true effective number of independent constraints.

• Improvement: Implement an active-learning AI that, given a new rotation curve measurement, automatically generates a suite of mock observations from competing models (flat, Keplerian, declining-but-not-Keplerian) and computes which future measurement would most strongly distinguish them.

• What it can do: Serve as an intelligent observation scheduler for telescopes, dynamically updating target lists based on which data points would most reduce model uncertainty. This accelerates the resolution of the dark matter debate by focusing resources on the most informative observations.