Reconciling the Diversity and Uniformity of Galactic Rotation Curves with Self-Interacting Dark Matter

Galactic rotation curves exhibit diverse behavior in the inner regions, while obeying an organizing principle, i.e., they can be approximately described by a radial acceleration relation or the Modified Newtonian Dynamics phenomenology. We analyze the rotation curve data from the SPARC sample, and explicitly demonstrate that both the diversity and uniformity are naturally reproduced in a hierarchical structure formation model with the addition of dark matter self-interactions. The required concentrations of the dark matter halos are fully consistent with the concentration-mass relation predicted by the Planck cosmological model. The inferred stellar mass-to-light ($3.6 \mu m$) ratios scatter around $0.5 M_\odot/L_\odot$, as expected from population synthesis models, leading to a tight radial acceleration relation and baryonic Tully-Fisher relation. The inferred stellar-halo mass relation is consistent with the expectations from abundance matching. These results indicate that the inner dark matter halos of galaxies are thermalized due to the self-interactions of dark matter particles.


I. INTRODUCTION
Galactic rotation curves of spiral galaxies show a variety of behavior in the inner parts even across systems with similar halo and stellar masses, which lacks a self-consistent explanation in the standard cold dark matter (CDM) model [1][2][3][4][5][6][7][8][9][10][11][12][13]. Along with this diversity, a long-standing observation is that many rotation curves can be understood in terms of Modified Newtonian Dynamics (MOND) phenomenology [14,15] (see [16] for a review), i.e., there exists a characteristic gravitational acceleration scale, g † ≈ 10 −10 m/s 2 ∼ cH 0 /7 with H 0 being the present Hubble expansion rate, below which the observed acceleration can be approximated as √ g † g bar with g bar being the baryonic acceleration (a.k.a. Milgrom's law). More recently, McGaugh et al. [17] analyzed the Spitzer Photometry and Accurate Rotation Curves (SPARC) dataset [18] and showed there is a tight relation between the total gravitational acceleration at any radius and the acceleration contributed by the baryons, assuming a constant stellar mass-to-light ratio Υ ,disk = 0.5M /L and Υ ,bulge = 0.7M /L in the 3.6 µm band. The scatter in this radial acceleration relation (RAR) is around 0.1 dex, and the tightness of this relation has been interpreted as a signature of MOND [19]. It has long been argued that the acceleration scale (including the cH 0 dependence) can emerge from hierarchical structure formation predicted in CDM [20,21]. Recent hydrodynamical simulations of galaxy formation with CDM have clearly shown that a RAR emerges [22][23][24]. However, these simulated galaxies do not represent the full range of the diversity in the SPARC dataset and they cannot yet explain the rotation curves of low and high surface brightness galaxies simultaneously.
In this paper, we show that self-interacting dark matter (SIDM) provides a unified way to understand the diverse rotation curves of spiral galaxies, while reproducing the RAR with a small scatter. We analyze the SPARC dataset based on the SIDM halo model proposed in [25,26] and demonstrate three key observations leading to this result.
• For cross section per unit mass σ/m 1 cm 2 /g, dark matter self-interactions only thermalize the inner regions at distances less than about 10% of the virial radius of galactic halos, while the outer regions remains unchanged. Thus, SIDM inherits essential features of the ΛCDM hierarchical structure formation model such as the halo concentration-mass relation, which sets the characteristic acceleration scale of halos.
• In the inner halo, thermalization ties dark matter and baryon distributions together [25,27,28], and the SIDM halo can naturally accommodate the diverse range of 'cored' and 'cusped' central density profiles, depending on how the baryons are distributed. Combined with the scatter in the concentration-mass relation, this provides the diversity required to explain the rotation curves [26,29,30]. We will demonstrate the SIDM fits are systematically superior to the MOND ones.
• For the same σ/m that addresses the diversity problem, the baryon content of the galaxies and the mass model of their host halos also lead to the RAR with a scatter as small as the one in [17]. In our SIDM fits, the inferred stellar Υ ,disk values for individual galaxies have a distribution peaked toward 0.5M /L , as expected from stellar population synthesis models [31].
The rest of the paper is organized as follows. In Sec. II, we present the SIDM fits to 135 galaxies from the SPARC sample, which exemplify the full range of the diversity. In Sec. III, we show the radial acceleration relation and the distribution of the stellar mass-to-light ratios from our SIDM fits, compared to the MOND fits. In Sec. IV, we discuss the host halo properties and the origin of the acceleration scale. In Sec. V, we show the predicted stellar -halo mass relation and the baryonic Tully-Fisher to the diverse rotation curves across a range of spiral galaxy masses, where we take σ/m = 3 cm 2 /g. The data points with error bars are from the SPARC dataset [18]. Each panel contains 14 galactic rotation curves that are selected to have similar flat rotation velocities at their furthest radial data points, and the corresponding V f bins are 79-91, 91-126, 139-172 and 239-315 km/s, spanning the mass range of the galaxies considered in this work. The galaxies are colored according to their relative surface brightness in each panel from low (red) to high (violet). relation (BTFR). We comment on future directions and conclude in Sec. VI. In the appendix, Methods, we provide detailed information about the model and the fitting procedure. In Supplementary Materials, we present SIDM and MOND fits to 135 individual galaxies from the SPARC sample and additional results that support the main text, including model fits to simulated halos.

II. THE DIVERSITY OF GALACTIC ROTATION CURVES
We select 135 out of 175 galaxies in the full SPARC sample based on the criteria that they must have a recorded value for the flat part of the rotation curve, V f . In our sample, 87, 42 and 6 galaxies have quality flags 1, 2 and 3, respectively. It spans a wide range of galaxy masses and inner shapes of rotation curves with V f ranging from 20 km/s to 300 km/s. In fitting to the data, we utilize the analytical SIDM halo model [26,29], where we assume the dark matter distribution in the inner halo follows the isothermal density profile, where ρ 0 is the central dark matter density, σ v0 is the one-dimensional dark matter velocity dispersion, Φ tot (R, z) is the total gravitational potential and R, z are cylindrical coordinates aligned with the stellar disk. We match this isothermal profile to a Navarro-Frenk-White (NFW) form [32,33] at r 1 , where a dark matter particle has scattered O(1) times over the age of the galaxy, assuming continuity in both the density and the enclosed mass at r 1 . In this way, the isothermal parameters (ρ 0 , σ v0 ) directly map on to the NFW parameters (r s , ρ s ) or (r max , V max ). This model provides an approximate way to calculate the SIDM distribution in a halo if its CDM counterpart is known, and vice versa. It correctly predicts the halo central density and its scalings with the outer halo properties, stellar profiles and cross section, as confirmed in both isolated and cosmological N-body simulations with and without baryons, see, e.g., [26,28,30,34,35]. See Methods and Supplementary Materials for a detailed description of the model and additional comparisons between model predictions and cosmological simulations.
We adopt two independent but complementary approaches to perform the analysis. In the controlled sampling (CS) approach, we demand that the host halos follow the concentration-mass relation within a 2σ range predicted in cosmological simulations [36]. We model the stellar distribution as an axisymmetric thin disk as in [29], which directly enters into the calculation of the density profile of SIDM through the gravitational potential Φ(R, z). In the CS fits, we start with the outer NFW halo and find the SIDM density profile that matches its mass and density at r 1 . In the second approach, we use the Markov Chain Monte Carlo (MCMC) sampling (MS) to explore the full likelihood. To save computational time, we assume spherical symmetry by spreading the mass within the disk at radius R into a sphere of the same radius [25,28]. The rotation curves generated from two approaches agree well and the differences in the fits are small (see Supplementary Materials). For our main results, we show inferences from both of the approaches.
In Fig. 1, we show the SIDM fits to the diverse rotation curves from the controlled sampling with σ/m = 3 cm 2 /g. In each panel, galaxies are selected to have similar flat rotation velocities at their outermost data points. The rise up to V f within their central regions displays a wide variety of slopes and the SIDM halo model provides equally good fits to the shallow and steeply rising rotation curves. The fits for the other galaxies in the sample are as good as those in Fig. 1 (see Supplementary Materials).
The success of the SIDM halo model stems from a combination of the following effects. First, SIDM thermalization ties the baryon and dark matter distributions together. For low surface brightness galaxies, thermalization leads to a shallow density core and a circular velocity profile that rises mildly with radius [37][38][39][40][41][42][43]. While, for high surface brightness ones, the core shrinks in response to the deeper baryonic potential and the central SIDM density increases accordingly [25,28,30,35]. The galaxies in our sample have a variety of central surface brightnesses, resulting in diverse central dark matter densities. Second, scatter in the cosmological halo concentration-mass relation leads to scatter in the characteristic SIDM core density and radius, which is reflected in the rotation curves [26]. Ref. [29] fitted 30 galaxies and illustrated the importance of these effects in explaining the diverse rotation curves. In this work, we fit a larger sample of galaxies and demonstrate that the observed galaxies are fully consistent with the SIDM predictions.
We have assumed a constant cross section to fit the SPARC sample because it is hard to pin down the cross section for individual galaxies. For low surface brightness galaxies with a large core, a large cross section, such as σ/m = 3 cm 2 /g is preferred [29]. However, since the central SIDM density varies mildly with the cross section in range of 1-10 cm 2 /g [44,45], a feature that is well-captured in our analytical model [26], an even larger cross section may work as well. For high surface brightness galaxies, to which most of galaxies with high V f belong, the fits are insensitive to the cross section because of the degeneracy between σ/m and Υ [29]. The effect in the SIDM fits induced by varying σ/m can be compensated by a minor change in the stellar mass-to-light ratio, and many of these systems are actually compatible with an NFW profile. The cross section may have a mild velocity dependence over the sample, as implied by the constraint from galaxy clusters [26,46,47], but it is impossible to extract it from the SPARC dataset given the reasons discussed above. In this work, we present the results for fixed σ/m = 3 cm 2 /g and they remain the same qualitatively for other values larger than ∼ 1 cm 2 /g on galaxy scales.
An important consequence of the large cross section is that the SIDM profile is driven quickly to be isothermal in the inner regions. This implies that the resultant SIDM fits will not depend sensitively on the formation history of individual galaxies [29], but the final stellar and gas distributions [25]. This has been explicitly confirmed in recent hydrodynamical SIDM simulations [34] and those with idealized disk growth [28]. Furthermore, in our fits r 1 is close to r s , which is well outside the stellar disk or budge in the galaxies. It is unlikely that a viable baryonic feedback process could change the halo mass profile significantly at that far distance. Thus, our analytical model takes into account the realistic baryon distribution for individual galaxies and encodes this effect on the SIDM halo profile through the matching procedure.

III. THE RADIAL ACCELERATION RELATION IN SIDM
In the RAR described in Ref. [17], the gravitational acceleration g tot at radius r is found to be related to the acceleration g bar at the same radius. This relation can be fit to a functional form with a single parameter g † : Their best-fit value of g † = 1.2 × 10 −10 m/s 2 is the oft-quoted MOND acceleration scale.
In the left panel of Fig. 2, we show the inferred total and baryonic acceleration values from the controlled sampling, where g mod tot and g mod bar are calculated from the SIDM fits, using the halo parameters and the best-fit Υ values for each galaxy. The intensity of color in Fig. 2 (left) reflects the density of points. After fitting the data with the empirical relation given in Eq. 2, we find the best-fit value of g † is 1.38 × 10 −10 cm 2 /g and the resulting dispersion in the residuals is 0.10 dex. Fig. 2 (middle) shows Υ ,disk distribution from the SIDM fits (solid). It is peaked toward Υ ,disk = 0.5M /L , in good agreement with predictions from stellar population synthesis models [31]. This is remarkable because no priors based on the stellar population synthesis models were used. We have also reproduced the analysis in Ref. [17] with Υ ,disk and Υ ,bulge were fixed to 0.5M /L and 0.7M /L , respectively. For this fixed Υ case, we obtained g † = 1.19 × 10 −10 m/s 2 and dispersion 0.12 dex, both in agreement with previous work [17]. rmax-Vmax distributions of the host halos in the SIDM fits with controlled (circles) and MCMC (squares) samplings. We also show the mean relation (black solid) and 2σ scatter (gray shaded) predicted in cosmological CDM simulations [36]. Middle: Halo virial mass vs galaxy stellar mass from the SIDM fits. The black solid line corresponds to the abundance matching inference from [49]. Right: total baryonic mass vs flat circular velocity for the 135 galaxies, where M bar is inferred from our SIDM fits (circles and squares). For comparison, we also show the case (triangles) when Υ ,disk and Υ ,bulge are fixed to 0.5M /L [50]. The black solid line is the mean baryonic Tully-Fisher relation from [50], derived from 118 SPARC galaxies with Υ ,disk = Υ ,bulge = 0.5M /L , at which the scatter is minimized.
For a more detailed comparison, we also fit the sample of 135 SPARC galaxies using the MOND relation in Eq. 2, where we fixed g † = 1.2 × 10 −10 m/s 2 , but varied Υ ,disk and Υ ,bulge using MCMC sampling (see also [48]). The results look similar if we set g † to 1.0 × 10 −10 or 1.4 × 10 −10 m/s 2 . The middle panel of Fig. 2 shows the Υ ,disk distribution from the MOND fits (dotted), which closely tracks the one from the SIDM fits. The right panel shows the distribution of minimum χ 2 /d.o.f. values for individual galaxies from the SIDM and MOND fits. The SIDM model provides a better fit than MOND for most of the galaxies (∼ 77%), while maintaining a tight RAR. In fact, 72% (45%) of them have χ 2 /d.o.f. ≤ 3 (1) in the controlled SIDM fits and those with a large χ 2 /d.o.f. value have either tiny errors or wiggles in the observed rotation curves that cannot be reproduced by MOND either.
We emphasize that the diversity in the inner rotation curves is also reflected in the g tot -g bar plane, as explicitly demonstrated in Supplementary Materials, where we show the g tot vs g bar plot, but now split the sample into two sets: radii outside and inside 2R d with R d being the scale radius of the stellar disk. The scatter is relatively large for radii < 2R d , and this is due to the different shapes in the inner rotation curves and not just the result of random errors (see also [51]). On the other hand, there is a clear ordered behavior of g tot vs g bar curves for radii > 2R disk , which is a reflection of the BTFR: the tight correlation between the flat circular velocity, V f , and the total baryonic mass, M bar for spiral galaxies [52]. In this regime, g tot ≈ √ g † g bar , where g tot ≈ V 2 f /r and g bar ≈ GM bar /r 2 , hence we have V 4 f /(GM bar ) ≈ g † . This is the success of MOND, i.e., if one assumes M bar ∝ V 4 f , then the normalization of the BTFR also predicts the rotation curve, which in many cases is a good fit to the observed one. Many studies do find M bar ∝ V s f with 3 < s < 4 [50,53,54], as we will also show in Sec. V; s = 4 is not forced upon us by the data, but it is not ruled out either. However, the MOND relation (Eq. 2) cannot explain the full range of the diversity in the inner rotation curves, while the success of SIDM is deeply rooted to hierarchical structure formation, as we discuss in the next section.

IV. THE CONCENTRATION-MASS RELATION AND ORIGIN OF THE CHARACTERISTIC ACCELERATION SCALE
We have demonstrated that SIDM explains both the diversity and the tight RAR exhibited in the rotation curves, as dark matter self-interactions thermalize the inner halo in the presence of the baryonic potential. Here, we show the host halos in the SIDM fits are consistent with predictions in the hierarchical structure formation model, see, e.g., [36,55,56]. Since the outer halo (r r 1 ) remains unchanged for σ/m = 3 cm 2 /g, we parameterize an SIDM halo using the concentration and mass or, equivalently, the maximal circular velocity (V max ) and the associated radius (r max ) of its CDM counterpart. Ideally, one would measure these halo parameters directly from the kinematics data and compare them with simulations. Unfortunately, most rotation curves do not have the radial extent needed to sufficiently constrain them. In this work, we impose the cosmological concentration-mass relation [36] as a prior similar to Ref. [57] and examine the consistency between its consequences and observations. In Fig. 3 (left), we show the r max -V max distributions from our controlled (circles) and MCMC (squares) samplings. For the former, we intend to seek the best SIDM fits to the rotation curves following the mean relation (solid) from simulations. For the sample we consider, 97% galaxies can be fitted within the 2σ band (gray shaded), calculated from the relation log 10 c 200 = 0.905 − 0.101 log 10 (M 200 /10 12 h −1 M ) with an intrinsic scatter of 0.11 (1σ) [36]. For the latter, we impose the c 200 -M 200 relation as a top-hat prior within the 3σ range in our MCMC sampling, together with an additional constrain on V max , 1/ The resulting inferences (median and 1σ error) are shown in the figure. The two results agree well with each other. It is remarkable that even with the stringent constraints on V max and r max (through the c 200 -M 200 relation), the SIDM halo model is able to fit the diverse rotation curves, as illustrated in Figs. 1 and 2 (left). Indeed, with the concentration-mass relation, we find the Υ ,disk distribution is peaked toward 0.5M /L in the fits, shown in Fig. 2

(middle).
To see the MOND acceleration scale emerging from the hierarchical structure formation model, we parametrize a CDM halo with its gravitational acceleration at r = 0 as g NFW (0) = GM/r 2 | r→0 ≈ 2πGρ s r s ≈ 2πV 2 max /(1.26r max ). Taking the mean cosmological V max -r max relation, r max = 27 kpc(V max /100 km/s) 1.4 , we have g NFW (0) ≈ 1.0 × 10 −10 m/s 2 (V max /240 km/s) 0.6 , which is close to the MOND acceleration parameter g † . This is the underlying reason why the empirical MOND relation captures the overall stellar kinematics of spiral galaxies well. In the presence of dark matter self-interactions and baryons, the actual central acceleration deviates from g NFW (0), but the general argument still holds. For example, we can characterize a halo with the acceleration at the scale radius r s , where the impact of dark matter self-interactions and influence of baryons tend to be small, g NFW (r s ) ≈ 0.39g NFW (0), slightly smaller than g NFW at the center. The characteristic halo acceleration has a mild dependence on V max , ranging from 20 to 300 km/s in the sample, and it also varies with the scatter in the cosmological relation. This variation is important, as shown in Fig. 3 (left). Since MOND does not have such flexibility, its overall fits are worse than the SIDM ones, as illustrated in Fig. 2 (right). We emphasize that g † = 1.38 × 10 −10 m/s 2 inferred from our SIDM fits in Sec. III is an average quantity over the sample after fitting to Eq. 2, not a universal value for all the galaxies as in MOND.
The calculation of the acceleration due to dark matter toward the center is more subtle. Inside a constant density core g SIDM (r) ∝ r, and we need to specify the radius where the acceleration is being computed. On average, the stellar halflight radius is empirically observed to track the virial radius as r 1/2 ≈ 0.015r vir [58], and we have r 1/2 ≈ 1.7R d for an exponential disk model. Without a significant contribution from baryons to the gravitational potential, SIDM predicts that g SIDM (r 1/2 ) = 10 −11 m/s 2 (V max /100 km/s) 0.2 for the median halo concentration, and its dependence on the halo mass is extremely mild. When baryons contribute, g tot does not increase linearly with g bar since both the central SIDM density and the core radius depend on the gravitational potential contributed by the baryons. The net result is a strong correlation between g tot and g bar , which is clearly evident in Fig. 2. The model predictions have a definite width in the g tot vs g bar plane and we have shown clearly that this scatter is required to fully explain the diversity in the rotation curve data.

V. THE CORRELATIONS BETWEEN THE TOTAL LUMINOUS AND DARK MATTER MASSES
We have seen the SIDM fits to the rotation curves require values for the halo concentration parameter that are completely in line with the expectations from the Planck experiment [59,60]. In addition, the stellar mass-to-light ratios are consistent with the results from stellar population models [31].
This leads to a natural question: what is the predicted halo mass for a given stellar mass in the SIDM model? Since we assume the primordial matter power spectrum is unchanged from the CDM one for the scales we are interested in, there should be a relation consistent with the abundance matching results in the literature. In the middle panel of Fig. 3, we show the stellar mass vs halo mass relation derived using the mass-to-light ratios from controlled (circles) and MCMC (squares) samplings. The error bars on the MCMC points denote the 1σ widths from the posteriors (16 th and 84 th percentiles). Our results are consistent with the overall trend in the relation from abundance matching (solid) [49] (see [51] for the CDM case).
We have already alluded to the importance of the BTFR in our discussion of the RAR. Lelli et al. [50] selected 118 SPARC galaxies and found that their V f -M bar inferences can be fitted with a simple relation: log(M bar ) = s log(V f ) + log(A), where s = 3.71 ± 0.08 and log(A) = 2.27 ± 0.18 for Υ ,disk = Υ ,bulge = 0.5M /L . The right panel of Fig. 3 shows the V f -M bar inferences with the Υ ,disk and Υ ,bulge values from the controlled (circles) and MCMC (squares) fits. The error bars in M bar on the MCMC points denote the 1σ widths in the stellar mass-to-light ratios from the posteriors, and the errors in V f are taken directly from the SPARC dataset [18]. We also show the fit from [50] as the solid line of Fig. 3 (right). Note that this fit used 118 galaxies and a few outliers at the low V f end were not included. For comparison, we plot the 135 galaxies in our sample as triangles by fixing Υ ,disk = Υ ,bulge = 0.5M /L . We see that their distribution in the V f -M bar plane is almost identical to the one from our SIDM fits. This is not surprising, as the Υ ,disk values inferred from the SIDM fits are peaked toward 0.5 M /L as shown in Fig. 2 (middle). Thus, we conclude that the SIDM fits also lead to a tight BTFR relation. For our fits, we find s ≈ 3.46 (CS), 3.27 (MS) and 3.58 (0.5M /L ), excluding six obvious outliers on the left side of the black line.
We note that there is no evidence in the data for s = 4, i.e., M bar ∝ V 4 f , the motivation for MOND, either in the constant Υ fits or in the SIDM fits. Many of the recent CDM simulations with efficient baryonic feedback seem to get something akin to the BTFR with s ≈ 3.6-3.8 [61][62][63][64], but it is fair to say that this is still not well understood theoretically, in particular, the smallness of the scatter in the BTFR, equivalent to the one seen in the RAR [65]. We expect that there will be interplay between dark matter self-interactions and baryonic feedback in changing the halo potential, and understanding how the BTFR emerges in SIDM is fertile territory for research in galaxy formation.

VI. DISCUSSION AND CONCLUSIONS
In this work, we have investigated SIDM as a solution to two puzzles that are present in galactic rotation curves: (1) the diversity of inner rotation curves in galaxies that have similar baryon content and similar flat circular velocities, and (2) the small scatter in the radial acceleration relation between the total gravitational acceleration and the one inferred from the baryonic mass content, i.e., uniformity.
We have fitted our SIDM halo model to the rotation curves of 135 SPARC galaxies, and found that it reproduces the observed diversity in the inner regions. The distribution of resulting 3.6 µm stellar disk mass-to-light ratios for the sample peaks at Υ ,disk ≈ 0.5 M /L , in good agreement with the stellar population models. Our fits lead to a radial acceleration relation described by the characteristic acceleration scale ∼ 10 −10 m/s, with tight scatter of 0.10 dex. The host halos are fully consistent with the Planck cosmology. The inferred stellar mass-halo mass relation agrees with the result from the abundance matching method, and the fits also predict a tight BTFR. These results provide compelling arguments in favor of the idea that the inner halos of galaxies are kinematically thermalized due to dark matter self-interactions.
The SIDM model automatically inherits all of the successes of the CDM model on large scales, as the predictions are indistinguishable at distances larger than about 10% of the virial radius of galactic halos. The required cross section is similar to the proton-neutron elastic scattering cross section and this may be a strong hint that the dark matter sector replicates some elements of the standard model. The large cross section keeps the inner halo isothermal and this makes the predictions for the central halo profile at later times insensitive to the star formation history, as confirmed in recent hydrodynamical N-body simulations [34,66]. This implies that a large variety of feedback models, e.g, [67][68][69][70][71][72], can be compatible with the SIDM model we have discussed here. The predictions are quantitatively the same for σ/m 1 cm 2 /g. This makes our results robust, but it makes hard to precisely determine the cross section from kinematic datasets on galaxy scales [29].
There are a number of promising directions that can further test SIDM and explore galaxy formation and evolution in this framework. Here, we highlight a few of them. SIDM simulations predict a correlation between the half-light radius of the stars and the dark matter core size in dwarf and low surface brightness galaxies [27], which should be further explored and may provide an observational test of SIDM. Similarly, the ultra-diffuse galaxies in the clusters could be a test laboratory [73]. A related issue is the origin of the large spread in the surface brightness of galaxies, which remains poorly understood. Interestingly, hydrodynamical simulations of galaxy clusters show that the stellar density profiles in SIDM are more diverse than in their CDM counterparts [34]. Is this a more general feature in SIDM due to the dynamical interplay between core formation and feedback? How does this interplay impact the emergence of the BTFR? Finally, at the lowest mass end, the dwarf spheroidal galaxies, including the so-called ultra-faint dwarfs, in the Local Group could provide a key test of SIDM (see [74,75]). Dedicated SIDM simulations with the baryons will be required to explore these exciting topics.
The predictive power of the SIDM model, the clear connection to cosmology, and its rich implications for other astrophysical observations and particle physics phenomenology [76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94][95], all taken together make a clear case that it should be treated on the same footing as the CDM model. The economical explanation, with the addition of just one parameter, for the diverse rotation curves across the entire range of observed galaxies argues in favor of the idea that the dark matter particles have a large affinity for the self-interactions.

Acknowledgments
We thank the authors of [30,39,44] for the simulation data, which are used for comparison with the analytical model, shown in Supplementary Materials. This work was supported by the National Science Foundation Grant No. PHY-1620638 (MK), the U. S. Department of Energy under Grant No. de-sc0008541 (HBY), the Hellman Fellows Fund (HBY), UCR Academic Senate Regents Faculty Development Award (HBY), and the National Science Foundation under Grant No. NSF PHY-1748958 as part of the KITP "The Small-Scale Structure of Cold (?) Dark Matter" and "High Energy Physics at the Sensitivity Frontier" workshops (MK and HBY). HBY also acknowledges Tsung-Dao Lee Institute, Shanghai, for its hospitality during the completion of this work.

Methods
We provide a detailed description of the analytical model developed previously [26,29] and the fitting procedure in this section. We divide the halo into an inner and an outer region [39] with the aim that the outer halo is not significantly changed by the self-scattering process. In the inner region, dark matter self-interactions thermalize the halo in the presence of the baryonic potential, and we model the dark matter distribution using the isothermal density profile, ρ iso ∝ exp(−Φ tot (R, z)/σ 2 v0 ). Poisson's equation relates Φ tot to the dark matter and baryon profiles as For the outer halo, where the self-scattering effect becomes negligible, we model the dark matter distribution with an NFW profile ρ NFW (r) = ρ s r 3 s /r(r + r s ) 2 . To construct the full SIDM halo profile, we define a radius r 1 , where dark matter particles had one interaction on average over the age of the galaxy. We join the spherically-averaged isothermal (ρ iso ) and spherical NFW (ρ NFW ) profiles at r = r 1 such that the mass and density are continuous at r 1 . Thus, the isothermal parameters (ρ 0 , σ v0 ) directly map on to the NFW parameters (r s , ρ s ) or (r max , V max ).
The value of r 1 is determined by the following condition, where σ is the self-scattering cross section, m is the dark matter particle mass, v rel is the dark matter relative velocity in the halo, ... denotes averaging over the Maxwellian velocity distribution, t age is the age of the galaxy, and N sc is a factor of order unity, to be determined by calibrating to simulations. In this work, we have set t age = 10 Gyr and N sc = 1, which reproduce simulation results well; see Supplementary Materials. In principle, we should use different ages for each galaxy, say between 10 Gyr and 13 Gyr. However, our model can only constrain the combination of the cross section and the age. More importantly, we have set σ/m to a large enough value that the SIDM density profile is insensitive to small changes in the cross section. We assume that this cross section is a constant over the SPARC sample, so σv rel = σ(4/ √ π)σ v0 . We take two independent but complementary approaches. In the first one, we assume a thin-disk profile for the stellar disk in solving Eq.
where Σ 0 is the central surface density and R d is the scale radius. For each galaxy, we reconstruct the Σ 0 and R d values by fitting the profile to the disk contribution of the rotation curve as in [29]. We neglect the baryonic influence on the SIDM halo from the gas and bulge potentials, but include all the mass components in modeling the total circular velocity. This is a reasonable approximation for the following reasons: (1) the gas is less centrally concentrated and so its impact on the SIDM density profile is smaller, (2) the bulge (when present) mainly affects the innermost region, while the disk contributes in this region as well as at farther radii. Ref. [29] solved Eq. (3) with the thin-disk approximation and created numerical templates for the isothermal density profile on the grid of a ≡ 8πGρ . When the stellar profile is known, the parameters a and b give the central density and dispersion of the isothermal dark matter halo, which completely specify the inner density profile. We interpolate the templates to generate rotation curves for any set of (ρ 0 , σ v0 , Σ 0 , R d ). The fixed value of the cross section allows us to match this density profile to the outer NFW density profile.
In fitting to the SPARC sample with the templates, we take a controlled sampling approach. For a given galaxy, we start with the mean r max -V max relation from cosmological ΛCDM simulations [36] and an NFW profile that matches the flat part of the rotation curve. Then, we choose an appropriate Υ ,disk (Υ ,bulge ) value to reproduce the inner rotation curve. We calculate a χ 2 /d.o.f. value for each fit and iterate this process manually by adjusting the parameters until a good fit is achieved. For most galaxies, the very first step provides decent fits, showcasing the simplicity of the model and its ability to fit the observed data simultaneously. For each galaxy, we demand the (r max , V max ) values to be within the ∼ 2σ band. In this way, we have good control over the halo parameters in the fits. The goal is to see to what degree are the galaxy halos of the SPARC sample consistent with predictions of the hierarchical structure formation scenario, and the extracted Υ values consistent with stellar population synthesis model results [31].
In our second approach, we perform a MCMC sampling of the SIDM model parameter space. Since it is computationally expensive to use the templates, we use a spherical approximation to model the baryon distribution [25,28]. We create a spherical baryonic mass profile from the stellar and gas masses, such that the baryonic mass within a sphere of radius r is M b (r) = (V 2 disk + V 2 bulge + V 2 gas )r/G, where V disk is the contribution to the rotation curve from the disk at radius r and similarly for the bulge and gas. We solve Eq. (3) in the spherical limit by taking r = √ R 2 + z 2 . We compute ρ iso (r) starting at a small radius assuming a core and integrate the equation to larger radii. The behavior of the baryonic density profile as r → 0 is chosen so that a core at small radii is physical [25]. We compared the isothermal halos from this spherical approximation to those from the axisymmetric case (templates) and found agreement within 10-20%. Thus, while we expect some variance in the inferred parameters between the two methods, the overall features should be very similar. This expectation is borne out by our final fits.
We match the isothermal density profile ρ iso , parameterized by (ρ 0 , σ v0 ), to the NFW density profile at r 1 , and this determines (V max , r max ). Thus, the spherical model has four parameters, two for the entire halo and two for the mass-to-light ratios: (ρ 0 , σ v0 , Υ ,disk , Υ ,bulge ). We use the emcee implementation of the Affine invariant MCMC ensemble sampler [96] to infer the posteriors of these four model parameters. To streamline the calculation of r 1 at each point in parameter space for matching onto the outer NFW radius, we use the rate of scatterings, Γ 0 = ρ 0 (σ/m)(4/ √ π)σ v0 , within the isothermal core as the MCMC parameter in lieu of the core density ρ 0 .
The prior distributions used for the halo parameters and the mass-to-light ratios in the MCMC scan are as follows: • Γ 0 : Uniform prior in the range of log 10 2 < log 10 (Γ 0 × 10 Gyr) < 5.
• Υ : Uniform prior in 0.1 < Υ ,disk , Υ ,bulge < 10 M /L . The parameter Υ ,bulge is only included for galaxies whose surface brightness profiles have a stellar bulge decomposition provided in the SPARC dataset. All galaxies have Υ ,disk as a parameter describing their stellar disk. Additionally, we also impose two regularization priors.
• We add 5% of V f in quadrature for calculating the likelihood function. This allows the code to disregard the points deep within the central regions and those with tiny errors. We have checked that it doesn't change the inference of cores/cusps. We do not include this regularization error when quoting χ 2 values.
• We impose a uniform regulation prior on V max : 1/ For most of the galaxies (∼ 80%), our MCMC program can find physical fits without this prior. However, in some cases, the MCMC sampler tends to pick up fits not consistent with hierarchical structure formation predictions -either the baryonic contribution dominates the total rotation velocity at all radii and the halo concentration is unreasonable low (for some high surface brightness galaxies), or the opposite (for some low surface brightness galaxies). This is typically due to the lack of an extended rotation curve to fully constrain the halo parameters. The additional regularization prior fixes this issue. We have also checked that the results are similar if we consider a more generous range 1/2 < V max /V f < 2 (see Supplementary Materials).

Supplementary Materials
We provide additional information and results to supplement the results in the main text.
• In Table S1, we list the galaxies that are shown in Fig. 1 of the main text.
• We show the total acceleration vs the baryonic acceleration for the inner and outer regions in Fig. S1.
• Fig. S2 shows r max -V max , M star -M halo , and M bar -V flat relations, similar to Fig. 3 of the main text, but we impose the top-hat prior on the concentration-mass relation with a wider V max regulation, 1/2 < V max /V f < 2. In addition, we show the results with a Gaussian prior on the concentration-mass relation and 1/ • Fig. S3 shows • Fig. S5 shows detailed SIDM and MOND fits to individual 135 SPARC galaxies.
•  FIG. S1: Upper: the total acceleration vs the baryonic acceleration (colored) for the inner (r ≤ 2R d , left) and outer (r > 2R d , right) regions, where R d is the scale radius of the stellar disk. Lower: The gtot-g bar relation with a different color scheme, where the intensity is proportional to the density of points. The scatter in the gtot-g bar relation of the inner regions is visibly larger (black solid).   Despite the fact that we impose the exact matching condition at r1, i.e., ρiso = ρNFW and Miso = MNFW, and the agreement is better than ∼ 5-20% for σ/m ≥ 1 cm 2 /g and the results change very mildly from tage = 10 Gyr to 13 Gyr. Sokolenko et al., 1806.11539, also showed the core sizes predicted in this analytical model are consistent with their simulations. The agreement can be further improved through tweaks to this model by including small halo mass or cross section dependence in the r1 definition or allowing freedom in the matching at the level of ∼ 5%. In the paper, we take σ/m = 3 cm 2 /g, tage = 10 Gyr and the exact matching condition. ■ CDM ■ 0.5 cm 2 /g ■ 1 cm 2 /g ■ 5 cm 2 /g ■ 10 cm 2 /g ■ CDM ■ 0.5 cm 2 /g ■ 1 cm 2 /g ■ 2 cm 2 /g ■ 5 cm 2 /g ■ 10 cm 2 /g    Fig. S3, with controlled sampling (first row associated with each galaxy) and MCMC sampling with the tophat prior (best-fit value, second row; medium with 1σ errors, third row). The α value is the logarithmic slope of the dark matter density profile at r = 1.5%rvir. Galaxies are listed alphabetically, corresponding to the order in Fig. S3.