# Cosmological implications of the Gaia Milky Way declining rotation curve

E. Coquery, Alain Blanchard

► To cite this version:

E. Coquery, Alain Blanchard. Cosmological implications of the Gaia Milky Way declining rotation curve. Astronomy & Astrophysics - A&A, 2025, 703, pp.A88. 10.1051/0004-6361/202556337. hal-05352818

HAL Id: hal-05352818

https://ut3-toulouseinp.hal.science/hal-05352818v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmological implications of the Gaia Milky Way declining rotation curve

E. Coquery<fnref n="1" />,*, and A. Blanchard<fnref n="2" />,†

Received 9 July 2025 / Accepted 19 September 2025

## ABSTRACT

Although the existence of dark matter has been widely acknowledged in the cosmology community, its nature remains unknown despite decades of research. This calls into questions its very existence. This never-ending search for dark matter has led us to consider alternatives. Since increasing the enclosed mass is the only way to explain the flat appearance of galaxies' rotation curves in a Newtonian framework, modified Newtonian dynamics (MOND) theory proposes a modification to Newton's dynamics when the acceleration of matter at given radius is around or below a threshold value, $a_0$. Observed rotation curves, generally flat at large distances, are then usually well reproduced by MOND with an $a_0$ of ~1.2 × 10⁻¹⁰ m/s<fnref n="2" />. However, recent Gaia evidence indicates a decline in the Milky Way rotation curve, which a distinct behavior. Therefore, we examine whether MOND can accommodate the Gaia declining rotation curve of the Milky Way. We first use a standard model to describe the Milky Way's baryonic components. Secondly, we show that a Navarro, Frenk, and White model is able to fit the decline, assuming a scale radius ($R_s$) on the order of 4 kpc. In a third step, we show that the usual MOND paradigm is not able to reproduce the declining part in a standard baryonic model. Finally, we examine whether the MOND theory can accommodate the declining part of the rotation curve when the characteristics of the baryonic components are relaxed. To do so, we used a Markov chain Monte Carlo method to determine the characteristics of the stellar and HI disks, including their masses, that would result in the best possible fit of the decline. We find that the stellar disk should be massive, on the order of 10¹<fnref n="1" /> M☉. The HI disk mass is capped at nearly 1.8 × 10¹<fnref n="1" /> M☉ but could also be negligible. Finally, $a_0$ is consistent with 0, with an upper limit of 0.53 × 10⁻¹⁰ m/s<fnref n="2" /> (95%), a value much lower than that usually advocated to explain standard flat rotation curves in MOND theory (~1.2 × 10⁻¹⁰ m/s<fnref n="2" />).

**Key words.** Galaxy: kinematics and dynamics – dark matter

## 1. Introduction

Studying galaxies' rotation curves (RCs) is of paramount importance in cosmology, as they hint at the existence of dark matter. Observations have shown that most galaxies' RCs are flat, which is in strong disagreement with Newtonian dynamics' predictions for known baryonic components. Historically, cosmologists' favorite answer to this crisis is to assume the existence of an invisible mass in order to accommodate the observations. This assumption is also supported by the observations of unrelated phenomena, including gravitational lensing (Refregier 2003) and the cosmological microwave background (Ade et al. 2016). Dark matter, whose quantity was estimated to be nearly 30% of the Universe's mass-energy content, has now become a pillar of cosmology's standard model, Λ cold dark matter (ΛCDM).

However, since dark matter has never been observed directly, we could alternatively envisage that Newtonian dynamics changes under certain conditions. Modified Newtonian dynamics (MOND) theory (Milgrom 1983) proposes a modification of Newtonian dynamics in regimes where the acceleration approaches or drops below a characteristic value of the acceleration, $a_0$. In the deep MOND limit, the rotation velocity of a circular orbit resumes to $v = (GMa_0)^{1/4}$, which explains the flat appearance of most RCs. Previous fits of

RCs using MOND agree on a value of about 1.2 × 10⁻¹⁰ m/s<fnref n="2" /> (Begeman et al. 1991), an acceleration threshold so low that MONDian effects could not be detected on Earth nor in the Solar System. This theory has been further developed, leading to more complex formulations such as AQUAL (Bekenstein & Milgrom 1984), QUMOND (Milgrom 2023), or the relativistic TeVeS (Bekenstein 2004). However, these theories face tensions in multiple fields, such as galaxy cluster dynamics (Sanders 2003), cosmological microwave background anisotropies, and difficulties to fit the matter power spectrum (Dodelson 2011). To alleviate these tensions, hybrid theories have been put forward (see for example Bruneton et al. 2009 or Skordis & Złośnik 2021 for ΛCDM cosmology with MONDian effects at galactic scales). Despite those tensions, MOND remains a simple theory for reproducing flat galaxies' RCs. However, according to a statistical analysis of galaxy RCs from Spitzer Photometry and Accurate Rotation Curves (SPARC), the agreement does not favor MOND (Khelashvili et al. 2024).

Although it is relatively easy to measure other galaxies' RCs (Corbelli & Salucci 2000; Roberts 1975), probing the Milky Way's from the inside is a challenge at high radii. Available data tend to confirm that the MW's RC is flat (Mroz et al. 2019), consistent with both dark matter and MOND. However, some moderate decline has been reported (Robin et al. 2022), in agreement with an earlier claim consistent with the radial acceleration relation (McGaugh 2019). A moderate decline was also inferred for the most luminous galaxies (Persic et al. 1996; Salucci et al.

<footnote n="1">¹ École Centrale de Lyon, 36 Avenue Guy de Collongue, 69134 Écully, France</footnote>

<footnote n="2">² Université de Toulouse, UPS-OMP, IRAP, CNRS, 14 Avenue Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote>* Corresponding authors: even.coquery@gmail.com; alain.blanchard@irap.omp.eu</footnote>
2007). However, Gaia's latest data release (Vallenari et al. 2022) sheds light on a velocity decline after ≈15 kpc, on the order of 3.5 km/s/kpc, which different analyses agree on (see for example Wang et al. 2022; Zhou et al. 2023; Jiao et al. 2023; Ou et al. 2024). This decline is also consistent with previous indications for the presence of a dip in the Milky Way's RC (Huang et al. 2016) with a higher significance level. Some other galaxies' RCs were also found to be declining: a sample of 22 galaxies was studied under the MOND paradigm by Zobnina & Zasov (2020). They conclude that some galaxies' RCs do not align with the usual MOND paradigm, needing a value of $a_0$ lower than values obtained from previous fits of the RCs. Studying Gaia's newfound decline might yield similar conclusions regarding the Milky Way. Milky Way RC measurements before the decline are nonetheless more scattered; different papers have found very different values (e.g. Iocco et al. 2015; Mroz et al. 2019; Zhou et al. 2023; Sylos Labini et al. 2023; Wang et al. 2022; McGaugh 2018; Robin et al. 2022). Such considerations led us to focus on the declining part of the RC without considering the inner part.

In this study we aim to shed light on the declining RC's implications regarding MOND as inferred from Gaia, without attempting to fit the inner part. We started by building a model for the baryonic components of the Milky Way in order to compute the RC at various radii from the Galactic center. To compare ΛCDM with MOND, we implemented this model both under a dark matter paradigm and under Milgrom's modified dynamics. We then used this model as a basis to find the optimal value of $a_0$ that fits the decline – if such a value exists. Finally, we compared MOND results with ΛCDM, and we detail why the MOND paradigm does not accommodate the Milky Way RC under reasonable assumptions. This interpretation of the data is based on the assumption of simple dynamics that are not violently disturbed in the outer parts of the Galaxy (Koop et al. 2024; Kroupa et al. 2024).

## 2. The Milky Way's rotation curve

### 2.1. Modeling the baryonic components of the Milky Way

In order to compute the RC, we first need to establish a model for the spatial distributions of the baryonic components of the Milky Way. We chose to work with the B2 model used in Jiao et al. (2023) and described in de Salas et al. (2019). This model consists of three components: (i) a spherical bulge modeled by a Hernquist potential,

$$ \Phi(r) = -\frac{GM_{\text{bulge}}}{r+r_b}, \qquad (1) $$

(ii) a thin stellar disk, and (iii) multiple gas disks. Given the values provided by de Salas et al. (2019) and references therein, we only considered the HI disk since the other gas masses are negligible.

Both the stellar disk and the HI disk densities are given by a double exponential,

$$ \rho^i(R, z) = \rho_0^i \exp\left(-\frac{R}{r_d^i} - \frac{|z|}{z_d^i}\right), \qquad (2) $$

where $i \in \{\text{st.}, \text{HI}\}$ and $\rho_0^i$ is a normalization constant:

$$ \left\{ \begin{aligned} \rho_0^{\text{st.}} &= \frac{M^{\text{st.}}}{4\pi(r_d^{\text{st.}})^2 z_d^{\text{st.}}} \\ \rho_0^{\text{HI}} &= \frac{M^{\text{HI}}}{4\pi r_d^{\text{HI}} z_d^{\text{HI}}} (R_t + r_d^{\text{HI}}) e^{-R_t/r_d^{\text{HI}}} \end{aligned} \right. $$

<table><caption>**Table 1.** B2 model parameters (de Salas et al. 2019).</caption><thead><tr><th>Parameter</th><th>Value (B2 Model)</th></tr></thead><tbody><tr><td>$M^{\text{st.}}$ (M<sub>☉</sub>)</td><td>3.65 × 10<sup>10</sup></td></tr><tr><td>$r_d^{\text{st.}}$ (kpc)</td><td>2.35</td></tr><tr><td>$z_d^{\text{st.}}$ (kpc)</td><td>0.14</td></tr><tr><td>$M^{\text{HI}}$ (M<sub>☉</sub>)</td><td>8.2 × 10<sup>9</sup></td></tr><tr><td>$r_d^{\text{HI}}$ (kpc)</td><td>18.24</td></tr><tr><td>$z_d^{\text{HI}}$ (kpc)</td><td>0.52</td></tr><tr><td>$R_t$ (kpc)</td><td>2.75</td></tr><tr><td>$M_{\text{bulge}}$ (M<sub>☉</sub>)</td><td>1.55 × 10<sup>10</sup></td></tr><tr><td>$r_b$ (kpc)</td><td>0.7</td></tr></tbody></table>

We used the $M^i$, $M_{\text{bulge}}$, $r_d^i$, $z_d^i$, $R_t$, and $r_b$ values from de Salas et al. (2019), which are listed in Table 1.

The stellar population of the Galaxy can be described by more complex models with several different components. The presence of a thick disk, for instance, impacts the vertical structure of the velocity dispersion as well as the possible flaring of the disk (López-Corredoira 2025). To determine what the presence of a thick disk would imply, we also computed the circular velocity in a model comprising two equal-mass components (Pouliasis et al. 2017), one in a thin disk and one in a thick disk, following the Miyamoto–Nagai profile (Miyamoto & Nagai 1975). From Figs. 1 and 2, one can see that the RCs are nearly identical. Concerning the flaring, we note that Sylos Labini (2024) conclude that its impact has only a marginal effect on the RC at large distances, i.e., at the scales relevant to our study. In our work the RC is examined at a distance far from the bulk of the stellar components.

### 2.2. Rotation curve under ΛCDM

In Newtonian dynamics, assuming a circular movement, each mass component produces an acceleration field that can be written in terms of a circular velocity associated with this component:

$$ \frac{(v^i)^2}{r} = a = K_r^i, \qquad (3) $$

where $a$ is the radial acceleration, $K_r^i$ is the total radial force per unit mass, and $i$ corresponds to any baryonic component described above (stellar disk, bulge, or gas disk) or to the added dark matter mass component: $i \in \{\text{st.}, \text{bulge}, \text{HI}, \text{dm}\}$.

This means that in order to compute the velocity, we first have to determine the radial force. We did so by integrating Poisson's equation with the density $\rho^i$. For the density given in Eq. (2), Poisson's equation can be solved in terms of Hessel transforms. We can thus find the radial force by considering the R derivative. Details of calculation can be found in Kuijken & Gilmore (1989). We can easily determine $v_{\text{bulge}}$ using Eq. (3) and the first derivative of Hernquist's potential, Eq. (1).

Finally, we need to specify a dark matter model to compute $v_{\text{dm}}$. In ΛCDM, dark matter potentials are expected to follow a Navarro, Frenk, and White (NFW) profile (Navarro et al. 1996). We therefore used a standard NFW model (Navarro et al. 1996):

$$ \rho_{\text{NFW}}(r) = \frac{\rho_{0,\text{NFW}}}{r/R_s(1+r/R_s)^2}, \qquad (4) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. RC decline example using NFW with $\rho_{0,NFW} = 1.8 \times 10^8$ M$_\odot$/kpc$^3$, $R_{\text{thin}}^2 = 3.97$ kpc, and $R_{\text{thin+thick}}^2 = 3.98$ kpc. The parameters for the baryonic components are those of the B2 model in Table 1 taken from Jiao et al. (2023). We also provide the RC (gray lines) with two disks, one thin and one thick (see Sect. 2.1).</figcaption></figure>

where $\rho_{0,NFW}$ and the scale radius ($R_s$) are free parameters. Then we only needed to integrate the density of Eq. (4) between 0 and a given radius ($r$) to get the enclosed mass in the corresponding sphere, and then used the circular movement assumption to obtain the velocity:

$$v_{\text{dm}}(r)^2 = \frac{GM_{\text{dm}}(r)}{r}.$$

The RC under a circular movement assumption can now be evaluated:

$$v^2 = v_{\text{st.}}^2 + v_{\text{gas}}^2 + v_{\text{bulge}}^2 + v_{\text{dm}}^2. \qquad (5)$$

Using Eq. (5), we can now compute the RC under the $\Lambda$CDM dark matter paradigm. We provide an example in Fig. 1, with parameters chosen to yield an acceptable fit to the declining part of the RC. This example shows that NFW is able to accommodate the declining part, with a $\chi^2$ of $\chi_{\text{thin}}^2 = 2.6$, $\chi_{\text{thin+thick}}^2 = 3.1$ (reduced $\chi$ of 0.2 and 0.24), which is low but acceptable.

This yields different results than Jiao et al. (2023). When computing the dynamical mass using the same critical density value as Jiao et al. (2023) using the parameters from our fit of the RC under NFW, we get $M_{\text{dyn}} = 4.28 \times 10^{11}$ M$_\odot$ for a virial radius of 153.5 kpc. These values are higher than the results from Jiao et al. (2023) and consistent with Sylos Labini (2024). They are also consistent with their newfound upper limit on the dynamical mass. In addition, extrapolation to a distance nearly ten times greater is highly uncertain: on galaxy scales, strong feedback is expected (Blanchard et al. 1992), which could seriously alter the dark matter distribution in the inner part of galaxies (Li et al. 2022). The concentration parameter ($c$) for the above characteristics is larger than found by Eilers et al. (2019), but our fit is performed only on $R > 13$ kpc. Our value of $c$ is much larger than the average value expected from $\Lambda$CDM simulations involving only dark matter (Bullock et al. 2001). The role of baryons may, however, considerably increase the concentration parameter on galaxy scales (Shao et al. 2023). We note that Lin & Li (2019) infer a determination of the Milky Way's RC up to 100 kpc, consistent with a NFW profile, with data roughly consistent with those used in Fig. 1, although not reproduced by the global fit that covers the RC in the range 4.6–98 kpc from Huang et al. (2016). Finally, we note that dynamical arguments from the stellar streams dynamics (Ibata et al. 2024) as well as from the local group kinematics (Benisty 2024) point consistently toward a Milky Way mass on the order of or above $10^{12}$ M$_\odot$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Decline fit under the MOND paradigm. $a_0^{\text{thin}} = 2.417 \times 10^{-10}$ m/s$^2$, and $a_0^{\text{thin+thick}} = 2.429 \times 10^{-10}$ m/s$^2$. The parameters for the baryonic components are the same as in Fig. 1.</figcaption></figure>

## 2.3. Rotation velocity under MOND

Using the standard MOND formulation instead of a more rigorous one in a RC computation leads to a difference of about 5% (López-Corredoira & Betancort-Rijo 2021). For the sake of simplicity, we thus stuck to Newton's second law as modified in Milgrom (2015):

$$\mu(a/a_0)a = K_r, \qquad (6)$$

where $a$ is the acceleration that now may differ from the total radial force per unit mass ($K_r$), and $\mu$ is a function chosen such that (McGaugh 2004)

$$\begin{cases} \lim_{a \gg a_0} \mu(a/a_0) = 1 \\ \lim_{a \ll a_0} \mu(a/a_0) = a/a_0. \end{cases}$$

We chose to work with $\mu(x) = \frac{x}{\sqrt{1+x^2}}$, as it is standard and was used by McGaugh (2004)$^1$. By only keeping the $a > 0$ solution, one can derive the total velocity from Eq. (6):

$$v(r) = \left( \frac{r^2}{2} \left( K_r^2 + \sqrt{K_r^2 (K_r^2 + 4a_0^2)} \right) \right)^{1/4}.$$

We then fit the declining part of the RC with this MOND paradigm under the baryonic models described above. The value of $a_0$ was determined by minimizing the $\chi^2$ and considering the Jiao et al. (2023) declining part of the RC. The best value of $a_0$ found here is $2.417 \times 10^{-10}$ m/s$^2$, which is not consistent with the value derived from the RC in other galaxies. The resulting RC is shown in Fig. 2. The corresponding $\chi^2$, 49.7 and 57.6 (reduced $\chi^2$ of 3.83 and 4.43), confirm the visual impression that the fits are not satisfying. We thus conclude that using the B2 model described in Jiao et al. (2023), both with thin and thick disks, the declining RC cannot be reproduced with Milgrom dynamics. The physical reason for the behavior of the RC can be understood: the bulk of the mass is within 13 kpc, and the mass beyond 13 kpc is so low that the dynamic is in the MONDian regime and the RC is nearly flat.

<footnote>$^1$ This choice is not important as we focus on the asymptotic behavior of the curve. We have verified that the flattening of the RC is similar to the so-called simple function (Famaey & Binney 2005).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** MCMC contours for $r_d = 2.35$ kpc and $r_d = 3.1$ kpc. The stellar disk mass and $a_0$ are left as free parameters. The red cross indicates the observed values of $M^{\text{stellar}}$ and the standard value of $a_0$ with their respective uncertainty (see the main text).</figcaption></figure>

# 3. Relaxing the baryonic components to fit the decline with MOND

## 3.1. Methodology

We have seen that MOND theory cannot accommodate the declining part of the RC of the Milky Way when the baryonic components are those of the B2 model as presented in Sect. 2.1. In this section we examine whether relaxing the properties of the baryonic components could alleviate this inconsistency. As the disk mass is the dominant baryonic component, it is clearly an important parameter in the Milky Way RC. We thus treated this quantity as a free parameter, as well as $a_0$. For the scale radius, we used two values: $r_d = 2.35$ kpc (which is the value measured by Misiriotis et al. 2006 and used by both Jiao et al. 2023 and de Salas et al. 2019) and $r_d = 3.1$ kpc, which are close to the bounds provided by Bland-Hawthorn & Gerhard (2016).

To determine our two parameters, the disk mass ($M^{\text{st.}}$) and $a_0$, we ran two Markov chain Monte Carlo (MCMC) analyses, one for each scale radius. We applied the MCMC ensemble slice sampler algorithm from ZEUS (Karamanis et al. 2021). $M^{\text{st.}}$ was allowed to vary between $3 \times 10^{10} M_{\odot}$ and $2 \times 10^{11} M_{\odot}$ and $a_0$ between 0 and $3 \times 10^{-10} m/s^2$.

## 3.2. Results

Applying the methodology described above yields constraints on the parameter space (see Fig. 3). The uncertainty on $M^{\text{st.}}$ is representative of the extreme values found in Binney & Tremaine (2011) and Bland-Hawthorn & Gerhard (2016), and the uncertainty on $a_0$ can be found in Milgrom (2015). One can see that there is a correlation between the stellar disk mass and $a_0$. Moreover, the value of $r_d$ does not have much of an impact on the results. From these contours we extracted the pair of values ($M^{\text{st.}}$, $a_0$) presented in Table 2 that minimizes the $\chi^2$: $\chi^2 = 5.53$ for $r_d = 2.35$ kpc, and $\chi^2 = 3.84$ for $r_d = 3.1$ kpc; these values

<table><caption>**Table 2.** Central values for $r_d = 2.35$ kpc and $r_d = 3.1$ kpc.</caption><thead><tr><th>r<sub>d</sub> (kpc)</th><th>M<sup>st.</sup> (M<sub>☉</sub>)</th><th>a<sub>0</sub> (m/s<sup>2</sup>)</th></tr></thead><tbody><tr><td>2.35</td><td>(11.01<sup>+0.63</sup><sub>-0.76</sub>) × 10<sup>10</sup></td><td>(0.65<sup>+0.12</sup><sub>-0.10</sub>) × 10<sup>-10</sup></td></tr><tr><td>3.1</td><td>(10.95<sup>+0.62</sup><sub>-0.66</sub>) × 10<sup>10</sup></td><td>(0.592<sup>+0.104</sup><sub>-0.092</sub>) × 10<sup>-10</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Fits using MCMC results (Table 2). The parameters for the other baryonic components are the same as in Fig. 2.</figcaption></figure>

<table><caption>**Table 3.** Prior used for the MCMC in Fig. 5.</caption><thead><tr><th>Parameter</th><th>Lower bound</th><th>Higher bound</th></tr></thead><tbody><tr><td><i>M</i><sup>st.</sup> (M<sub>☉</sub>)</td><td>3 × 10<sup>10</sup></td><td>2 × 10<sup>11</sup></td></tr><tr><td><i>a</i><sub>0</sub> (m/s<sup>2</sup>)</td><td>0</td><td>3 × 10<sup>-8</sup></td></tr><tr><td><i>M</i><sup>gas</sup> (M<sub>☉</sub>)</td><td>0</td><td>3 × 10<sup>11</sup></td></tr><tr><td><i>r</i><sub>d</sub> (kpc)</td><td>0</td><td>20</td></tr><tr><td><i>z</i><sub>d</sub> (kpc)</td><td>0</td><td>2</td></tr></tbody></table>

are slightly higher than with NFW but still acceptable considering we only have two free parameters (see Fig. 4).

The obtained value of $a_0$ is smaller than the standard $1.2 \times 10^{-10} m/s^2$ and is not consistent with previous results, namely from Begeman et al. (1991). In other words, according to our study, there is no way for MOND to explain both the Milky Way RC and other galaxies' RCs with the same value of $a_0$. A similar conclusion is obtained by Chan & Law (2023).

Moreover, the stellar disk mass used to find such a value is about $10 \times 10^{10}$ solar masses, which is not consistent with the value $(3.5 \pm 1) \times 10^{10} M_{\odot}$ from Bland-Hawthorn & Gerhard (2016). The physical reason is that in order to produce a decreasing RC at $r > 15$ kpc, the dynamic needs to be close to the Newtonian regime, which in turn requires a high mass (more than 95% of the mass lies within a radius of 15 kpc). The MONDian regime starts at a larger radius with a value of $a_0$ lower than the standard value.

## 3.3. More freedom on the baryonic components

Since relieving constraints on the stellar disk mass alone does not yield satisfying results, in this section we explore possibilities by increasing the number of free parameters on the baryonic components. We launched another MCMC run with $M^{\text{st.}}$, $a_0$, $M^{\text{HI}}$, $r_d$, and $z_d$ as free parameters. Each prior used for this MCMC can be found in Table 3.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** MCMC contours with a free stellar disk mass, free HI disk mass, free $a_0$, free scale radius, and free scale height.</figcaption></figure>

A number of interesting properties can be inferred from Fig. 5: it is possible to adequately fit the RC with the baryonic components provided that they are allowed to take values well above the standard values: for the models taken in the $1\sigma$ domain, the $\chi^2$ typically lies between 3 and 6, which is acceptable and similar to the $\Lambda$CDM case. The scale height of the disk has no correlation with other parameters (except a weak correlation with the disk mass) and has no preferred value. The HI disk mass ($M^{HI}$) has little to no impact on the results; a wide range of values allowed us to fit the decline. High HI disk masses are allowed, far above the standard value but still below the preferred stellar disk mass. The $1\sigma$ intervals for the scale radius of the stellar disk and for its mass are higher than found in Bland-Hawthorn & Gerhard (2016) and show a correlation with the stellar disk mass. Relieving the constraints on the baryonic matter distribution yields a lower value for $a_0$ than found in the previous MCMC analysis, and $a_0$ does not appear to be correlated with the baryonic distribution parameters. One can see that $a_0 = 0 \text{ m/s}^2$ is in the $1\sigma$ interval, assuming a heavy stellar disk and a scale radius ($r_d$) > 3 kpc. The fact that the stellar disk is heavier than usual observations (Bland-Hawthorn & Gerhard 2016) can be explained by the close-to-zero value of $a_0$: since little to no MONDian effect is required and no dark halo is

assumed, one needs extra matter in the disk to reach an acceptable value of the circular velocity. Essentially, this consists in a dark matter disk instead of the halo used in Fig. 1, and models using a lighter disk and a nonvanishing $a_0$ are not significantly preferred over the $a_0 = 0$ paradigm.

**4. Discussion and conclusion**

In this work we examined the two major solutions to the missing mass problem in galaxies, applied to the Milky Way. More precisely, we compared the ability of $\Lambda$CDM and MOND to fit the Milky Way's declining RC measured by Gaia, assuming a simple dynamic that is not significantly disturbed in the outer parts of the Galaxy (Koop et al. 2024; Kroupa et al. 2024). Using a standard model for the baryonic components of the Milky Way, we show that a simple dark matter distribution model like NFW, expected in $\Lambda$CDM, is able to explain Gaia's decline with ease, albeit with an extreme value for the concentration parameter. On the other hand, the MOND formulation cannot accommodate the decline under the B2 model, even when allowing $a_0$ to be a free parameter. We then considered relieving constraints on baryonic parameters as well as the value of $a_0$ in order to examine whether MOND could accommodate the decline. For this
we performed an MCMC analysis on *a*₀ and the parameters of the baryonic components using the data of the declining part of the Milky Way RC. Preferred models have disk masses at odds with values inferred from observations, with no significant preference for a nonvanishing *a*₀. We get an upper limit on *a*₀ of 0.53 × 10⁻¹⁰ m/s⁻² (95%), which is significantly lower (by 5σ) than what has been found necessary to fit flat RCs in other galaxies with MOND.

We conclude that the declining RC of the Milky Way as recently inferred from Gaia data can be interpreted as being due to the presence of an NFW-type dark matter halo. It is not easy to reproduce in the MOND alternative.

**Acknowledgements.** This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. This work was supported by CNES.

## References

Ade, P. A. R., Aghanim, N., Arnaud, M., et al. 2016, A&A, 594, A13

Begeman, K. G., Broeils, A. H., & Sanders, R. H. 1991, MNRAS, 249, 523

Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509

Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7

Benisty, D. 2024, A&A, 689, L1

Binney, J., & Tremaine, S. 2011, Galactic Dynamics Second Edition (Princeton University Press)

Blanchard, A., Valls-Gabaud, D., & Mamon, G. A. 1992, A&A, 264, 365

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Bruneton, J.-P., Liberati, S., Sindoni, L., & Famaey, B. 2009, JCAP, 2009, 021

Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559

Chan, M. H., & Law, K. C. 2023, ApJ, 957, 24

Corbelli, E., & Salucci, P. 2000, MNRAS, 311, 441

de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, JCAP, 2019, 037

Dodelson, S. 2011, Int. J. Mod. Phys. D, 20, 2749

Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120

Famaey, B., & Binney, J. 2005, MNRAS, 363, 603

Huang, Y., Liu, X. W., Yuan, H. B., et al. 2016, MNRAS, 463, 2623

Ibata, R., Malhan, K., Tenachi, W., et al. 2024, ApJ, 967, 89

Iocco, F., Pato, M., & Bertone, G. 2015, Nat. Phys., 11, 245

Jiao, Y., Hammer, F., Wang, H., et al. 2023, A&A, 678, A208

Karamanis, M., Beutler, F., & Peacock, J. A. 2021, MNRAS, 508, 3589

Khelashvili, M., Rudakovskyi, A., & Hossenfelder, S. 2024, ArXiv e-prints [arXiv:2401.10202]

Koop, O., Antoja, T., Helmi, A., Callingham, T. M., & Laporte, C. F. P. 2024, A&A, 692, A50

Kroupa, P., Pflamm-Altenburg, J., Mazurenko, S., et al. 2024, ApJ, 970, 94

Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571

Li, P., McGaugh, S. S., Lelli, F., Schombert, J. M., & Pawlowski, M. S. 2022, A&A, 665, A143

Lin, H.-N., & Li, X. 2019, MNRAS, 487, 5679

López-Corredoira, M. 2025, ApJ, 978, 45

López-Corredoira, M., & Betancort-Rijo, J. E. 2021, ApJ, 909, 137

McGaugh, S. S. 2004, ApJ, 609, 652

McGaugh, S. S. 2018, RNAAS, 2, 156

McGaugh, S. S. 2019, ApJ, 885, 87

Milgrom, M. 1983, ApJ, 270, 365

Milgrom, M. 2015, Can. J. Phys., 93, 107

Milgrom, M. 2023, Phys. Rev. D, 108, 084005

Misiriotis, A., Xilouris, E. M., Papamastorakis, J., Boumis, P., & Goudis, C. D. 2006, A&A, 459, 113

Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533

Mroz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Ou, X., Eilers, A.-C., Necib, L., & Frebel, A. 2024, MNRAS, 528, 693

Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27

Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66

Refregier, A. 2003, ARA&A, 41, 645

Roberts, M. S. 1975, in Dynamics of the Solar Systems, ed. A. Hayli, IAU Symp., 69, 331

Robin, A. C., Bienaymé, O., Salomon, J. B., et al. 2022, A&A, 667, A98

Salucci, P., Lapi, A., Tonini, C., et al. 2007, MNRAS, 378, 41

Sanders, R. H. 2003, MNRAS, 342, 901

Shao, M. J., Anbajagane, D., & Chang, C. 2023, MNRAS, 523, 3258

Skordis, C., & Złosnik, T. 2021, Phys. Rev. Lett., 127, 161302

Sylos Labini, F. 2024, ApJ, 976, 185

Sylos Labini, F., Chrobáková, Ž., Capuzzo-Dolcetta, R., & López-Corredoira, M. 2023, ApJ, 945, 3

Vallenari, A., Brown, A., & Prusti, T. 2022, A&A, 674, A1

Wang, H.-F., Chrobáková, Ž., López-Corredoira, M., & Sylos Labini, F. 2022, ApJ, 942, 12

Zhou, Y., Li, X., Huang, Y., & Zhang, H. 2023, ApJ, 946, 73

Zobnina, D. I., & Zasov, A. V. 2020, Astron. Rep., 64, 295