| Issue |
A&A
Volume 703, November 2025
|
|
|---|---|---|
| Article Number | L17 | |
| Number of page(s) | 10 | |
| Section | Letters to the Editor | |
| DOI | https://doi.org/10.1051/0004-6361/202557021 | |
| Published online | 25 November 2025 | |
Letter to the Editor
The odd primordial halo of the Milky Way implied by Gaia: A shallow core, but a steep decline
1
School of Astronomy and Space Science, Nanjing University, Nanjing, Jiangsu 210023, China
2
Department of Astronomy, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA
3
Leibniz-Institute for Astrophysics, An der Sternwarte 16, 14482 Potsdam, Germany
4
Observatoire de Paris, Université PSL, CNRS, Place Jules Janssen, Meudon 92195, France
5
Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85722, USA
⋆ Corresponding author: pli@nju.edu.cn
Received:
28
August
2025
Accepted:
3
November
2025
Primordial dark matter halos are well understood from cold dark matter-only simulations. Since they can contract significantly as baryons settle into their centers, direct comparisons with observed galaxies are complicated. We present an approach to reversing the halo contraction by numerically calculating the halo response to baryonic infall and iterating the initial condition. This allowed us to derive spherically averaged primordial dark matter halos for observed galaxies. We applied this approach to the Milky Way and found that the latest Gaia measurements for the rotation velocities imply an odd primordial Galactic halo: Its concentration and total mass differ by more than 3σ from the predictions, and the density profile presents an inner core that is too shallow and an outer decline that is too steep to be compatible with the cold dark matter paradigm.
Key words: Galaxy: general / Galaxy: kinematics and dynamics / Galaxy: structure / dark matter
© The Authors 2025
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Disk galaxies have often been used to study the nature of dark matter because their rotation velocities are excellent indicators of their mass distributions (e.g. see Li et al. 2020). With the advent of the European Space Agency cornerstone project Gaia, astronomers can measure the parallaxes and proper motions for millions of nearby stars with high accuracy (Gaia Collaboration 2023a). In this way, it becomes possible to measure the circular Galactic velocity using the 6D phase-space information of individual stars, which cannot be achieved for external galaxies. This makes our Galaxy a unique laboratory for studying dark matter.
The rotational velocity profile of the Milky Way (MW) has been known to decline at large radii (Eilers et al. 2019; Ou et al. 2024), and the latest measurements using Gaia DR3 even reported that the decline is Keplerian (Jiao et al. 2023; Hammer et al. 2024). This implies that the density of its dark matter halo is extremely low at large radii. The halo density profile is not the same as that of the primordial halo, however, because the latter evolves significantly due to baryonic feedback and compression. It is important to note that the nature of dark matter determines the structure of the pristine halo, which is free of baryonic effects and has been well established with dark matter-only simulations (Navarro et al. 1996, 2004). For massive galaxies such as the Milky Way, studies have shown that stellar feedback has a negligible effect (Di Cintio et al. 2014; Read et al. 2016), while baryonic compression leads to significant halo contraction (Li et al. 2022a,b). The halo contraction can be well modeled and numerically calculated via an adiabatic process (Sellwood & McGaugh 2005; Gnedin et al. 2011; Piffl et al. 2015; Binney & Piffl 2015).
We present an approach for inferring the primordial halos by implementing baryonic compression in the Markov chain Monte Carlo (MCMC) approach (Foreman-Mackey et al. 2013a). We apply this approach to the Milky Way in this paper.
2. Numerical modeling and Galactic data
We used the algorithm proposed by Young (1980) to model the adiabatic contraction of dark matter halos. The main idea was to make use of the conservation of three adiabatic actions (the third action vanishes due to the spherical symmetry, so that only two actions remain, one radial and one azimuthal), and write the distribution function of dark matter particles in terms of the actions, so that they were invariant as the halo contracts. The baryonic potential was then gradually added to the total potential, which led to a gradual change in the density profile. A realization of this algorithm was developed originally as part of the N-body simulation code GALAXY (Sellwood 2014), and it was applied to observational data for the first time by Sellwood & McGaugh (2005). We further developed it to be suitable for diverse observational data and iterable for an MCMC implementation. The MCMC implementation is critical because it allowed us to sample and choose primordial halos based on observed rotation curves. Therefore, we were able to directly link galaxy rotation curves to primordial dark matter halos, which are more comparable to dark matter-only simulations. The code is now made independent as COMPRESS and is released at GitHub1. The MCMC implementation is also included as an example to show its usage. We provide an in-depth description of the code in Appendix A.1.
As an application, we study the primordial halo of the Milky Way in this work. Unlike external galaxies, for which modeling baryonic mass distributions is hampered by the uncertainties on the stellar mass-to-light ratios, our Galaxy has more constraints on the stellar mass, such as the vertical kinematics of local stars (Bovy & Rix 2013) and microlensing (Wegg et al. 2016). We therefore did not treat baryonic mass distributions as free parameters. Instead, we adopted a wide variety of popular baryonic models from the literature to ensure that our results are robust and independent of the choice of models. In total, we adopted 12 baryonic models (for more details, see Appendix A.2).
The circular Galactic velocity has been extensively measured using the data from Gaia DR3 with improved parallaxes and proper motions (e.g. Gaia Collaboration 2023b; Wang et al. 2023; Zhou et al. 2023; Jiao et al. 2023; Ou et al. 2024). The measured velocity profiles from different groups are reasonably consistent with each other, with only small differences at small and large radii. In particular, the measurement by Jiao et al. (2023) considered systematic uncertainties due to the neglected cross-term in the Jeans equation, the uncertainty on the disk scale length, variations in the stellar density profile, and azimuthally varying stellar samples. Their error estimates are therefore relatively more reasonable. We therefore adopted their measurements. The interpretation of these data remains open to debate (Koop et al. 2024; Hammer et al. 2024; Ou et al. 2025). We investigate here the primordial halo that is implied by the current Gaia data when baryonic compression is considered.
To parameterize the primordial halo, we initially used the NFW model (Navarro et al. 1996) because it is a prediction of the standard CDM model. However, the NFW model does not fit the Galactic rotation for any baryonic model because the Gaia velocity profile declines fast at large radii and is steeper than the outer density slope of the NFW model (see the fits in Fig. A.3). We therefore employed the more flexible Einasto profile (Einasto 1965), which has an additional shape parameter α. Navarro et al. (2004) found that the Einasto profile with a shape parameter α ∼ 0.17 fits their simulated dark matter halos better than the NFW model. This shows that parameterizing primordial halos using the Einasto model can recover simulated halos if the data agree with the CDM model. We used the code EMCEE (Foreman-Mackey et al. 2013b) to map the posterior distributions of the primordial halo parameters. The detailed setup is given in Appendix A.3.
3. Results
Figure 1 shows the rotation curve fit using the Einasto model (Einasto 1965) for primordial halos and the McGaugh19 model (McGaugh 2019) for baryonic mass distributions. For all the adopted baryonic models, the Gaia rotation velocity profile fit well with compressed Einasto halos (see Fig. 1). The results are summarized in Table A.1. Since the rotation curve data by Jiao et al. (2023) only extend to ≥10 kpc, we also plot the inner measurements from Gaia DR2 (Eilers et al. 2019) and the VVV survey (Portail et al. 2017) as a reference, although these data were not fit. The extrapolations of the best Einasto fit can also describe the inner rotation velocity when the McGaugch19 model for baryons is used. This is not always true for other models, however. The primordial and compressed halos contribute rather different rotation velocities. The difference is more profound from 5 to 20 kpc. This shows that the halo structure within 20 kpc changes dramatically after baryonic compression.
![]() |
Fig. 1. Example of a circular velocity fit using the McGaugh19 model for baryonic mass distributions. The purple, blue, and green lines represent the contributions of the bar, disk, and gas components, respectively. The solid and dashed black lines show the current and primordial dark matter halos, respectively. The solid red line indicates the total velocity profile. The black points show the latest Gaia measurements (Jiao et al. 2023), and the gray upward triangles and squares show the terminal velocities from (McClure-Griffiths & Dickey 2007, 2016), and Portail et al. (2017), respectively. The data marked with open symbols were not fit because they do not consider the systematic uncertainties. The fits with 12 baryonic models are shown in Fig. A.2. |
Figure 2 shows the total halo masses and concentrations of the primordial halos using different baryonic models (also see Table A.1). For comparison, we also present the halos from the direct fits that did not implement adiabatic contraction. These are simply the current halos and are same as the compressed halos when adiabatic contraction is implemented in our fitting procedure because we fit the same data. Fits without baryonic compression were performed before, and our results are roughly consistent with the results by Ou et al. (2024). Regardless of whether baryonic compression is included, the results assuming different baryonic models are fairly consistent, which ensures that our results are model independent. Baryonic compression does not change the inferred total halo mass because it only redistributes dark matter particles. The total halo mass ranges from 1.09 to 1.42 × 1011 M⊙, which is consistent with previous findings (Jiao et al. 2023). They are lower by more than 3σ than the predicted values from the abundance-matching relation, however (Moster et al. 2013). This correlation links the predicted halo mass function to observed stellar mass functions. Baryonic compression has a considerable effect on the halo concentration. Neglecting it leads to significantly overestimated halo concentrations. The inferred primordial halos therefore have rather low concentrations, which strongly contradicts the predictions from N-body simulations with cold dark matter (Dutton & Macciò 2014).
![]() |
Fig. 2. Halo masses and concentrations of the primordial Galactic halos derived from the Gaia circular velocity fits using 12 baryonic models. The red and blue stars with errors represent the halos with and without adiabatic contraction, respectively. The predicted halo mass-concentration relation within 1σ from simulations (Dutton & Macciò 2014) is shown as the declining band. The vertical band shows the expected range of the MW halo mass according to the abundance-matching relation (Moster et al. 2013). The upper and lower limits are set by the highest stellar mass (model A&S) plus 1σ and the lowest stellar mass (model I) minus 1σ, respectively. |
The structures of primordial and current halos are quite different (Fig. 3). Baryonic compression significantly increases the density of the inner halo. Thus, the inner structure of the inferred primordial halos is rather flat. As a trade-off, the dark matter density decreases more quickly at large radii, which slightly helps to explain the quickly declining circular velocity from Gaia (Jiao et al. 2023). CDM-only simulations predict a universal cuspy profile, which can be described by the Einasto model with a shape parameter of α ∼ 0.17 (Navarro et al. 2004). In logarithmic space, the density profile of CDM halos is almost linear in the considered radial range (Fig. 3). In contrast, the structures of primordial and current halos are apparently curved. Their inner halo structures have a shallow core, which suggests a core-cusp problem. This is a classical problem in dwarf galaxies (Oh et al. 2011). It is usually thought that massive galaxies like the Milky Way host cuspy halos, which is consistent with CDM predictions. The Gaia rotation velocities suggest the opposite, however. When baryonic compression is considered, the core of the inferred primordial halo is significantly shallower, which exacerbates the core-cusp problem. The outer density profiles decline steeply as a result from the declining circular velocity profiles. Although including baryonic compression helps to reduce the dark matter density in the outer region, the outer decline of the inferred primordial halo is still too steep (Gaussian-like) compared to that of CDM halos.
![]() |
Fig. 3. Structure of the inferred primordial and current Galactic halos, along with predictions for the cold and warm dark matter. The density profiles are scaled so that there is no need to assume or consider the masses or concentrations for these halos. The gray band indicates the range of the current halos derived from the Gaia velocity fits using the 12 baryonic models, and the red band shows their corresponding primordial halos within 1σ. The blue band presents the simulated halos with cold dark matter only (Dutton & Macciò 2014). The purple band shows the warm dark matter halos (normalized to match the primordial Galactic halo) with a core size spanning from 4.56 kpc (WDM5 in Macciò et al. 2012) to 7.0 kpc, corresponding to a particle mass of 0.05 keV and lower. |
A popular solution for the core-cusp problem is to introduce baryonic feedback, i.e. some physical processes occurring during galaxy formation that drive outflows and help flatten the inner structure of the CDM halo. Cole & Binney (2017) showed that a centrally heated halo can have a shallow core in their numerical approach based on distribution functions. The question is whether the amount of energy is sufficient to heat the central dark matter particles. Unlike dwarf galaxies, in which stellar feedback (e.g., supernova explosion and stellar winds) can efficiently transform primordial cusps to cores, our Milky Way is such a massive galaxy that more powerful feedback is required. The only known powerful feedback comes from the accretion of matter onto black holes in galaxy nuclei, known as active galactic nuclei (AGNs; Magorrian et al. 1998; Silk & Rees 1998). AGN feedback has been introduced in the Illustris TNG simulations (Weinberger et al. 2017; Pillepich et al. 2018). It was found, however, that the CDM halos hosted by massive galaxies in TNG have an inner density slope of ∼2 (Wang et al. 2020; Pillepich et al. 2024), which is steeper than that of their primordial counterparts. This is roughly consistent with the finding that only considered baryonic compression (Li et al. 2022a). This shows that AGN feedback as implemented in the TNG simulations is not sufficiently powerful to counteract the effect of baryonic compression in the inner region. In addition, the outer halo structure is insensitive to baryonic feedback. Baryon-driven outflows can cause the outer decline to be slightly shallower at most (Di Cintio et al. 2014), and this trend is opposite to what is required by the quickly declining rotation velocities.
Lighter dark matter has also been proposed to generate cored halos, such as warm dark matter (WDM). In Fig. 3 we compare our results with the WDM simulations of the MW mass (Macciò et al. 2012). Since the primordial Galactic halo we inferred has a large core, we only show the halo with the lowest particle mass (0.05 keV; i.e., WMD5 in Macciò et al. 2012) corresponding to a core size of 4.56 kpc (the lower limit of the purple band), which is the largest core in their simulations, but is still smaller than that of the MW halo. To reproduce the inner structure of the Galactic halo, a core size of ∼7 kpc is required (the upper limit of the purple band), but this would also extend the halo more, which contradicts the steep decline suggested by Gaia. Moreover, to generate a larger core, the WDM mass would further need to be reduced to below 0.05 keV. As pointed out by Macciò et al. (2012), a particle mass lower than 0.1 keV would prevent the formation of dwarf galaxies, causing a catch-22 problem.
Fuzzy dark matter (Hu et al. 2000) avoids this problem by its nonrelativistic Bose-Einstein condensation (Elgamal et al. 2024), and its quantum stress helps to generate large soliton cores. The inner halo structure might thus be roughly consistent with that of the Milky Way, although we cannot make a robust comparison because our approach only considers gravity, while the halo center is subject to strong quantum stress. In the outer region, however, the halo structure is clearly granular, which is indistinguishable from the CDM halos (Schive et al. 2014, 2014). Fuzzy dark matter can therefore hardly reproduce the steep decline of the Galactic halo.
The discrepancies in the inner and outer regions therefore cannot be simultaneously reconciled by simply twisting the nature of dark matter because the structure of the primordial Galactic halo is rather odd: The core is too shallow, and the decline is too steep. Generally speaking, a steep decline suggests that the halo is compact, so that a central cusp is expected, while a shallow core implies that dark matter particles are difficult to aggregate, so that the halo is expected to be extended, corresponding to a slow decline. The shallow core and the steep decline are not expected to be present simultaneously, however, because they imply opposite natures of dark matter. The fact that they do coexist might suggest that some exotic nature of dark matter and more powerful baryonic feedback need to be combined: A steep decline can appear if dark matter particles can gather more closely as a result of some unknown nature; and a shallow core can be generated if the feedback is more powerful. More theoretical and simulation work is required to determine whether it is possible to reproduce the too-shallow-then-too-steep structure.
4. Discussion and conclusion
We demonstrated the importance of considering baryonic compression when studying the nature of dark matter with galaxy dynamics. Using an innovative technique, we inferred the primordial dark matter halo of the Milky Way from the Gaia rotation curve. Jiao et al. (2023) showed that the declining Galactic rotation curve leads to a significantly lower dark matter halo mass than is expected from the CDM paradigm, more specifically, the abundance-matching relation (e.g. Moster et al. 2013; Behroozi et al. 2013; Kravtsov et al. 2018). We focused on the detailed structure of the primordial halo, which is directly related to the nature of dark matter. We found that the structure of the primordial Galactic halo is too shallow and then too steep, which is inconsistent with the standard CDM paradigm built from external galaxies and the large-scale structure. If our Milky Way is not special, the nature of dark matter and the model of galaxy formation might need to be adjusted.
The discrepancy at small radii was first indicated by Binney & Piffl (2015), who found that an adiabatically compressed NFW halo cannot reconcile various observational constraints, such as microlensing (Popowski et al. 2005). Their focus was on the current configuration of the Milky Way, while we advanced the technique to reverse the adiabatic contraction and thereby infer the primordial halo. Cautun et al. (2020) also investigated the contraction of the Galactic dark matter halo. They selected MW-mass halos from various hydrodynamical and dark matter-only (DMO) simulations. By comparing these halos from the two sets of simulations, they were able to construct a mean mass profile for the compressed halo that was then used to fit the Gaia rotation curve (Eilers et al. 2019). Since their simulations assumed cold dark matter, the halos from the DMO simulations are mostly NFW. Initially, we assumed that primordial halos are NFW as well, but we found that the slowly decreasing density profile of NFW halos did not fit the quickly declining Gaia rotation curve at large radii. This is also evident from Fig. 10 in Cautun et al. (2020) and similarly in Binney & Vasiliev (2023, see Fig. 15). We chose a more flexible Einasto profile for primordial halos. For this fundamental reason, they did not find the discrepancy we identified.
It is important to stress that we relied upon the rotation curve deduced from the Gaia data. Since different groups found consistent declining Galactic rotation curves using the Gaia data, the problem we identified is not specific to the work of Jiao et al. (2023), but typical to the Gaia data. We therefore conclude that these results from Gaia differ considerably from the ΛCDM expectation for our Milky Way.
Acknowledgments
We thank Jo Bovy for assistant in using Galpy. P. L. is supported by the startup funding from Nanjing University. M. S. P. acknowledges funding of a Leibniz-Junior Research Group (project number J94/2020).
References
- Allen, C., & Santillan, A. 1991, Rev. Mex. Astron. Astrofis., 22, 255 [NASA ADS] [Google Scholar]
- Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 [NASA ADS] [CrossRef] [Google Scholar]
- Binney, J., & Piffl, T. 2015, MNRAS, 454, 3653 [NASA ADS] [CrossRef] [Google Scholar]
- Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, USA: Princeton University Press) [Google Scholar]
- Binney, J., & Vasiliev, E. 2023, MNRAS, 520, 1832 [NASA ADS] [CrossRef] [Google Scholar]
- Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27 [Google Scholar]
- Bovy, J. 2015, ApJS, 216, 29 [NASA ADS] [CrossRef] [Google Scholar]
- Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115 [Google Scholar]
- Calchi Novati, S., & Mancini, L. 2011, MNRAS, 416, 1292 [NASA ADS] [CrossRef] [Google Scholar]
- Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291 [Google Scholar]
- Cole, D. R., & Binney, J. 2017, MNRAS, 465, 798 [Google Scholar]
- de Jong, J. T. A., Yanny, B., Rix, H.-W., et al. 2010, ApJ, 714, 663 [NASA ADS] [CrossRef] [Google Scholar]
- de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, JCAP, 2019, 037 [Google Scholar]
- Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014, MNRAS, 441, 2986 [Google Scholar]
- Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359 [Google Scholar]
- Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120 [Google Scholar]
- Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87 [NASA ADS] [Google Scholar]
- Elgamal, S., Nori, M., Macciò, A. V., Baldi, M., & Waterval, S. 2024, MNRAS, 532, 4050 [Google Scholar]
- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013a, PASP, 125, 306 [Google Scholar]
- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013b, PASP, 125, 306 [Google Scholar]
- Gaia Collaboration(Vallenari, A., et al.) 2023a, A&A, 674, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Gaia Collaboration (Drimmel, R., et al.) 2023b, A&A, 674, A37 [CrossRef] [EDP Sciences] [Google Scholar]
- Gnedin, O. Y., Ceverino, D., Gnedin, N. Y., et al. 2011, ArXiv e-prints [arXiv:1108.5736] [Google Scholar]
- Hammer, F., Jiao, Y. J., Mamon, G. A., et al. 2024, A&A, 692, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Hernquist, L. 1990, ApJ, 356, 359 [Google Scholar]
- Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [Google Scholar]
- Iocco, F., Pato, M., & Bertone, G. 2015, Nat. Phys., 11, 245 [NASA ADS] [CrossRef] [Google Scholar]
- Jiao, Y., Hammer, F., Wang, J. L., & Yang, Y. B. 2021, A&A, 654, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Jiao, Y., Hammer, F., Wang, H., et al. 2023, A&A, 678, A208 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [Google Scholar]
- Klypin, A., Zhao, H., & Somerville, R. S. 2002, ApJ, 573, 597 [Google Scholar]
- Koop, O., Antoja, T., Helmi, A., Callingham, T. M., & Laporte, C. F. P. 2024, A&A, 692, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astron. Lett., 44, 8 [Google Scholar]
- Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
- Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, ApJS, 247, 31 [Google Scholar]
- Li, P., McGaugh, S. S., Lelli, F., Schombert, J. M., & Pawlowski, M. S. 2022a, A&A, 665, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Li, P., McGaugh, S. S., Lelli, F., et al. 2022b, ApJ, 927, 198 [NASA ADS] [CrossRef] [Google Scholar]
- Macciò, A. V., Paduroiu, S., Anderhalden, D., Schneider, A., & Moore, B. 2012, MNRAS, 424, 1105 [CrossRef] [Google Scholar]
- Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285 [Google Scholar]
- McClure-Griffiths, N. M., & Dickey, J. M. 2007, ApJ, 671, 427 [CrossRef] [Google Scholar]
- McClure-Griffiths, N. M., & Dickey, J. M. 2016, ApJ, 831, 124 [NASA ADS] [CrossRef] [Google Scholar]
- McGaugh, S. S. 2019, ApJ, 885, 87 [Google Scholar]
- McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
- McMillan, P. J. 2017, MNRAS, 465, 76 [NASA ADS] [CrossRef] [Google Scholar]
- Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
- Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121 [Google Scholar]
- Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996, MNRAS, 283, L72 [NASA ADS] [Google Scholar]
- Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039 [Google Scholar]
- Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, R. C., Jr 2011, AJ, 141, 193 [NASA ADS] [CrossRef] [Google Scholar]
- Ou, X., Eilers, A.-C., Necib, L., & Frebel, A. 2024, MNRAS, 528, 693 [NASA ADS] [CrossRef] [Google Scholar]
- Ou, X., Necib, L., Wetzel, A., et al. 2025, ArXiv e-prints [arXiv:2503.05877] [Google Scholar]
- Piffl, T., Penoyre, Z., & Binney, J. 2015, MNRAS, 451, 639 [NASA ADS] [CrossRef] [Google Scholar]
- Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077 [Google Scholar]
- Pillepich, A., Sotillo-Ramos, D., Ramesh, R., et al. 2024, MNRAS, 535, 1721 [NASA ADS] [CrossRef] [Google Scholar]
- Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
- Popowski, P., Griest, K., Thomas, C. L., et al. 2005, ApJ, 631, 879 [NASA ADS] [CrossRef] [Google Scholar]
- Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017, MNRAS, 465, 1621 [NASA ADS] [CrossRef] [Google Scholar]
- Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573 [NASA ADS] [CrossRef] [Google Scholar]
- Schive, H.-Y., Chiueh, T., & Broadhurst, T. 2014, Nat. Phys., 10, 496 [NASA ADS] [CrossRef] [Google Scholar]
- Schive, H.-Y., Liao, M.-H., Woo, T.-P., et al. 2014, Phys. Rev. Lett., 113, 261302 [NASA ADS] [CrossRef] [Google Scholar]
- Sellwood, J. A.. 2014, ArXiv e-prints [arXiv:1406.6606] [Google Scholar]
- Sellwood, J. A., & McGaugh, S. S. 2005, ApJ, 634, 70 [NASA ADS] [CrossRef] [Google Scholar]
- Silk, J., & Rees, M. J. 1998, A&A, 331, L1 [NASA ADS] [Google Scholar]
- Stanek, K. Z., Udalski, A., Szymański, M., et al. 1997, ApJ, 477, 163 [NASA ADS] [CrossRef] [Google Scholar]
- Wang, Y., Vogelsberger, M., Xu, D., et al. 2020, MNRAS, 491, 5188 [Google Scholar]
- Wang, H.-F., Chrobáková, Ž., López-Corredoira, M., & Sylos Labini, F. 2023, ApJ, 942, 12 [NASA ADS] [CrossRef] [Google Scholar]
- Wegg, C., Gerhard, O., & Portail, M. 2016, MNRAS, 463, 557 [CrossRef] [Google Scholar]
- Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291 [Google Scholar]
- Young, P. 1980, ApJ, 242, 1232 [NASA ADS] [CrossRef] [Google Scholar]
- Zhou, Y., Li, X., Huang, Y., & Zhang, H. 2023, ApJ, 946, 73 [CrossRef] [Google Scholar]
Appendix A: Method and data
A.1. Modeling halo contraction
Primordial dark matter halos are self-supported and free of baryonic effects. They start to contract as baryonic gas is accreted into their centers. We modeled the process with the code COMPRESS. The code starts with setting up an isolated, extensive and massive spherical initial halo, which must have either a known equilibrium distribution function (DF) or an isotropic one that can be derived by Eddington inversion (Binney & Tremaine 2008). Since the central attraction is changed by embedding a disk and/or bulge at the center of the halo, we need to derive a revised mass profile for the halo in the composite model. Following Young (1980), we assume that the masses of the disk and bulge were increased adiabatically from zero to their present values, an assumption that does not change the actions of the halo particles (Binney & Tremaine 2008). It was previously assumed angular momentum, one of the actions of an orbit (Jϕ = 2πL), was alone conserved (Blumenthal et al. 1986; Klypin et al. 2002), so that the formula would apply only if the orbits of all particles were initially, and remained, precisely circular. The orbits of particles in all reasonable spherical models librate radially, and therefore one must take conservation of radial action (Jr = 2∫rminrmaxvrdr with vr being the radial velocity) into account when computing the density response to adiabatic changes to the potential. The pressure of radial motions makes a realistic halo more resistant to compression than the naive model with no radial action would predict. Given the conservation of the number of particles and the invariance of adiabatic actions, the DF, when written in terms of actions, is the same after the adiabatic change as before, i.e. f0(Jr, Jϕ) = fn(Jr, Jϕ).
As described by Young (1980) and by Sellwood & McGaugh (2005), a first revised density profile of the halo can be obtained from the original DF by assuming it is an unchanged function of the actions in the new combined potential of the disk and original halo,
where the distribution function f(E, L) is converted from the invariant f(Jr, Jϕ = 2πL) using the function Jr(E, L), and the total potential Φ(r) is updated with the addition of baryonic contribution, i.e. Φ = Φ0 + Φbar. This is the driven force that changes the dark matter density function ρ(r). We then use the first revised density profile of the halo to compute a second revised halo density profile from the DF by again assuming it is an unchanged function of the actions in the new combined potential of the disk and the first revised halo. We iterate in this manner until the change to the halo density between iterations is acceptably small. An initially isotropic DF may become mildly radially biased, which can still be represented by the unchanged actions. Unlike the approximate power-law fit (Gnedin et al. 2011), our procedure gives the detailed mass profile of the compressed halo that results from the adopted initial halo and the observed disk and bulge, which we need to fit the extensive rotation curve.
Note that our procedure assumes the potential remains spherically symmetric, so we must approximate the disk potential by the monopole only term, i.e. using only the disk mass enclosed in a sphere of radius r. From a comparison with a simulation in which a disk was grown slowly inside a spherical halo, the aspherical part of the disk potential caused negligibly small changes to the spherically averaged halo potential. The adiabatic modeling of halo contraction has been tested using N-body simulations in (Sellwood & McGaugh 2005). The spherically averaged density profiles agree remarkably well (see their Fig. 1).
It is important to note that our procedure does take account of adiabatic baryonic feedback, that may rearrange the disk mass through galactic fountains, etc. Of course, explosive feedback that impulsively ejects a significant disk mass fraction is a non-adiabatic change that contradicts our working assumption; such behavior might occur during the assembly of dwarf galaxies, but not those of the mass of the Milky Way.
Note that the adiabatic compression procedure is efficient and runs quite fast, so that we can implement it into the MCMC approach. This way, we randomly sample a large pool of primordial halos with Monte Carlo, and calculate their contracted halos with COMPRESS. The contracted halos are then compared to the Gaia rotation velocities to build the likelihood function and calculate their probabilities. We then run the Markov chain to generate new pools of primordial halos with larger probabilities, until we obtain a sufficient number of stable pools of primordial halos. Eventually, we select the primordial halo with the highest probability, which will match the Gaia circular velocities the best. This implementation improves over the previous approach (Li et al. 2022a), which essentially uses a single walker to find the best initial halo. The latter is faster, but it cannot estimate the uncertainties of the fitting parameters. Combining the numerical modeling and the MCMC approach, we can infer the primordial halos from observations.
A.2. Baryonic models of the Milky Way
We adopted 12 baryonic models in this paper, including Model A&S (Allen & Santillan 1991), Model I (Pouliasis et al. 2017), Bovy15 (Bovy 2015), McMillan17 (McMillan 2017), B2 (de Salas et al. 2019), and McGaugh19 (McGaugh 2019) as detailed below.
Model A&S The classical Model A&S models a single disk using the Miyamoto–Nagai potential (Miyamoto & Nagai 1975), given by
where the parameters are set as M = 8.561 × 1010 M⊙, a = 5.3178 kpc, b = 0.25 kpc. The bar is approximated as a Plummer sphere (Plummer 1911),
where the total mass Mbar = 1.406 × 1010 M⊙ and b = 0.3873 kpc.
Bovy15 The Bovy15 model also adopts the Miyamoto–Nagai disk (Eq. A.2) with the total mass M = 6.8 × 1010 M⊙, a = 3.0 kpc, b = 0.28 kpc. It models the bar using a power-law density profile with an exponential cutoff,
where the scale density ρ0 = 0.005858 M⊙ pc−3 corresponding to the total bar mass of 0.5 × 1010 M⊙, α = 1.8, rc = 1.9 kpc, and r1 = 8 kpc.
Model I Model I improves over Model A&S by including two Miyamoto–Nagai disks (Eq. A.2), a thin disk with athin = 5.3 kpc, bthin = 0.25 kpc, and a thick disk with athick = 2.6 kpc, bthick = 0.8 kpc. Both have the same total stellar mass M = 3.944 × 1010 M⊙. The Plummer model (Eq. A.3) is also used to approximately describe the bar: Mbar = 1.067 × 1010 M⊙ and b = 0.3 kpc.
McMillan17 The MCMillan17 model includes thin and thick disks and models their mass density with a double exponential function,
where M, L, and H are the total stellar mass, the scale length, and the scale height, respectively. The values that match various observational constraints are: Mthin = 3.5 × 1010 M⊙, Lthin = 2.50 kpc, and Hthin = 0.3 kpc for the thin disk; Mthick = 1.0 × 1010 M⊙, Lthick = 3.02 kpc, and Hthick = 0.9 kpc for the thick disk. The bar is modeled with an axisymmetric function,
where α = 1.8, r0 = 0.075 kpc, rcut = 2.1 kpc, the axis ratio q = 0.5, and ρ0 = 98.4 M⊙ pc−3. McMillan17 also models gas distributions using an exponential function with a central hole,
For the H I disk, Σ0 = 53.1 M⊙ pc−2, Rd = 7 kpc, Rm = 4 kpc, zd = 0.085 kpc; while for the molecular gas disk, Σ0 = 2180 M⊙ pc−2, Rd = 1.5 kpc, Rm = 12 kpc, zd = 0.045 kpc. The contributions of each component to circular velocity and surface brightness are calculated using the open software GalPot2.
B2 The B2 model adopts a single double exponential disk (Eq. A.5) with M = 3.65 × 1010 M⊙, L = 2.35 kpc, H = 0.14 kpc. It also includes components for warm and cold dust, H I and molecular gas. These components are modeled with the double exponential model as well, and the parameters are given by
The units for mass and scale height/lenth are M⊙ and kpc, respectively. The bar is modeled using the Hernquist potential (Hernquist 1990),
where Mbar = 1.55 × 1010 M⊙ and rb = 0.7 kpc.
McGaugh19 Unlike previous models, the McGaugh19 model does not assume smooth functions, but numerically maps the azimuthally averaged surface mass densities and the contributions to circular velocity of the stellar disk, the gas disk, and the bar by applying the radial acceleration relation (RAR; McGaugh et al. 2016; Lelli et al. 2017) to measured circular velocities. This approach is based on the fact that the RAR is shown to work in individual disk galaxies (Li et al. 2018). As such, the resulting mass distributions have bumps and wiggles that are typically present in external galaxies.
Combinations of bar and disk models Following the literature (Iocco et al. 2015; Jiao et al. 2021), we also consider baryonic models from the combinations of bar models and disk models. Two triaxial bar models are considered: E2 (Stanek et al. 1997) defined as
where xb = 0.899 kpc, yb = 0.386 kpc, zb = 0.250 kpc, and the scale density is determined from the total mass Mbar = 2.41 × 1010 M⊙ using Mbar = 8πxbybzbρ0; G2 (Stanek et al. 1997) given by
where xb = 1.239 kpc, yb = 0.609 kpc, zb = 0.438 kpc, and ρ0 can be calculated from the total mass Mbar = 2.12 × 1010 M⊙ via Mbar = 6.57252 πxbybzbρ0. Although we adopt triaxial models for the bar, we calculate their azimuthally averaged surface mass density and contributions to circular velocity. This is because the bar is rotating. Averaging azimuthally is equivalent to averaging over time. The disk models are CM (Calchi Novati & Mancini 2011), J (Jurić et al. 2008) and dJ (de Jong et al. 2010). They all model thin and thick disks using the double exponential function (Eq. A.5), and the parameters are given as
The units for the total stellar mass (M) and the scale length/height (L/H) are 1010 M⊙ and kpc, respectively. Combining these bar and disk models, we obtain 6 additional baryonic models: E2CM (E2 bar + CM disk), E2J (E2 bar + J disk), E2dJ (E2 bar + dJ disk), G2CM (G2 bar + CM disk), G2J (G2 bar + J disk), and G2dJ (G2 bar + dJ disk). In total, we adopt 12 baryonic models in this paper. Except for McMillan17 and McGaugh19, we calculate their surface mass density and their contributions to circular velocity using the open Python package galpy3 (Bovy 2015). Those data are publicly available on GitHub4.
A.3. Halo model and setup
Inferred properties of the primordial halos parameterized with the Einasto profile (see Appendix A.3) before adiabatic contraction assuming a variety of baryonic models (see Appendix A.2).
We adopt the Einasto model (Einasto 1965; Navarro et al. 2004) for primordial halos. Its density profile is given by
where the shape parameter α, scale density ρs and scale radius rs are free parameters. Its cumulative mass is hence
where x = r/rs, Γ(a, x) = ∫0xta − 1e−tdt is the incomplete Gamma function. The total halo mass M200 and concentration C200 are defined as
where r200 is the radius that encloses a mean halo density of 200 times the critical density of the Universe, and V200 is the rotation velocity due to the dark matter halo at r200. The velocity profile is then given by
The total velocity is the summation of baryonic contributions and the dark matter contribution,
where Vbar, Vdisk, and Vgas are the rotation velocity contributed by the bar, the disk, and the gas (including dusts) components, respectively; VEIN, compress is the contribution from the compressed Einasto halo. Those baryonic components are fixed for each given baryonic model, and we consider 12 baryonic models from the literature as described in A.2. Therefore, the free parameters are only on the halo, i.e. V200, C200, and α.
We use EMCEE (Foreman-Mackey et al. 2013a) to map the posterior distributions of the fitting parameters and impose flat priors within the following hard boundaries: 0.01 < α < 4.0, 5.0 < C200 < 30.0 and 50 < V200 < 200 km/s. We checked and verified these ranges are sufficiently wide. For comparison, we fit the circular velocities without implementing adiabatic contraction of dark matter halos. The fit results and posterior distributions of the free parameters are shown in Fig. A.1. When implementing baryonic compression, the computational time increased dramatically. For the Einasto model, one has to generate different input tables for different values of the shape parameter α. This would require significantly more computational time. For simplicity, we fix the value of α for each baryonic model using the fitting results without implementing halo contraction. This simplification is supported by the fact that the declining part of the circular velocity profile from Gaia strongly constrains the outer density slope of the current halo, which is insensitive to baryonic compression. As such, including halo contraction or not does not significantly influence the value of α. The decent fits in Fig. A.2 further justify the simplification.
![]() |
Fig. A.1. Fits of Galactic circular velocities using the Einasto model and the posterior distributions of fitting parameters without implementing adiabatic halo contraction using 12 baryonic models. |
![]() |
Fig. A.2. Fits of Galactic circular velocities using the Einasto model and the posterior distributions of fitting parameters implementing adiabatic halo contraction using 12 baryonic models. |
![]() |
Fig. A.3. Fits of Galactic circular velocities using the NFW model implementing adiabatic halo contraction using 12 baryonic models. Data points with errors are the rotation velocities from Jiao et al. (2023), while open triangles show the data from Eilers et al. (2019), which are not fitted. Blue, purple, green and black solid lines correspond to the contributions by the stellar disk, central bar, gas (and dust if any), and compressed dark matter halo, respectively. The total contributions are shown using red solid lines. Black dashed lines are the inferred primordial halos. |
All Tables
Inferred properties of the primordial halos parameterized with the Einasto profile (see Appendix A.3) before adiabatic contraction assuming a variety of baryonic models (see Appendix A.2).
All Figures
![]() |
Fig. 1. Example of a circular velocity fit using the McGaugh19 model for baryonic mass distributions. The purple, blue, and green lines represent the contributions of the bar, disk, and gas components, respectively. The solid and dashed black lines show the current and primordial dark matter halos, respectively. The solid red line indicates the total velocity profile. The black points show the latest Gaia measurements (Jiao et al. 2023), and the gray upward triangles and squares show the terminal velocities from (McClure-Griffiths & Dickey 2007, 2016), and Portail et al. (2017), respectively. The data marked with open symbols were not fit because they do not consider the systematic uncertainties. The fits with 12 baryonic models are shown in Fig. A.2. |
| In the text | |
![]() |
Fig. 2. Halo masses and concentrations of the primordial Galactic halos derived from the Gaia circular velocity fits using 12 baryonic models. The red and blue stars with errors represent the halos with and without adiabatic contraction, respectively. The predicted halo mass-concentration relation within 1σ from simulations (Dutton & Macciò 2014) is shown as the declining band. The vertical band shows the expected range of the MW halo mass according to the abundance-matching relation (Moster et al. 2013). The upper and lower limits are set by the highest stellar mass (model A&S) plus 1σ and the lowest stellar mass (model I) minus 1σ, respectively. |
| In the text | |
![]() |
Fig. 3. Structure of the inferred primordial and current Galactic halos, along with predictions for the cold and warm dark matter. The density profiles are scaled so that there is no need to assume or consider the masses or concentrations for these halos. The gray band indicates the range of the current halos derived from the Gaia velocity fits using the 12 baryonic models, and the red band shows their corresponding primordial halos within 1σ. The blue band presents the simulated halos with cold dark matter only (Dutton & Macciò 2014). The purple band shows the warm dark matter halos (normalized to match the primordial Galactic halo) with a core size spanning from 4.56 kpc (WDM5 in Macciò et al. 2012) to 7.0 kpc, corresponding to a particle mass of 0.05 keV and lower. |
| In the text | |
![]() |
Fig. A.1. Fits of Galactic circular velocities using the Einasto model and the posterior distributions of fitting parameters without implementing adiabatic halo contraction using 12 baryonic models. |
| In the text | |
![]() |
Fig. A.2. Fits of Galactic circular velocities using the Einasto model and the posterior distributions of fitting parameters implementing adiabatic halo contraction using 12 baryonic models. |
| In the text | |
![]() |
Fig. A.3. Fits of Galactic circular velocities using the NFW model implementing adiabatic halo contraction using 12 baryonic models. Data points with errors are the rotation velocities from Jiao et al. (2023), while open triangles show the data from Eilers et al. (2019), which are not fitted. Blue, purple, green and black solid lines correspond to the contributions by the stellar disk, central bar, gas (and dust if any), and compressed dark matter halo, respectively. The total contributions are shown using red solid lines. Black dashed lines are the inferred primordial halos. |
| In the text | |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.








![$$ \begin{aligned} \rho _{\rm bar} = \frac{\rho _0}{(1+r^{\prime }/r_0)^\alpha }\exp {\Big [-(r^{\prime }/r_{\rm cut})^2\Big ]},\ r^{\prime }=\sqrt{R^2+(z/q)^2}, \end{aligned} $$](/articles/aa/full_html/2025/11/aa57021-25/aa57021-25-eq6.gif)






![$$ \begin{aligned} \rho _{\rm EIN}(r) = \rho _s\exp \Big \{-\frac{2}{\alpha }\Big [\Big (\frac{r}{r_s}\Big )^\alpha -1\Big ]\Big \}, \end{aligned} $$](/articles/aa/full_html/2025/11/aa57021-25/aa57021-25-eq49.gif)






