An overlooked source of uncertainty in the mass of the Milky Way

In the conventional approach to decomposing a rotation curve into a set of contributions from mass model components, the measurements of the rotation curve at different radii are taken to be independent. It is clear, however, that radial correlations are present in such data, for instance (but not only) because the orbital speed depends on the mass distribution at all (or, minimally, inner) radii. We adopt a very simple parametric form for a covariance matrix and constrain its parameters using Gaussian process regression. Applied to the rotation curve of the Milky Way, this suggests the presence of correlations between neighbouring rotation curve points with amplitudes $<10\,\mathrm{km}\,\mathrm{s}^{-1}$ over length scales of $1.5$-$2.5\,\mathrm{kpc}$ regardless of the assumed dark halo component. We show that accounting for such covariance can result in a $\sim 50$ per cent lower total mass estimate for the Milky Way than when it is neglected, and that the statistical uncertainty associated with the covariance is comparable to or exceeds the total systematic uncertainty budget. Our findings motivate including more detailed treatment of rotation curve covariance in future analyses.

Rotation curves are a widely-used dynamical tracer of the mass content of and distribution within galaxies including the Milky Way (e.g.Carignan et al. 2006;de Blok et al. 2008;Kuzio de Naray et al. 2008;Swaters et al. 2009;McMillan 2011;Bovy et al. 2012;Sofue 2012;Adams et al. 2014;Pato et al. 2015;Lelli et al. 2016;Eilers et al. 2019;Cautun et al. 2020;Ou et al. 2024, amongst many others).The atomic measurements comprising a rotation curve -the orbital speed at fixed radiusare usually taken to be independent from and uncorrelated with their counterparts at other radii (including in all the references above; see also Põder et al. 2023, for some discussion of concerns around such correlations).However, it is easy to see that correlations across radii must exist.One unambiguous source is the connection between gravitationally-driven kinematics and integrals of the mass distribution -for example, the mass within some central aperture contributes to the kinematics at all larger radii, introducing a correlation.Correlations may also arise due to instrumental effects such as beam smearing (Swaters et al. ⋆ kyle.a.oman@durham.ac.uk 2009, and references therein), modelling effects such as a geometrically thick disc being imperfectly separated into rings, or physical effects such as a spiral arm coherently perturbing the kinematics over a range in radius.The presence of correlations can further be inferred by noticing that the scatter implied by the statistical uncertainties on rotation curve measurements often exceeds the point-to-point scatter measured (but misestimates of the uncertainties could also contribute, including incorrectly assuming that uncertainties are Gaussian distributed).
The heterogenous origins of radial correlations in rotation curve measurements makes them challenging to model explicitly.Posti (2022) proposed the pragmatic approach of assuming a parametric form for the covariance matrix: This describes a correlation with amplitude 1 a k between points of a given radial separation Ri − Rj that decays exponentially with a length scale s k ; the uncertainties on the individual measurements are σi.There is freedom in the choice of 'kernel func- Posti (2022).
tion' (the first term in Eq. ( 1)) but Posti (2022) reports that reasonable variations in the choice do not lead to large differences in results, which is enough for our illustrative purposes in this work.Posti (2022) argued that marginalizing over the possibility of such correlations across radii, even in such a simplistic manner, leads to more realistic and less biased confidence intervals on model parameters of interest, such as those describing the dark halo component of a galaxy.Analysis of the rotation curve of the Milky Way is distinct from other galaxies in two important ways.First, the systematic uncertainties affecting the measurement are distinct from those relevant in other galaxies due to our unique vantage point, especially for recent measurements incorporating high-precision proper motion measurements of stars.Second, we have more constraints on the visible matter content and distribution of our Galaxy than for external galaxies (e.g.Binney & Tremaine 2008, sec. 2.7, McMillan 2017;Sofue 2020, and references therein), making mass models of the Milky Way more tightly constrained.These considerations motivate us to apply the approach of Posti (2022), who illustrated it using measurements of NGC 2403, to the Milky Way in order to assess whether accounting for correlations in the measurements make up a significant portion of the uncertainty budget in mass models of the Milky Way.We also explore whether failing to account for such correlations is likely to lead to biases in rotation curve-based measurements of the mass of the Milky Way.

MASS MODEL COMPONENTS AND FITTING METHODOLOGY
We use the rotation curve reported by Ou et al. (2024, see their table 1) as a representative example of recent measurements (see also Wang et al. 2023;Zhou et al. 2023;Jiao et al. 2023) incorporating data from the Gaia mission (Gaia Collaboration et al. 2023).We also adopt the structural parameters for the baryonic components of the mass model of Ou et al. (2024, see their table 2).The circular velocity curve specified by our model given its parameters is: We have grouped together the contribution of the stellar bulge and disc (both held fixed in model optimisation) in vstars, and likewise the contributions of H I gas, H2 gas, warm dust and cold dust in vgas (also all held fixed).Whereas Ou et al. (2024) derived the circular velocity of each component from the enclosed mass at each radius, for disc-like components we instead use the expression in terms of the gradient of the potential Φ in the disc midplane: A derivation of the potential of a thick exponential disc can be found in Binney & Tremaine (2008, sec. 2.6.1c);we evaluate the relevant integrals numerically and use a cubic spline approximation to measure its radial derivative.For sphericallysymmetric components we use the usual expression vcirc = GM enclosed /r.We use two models for v dark matter : the Navarro et al. (1996, NFW) profile and pseudo-isothermal (Gunn & Gott 1972) profile.Both can be expressed as: where 1 3 and ∆crit = 200.The mass enclosed within a sphere within which the average density is ∆crit times the critical density for closure, M200c, is a free parameter.The models differ in the definition of the second free parameter -the 'concentration paramter' c.For the NFW profile we adopt the usual definition cNFW = R200c/Rs, where Rs is the 'scale radius' in the density profile ρ(R) ∝ (R(1 + R/Rs)) −2 .For the pseudo-isothermal profile we define cPI = R200c/Rc where Rc is the 'core radius' in the density profile ρ(R) ∝ (1 + (R/Rc) 2 ) −1 .Finally, the function fc is defined for the two models respectively as 2 : and We stress that the choice of these two models for v dark matter is not motivated by their suitability to describe the Milky Way rotation curve data.Instead, they are chosen for their similarity in terms of simplicity of mathematical form (e.g. two free parameters each, analytic expressions for all needed quantities) but stark dissimilarity in structure (central ρ ∝ r −1 and outer ρ ∝ r −3 density profile for the NFW model vs. central ρ ∝ r 0 and outer ρ ∝ r −2 density profile for the pseudoisothermal model).We will show below that some outcomes of allowing for correlations in the rotation curve data are common to both models, suggesting that the lessons learned are quite general.
We optimize the two free parameters of the model following exactly the same methodology as Posti (2022) -in fact we use their software implementation, with some modification to include the stellar bulge component and the pseudo-isothermal model in addition to the NFW model, in both the case assuming independent measurements of the rotation curve at each radius and that where correlations across radii are modelled with Gaussian process (GP) regression.

IMPACT OF CORRELATION MODELING ON INFERRED MILKY WAY MASS PROFILES
We show the mass models resulting from our modelling for the two halo models in the cases without (upper panels) and with (lower panels) the GP regression model for radial correlations in the Milky Way rotation curve in Fig. 1.Corresponding best-fitting parameter values and uncertainties are reported in Table 1, and one-and two-dimensional marginalised posterior probability distributions are shown in the Appendix.In the cases without GP regression the confidence intervals on the models' total circular velocity curves are unrealistically small, and the models' inability to adequately capture the data is reflected by a reduced chi-squared (χ 2 ν ) of about 10 (NFW halo)  2024).In all cases the stellar disc (dashed yellow) and bulge (dotted yellow; combined stellar components shown with solid yellow), and 'gas' (green; including H I gas, H 2 gas, warm dust and cold dust) are kept fixed to the model proposed by Ou et al. (2024); see Sec. 2 for details.In the left panels the halo component (solid purple) is an NFW model, while in the right panels it is a pseudo-isothermal sphere.In the upper panels no correlation ('no GP') between the observed rotation speeds is assumed, while in the lower panels the covariance matrix for the observations is estimated using Gaussian process ('GP') regression as described by Posti (2022).Shaded bands mark the regions enclosing 95 per cent of model curves at each radius for the halo component and model total (in the upper panels these bands are very thin).
or 5 (pseudo-isothermal halo).Qualitatively, these fits are being driven by the measurements near 10 < R/kpc < 15 that have very small uncertainties and therefore miss the measurements near 20 < R/kpc < 25 where the uncertainties are larger.
In the cases with GP regression the models achieve a more balanced 'compromise' fit because a moderate deviation from the points with small uncertainties (10 < R/kpc < 15) is mitigated by the assumption that these measurements are correlated.The correlation amplitudes and scales preferred by the models seem intuitively plausible for a galaxy like the Milky Way: (a k , s k ) ∼ (5 km s −1 , 2.4 kpc) for an NFW halo, or (3 km s −1 , 2.0 kpc) for a pseudo-isothermal halo.The confidence intervals on the total circular velocity curves and dark halo components are larger and, we feel, more representative of the statistical uncertainty in the measurements.The fits are for-Table 1.The first four columns show the best-fitting and marginalized 16 th -84 th percentile confidence intervals for the free parameters of our mass models for the two halo models -NFW and pseudo-isothermal (P-I) -and the case where correlations between rotation curve points are ignored (no GP) or accounted for through Gaussian process regression (GP).The last column shows the inference on the mass of the dark matter halo component of the model within a 20 kpc spherical aperture.Values on a linear scale are given here for ease of reference, but we note that the probability distributions are either close to log-normal (M 200c , c, M DM (r < 20 kpc)) or asymmetric (a k , s k ) -see figures in the Appendix.mally much better, with χ 2 ν ∼ 1 for both halo models, justifying the addition of the two additional free parameters.
In addition to wider confidence intervals, the models with GP regression are also biased towards having somewhat more centrally-concentrated dark halo components, visible in Fig. 1 as an elevated circular velocity for the halo component near the centre (R < ∼ 10 kpc).This is compensated by lower halo masses such that the total dark matter mass within 16 and 25 kpc is the same as in the models without GP regression for the NFW and pseudo-isothermal halo cases, respectively (the dark matter mass within 20 kpc for each model is tabulated in Table 1).The marginalised posterior probability distributions for the halo mass for the 4 models plotted in Fig. 1 are shown in Fig. 2 (see also the Appendix).The variants with GP regression have sys-tematically lower M200c, by ∼ 0.2 dex (58 per cent) for the NFW halo or ∼ 0.1 dex (26 per cent) for the pseudo-isothermal halo.

Appropriateness of the halo models
We briefly discuss whether the models shown in Fig. 1 are realistic.The NFW halo model, in particular, comes with a strong implied prior on the concentration cNFW given the mass M200c (e.g.Ludlow et al. 2014).We have chosen not to impose this prior in our modelling, primarily to enable a fair comparison with the pseudo-isothermal halo model where no similarly strong prior exists (but see Kormendy & Freeman 2004;de Blok et al. 2008).Given the existence of a mass-concentration relation for the NFW model, it makes sense to ask whether our models are consistent with it.We have deliberately not shown the relation in the c versus M200c panel of Fig. A1 -in fact it lies largely outside of the axes.Both of our models with an NFW halo (with and without GP regression) prefer a region of the parameter space many standard deviations above the locus of the mass-concentration relation.This can reasonably be attributed to neglecting any response of the halo to the assembly of the Galaxy (see e.g.Cautun et al. 2020, and references therein), possible mismodelling of its baryonic components, or, likely, a combination of these.
We stress, however, that our objective in this work is not to create a realistic mass model for the Milky Way, but to highlight the kinds of systematic biases that can arise when correlations between points in the rotation curve measured at different radii are ignored.With this end in mind, the choice of halo models and how realistic they are (within reason) is irrelevant: our modelling clearly shows that biases of the same sign and similar amplitude arise in both models that we explore, despite their dissimilarity.We therefore expect that more realistic models, which likely have a central density profile somewhere between a steep NFW cusp and the flat core of the pseudo-isothermal model, suffer from the same sorts of systematic biases (e.g. in halo mass; Fig. 2) as our models.

Importance of statistical uncertainty in the uncertainty budget
Another useful question to consider is whether the statistical uncertainty associated with correlations between circular velocity curve measurements at different radii is comparable to other leading sources of systematic uncertainty in the determination of the circular velocity curve.Ou et al. (2024) provide estimates for the systematic uncertainties associated with effects such as varying the assumed density profile of the kinematic tracer population, the choice to neglect a certain higher-order term in the Jeans' equations, the uncertainty in the galactocentric Solar radius, and others in their fig.6.We reproduce the curve showing the sum in quadrature of all of the systematic uncertainties considered in that work as a function of radius with the black line in Fig. 3.The systematic uncertainty budget is about 2-3 per cent, with a gradual increase between about 5 and 20 kpc, and a sharp increase to > 10 per cent at larger radii.
On the same figure we show the statistical uncertainty in the total model circular velocity curve as a function of radius for our models with GP regression.Comparing to the models without GP regression (e.g.shaded bands in Fig. 1), it is clear that the statistical uncertainty is severely underestimated if correlations between rotation curve measurements at different radii are neglected.When correlations are accounted for, the statistical uncertainty is about 2 to 5 per cent (purple solid and dashed lines for the NFW and pseudo-isothermal halo models, respectively).Whereas in the case without GP regression we would conclude that the total uncertainty budget for the model circular velocity curve is strongly systematics-dominated, in the case with GP regression it is clear that the statistical uncertainty cannot be neglected.

Possible pitfalls in interpretation
The approach of modelling correlations in rotation curve measurements across radii by adding degrees of freedom (a k and s k ) to the model and marginalizing over them requires care in interpretation.We have experimented with applying the same methodology to external galaxies and have encountered cases where the proposed mass model is unable to describe the data (an NFW dark halo being fit to a linearly rising rotation curve, for example).In such cases the parameter search responds by moving to very large values of a k (e.g.> 1000 km s −1 ) and s k (e.g.> 100 kpc), bounded above only by the boundary imposed on the (flat) prior.Such cases are clear model failures, but serve to illustrate that the parameters a k and s k will respond to essentially any discrepancy between the model and data by proposing a stronger correlation.In our exploration of the Milky Way, it is clear that something is missing from the models because χ 2 ν ≫ 1 (see Figs. A1 & A2).Allowing for correlations in the measurements provides a plausible extension to the models: the goodness of fit improves and the additional parameters converge to plausible values corresponding to a fraction of the rotation speed and linear size of the disc.The question of whether this provides a more compelling explanation than other ways of accounting for the discrepancy between the data and models when correlations are neglected remains.
It is very challenging to account for all possible sources of correlation.Attempting to write down the full covariance matrix for a rotation curve measurement is currently infeasibleindeed this provides the motivation for the approach of Posti (2022).However, taking some burden of capturing correlations off of the model by describing them in the data and its associated uncertainties would undoubtedly help to mitigate the possible pitfalls described above.

SUMMARY
We have shown that accounting for the possibility that rotation curve measurements of the Milky Way at different radii are statistically correlated has significant implications for inference of the Milky Way's total mass and the structure of its dark halo.In particular: • accounting for such correlations can lead to a difference in the inferred mass of the Milky Way (lower by about 50 per cent); • it also results in larger uncertainty estimates for mass model parameters that we feel are more representative of the constraining power of the data; • the uncertainty budget in mass models of the Milky Way is likely not dominated by systematic uncertainties once the statistical uncertainty associated to correlations in the measurements is accounted for.
The above conclusions hold whether we assume an NFW or a pseudo-isothermal form for the dark halo component in our mass models, suggesting that they are probably generic for any broadly plausible choice of dark halo model.This strongly motivates including an allowance for correlations in the rotation curve measurements in future efforts to decompose the rotation of the Milky Way into components, but we caution that our assumed form for the covariance matrix (Eq. 1) is likely too simple to fully capture the correlations likely to be present in the measurements.

Figure 1 .
Figure1.Illustrative mass models for the Milky Way, with and without modeling radial correlations.The measurements (points with 1σ error bars) are as presented inOu et al. (2024).In all cases the stellar disc (dashed yellow) and bulge (dotted yellow; combined stellar components shown with solid yellow), and 'gas' (green; including H I gas, H 2 gas, warm dust and cold dust) are kept fixed to the model proposed byOu et al. (2024); see Sec. 2 for details.In the left panels the halo component (solid purple) is an NFW model, while in the right panels it is a pseudo-isothermal sphere.In the upper panels no correlation ('no GP') between the observed rotation speeds is assumed, while in the lower panels the covariance matrix for the observations is estimated using Gaussian process ('GP') regression as described byPosti (2022).Shaded bands mark the regions enclosing 95 per cent of model curves at each radius for the halo component and model total (in the upper panels these bands are very thin).

Figure 2 .
Figure 2. Marginalized posterior probability distributions for the model halo mass.The four cases shown are for an NFW halo (filled histograms) and a pseudo-isothermal halo (open histograms) for the nocorrelations ('no GP'; grey) and covariance estimated with Gaussian process regression ('GP', purple) models.For both halo models, including an estimate of the coviariance results in a significant bias to lower halo mass; the posterior probability distribution also becomes wider.One-and two-dimensional marginalised posterior probability distributions for all model parameters ('corner plots') can be found in Figs.A1 & A2.

Figure 3 .
Figure 3.Comparison of the statistical uncertainty as a function of radius for our mass models including Gaussian process regression estimates of the covariance matrix.Measurements in the case with the NFW halo model (solid purple) and pseudo-isothermal halo model (dashed purple) are shown.The curves represent the width of the interval enclosing 68 per cent of model curves at each radius, normalised by the model circular velocity at that radius.These confidence intervals are comparable in width to the estimate of the total (1σ) systematic uncertainty estimate of Ou et al. (2024, black, reproduced from their fig.6).

Figure A1 .
Figure A1.One-and two-dimensional marginalised posterior probability distributions for the parameters of our mass models with an NFW halo model in the 'no GP' (grey) and 'GP' (purple) cases.Model parameters in the 'no GP' case are the halo mass M 200c and concentration c NFW ≡ R 200c /Rs.The GP case supplements these with a correlation amplitude a k and length scale s k .The upper right panel shows the distribution of χ 2 ν values for samples in the Markov chains.