Rotation curve decompositions with Gaussian Processes: taking into account data correlations leads to unbiased results

Correlations between velocity measurements in disk galaxy rotation curves are usually neglected when fitting dynamical models. Here I show how data correlations can be taken into account in rotation curve decompositions using Gaussian Processes. I find that marginalizing over correlation parameters proves critical to obtain unbiased estimates of the luminous and dark matter distributions in galaxies.


INTRODUCTION
Rotation curves provide robust constraints on the gravitating mass of disk galaxies, have been historically instrumental to characterize the "missing mass problem", and are useful today to link the luminous and dark masses in galaxies (Posti et al. 2019). The distribution of dark matter (DM) in galaxies is often constrained with rotation curve decompositions employing Bayesian analyses (Posti et al. 2019;Li et al. 2020). Usually such models assume that each point in the curve is independent, since measurements weigh equally in standard χ 2 fits that do not include correlation terms. When extracting rotation curves from observations, correlations between velocity measurements are typically neglected and each annulus is treated independently. This, at best, is a first-order approximation, since projection effects, resolution effects, and physical phenomena e.g. non-circular motions, can correlate velocities in adjacent annuli.
Here, I demonstrate how to take correlations between rotation curves' datapoints into account using Gaussian Processes (GPs) and I show that doing so is critical to get unbiased estimates of DM halo parameters. A python code is made publicly available (Posti 2022).

ROTATION CURVE DECOMPOSITIONS
Data I consider the composite Hα-HI rotation curve of the galaxy NGC2403 obtained with Fabry-Perot and radio interferometry, respectively by Daigle et al. (2006) and Fraternali et al. (2002). I selected this curve from the SPARC database (Lelli et al. 2016) because of its fine sampling and extension, so that my experiment is safe from artifacts due poor data quality.
The rotation curve of a galaxy is where V rot (R) is the measured circular velocity, while V gas (R), V (R), and V DM (R) are the contributions from gas, stars, and DM. V gas and V are determined directly from observations of the gas and stellar surface brightness, with the unknown scaling of the mass-to-light ratio M/L. V DM is the objective of the fit and here is parametrized assuming Navarro et al. (1996, hereafter NFW) profiles for the DM, i.e. with mass M h and concentration c as free parameters. I use Eq.
(1) to fit for the stellar mass-to-light ratio M/L, halo mass M h and concentration c, and I take at face value the data from the SPARC catalog: V rot,i and uncertainties σ Vrot,i , V gas,i , and V ,i for each radius R i .

Model without GPs
Priors are standard in this context (e.g. Li et al. 2020), i.e. uninformative for log M h , Gaussian for log c, with mean following the c − M h relation from cosmological simulations, and Gaussian for log (M/L), centered on 1 M/L = 0.5. This likelihood implicitly assumes that the datapoints are independent, because each term in the sum depends only measurements at a specific radius R i . The results are presented in Figure 1, where the left-hand panel shows the corner plot of the marginalized posterior in blue. The posterior is well sampled and narrow, implying a tiny uncertainty i.e. log M h /M = 11.33 ± 0.02. The top-right panel of Fig. 1 shows the decomposed rotation curve where the blue curve is the prediction of the median model. Note that I also plot a light blue band representing the 2 − σ uncertainty on the predicted V rot (R); however, the band is small compared to width of the blue curve and is barely visible at R 20 kpc.
These results are compatible with those of previous works (Posti et al. 2019;Li et al. 2020).

Model with GPs: allowing for correlated data
To relax the assumption of independent V rot,i datapoints I introduce the covariance matrix K, whose values K ij are given by the kernel function k as where δ ij is the Kronecker delta and the parameters A k and s k are a characteristic amplitude and scale. k is maximal, implying strong correlation, for points at small distances, while it smoothly approaches zero, implying no correlation, for points where |R i − R j | s k . The functional form of k(R i , R j ) is arbitrarily chosen as it is common in GP studies (e.g. Aigrain & Foreman-Mackey 2022). Fortunately, in my testing, varying this function does not make an impactful difference.
The new log-likelihood is where V res is the vector with components V res,i = V rot,i − V model (R i |θ V ), while K −1 and |K| are the inverse and determinant of K. With this notation the problem now becomes a GP regression, for which I use the library tinygp . I assume uninformative priors on the two additional free parameters θ k = (log A k , log s k ) and I employ the same setup of the previous run. The resulting the marginalized posterior is shown in the left-hand panel of Fig. 1 as the orange model. The blue and the orange models have an incompatible mean at 2 − σ, which suggests that the results of the less flexible model (blue) may be biased.
The bottom-right panel of Fig. 1 shows the curve decomposition of the median orange model. When including GPs, the model predicts larger uncertainties in V rot (light orange band), providing a better fit to observations, and a more massive stellar disk, which is dominant over DM within R 4 kpc, compared to the model without GPs. The posterior also becomes wider when including GPs, implying larger uncertainties (e.g. log M h = 11.50±0.09), because accounting for data correlations introduces additional degrees of freedom, which increase the model's capacity and reduce biases. Such θ k parameters are actually nuisance parameters that when marginalized over yield more realistic and less biased predictions on θ V .
The inset in the bottom-right panel of Fig. 1 shows the resulting kernel k(R i , R j ). The parameters θ k are well constrained: log A k = 1.10 ± 0.26 and log s k /kpc = 0 ± 0.1, i.e. points closer than ∼ 1 kpc have covariance of ∼ 10 km 2 /s 2 (to be compared with the typical σ 2 Vrot 7 km 2 /s 2 ). The two-block shape is due to the composite nature of the curve: optical in the center with 0.1 kpc spacing, radio in the outskirts with 0.5 kpc spacing.

CONCLUSIONS
In this paper I used state-of-the-art observations to show that assuming independent velocity datapoints in rotation curve fits leads to biased parameter constraints and unrealistically small uncertainties on the predicted velocity profile. GPs can be used to account for unknown data correlations with minimal modifications to existing pipelines.
Deriving rotation curves from optical or radio observations is a complex task, unlikely to result in a series of uncorrelated velocity measurements. GPs may provide a simple and computationally cheap solution to use the wealth of literature rotation curve data while accounting for possible underlying correlations.