# A Critical Assessment of Solutions to the Galaxy Diversity Problem

Aidan Zentner  Department of Physics, Princeton University, Princeton, NJ
08544, USA<EMAIL_ADDRESS>Siddharth Dandavate  Department of Physics,
Princeton University, Princeton, NJ 08544, USA  Oren Slone Department of
Physics, Princeton University, Princeton, NJ 08544, USA  Center for Cosmology
and Particle Physics, Department of Physics,
New York University, New York, NY 10003, USA  Mariangela Lisanti Department of
Physics, Princeton University, Princeton, NJ 08544, USA  Center for
Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY
10010, USA

###### Abstract

Galactic rotation curves exhibit a diverse range of inner slopes.
Observational data indicates that explaining this diversity may require a
mechanism that correlates a galaxy’s surface brightness with the central-most
region of its dark matter halo. In this work, we compare several concrete
models that capture the relevant physics required to explain the galaxy
diversity problem. We focus specifically on a Self-Interacting Dark Matter
(SIDM) model with an isothermal core and two Cold Dark Matter (CDM) models
with/without baryonic feedback. In contrast to the CDM case, the SIDM model
can lead to the formation of an isothermal core in the halo, and is also
mostly insensitive to baryonic feedback processes, which act on longer time-
scales. Using rotation curves from 90 galaxies in the _Spitzer_ Photometry &
Accurate Rotation Curves (SPARC) catalog, we perform a comprehensive model
comparison that addresses issues of statistical methodology from prior works.
The best-fit halo models that we recover are consistent with standard CDM
concentration-mass and abundance matching relations. We find that both the
SIDM and feedback-affected CDM models are better than a CDM model with no
feedback in explaining the rotation curves of low and high surface brightness
galaxies in the sample. However, when compared to each other, there is no
strong statistical preference for either the SIDM or the feedback-affected CDM
halo model as the source of galaxy diversity in the SPARC catalog.

##  1 Introduction

Galactic rotation curves have historically served as a powerful tool for
understanding the distribution of dark and visible matter in the Universe.
With recent catalogs of rotation curves that span a wide range of galaxy
masses, we now have the unique opportunity to test the consistency of dark
matter (DM) models across a broad population of galaxies. One important
observation that must be accounted for by any DM model is the diversity of
galactic rotation curves and their apparent correlation with the baryonic
properties of a galaxy (Flores & Primack, 1994; de Blok et al., 2001; Kuzio de
Naray et al., 2010; Oh et al., 2011; Oman et al., 2015). In this work, we use
the _Spitzer_ Photometry & Accurate Rotation Curves (SPARC) catalog (Lelli et
al., 2016a) to critically assess the preference for models of collisionless
Cold DM (CDM) with and without feedback, and for self-interacting DM (SIDM).

The diversity problem is the simple observation that the rotation curves of
spiral galaxies exhibit a diverse range of inner slopes. Specifically, high
surface brightness galaxies typically have more steeply rising rotation curves
than do their low surface brightness counterparts, even if both have similar
halo masses. This behavior cannot easily be explained if halo formation is
self-similar, for example if the halo is modeled using the standard Navarro-
Frenk-White (NFW) profile (Navarro et al., 1997). Within the CDM framework,
two galaxies with similar halo masses will have comparable total stellar mass
and concentration. If no feedback mechanism is invoked, this sets the
properties of the DM profile, typically modeled as an NFW profile, such as its
scale radius and overall normalization. The inner slope of the halo density is
insensitive to baryonic feedback and is fixed at a constant value. In this
scenario, therefore, the halo model has little flexibility to account for
correlations with the surface brightness of the galaxy (Oman et al., 2015;
Katz et al., 2016; Li et al., 2019, 2020).

One way to address galaxy diversity within the CDM framework is to use
baryonic feedback to inject energy into the halo and redistribute DM to form
larger cores near the halo’s center (Navarro et al., 1996; Gnedin & Zhao,
2002; Read & Gilmore, 2005; Governato et al., 2010; Pontzen & Governato, 2012;
Di Cintio et al., 2014; Madau et al., 2014; Di Cintio et al., 2014; Chan et
al., 2015; Tollet et al., 2016; Lazar et al., 2020; Macciò et al., 2020). For
example, repetitive star bursts have been demonstrated to effectively produce
cores in the brightest dwarf galaxies (Pontzen & Governato, 2012). For less
massive dwarfs, such feedback is not strong enough to efficiently remove a
galaxy’s inner-most mass. For larger systems, the galaxies are so massive that
feedback is insufficiently powerful to create coring and profiles similar to
NFW are preserved. Baryonic feedback processes therefore link the properties
of a galaxy’s DM density profile to its baryonic content, as needed to explain
the galaxy diversity, and provide improved fits to rotation curve data (Katz
et al., 2016; Allaert et al., 2017; Li et al., 2019, 2020; Santos-Santos et
al., 2020).

An alternate approach to address the galaxy diversity problem is to change the
properties of the DM model itself. A well-known example is that of SIDM
(Spergel & Steinhardt, 2000). If DM can self interact, heat transfer in the
central regions of the halo is efficient, creating an isothermal region within
the density profile (Vogelsberger et al., 2012; Zavala et al., 2013; Rocha et
al., 2013b; Peter et al., 2013). This central, isothermal profile is
exponentially dependent on the baryonic potential of the galaxy. While a DM-
dominated system will form a core, a baryon-dominated system will retain a
more cuspy center, similar to that of an NFW profile (Kaplinghat et al.,
2014a). Thus, in the SIDM framework, the size of the core naturally correlates
with the baryonic surface brightness of the galaxy, allowing for a mechanism
that can account for galaxy diversity (Kamada et al., 2017; Creasey et al.,
2017; Ren et al., 2019; Kaplinghat et al., 2020).

When the time-scale for forming the isothermal core is shorter than the time-
scale for baryonic feedback, SIDM halos are insensitive to the details of
feedback processes (Ren et al., 2019). In this case, an SIDM halo thermalizes
fast enough that the resulting density distribution is only sensitive to the
present-day baryonic distribution, and not the formation history of the
galaxy. These results have been confirmed with hydrodynamic simulations,
which—for the self-interaction cross sections studied—have found that SIDM
density profiles are more robust to stellar feedback than their CDM
counterparts (Vogelsberger et al., 2014; Fry et al., 2015; Robles et al.,
2017; Robertson et al., 2018; Sameie et al., 2021; Vargya et al., 2021; Burger
et al., 2022).

This paper presents a systematic comparison of three models with equal numbers
of degrees-of-freedom using SPARC galaxies. These models attempt to capture
the important physics of: 1) the isothermal core in SIDM, 2) baryonic feedback
in CDM, and 3) CDM with no baryonic feedback. Our analysis highlights aspects
of the SPARC dataset that can potentially bias the model comparison and also
addresses several issues with the statistical methodology of previous SPARC
rotation curve studies (Katz et al., 2016; Li et al., 2019; Kaplinghat et al.,
2020; Li et al., 2020). We find that the CDM model with no baryonic feedback
is strongly disfavored for low surface brightness galaxies, although less-so
for high surface brightness galaxies. Additionally, we find no strong
statistical preference for the SIDM model or the CDM model with baryonic
feedback, for either the low or high surface brightness systems.

The paper is organized as follows. In Sec. 2, we discuss the SPARC dataset,
emphasizing several aspects of the publically available data that can
systematically affect the results of this study. Sec. 3 describes the data
selection, theory modeling, and likelihood analysis. Sec. 4 presents the
rotation curve fits to a few example galaxies and discusses over-fitting
issues to the data. Sec. 5 presents the statistical model comparison for the
scenarios considered in this work. We then conclude in Sec. 6. App. A provides
rotation curve fits and corner plots for example galaxies spotlighted in the
text. App. B focuses on several anomalous galaxies in the SPARC dataset.

Rotation curve fits and corner plots for all galaxies used in this study are
publicly available on github.111https://github.com/siddanda/DarkModels

##  2 The SPARC Dataset

The SPARC database includes 3.6 $\mu$m surface photometry measurements and
H1/H$\alpha$ rotation curves for 175 late-type galaxies (Lelli et al., 2016a).
This comprehensive sample includes galaxies with stellar masses $3\times
10^{7}\lesssim M_{\star}\lesssim 3\times 10^{11}~{}M_{\odot}$ and a wide range
of luminosities, morphologies, and surface brightnesses. While it is
representative of nearby disk galaxies, it does not form a statistically
complete sample (Katz et al., 2016; Lelli et al., 2016a).

The SPARC catalog is publicly available and provides the observed rotation
curve measurements, $V_{i,{\rm obs}}$, and their errors, $\delta V_{i,{\rm
obs}}$, for each galaxy, where the subscript $i$ corresponds to the radius at
which the measurement is reported. The baryonic mass models are given as the
circular velocities for the disk ($V_{i,{\rm disk}}$), bulge ($V_{i,{\rm
bul}}$), and gas ($V_{i,{\rm gas}}$) components. These velocities were found
by solving the Poisson equation for the observed surface brightness at 3.6
$\mu$m for disk and bulge contributions and from the HI surface brightness for
the gas contribution, assuming mass-to-light ratios of unity for all,
$\Upsilon_{\star,\rm{disk}}=\Upsilon_{\star,\rm{bul}}=1\,M_{\odot}/L_{\odot}$
(Lelli et al., 2016a).

If a model predicts a galaxy’s rotation curve contribution from DM at a
certain radius, $V_{\rm DM}$, these can be used to infer the galaxy’s
predicted total rotation curve through,

$\displaystyle V_{\rm tot}^{2}=$ $\displaystyle V_{\rm DM}\,|V_{\rm
DM}|+V_{\rm gas}\,|V_{\rm gas}|+\Upsilon_{\star,\rm{disk}}V_{\rm
disk}\,|V_{\rm disk}|$ (1) $\displaystyle\hskip
7.22743pt+\,\Upsilon_{\star,\rm{bul}}V_{\rm bul}\,|V_{\rm bul}|\,,$

where the $i$ subscripts have been suppressed for notational
simplicity.222Absolute values are needed to account for negative velocities
that are provided in the SPARC catalog—see Lelli et al. (2016a). The model can
then be compared to data to quantify its goodness of fit, as was done in Katz
et al. (2016); Kamada et al. (2017); Li et al. (2019); Ren et al. (2019);
Kaplinghat et al. (2020); Li et al. (2020).

There are three important aspects of the SPARC data that should be carefully
considered when performing such statistical model comparisons. First,
uncertainties on the baryonic contribution to rotation curves are not included
in the SPARC catalog. Technically, such uncertainties should be directly
propagated to the predicted $V_{\rm tot}$ through $V_{\rm gas}$, $V_{\rm
disk}$ and $V_{\rm bul}$ in Eq. (1). In line with other studies performed to
date, we are unable to properly account for these uncertainties given the
information that is publicly available in the SPARC catalog. Moreover, in the
case of the SIDM and feedback-affected CDM model considered here, the
predicted value of $V_{\rm DM}$ itself depends on the galaxy’s baryonic
distribution and the uncertainties in baryonic velocities should in principle
be propagated into the DM model. Unfortunately, it is not possible to estimate
how this affects the conclusions of this study without knowing the magnitude
of the baryonic uncertainties.

Second, some models may over-fit the data. Many SPARC galaxies have data
points whose spread around the best-fit model predictions is smaller than
expected from the uncertainties on those data points. This could be indicative
of correlations between data points that are not accounted for in the model,
or of over-estimation of errors for those data points. To address this issue,
potentially problematic galaxies will be flagged when presenting results. The
selection tests for over-fit galaxies are discussed in detail in Sec. 4.

Lastly, due to the possibility of unknown systematic biases in the selection
of SPARC galaxies, one might be concerned that the best-fit parameters of a
Bayesian analysis will be strongly prior-dependent. For this reason, the
benchmark analysis of this study is designed with reasonably loose priors.
Sec. 5 includes a summary of how variations on these priors affect the primary
conclusions.

##  3 Methodology

From the SPARC catalog, 98 galaxies with at least 10 data points are selected
that have a well-defined velocity for the flat part of the rotation curve
(Lelli et al., 2016b), $V_{\text{f}}$, with inclination $i\geq 30^{\circ}$, as
well as those suitable for dynamical studies with quality flag $Q\neq 3$. For
each galaxy, our analysis uses the observed circular velocity and associated
uncertainty as a function of galactocentric radius, as well as the derived
contributions to the total velocity from the gas, galactic disk, and galactic
bulge, if present.

The following subsections describe how we model mass distributions for CDM
(Sec. 3.1) and SIDM (Sec. 3.2) halos and how the predicted velocity rotation
curves are then fit to data (Sec. 3.3). Throughout, we work with virial
quantities and use the Planck 2015 cosmology with $H_{0}=67.8$ km s-1 Mpc-1,
$\Omega_{m}=0.308$, and $\Omega_{\Lambda}=0.692$ (Planck Collaboration, 2016).
The virial mass, $M_{\rm vir}$, is defined as the enclosed halo mass at the
virial radius, $r_{\rm vir}$, where the average enclosed density is a factor
$\Delta_{c}$ greater than the critical density. Following Bryan & Norman
(1998) and assuming the Planck cosmology, $\Delta_{c}\approx 102$.

###  3.1  Mass Models for Cold Dark Matter

Two spherically-symmetric halo models for CDM are considered, one with and one
without baryonic feedback processes accounted for. Both CDM halos are modeled
using the $(\alpha,\beta,\gamma)$ density profile,

$\rho_{\scriptscriptstyle{\rm
CDM}}(r)=\frac{\rho_{s}}{\left(r/r_{s}\right)^{\gamma}\left[1+\left(r/r_{s}\right)^{\alpha}\right]^{(\beta-\gamma)/\alpha}}\,,$
(2)

where $\rho_{s}$ is the scale density and $r_{s}$ is the scale radius
(Hernquist, 1990). The logarithmic slope of the density profile goes as
$d\log\rho/d\log r=-2$ at a galactocentric distance of

$r_{\scriptscriptstyle{-2}}=\left(\frac{2-\gamma}{\beta-2}\right)^{1/\alpha}r_{s}\,.$
(3)

The enclosed mass, $M(r)$, takes an analytic form for this profile and the
prediction for the DM contribution to the circular velocity at any radius is

$V_{\rm DM}^{2}(r)=\frac{GM(r)}{r}\,,$ (4)

where $G$ is Newton’s gravitational constant. Once the $\alpha,\beta,\gamma$
exponents are specified, there are two free parameters remaining. Without loss
of generality, they are chosen to be the concentration, $c_{\rm vir}=r_{\rm
vir}/r_{\scriptscriptstyle{-2}}$, and the virial velocity,

$V_{\rm vir}=r_{s}H_{0}c_{\rm vir}\sqrt{\frac{\Delta_{c}}{2}}\,.$ (5)

To predict the total rotational velocity $V_{\rm tot}$ and compare to data,
the mass-to-light ratio parameters, $\Upsilon_{\star,\rm{disk}}$ and
$\Upsilon_{\star,\rm{bul}}$, must also be specified. Thus, the CDM models
considered here have a total of four free parameters.

The baseline CDM scenario considered in this work assumes no baryonic
feedback. In this case, the density distribution is given by the NFW profile
with $(\alpha,\beta,\gamma)=(1,3,1)$ and $r_{\scriptscriptstyle{-2}}=r_{s}$.
This case is also a good proxy for CDM scenarios with weak baryonic feedback.

As an example of a feedback-affected CDM halo, we use the DC14 model (Di
Cintio et al., 2014; Di Cintio et al., 2014), which is based on galaxies of
dwarf to Milky Way masses from the Making Galaxies In a Cosmological Context
(MAGICC) project (Stinson et al., 2012). These hydrodynamic MAGICC simulations
account for stellar feedback, which can lead to a cored density profile
depending on the galaxy’s stellar-to-halo mass, $M_{\star}/M_{\rm vir}$. In
this case, the $(\alpha,\beta,\gamma)$ exponents are:

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle
2.94-\log_{10}\left[\left(10^{X+2.33}\right)^{-1.08}+\left(10^{X+2.33}\right)^{2.29}\right]$
$\displaystyle\beta$ $\displaystyle=$ $\displaystyle 4.23+1.34X+0.26X^{2}$
$\displaystyle\gamma$ $\displaystyle=$
$\displaystyle-0.06+\log_{10}\left[\left(10^{X+2.56}\right)^{-0.68}+10^{X+2.56}\right]\,,$

with $X=\log_{10}\left(M_{\star}/M_{\text{vir}}\right)$. The DC14 model is
valid in the mass range $-4.1<\log_{10}(M_{\star}/M_{\rm vir})<-1.3$. Because
the model does not capture potential feedback from active galactic nuclei
(AGN) at masses above this range, these galaxies are excluded from our
analysis after the fitting procedure is completed, further reducing the sample
to 90 galaxies from 98. $M_{\star}$ is estimated by using a combination of the
total and bulge luminosities reported in the SPARC catalog, $L_{\rm tot}$ and
$L_{\rm bul}$, and the mass-to-light ratio parameters,

$M_{\star}\approx\Upsilon_{\star,\rm{disk}}\left(L_{\rm tot}-L_{\rm
bul}\right)+\Upsilon_{\star,\rm{bul}}L_{\rm bul}\,,$ (6)

under the assumption that the gas is a negligible contribution to $L_{\rm
tot}$.

The overall properties of the DC14 model have been confirmed with other
hydrodynamic simulation codes that include stellar feedback (Chan et al.,
2015; Tollet et al., 2016; Lazar et al., 2020; Macciò et al., 2020). In
general, core-formation is suppressed for $M_{\star}/M_{\rm vir}\lesssim
10^{-4}$ and becomes maximal when approaching $M_{\star}/M_{\rm vir}\sim
5\times 10^{-3}$. For $M_{\star}/M_{\rm vir}\gtrsim 10^{-2}$, the halo density
profile becomes preferentially cuspy once again, though the degree to which
this occurs depends on AGN feedback (Macciò et al., 2020). The details on how
rapidly the transitions between cuspy and cored occur as a function of
$M_{\star}/M_{\rm vir}$, as well as the expected degree of scatter about the
median distribution, does vary between simulation codes.

Bose et al. (2019); Dutton et al. (2019); Benítez-Llambay et al. (2019);
Dutton et al. (2020) emphasize that bursty star formation histories are not
the only requirement for core formation. Additionally, one needs a star
formation threshold that is large enough to allow enough gas to build up near
the center of the galaxy, so that when it is blown out by supernova
explosions, it can effectively perturb the DM orbits away from the center of
the galaxy. As a result, simulation codes with lower star formation thresholds
do not see the same level of coring in dwarf systems as captured by the DC14
model.

###  3.2  Mass Model for SIDM

To model the SIDM density distribution, we follow the analytic approach of
Rocha et al. (2013a); Kaplinghat et al. (2014b, 2016), which has been verified
with numerical simulations of isolated SIDM halos (Creasey et al., 2017;
Elbert et al., 2018; Sameie et al., 2018). This procedure requires dividing
the halo into an inner and an outer region, separated at a galactocentric
radius $r_{1}$ where DM particles average one interaction over the age of the
galaxy. Beyond $r_{1}$, SIDM interactions are negligible, and the DM profile
is expected to resemble the CDM result, modeled using the NFW profile. Below
$r_{1}$, the self interactions thermalize the DM distribution, and the density
distribution is well-modeled as an isothermal profile,

$\rho_{\rm
iso}(r)=\rho_{\text{0}}\exp(-\Delta\Phi_{\text{tot}}(r)/\sigma_{\text{v0}}^{2})\,,$
(7)

where $\Delta\Phi_{\text{tot}}\equiv\Phi_{\rm tot}(r)-\Phi_{\rm tot}(0)$ is
the total gravitational potential, $\rho_{0}$ is the central DM density, and
$\sigma_{\text{v0}}$ is the one-dimensional velocity dispersion. Thus, the
SIDM density profile is given by

$\rho_{\rm SIDM}(r)=\begin{cases}\rho_{\rm iso}(r)&r<r_{1}\\\
\rho_{\scriptscriptstyle{\rm NFW}}(r)&r\geq r_{1}\,.\end{cases}$ (8)

The NFW and isothermal distributions are “stitched” together at $r_{1}$ such
that the enclosed mass and density are continuous. Additionally, the value of
$r_{1}$ is set by the relation

$\rho_{\rm
iso}(r_{1})\,t_{\text{age}}\braket{\sigma\,v_{\text{rel}}}/m_{\chi}=1\,,$ (9)

where $\sigma/m_{\chi}$ is the self-scattering cross section per DM particle
mass, $v_{\text{rel}}$ is the relative velocity of the scattering particles,
and $t_{\text{age}}$ is the age of the galaxy. For an isothermal density
distribution and Maxwellian velocity distribution, $\braket{\sigma
v_{\text{rel}}}=\sigma(4/\sqrt{\pi})\sigma_{\text{v0}}$.

|  Parameter |  Units |  Prior Type |  Prior Range
---|---|---|---|---
NFW, DC14 |  $V_{\rm vir}$ |  km/s |  uniform |  [0, 500]
SIDM |  $\Gamma_{0}$ |  $0.1~{}\text{Gyr}^{-1}$ |  log-uniform |  [2, 105]
|  $\sigma_{\text{v0}}$ |  km/s |  uniform |  [2, 500]
All Models |  $\log_{10}c_{\rm vir}$ |  – |  normal |  Eqs. (15), (16)
|  $\Upsilon_{\star,\text{disk}}$ |  $M_{\odot}/L_{\odot}$ |  log-normal |  $\mu=0.5,\sigma=0.1~{}\text{dex}$
|  $\Upsilon_{\star,\text{bul}}$ |  $M_{\odot}/L_{\odot}$ |  log-normal |  $\mu=0.7,\sigma=0.1$ dex
|  $V_{\text{max}}/V_{\text{f}}$ |  – |  uniform |  [2-0.5, 20.5]
|  $M_{\rm b,tot}/M_{\rm vir}$ |  – |  uniform |  $[0,\ 0.2]$
|  |  |  |
Table 1: Priors on the free parameters used to model the dark matter and
baryonic contributions to the rotation curve data from the SPARC dataset. The
free parameters for the NFW and DC14 models are the virial velocity, $V_{\rm
vir}$, and concentration, $c_{\rm vir}$. The prior on $c_{\rm vir}$ follows
the best-fit concentration-mass relation of Dutton & Macciò (2014). The free
parameters for the SIDM model are $\Gamma_{0}$ and $\sigma_{\text{v0}}$, as
defined in Eqs. (7) and (13); the prior on concentration is also enforced for
SIDM halos. In all cases, we place a prior motivated by Li et al. (2020), and
references therein, on the mass-to-light ratios of the stellar disk and bulge,
$\Upsilon_{\star,\text{disk}}$ and $\Upsilon_{\star,\text{bul}}$, as well as a
regularization prior on $V_{\rm max}/V_{\text{f}}$. The benchmark analysis
also requires that the ratio of total baryonic to DM mass, $M_{\rm
b,tot}/M_{\rm vir}$, be less than $0.2$. Sec. 5 discusses how the conclusions
of the study are affected by variations on these benchmark priors.

As a concrete example, we set $t_{\text{age}}=10\ \text{Gyr}$ and
$\sigma/m_{\chi}=3\ \text{cm}^{2}/\text{g}$, but note that the resulting
profile is only weakly dependent on these choices (Elbert et al., 2015;
Sokolenko et al., 2018; Ren et al., 2019). This can be understood intuitively
because typically $r_{1}\approx r_{-2}$ of the corresponding NFW profile and
so there is an approximate scaling $r_{1}\propto\sqrt{t_{\rm
age}\times\sigma/m_{\chi}}$. Therefore, $r_{1}$ is only weakly dependent on
$t_{\rm age}\times\sigma/m_{\chi}$. Because $r_{1}$ depends only on the
product of $t_{\rm age}$ and $\sigma/m_{\chi}$, there is an inherent
degeneracy between the two values.

Given a baryonic density profile, $\rho_{\text{b}}(r)$, the isothermal SIDM
profile can be found via Eqs. (7–9) together with the Poisson equation,

$\nabla^{2}\Phi_{\text{tot}}(r)=4\pi G\left[\rho_{\rm
iso}(r)+\rho_{\text{b}}(r)\right]\,.$ (10)

To obtain an estimate for $\rho_{b}(r)$, the baryonic mass profile is taken to
be spherically symmetric with total mass at any radius $r_{i}$ (where a
measurement is reported) given by,

$\displaystyle M_{\text{b}}(r_{i})=$
$\displaystyle(\Upsilon_{\star,\text{disk}}V_{\text{disk}}\,|V_{\text{disk}}|+\Upsilon_{\star,\text{bul}}V_{\text{bul}}\,|V_{\text{bul}}|$
(11) $\displaystyle\hskip
7.22743pt+\,V_{\text{gas}}\,|V_{\text{gas}}|)r_{i}/G\,.$

An approximate density profile is then created as follows. A function for
$M_{\rm b}(r)$ is created by quadratic interpolation between the smallest and
largest radii, $r_{\rm min}$ and $r_{\rm max}$, at which measurements exist.
The density profile is then taken to be,

$\rho_{\text{b}}(r)=\begin{cases}M_{\text{b}}(r_{\rm min})/(4\pi r_{\rm
min}^{3}/3)&r<r_{\rm min}\\\ \left(dM_{\rm b}(r)/dr\right)/\left(4\pi
r^{2}\right)&r_{\rm min}\leq r\leq r_{\rm max}\\\ \rho_{\text{b}}(r_{\rm
max})\left(r/r_{\rm max}\right)^{-4}&r_{\rm max}<r\,.\end{cases}$ (12)

The rotational velocity associated with the SIDM halo follows from $\rho_{\rm
SIDM}$ through Eq. (4).

There are four free parameters for the SIDM model, which are chosen to be
$\Gamma_{0}$, $\sigma_{\text{v0}}$, $\Upsilon_{\star,\rm{disk}}$ and
$\Upsilon_{\star,\rm{bul}}$. As in Ren et al. (2019), the scattering rate is
defined as,

$\Gamma_{0}=\rho_{0}(\sigma/m_{\chi})(4/\sqrt{\pi})\sigma_{\text{v0}}.$ (13)

We do not consider a feedback-affected SIDM halo, having confirmed that the
best-fit thermalization timescale is typically less than 0.25 Gyr for most of
the SPARC galaxies in our sample. Feedback could potentially change our
results if it acts on a timescale which is much shorter than this. However, as
mentioned above, simulations show that this does not typically occur, even for
larger cross sections than $3$ cm2/g.

###  3.3  Likelihood Procedure

For a model with free parameters $\boldsymbol{\theta}$, the log-likelihood,
$\ln{\mathcal{L}}(\boldsymbol{\theta})$, is defined to be,

$\ln{\mathcal{L}}\left(\boldsymbol{\theta}\right)=-\frac{1}{2}\sum_{i}\left(\frac{V_{i,\text{obs}}-V_{\text{tot}}(\boldsymbol{\theta},r_{i})}{\delta
V_{i,\text{obs}}}\right)^{2}\,,$ (14)

where $V_{i,\text{obs}}$ and $\delta V_{i,\text{obs}}$ are the velocity and
associated uncertainty of the $i^{\text{th}}$ data point in a galaxy’s
rotation curve, and $V_{\text{tot}}(\boldsymbol{\theta},r_{i})$ is the
corresponding model prediction at radius $r_{i}$. The parameter space is
scanned using the affine-invariant Markov Chain Monte Carlo (MCMC) sampler,
emcee (Goodman & Weare, 2010; Foreman-Mackey et al., 2013). We use 48 walkers
and 20,000 steps, with a burn-in period of 5,000 steps. The convergence of the
post burn-in chain for each parameter is verified by ensuring that its length
is at least 50 times the estimated autocorrelation time (Goodman & Weare,
2010; Foreman-Mackey et al., 2013). A limited number of galaxies fail this
test, oftentimes if a parameter’s posterior distribution is not unimodal.
These galaxies are run for an additional 10,000 steps. Those that are still
not converged after this step are NGC6015 for NFW; none for DC14; and NGC0247
for SIDM.

We perform this sampling in total for three different models: NFW, DC14, and
SIDM. The priors used for the benchmark analysis in this study are summarized
in Tab. 1. All the models considered share the same prior on the disk and
bulge mass-to-light ratios and have total baryonic mass to DM mass $M_{\rm
b,tot}/M_{\rm vir}\leq 0.2$. Additionally, as in Ren et al. (2019), we
restrict $2^{-0.5}\leq V_{\rm max}/V_{\rm f}\leq 2^{0.5}$ to enforce that the
predicted flat rotation curve corresponds to the measured $V_{\rm f}$.

Figure 1: Rotation curve data for six example galaxies in the SPARC catalog.
Curves for the NFW, DC14, and SIDM models are indicated by the solid yellow,
red, and aqua lines, respectively. In each case, the central curve and colored
band correspond to the 16–50–84th percentile of the posterior distribution.
Note the change in axis range between panels. _(Top Row)_ NGC0055 is a low
surface brightness galaxy in the SPARC sample where the SIDM and DC14 model do
equally well and are both better fits than the NFW model. All three models fit
the data well for NGC3741 and NGC3769, which are examples of low and high
surface brightness galaxies, respectively. _(Bottom Row)_ When performing a
model comparison using rotation curve data, there can be considerable galaxy-
to-galaxy scatter. It is important to keep in mind that the variability of the
data between galaxies can contribute to this and to therefore look for trends
across an ensemble of galaxies. Provided here are three illustrative examples:
the UGC06787 curve is highly non-monotonic, the IC2574 curve does not include
much data in the region where it is expected to flatten, and the NGC4013 curve
does not extend down to low radii. The relevant rotation curve fits and corner
plots for these six galaxies are provided in App. A; other relevant summary
plots can be found for the remaining SPARC galaxies on the github.

In addition, for all models a prior is placed on the halo concentration-mass
relation, assuming that it follows Dutton & Macciò (2014) with

$\displaystyle\log_{10}c_{\rm vir}=1.025-0.097\log_{10}\left(M_{\rm
vir}/\left[10^{12}h^{-1}M_{\odot}\right]\right)\,,$ (15)

and intrinsic scatter $0.11$ dex. The concentration of a DC14 halo, $c_{\rm
vir}^{\scriptscriptstyle{\rm DC14}}$, is defined in terms of the NFW
concentration using

$c_{\rm vir}^{\scriptscriptstyle{\rm DC14}}=c_{\rm
vir}\,\left(1+0.00001e^{3.4\left(X+4.5\right)}\right)\,,$ (16)

as provided in Di Cintio et al. (2014), but since updated.333Private
communication with A. DiCintio.

Sec. 5 explores several variations on these benchmark priors and their effects
on the conclusions of the study.

##  4 Best-Fit Rotation Curve Examples

The procedure described in Sec. 3 has been performed for all 90 galaxies
within the SPARC catalog that pass the quality cuts. For each galaxy, and for
each of the three models considered in this study, we recover posteriors of
the free parameters and a best-fit rotation curve. Fig. 1 provides examples of
six galaxies from the SPARC catalog together with their best-fit rotation
curves. These examples are not meant to represent average trends, but rather
have been chosen to highlight some relevant aspects of our analysis. The full
set of rotation curves are publicly available (see Footnote 1).

Many galaxies in the dataset prefer the DC14 or SIDM models over NFW, while
other galaxies show either no preference or some preference for NFW. These
results may depend on a galaxy’s surface brightness, defined as
$\Sigma_{0}=\Upsilon^{-1}M_{\star}/(2\pi R_{\text{d}}^{2})$, where
$R_{\text{d}}$ is the disk scale radius and the mass-to-light ratio is taken
to be $\Upsilon=1M_{\odot}/L_{\odot}$. Throughout, we consider low surface
brightness galaxies to be those with
$\log_{10}\left(\Sigma_{0}/\left(L_{\odot}/\text{pc${}^{2}$}\right)\right)<2.5$
and high surface brightness galaxies to be those with values above this range.

In the top row of Fig. 1, NGC0055 is an example of a low surface brightness
galaxy where clearly the NFW model struggles to reproduce the measured
rotation curve at small galactic radii. However, NGC3741 is a counter example
where all three models fit the data well. The galaxy NGC3769 is a high surface
brightness galaxy where again there is no strong preference for any of the
three models considered.

In general, we observe galaxy-to-galaxy variation in the recovered model
preferences, which underscores the importance of studying an ensemble of
galaxies. One contributing factor to this spread is the fact that some
galaxies have more/less rotation curve data available and—in some
cases—unusually shaped curves. The second row of Fig. 1 provides some
illustrative examples. UGC06787 is an example of a galaxy whose rotation curve
is highly non-monotonic, while the reported uncertainties for the rotation
values are remarkably small. Clearly, none of the models considered in this
study are able to capture this behavior. IC2574 is an example of a galaxy that
passes the cut of a well-defined $V_{\text{f}}$, however clearly is not
measured out to the radius at which the rotation becomes constant. Finally,
NGC4013 is an example of a galaxy whose inner-most measurement is only at
$\sim 5$ kpc and thus does not help distinguish any model where differences
become apparent closer to the galactic center.

Comparing the various models requires specifying a test statistic. Prior
rotation curve studies have performed model comparisons using the reduced chi-
squared test (Katz et al., 2016; Li et al., 2019; Ren et al., 2019; Kaplinghat
et al., 2020; Li et al., 2020). In particular, the chi-squared per degree of
freedom, $\chi^{2}_{\nu}$, was evaluated for each galaxy in the SPARC sample,
for a model of interest. In some cases, the cumulative distribution function
(CDF) of the reduced chi-squared values was then plotted. A model was
considered to be preferred over another if it resulted in more galaxies with
lower $\chi^{2}_{\nu}$. There are, however, several fundamental issues with
comparing models in this fashion.

The most important of these issues is that the reduced chi-squared is not an
accurate measure of the goodness-of-fit for the rotation curve models under
consideration. The problem lies in how the degrees of freedom, $\nu$, are
determined (Andrae et al., 2010). By naive degree-of-freedom counting, one
usually defines $\nu$ to be the difference between the number of data points
in the sample, $N$, and the number of fit parameters, $P$. However, as
highlighted in Andrae et al. (2010), this is only correct for a linear model,
$f\left(\boldsymbol{\theta},\boldsymbol{x}\right)$, defined as

$f\left(\boldsymbol{\theta},\boldsymbol{x}\right)=\sum_{p=1}^{P}\theta_{p}B_{p}\left(\boldsymbol{x}\right)\,,$
(17)

where $\theta_{p}$ are the free parameters and $B_{p}\left(\mathbf{x}\right)$
are basis functions. If the basis functions are linearly independent, then the
number of degrees of freedom is given by $\nu=N-P$. If this requirement fails,
then determining $\nu$ is highly non-trivial. Indeed, for non-linear models,
$\nu\in[0,N-1]$ and may not even remain constant during the fitting procedure
(Andrae et al., 2010). The NFW, DC14, and SIDM models defined in Secs. 3.1 and
3.2 are non-linear as they clearly do not take the form of Eq. (17). While the
chi-squared can still be robustly evaluated—and thus the likelihood procedure
remains intact—the _reduced_ chi-squared is not an adequate measure of the
goodness-of-fit.

Figure 2: Rotation curve data for the SPARC galaxy F583-1 (black points). The
curves/bands for each halo mass model indicate the 16–50–84th percentiles of
the posterior distribution. The NFW, DC14, and SIDM models are indicated by
the yellow, red, and aqua lines/bands, respectively. _(Right)_ The normalized
residuals for the SIDM, DC14, and NFW models. The solid black curve
corresponds to a Gaussian distribution with mean $\mu=0$ and standard
deviation $\sigma^{2}=1$. If the distribution of residuals has $\sigma^{2}<1$,
as is the case for the SIDM and DC14 fits to this galaxy, then it indicates an
issue with over-fitting the data.

On its own, this observation makes it impossible to properly interpret the CDF
of $\chi^{2}_{\nu}$. We should also briefly note, however, that even if the
use of reduced chi-squared were appropriate, there would be further issues in
the statistical interpretations made in Katz et al. (2016); Li et al. (2019);
Kaplinghat et al. (2020); Li et al. (2020) that rely on $\chi^{2}_{\nu}$.
First, models with $\chi^{2}_{\nu}\ll 1$ are indicative of over-fitting and
thus should not be interpreted as good fits to the data. Second, model
comparisons should be evaluated on a galaxy-by-galaxy basis by considering the
change in chi-squared, or $\Delta\chi^{2}$, between the different model fits.
However, such comparisons are only valid when comparing nested models,444Two
models are considered nested if one includes all the fit parameters of the
other. which is not the case for the NFW, SIDM, and DC14 cases.

To interpret the goodness-of-fit of each model to the rotation curve data, we
use the distribution of normalized residuals,

$R_{i}(\boldsymbol{\theta})=\frac{V_{i,\text{obs}}-V_{\rm
tot}(\boldsymbol{\theta},r_{i})}{\delta V_{i,\text{obs}}}\,,$ (18)

evaluated for all $i$ data points. The distribution of $R_{i}$ for the true
model will follow a Gaussian distribution with mean $\mu=0$ and standard
deviation $\sigma^{2}=1$. Just because a model results in residuals that
follow this Gaussian distribution does not mean that it is necessarily the
true model, only that it is a viable candidate.

Figure 3: _(Top Row)_ The best-fit concentration, $c_{\rm vir}$, for each
SPARC galaxy plotted as a function of its best-fit halo virial mass for the
NFW (left), DC14 (middle), and SIDM (right) models. The solid black line
corresponds to the best-fit concentration-mass relation from Dutton & Macciò
(2014); the gray bands correspond to $\pm 1$ and $2\sigma$ spread with
$\sigma=0.11$ dex. Note that our benchmark analysis enforces this
concentration-mass relation at the prior level. _(Bottom Row)_ The best-fit
stellar mass for each SPARC galaxy plotted as a function of its best-fit halo
virial mass for the NFW (left), DC14 (middle), and SIDM (right) models. The
solid black line denotes the abundance matching relation from Behroozi et al.
(2019); the gray bands correspond to $\pm 1$ and $2\sigma$ spread with
$\sigma=0.3$ dex. Importantly, the abundance-matching relation is not enforced
at the prior level for the benchmark analysis, but it is clearly recovered by
the fitting procedure. For all panels, each point corresponds to the median
parameter value, with the lines spanning the 16–84${}^{\text{th}}$ percentiles
of the posterior distribution. A galaxy with a suspected over (under)-fitting
issue is denoted by an upwards (downwards)-pointing triangle instead of a
circle. If a galaxy fails the autocorrelation test, it is denoted as a filled
marker. For each CDM panel, the aqua points denote galaxies that favor SIDM,
the yellow (red) points denote galaxies that favor the NFW (DC14) model, and
the gray points denote galaxies with no strong model preference
($|\Delta\text{BIC}|<6$).  Figure 4: The inner slope, $\alpha$, of the best-
fit DM density profile as a function of $M_{\star}/M_{\text{vir}}$ (top row)
and surface brightness $\log_{10}\Sigma_{0}$ (bottom row). The inner slope is
evaluated at 1.5% of the virial radius, $r_{\text{vir}}$. Results are provided
for the NFW (left), DC14 (middle), and SIDM (right) models. When the NFW model
is used, the inner slope is consistently cuspy. In contrast, the DC14 and SIDM
models have enough flexibility to yield cored profiles for low-mass galaxies
and cuspier profiles for higher-mass systems. A similar transition from cored
to cuspy profiles is observed moving from low to high surface brightness
systems. Compared to DC14, the SIDM model yields a larger spread in inner
slopes for the high surface brightness (also, high-mass) galaxies. The plot
formatting is the same as in Fig. 3.

For the rotation curve analysis, it is important to identify galaxies where a
particular model severely over- or under-fits the data. These cases can be
identified if the distribution of residuals is markedly different from a
Gaussian with $\mu=0$ and $\sigma^{2}=1$. To flag these instances, we use the
Kolmogorov-Smirnov (KS) test, which allows one to determine whether an
observed distribution is drawn from some reference probability distribution
function. We also use the Shapiro-Wilk (SW) test, which is sensitive to the
normality of the observed distribution. Whether the KS or SW test is better
able to distinguish the observed distribution of normalized residuals from the
Gaussian distribution depends on a variety of factors, including the number of
measured data points and the dispersion of $R_{i}$ values. As a result, we
identify any galaxy with a potential over/under-fitting issue if it has a KS
or SW p-value of $<0.05$. We require that all galaxies have at least 10 data
points in their rotation curve to ensure that these p-values are robust.

Fig. 2 provides a concrete example of these considerations. The left panel
shows the rotation curve data for F583-1. The corresponding best-fit SIDM,
DC14 and NFW models are also provided. The right panel shows the corresponding
residuals, with the black curve indicating the target normal distribution.
Both the DC14 and SIDM models over-fit the data as the distribution of
residuals has $\sigma^{2}<1$; the NFW residuals do not suffer this issue. This
is also visually apparent from the rotation curve in the left panel, where it
is clear that the aqua and red lines pass almost perfectly through each data
point, despite the size of the error bars on the data.

Of the 90 galaxies studied, 7 are flagged as over-fit for the NFW model, 16
for DC14, and 13 for SIDM. Another 20 are flagged as under-fit for the NFW
model, 16 for DC14 and 16 for SIDM. Between all three models, there are 4 (11)
galaxies that are jointly classified as over-fit (under-fit). Between the DC14
and SIDM models alone, there are 10 (13) galaxies that are jointly classified
as over-fit (under-fit). In general, we find that the under-fitting issues for
the DC14 and SIDM models are overwhelmingly with the high surface brightness
galaxies, and the over-fitting issues are with the low surface brightness
galaxies. For the NFW scenario, we do not find as clear a pattern for the
under-fit and over-fit galaxies.

##  5 Model Comparison

Before comparing the different DM models, one should verify that all three
scenarios considered (NFW, DC14, and SIDM) recover well-understood
concentration-mass and stellar mass-halo mass (SMHM) relations. Fig. 3
summarizes our results for the NFW, DC14 and SIDM scenarios, shown in the
left, middle, and right columns, respectively. If a galaxy is indicated by a
triangle, then it failed the goodness-of-fit criteria; an upwards-pointing
triangle corresponds to an over-fit model and a downwards-pointing triangle
corresponds to an under-fit model. If a galaxy is denoted by a filled marker,
then it failed the autocorrelation test.

In the top row, the concentration, $c_{\rm vir}$, is plotted as a function of
halo virial mass, $M_{\rm vir}$. The best-fit Planck concentration-mass
relation from Dutton & Macciò (2014) is indicated by the black line, with the
gray bands showing the $\pm 1$ and $2\sigma$ spread with $\sigma=0.11$ dex.
This relation was enforced with a log-normal prior (see Tab. 1). The recovered
concentrations correspond reasonably well with the theoretical expectation for
all three cases, although with larger scatter. In particular, we find that
approximately 37% (43%) (36%) of the galaxies fall inside the $\pm 1\sigma$
band for NFW (DC14) (SIDM); the corresponding fraction of galaxies within the
$\pm 2\sigma$ band is 66% (74%) (74%).

The recovered SMHM for each model is provided in the bottom row of Fig. 3. The
best-fit distribution from Behroozi et al. (2019) is also shown for
comparison. The gray bands correspond to an approximate $\pm 1$ and $2\sigma$
spread with $\sigma=0.3$ dex.555While the best-fit scatter in stellar mass
depends on the halo mass (Behroozi et al., 2019), this approximation is a good
proxy for the range of halo masses relevant to the SPARC catalog. It should be
stressed that the SMHM relation has not been included at the prior level.
Indeed, the only requirement on the stellar mass in our benchmark analysis
comes from the conservative prior $M_{\rm b,tot}/M_{\rm vir}\leq 0.2$. The
recovered SMHM relation roughly tracks the expectation for all three models,
although the DC14 case yields a dispersion that is most closely aligned with
the expectation. In this case, we find that approximately 36% (52%) (51%) of
the galaxies fall inside the $\pm 1\sigma$ band for NFW (DC14) (SIDM); the
corresponding fraction of galaxies within the $\pm 2\sigma$ band is 70% (81%)
(69%).

Figure 5: Distribution of $\Delta$BIC values for all SPARC galaxies with
$\log_{10}\left(\Sigma_{\rm
eff}/\left(L_{\odot}/\text{pc}^{2}\right)\right)<2.5$ (left) and $>2.5$
(right). Positive values indicate a preference for SIDM over NFW (orange) and
DC14 (red). The vertical gray band indicates regions of $\Delta$BIC where the
preference for either model is not strong. For each panel, the colored box and
whisker plots summarize the spread of the corresponding histogram. The box
extends from the first to the third quartile of $\Delta$BIC values, with the
median indicated by the vertical black line. The whiskers extend a factor of
1.5 times the interquartile range. In general, the SIDM model provides a
better fit to the rotation curve data than the NFW model, especially for low
effective surface brightness galaxies. Additionally, there is no significant
preference for either the SIDM or DC14 models for either the low or high
effective surface brightness galaxies.

For each galaxy in the sample, we compare the Bayesian Information Criterion
(BIC) obtained for each model and consider

$\Delta\text{BIC}=\text{BIC(model 1)}-\text{BIC(model 2)}\,.$ (19)

When $\Delta$BIC $>6$, this indicates a strong preference for Model 2 over
Model 1 and vice versa when $\Delta$BIC $<-6$ (Kass & Raftery, 1995). For all
other cases, we consider the two models statistically indistinguishable. In
Fig. 3, we illustrate the preference for SIDM against either NFW or DC14. Each
point in the left-most and central panels is colored by its $\Delta\text{BIC}$
value with aqua indicating a strong preference for SIDM, yellow indicating a
strong preference for NFW, and red indicating a strong preference for DC14.
Gray indicates no preference between SIDM and the CDM model being considered.

The SIDM model is preferred over NFW for many of the galaxies over the entire
range of masses. In contrast, most galaxies in the low-mass range show no
strong preference between SIDM and DC14. At high masses, there are
approximately an equal number of galaxies with preference for either DC14 or
SIDM. This suggests that both the SIDM and DC14 models provide enough
flexibility to capture the diversity of the rotation curves in the lowest-mass
SPARC systems, especially compared to NFW. Of galaxies with larger virial
masses, some exhibit a strong preference for one model over another, but there
is no systematic preference for any of the models when the galaxies are
considered in aggregate.

Fig. 4 illustrates how the inner slope of the best-fit density profile,
$\alpha$, depends on $M_{\star}/M_{\text{vir}}$ (top row) and on surface
brightness $\log_{10}\Sigma_{0}$ (bottom row). As in Fig. 3, the $\Delta$BIC
value for each galaxy is indicated by the color of the point; if a data point
is marked as a triangle, it indicates that the galaxy failed the goodness-of-
fit criteria. The inner slope is evaluated at 1.5% of the virial radius. In
the NFW case, the best-fit density profiles are consistently cuspy, as is
enforced by the shape of the density profile, with $\alpha\in[-1.5,-1]$ across
the complete galaxy sample. For the DC14 case, the low surface brightness
galaxies exhibit coring, with best-fit inner slopes closer to $\alpha\sim 0$.
Towards higher-surface brightness galaxies, the inner-slopes decrease in value
and the density profiles become cuspier, similar to the expectation for NFW
profiles. As for the DC14 model, the SIDM model results in cored low surface
brightness galaxies. At the high surface brightness range, the SIDM fits
result in a wide range of inner slopes, spanning $\alpha\in[-3,-0.5]$.666Note
that our analytic model for the SIDM density profile does not account for
adiabatic contraction of the halo (Blumenthal et al., 1986; Gnedin et al.,
2004), which can have an effect for $\log_{10}M_{\star}/M_{\rm vir}\gtrsim-2$
and increase the inner density of the DM profile. The inclusion of adiabatic
contraction effects in SIDM halos will be explored in future work (see Jiang
(2021)).

Fig. 5 shows the histogram of $\Delta$BIC values in the low and high
_effective_ surface brightness regimes, comparing either NFW or DC14 to SIDM
(the effective surface brightness is reported in the SPARC catalog). These
histograms underscore the need of performing such an analysis over a large
ensemble of galaxies. Indeed, the galaxy-to-galaxy scatter in the $\Delta$BIC
values can be quite substantial, and it is thus important to look at the
trends of the ensemble.

In the low effective surface brightness regime with
$\log_{10}\left(\Sigma_{\rm
eff}/\left(L_{\odot}/\text{pc}^{2}\right)\right)<2.5$, SIDM is preferred over
NFW with the 25-50-75th percentiles of the $\Delta$BIC values being 1.39,
6.81, and 23.24. In contrast, there is no significant preference for SIDM over
DC14, with the 25-50-75th percentiles of the $\Delta$BIC values being
$-1.63,-0.22,$ and $1.10$. In the high effective surface brightness regime,
the preference for the SIDM model over NFW continues, although it is much
weaker. In this case, the 25-50-75th percentiles of the $\Delta$BIC values are
$-2.69,0.71,$ and $11.79$. There continues to be no significant preference for
SIDM over DC14 in this regime, with the 25-50-75th percentiles of the
$\Delta$BIC values being $-6.46,-0.37,$ and $8.23$. These results underscore
that the SIDM model is more successful at fitting the rotation curves than the
NFW model, but that this behavior is ameliorated when using the DC14 model,
which allows for coring due to baryonic feedback.

To better understand how our conclusions depend on some key assumptions, we
perform several variations of the benchmark analysis. For one test, we
consider the effect of imposing the abundance-matching relation from Behroozi
et al. (2019) at the prior level, as opposed to the prior on $M_{\rm
b,tot}/M_{\rm vir}$. Specifically, we use the best-fit parameters associated
with all star-forming galaxies (row 6 in their Table J1) and assume a constant
scatter on the stellar mass of $\sigma\sim 0.3$ as a function of $M_{\rm
vir}$. For a separate test, we follow Tab. 1, but use a linear-flat prior on
the disk and bulge mass-to-light ratios:
$\Upsilon_{\star,\text{disk}},\Upsilon_{\star,\text{bul}}\in[0.1,10]$. For
both tests, the results are highly consistent with those presented in Figs.
3–5. All relevant plots for these tests are available at the web address in
Footnote 1.

As an outcome of our study, we are able to scrutinize individual rotation
curves and their best-fit results for each model. A recent study by Santos-
Santos et al. (2018) found that NIHAO simulated galaxies were unable to
reproduce SPARC rotation curves that have either very small or order-unity
values for the ratio of circular velocity at 2 kpc to circular velocity at
last measured radius. In Appendix B, we present our results for five galaxies
that are both in our sample and are marked as outliers in Santos-Santos et al.
(2018). We find that only one out of the five has a substantial preference for
DC14 over SIDM, with the others showing no strong preference between these two
models. We also find that many of these galaxies are outliers in terms of
their concentrations. The latter result could potentially explain why NIHAO
was unable to reproduce galaxies of this type.

A final point worth discussion is that in this study we consider a constant
cross section per unit mass of $3$ cm2/g for the SIDM scenario. That said,
constraints on SIDM parameter space point to a velocity-dependent cross
section and therefore one should think of the constant value chosen here as an
approximation to the cross section at the typical velocities of SPARC galaxies
in our sample. However, we find that these galaxies actually exhibit an order
of magnitude spread in $\sigma_{\rm v0}$ of around $30$–$250$ km/s, and
therefore any velocity-dependent cross section cannot vary too much over that
range. This requirement is likely to be in conflict with the combination of
constraints arising from cluster and galaxy group relaxation (Sagunski et al.,
2021) and from central densities of satellite galaxies (Slone et al., 2021).
Thus, a more dedicated study of this possibility is warranted.

##  6 Conclusions

In this paper, we use the SPARC dataset to study several solutions to the
galaxy diversity problem. We focus on three specific models that capture the
relevant physics of SIDM and CDM with/without baryonic feedback. For the
latter, we use the DC14 model (Di Cintio et al., 2014) as an example of a
feedback-affected halo and an NFW model as an example of a CDM halo with no
(or weak) feedback. Improving upon the statistical methodology used in
previous studies of SPARC galaxies, we perform a comprehensive model
comparison to better understand the preference for the three models we
consider.

Our benchmark analysis, which takes the Planck concentration-mass relation of
Dutton & Macciò (2014) as a prior, roughly recovers the SMHM relation of
Behroozi et al. (2019) for all three scenarios. Both the SIDM and DC14 cases
return profiles that are more cored at low surface brightness than the NFW
expectation. For high surface brightness galaxies, all three models predict
cuspier profiles, although the SIDM model results in the largest spread in
best-fit inner slopes.

Because of galaxy-to-galaxy variance in the model preferences, it is important
to consider the trends of the sample in aggregate. Overall, we find that the
high surface brightness galaxies in the SPARC catalog have no strong
statistical preference for SIDM over DC14 model and only a weak preference for
SIDM over NFW. For the low surface brightness galaxies, the SIDM model is
strongly preferred over NFW, but is as good a fit to the data as the DC14
model. We thus conclude that both SIDM as well as CDM with baryonic feedback
provide adequate flexibility to explain the diversity of rotation curves in
the SPARC catalog, and that the current data is not adequate to differentiate
between these two scenarios. The underlying reason for this is that the
specific range of $M_{\star}/M_{\rm vir}$ for which coring enforced by SIDM or
by baryonic feedback is most similar is also the approximate range probed by
the SPARC catalog. In contrast, lower mass halos would demonstrate inherent
differences between the models because $M_{\star}/M_{\rm vir}$ decreases with
the virial mass, so coring from SIDM remains efficient even when coring from
baryonic feedback terminates (Robles et al., 2017; Lazar et al., 2020). Thus,
one of the main conclusions of this study is that distinguishing between these
types of physical scenarios requires rotation curve data for more low-mass
objects than currently available in SPARC.

##  7 Acknowledgements

The authors gratefully acknowledge P. Behroozi, A. Brooks, A. Di Cintio, F.
Jiang, M. Kaplinghat, H. Liu, S. Mishra-Sharma, M. Moschella, and L. Necib for
fruitful conversations. SD and AZ were supported in part by the Princeton
Physics Department’s undergraduate summer research program. ML and OS are
supported by the DOE under Award Number DE-SC0007968 and the Binational
Science Foundation (grant No. 2018140). This work was performed in part at the
Aspen Center for Physics, which is supported by NSF grant PHY-1607611. The
work presented in this paper was performed on computatational resources
managed and supported by Princeton Research Computing. This research made
extensive use of the following publicly available codes: IPython (Perez &
Granger, 2007), Jupyter (Kluyver et al., 2016), matplotlib (Hunter, 2007),
NumPy (van der Walt et al., 2011), SciPy (Jones et al., 2001), emcee (Foreman-
Mackey et al., 2013), and Pandas (Wes McKinney, 2010).

##  Appendix A Corner Plots

This appendix provides rotation curves and corresponding corner plots for the
NFW, DC14, and SIDM rotation curve fits performed for NGC0055, NGC3741,
NGC3769, UGC06787, IC2574, and NGC4013—the six example galaxies spotlighted in
Fig. 1. The posterior distributions are well-converged in all six cases. Each
galaxy also passes the autocorrelation test described in the main text (see
figures provided in the github).

Figure A1: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the NGC0055 rotation curve with corresponding corner plots.

Figure A2: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the NGC3741 rotation curve with corresponding corner plots.

Figure A3: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the NGC3769 rotation curve with corresponding corner plots.

Figure A4: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the UGC06787 rotation curve with corresponding corner plots.

Figure A5: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the IC2574 rotation curve with corresponding corner plots.

Figure A6: NFW (left column), DC14 (middle column), and SIDM (right column)
fits to the NGC4013 rotation curve with corresponding corner plots.

##  Appendix B Outlier Galaxies

This appendix provides the best-fit rotation curves (Fig. A7) and
concentration values (Fig. A8) for galaxies that are both in our sample and
are marked as outliers in Santos-Santos et al. (2018). Specifically, these are
galaxies which seem to have values of the ratio of $V_{\rm tot}$ at 2 kpc to
its value at the last measured radius that are not reproduced in NIHAO
simulations. The inset within each rotation curve provide values for
$\Delta$BIC between NFW and SIDM, and between DC14 and SIDM. We find that only
one out of the 5 galaxies (IC2574) shows a strong preference for DC14 over
SIDM, with the other showing no strong preference between these two cases. The
other galaxies have $\Delta$BIC values too small to make any conclusive
claims. Additionally, as shown in Fig. A8, we find that these galaxies tend to
be outliers in terms of their concentration relations. This could potentially
explain why simulations find it challenging to reproduce galaxies of this
type.

Figure A7: Best-fit rotation curves for galaxies which are both in our sample
and are marked as outliers in Santos-Santos et al. (2018). The insets within
each of the panels provide values for $\Delta$BIC between NFW and SIDM, and
between DC14 and SIDM. Positive values correspond to a preference for SIDM
over the CDM model. Figure A8: Concentration relations for galaxies in our
sample with those marked as outliers by Santos-Santos et al. (2018) shown in
green.

## References

  * Allaert et al. (2017) Allaert, F., Gentile, G., & Baes, M. 2017, Astron. Astrophys., 605, A55, doi: 10.1051/0004-6361/201730402
  * Andrae et al. (2010) Andrae, R., Schulze-Hartung, T., & Melchior, P. 2010, arXiv e-prints, arXiv:1012.3754.  https://arxiv.org/abs/1012.3754
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Benítez-Llambay et al. (2019) Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387, doi: 10.1093/mnras/stz1890
  * Blumenthal et al. (1986) Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27, doi: 10.1086/163867
  * Bose et al. (2019) Bose, S., et al. 2019, MNRAS, 486, 4790, doi: 10.1093/mnras/stz1168
  * Bryan & Norman (1998) Bryan, G. L., & Norman, M. L. 1998, ApJ, 495, 80, doi: 10.1086/305262
  * Burger et al. (2022) Burger, J. D., Zavala, J., Sales, L. V., et al. 2022, MNRAS, doi: 10.1093/mnras/stac994
  * Chan et al. (2015) Chan, T. K., Kereš, D., Oñorbe, J., et al. 2015, MNRAS, 454, 2981, doi: 10.1093/mnras/stv2165
  * Creasey et al. (2017) Creasey, P., Sameie, O., Sales, L. V., et al. 2017, MNRAS, 468, 2283, doi: 10.1093/mnras/stx522
  * de Blok et al. (2001) de Blok, W. J. G., McGaugh, S. S., Bosma, A., & Rubin, V. C. 2001, ApJ, 552, L23–L26, doi: 10.1086/320262
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014, MNRAS, 441, 2986, doi: 10.1093/mnras/stu729
  * Di Cintio et al. (2014) Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415, doi: 10.1093/mnras/stt1891
  * Dutton et al. (2020) Dutton, A. A., Buck, T., Macciò, A. V., et al. 2020, MNRAS, 499, 2648, doi: 10.1093/mnras/staa3028
  * Dutton et al. (2019) Dutton, A. A., Macciò, A. V., Buck, T., et al. 2019, MNRAS, 486, 655, doi: 10.1093/mnras/stz889
  * Dutton & Macciò (2014) Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359–3374, doi: 10.1093/mnras/stu742
  * Elbert et al. (2015) Elbert, O. D., Bullock, J. S., Garrison-Kimmel, S., et al. 2015, MNRAS, 453, 29, doi: 10.1093/mnras/stv1470
  * Elbert et al. (2018) Elbert, O. D., Bullock, J. S., Kaplinghat, M., et al. 2018, ApJ, 853, 109, doi: 10.3847/1538-4357/aa9710
  * Flores & Primack (1994) Flores, R. A., & Primack, J. R. 1994, ApJ, 427, L1, doi: 10.1086/187350
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306–312, doi: 10.1086/670067
  * Fry et al. (2015) Fry, A. B., Governato, F., Pontzen, A., et al. 2015, MNRAS, 452, 1468, doi: 10.1093/mnras/stv1330
  * Gnedin et al. (2004) Gnedin, O. Y., Kravtsov, A. V., Klypin, A. A., & Nagai, D. 2004, ApJ, 616, 16, doi: 10.1086/424914
  * Gnedin & Zhao (2002) Gnedin, O. Y., & Zhao, H. 2002, MNRAS, 333, 299, doi: 10.1046/j.1365-8711.2002.05361.x
  * Goodman & Weare (2010) Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65 , doi: 10.2140/camcos.2010.5.65
  * Governato et al. (2010) Governato, F., Brook, C., Mayer, L., et al. 2010, Nature, 463, 203, doi: 10.1038/nature08640
  * Hernquist (1990) Hernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Jiang (2021) Jiang, F. 2021, Private Communication
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python.  http://www.scipy.org/
  * Kamada et al. (2017) Kamada, A., Kaplinghat, M., Pace, A. B., & Yu, H.-B. 2017, Phys. Rev. Lett., 119, doi: 10.1103/physrevlett.119.111102
  * Kaplinghat et al. (2014a) Kaplinghat, M., Keeley, R. E., Linden, T., & Yu, H.-B. 2014a, Phys. Rev. Lett., 113, 021302, doi: 10.1103/PhysRevLett.113.021302
  * Kaplinghat et al. (2014b) —. 2014b, Phys. Rev. Lett., 113, 021302, doi: 10.1103/PhysRevLett.113.021302
  * Kaplinghat et al. (2020) Kaplinghat, M., Ren, T., & Yu, H.-B. 2020, J. Cosmology Astropart. Phys, 2020, 027, doi: 10.1088/1475-7516/2020/06/027
  * Kaplinghat et al. (2016) Kaplinghat, M., Tulin, S., & Yu, H.-B. 2016, Phys. Rev. Lett., 116, doi: 10.1103/physrevlett.116.041302
  * Kass & Raftery (1995) Kass, R., & Raftery, A. 1995, Journal of the American Statistical Association, 773\.  http://www.stat.washington.edu/people/raftery/Research/PDF/kass1995.pdf
  * Katz et al. (2016) Katz, H., Lelli, F., McGaugh, S. S., et al. 2016, MNRAS, 466, 1648, doi: 10.1093/mnras/stw3101
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in ELPUB
  * Kuzio de Naray et al. (2010) Kuzio de Naray, R., Martinez, G. D., Bullock, J. S., & Kaplinghat, M. 2010, ApJ, 710, L161–L166, doi: 10.1088/2041-8205/710/2/l161
  * Lazar et al. (2020) Lazar, A., Bullock, J. S., Boylan-Kolchin, M., et al. 2020, MNRAS, 497, 2393, doi: 10.1093/mnras/staa2101
  * Lelli et al. (2016a) Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016a, AJ, 152, 157, doi: 10.3847/0004-6256/152/6/157
  * Lelli et al. (2016b) —. 2016b, ApJ, 816, L14, doi: 10.3847/2041-8205/816/1/L14
  * Li et al. (2020) Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, Astrophys. J. Suppl., 247, 31, doi: 10.3847/1538-4365/ab700e
  * Li et al. (2019) Li, P., Lelli, F., McGaugh, S. S., Starkman, N., & Schombert, J. M. 2019, MNRAS, 482, 5106, doi: 10.1093/mnras/sty2968
  * Macciò et al. (2020) Macciò, A. V., Crespi, S., Blank, M., & Kang, X. 2020, MNRAS, 495, L46, doi: 10.1093/mnrasl/slaa058
  * Madau et al. (2014) Madau, P., Shen, S., & Governato, F. 2014, ApJ, 789, L17, doi: 10.1088/2041-8205/789/1/L17
  * Navarro et al. (1996) Navarro, J. F., Eke, V. R., & Frenk, C. S. 1996, MNRAS, 283, L72, doi: 10.1093/mnras/283.3.L72
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Oh et al. (2011) Oh, S.-H., de Blok, W. J. G., Brinks, E., Walter, F., & Kennicutt, R. C. 2011, AJ, 141, 193, doi: 10.1088/0004-6256/141/6/193
  * Oman et al. (2015) Oman, K. A., Navarro, J. F., Fattahi, A., et al. 2015, MNRAS, 452, 3650–3665, doi: 10.1093/mnras/stv1504
  * Perez & Granger (2007) Perez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Peter et al. (2013) Peter, A. H. G., Rocha, M., Bullock, J. S., & Kaplinghat, M. 2013, MNRAS, 430, 105–120, doi: 10.1093/mnras/sts535
  * Planck Collaboration (2016) Planck Collaboration. 2016, Astron. and Astrophys., 594, A13, doi: 10.1051/0004-6361/201525830
  * Pontzen & Governato (2012) Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464, doi: 10.1111/j.1365-2966.2012.20571.x
  * Read & Gilmore (2005) Read, J. I., & Gilmore, G. 2005, MNRAS, 356, 107, doi: 10.1111/j.1365-2966.2004.08424.x
  * Ren et al. (2019) Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Phys. Rev. X, 9, 031020, doi: 10.1103/PhysRevX.9.031020
  * Robertson et al. (2018) Robertson, A., et al. 2018, MNRAS, 476, L20, doi: 10.1093/mnrasl/sly024
  * Robles et al. (2017) Robles, V. H., Bullock, J. S., Elbert, O. D., et al. 2017, MNRAS, 472, 2945, doi: 10.1093/mnras/stx2253
  * Rocha et al. (2013a) Rocha, M., Peter, A. H. G., Bullock, J. S., et al. 2013a, MNRAS, 430, 81, doi: 10.1093/mnras/sts514
  * Rocha et al. (2013b) —. 2013b, MNRAS, 430, 81–104, doi: 10.1093/mnras/sts514
  * Sagunski et al. (2021) Sagunski, L., Gad-Nasr, S., Colquhoun, B., Robertson, A., & Tulin, S. 2021, JCAP, 01, 024, doi: 10.1088/1475-7516/2021/01/024
  * Sameie et al. (2018) Sameie, O., Creasey, P., Yu, H.-B., et al. 2018, MNRAS, 479, 359, doi: 10.1093/mnras/sty1516
  * Sameie et al. (2021) Sameie, O., Boylan-Kolchin, M., Sanderson, R., et al. 2021, MNRAS, 507, 720, doi: 10.1093/mnras/stab2173
  * Santos-Santos et al. (2018) Santos-Santos, I. M., Di Cintio, A., Brook, C. B., et al. 2018, MNRAS, 473, 4392, doi: 10.1093/mnras/stx2660
  * Santos-Santos et al. (2020) Santos-Santos, I. M. E., Navarro, J. F., Robertson, A., et al. 2020, MNRAS, 495, 58, doi: 10.1093/mnras/staa1072
  * Slone et al. (2021) Slone, O., Jiang, F., Lisanti, M., & Kaplinghat, M. 2021.  https://arxiv.org/abs/2108.03243
  * Sokolenko et al. (2018) Sokolenko, A., Bondarenko, K., Brinckmann, T., et al. 2018, JCAP, 12, 038, doi: 10.1088/1475-7516/2018/12/038
  * Spergel & Steinhardt (2000) Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760, doi: 10.1103/PhysRevLett.84.3760
  * Stinson et al. (2012) Stinson, G. S., Brook, C., Macciò, A. V., et al. 2012, MNRAS, 428, 129, doi: 10.1093/mnras/sts028
  * Tollet et al. (2016) Tollet, E., Macciò, A. V., Dutton, A. A., et al. 2016, MNRAS, 456, 3542, doi: 10.1093/mnras/stv2856
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vargya et al. (2021) Vargya, D., Sanderson, R., Sameie, O., et al. 2021, arXiv e-prints, arXiv:2104.14069.  https://arxiv.org/abs/2104.14069
  * Vogelsberger et al. (2012) Vogelsberger, M., Zavala, J., & Loeb, A. 2012, MNRAS, 423, 3740–3752, doi: 10.1111/j.1365-2966.2012.21182.x
  * Vogelsberger et al. (2014) Vogelsberger, M., Zavala, J., Simpson, C., & Jenkins, A. 2014, MNRAS, 444, 3684, doi: 10.1093/mnras/stu1713
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Zavala et al. (2013) Zavala, J., Vogelsberger, M., & Walker, M. G. 2013, MNRAS Lett., 431, L20–L24, doi: 10.1093/mnrasl/sls053