A Self-interacting Dark Matter Solution to the Extreme Diversity of Low-mass Halo Properties

The dark matter (DM) distribution on galactic scales and below is a powerful probe of DM particle properties (e.g., Bechtol et al. 2022). One key measure of this distribution is its diversity, i.e., the range of DM halo properties at fixed mass. Studies have long shown that the inner circular velocities of spiral and dwarf irregular galaxies in the field span a wide range—in excess of cold, collisionless DM (CDM)-only predictions—even though their total halo masses are similar (Kuzio de Naray et al. 2010; Oman et al. 2015; Ren et al. 2019; Santos-Santos et al. 2020). Although baryonic feedback within galaxies likely explains some of this diversity (e.g., see Sales et al. 2022 for a review), more recent observations challenge this as a complete explanation.

In particular, Minor et al. (2021, hereafter M21) showed that the strong lens SDSSJ0946+1006 is perturbed by a DM (sub)structure (Vegetti et al. 2010) with an extremely high density, ≳3σ (≳5σ) away from the median CDM expectation when interpreted as a subhalo (line-of-sight halo), and Ballard et al. (2023) reported a comparable (though slightly weaker) tension. Meanwhile, a class of large, gas-rich isolated galaxies has recently been discovered (Leisman et al. 2017; Mancera Piña et al. 2019, 2020, 2022; Guo et al. 2020; Shi et al. 2021). These ultradiffuse galaxies (UDGs) contain less DM than expected based on their baryonic content; Kong et al. (2022, hereafter K22) showed that these UDGs must inhabit extremely low-concentration halos, ≳3σ to 5σ away from the median CDM expectation (also see Mancera Piña et al. 2019, 2022). However, K22 find that CDM halos from the IllustrisTNG hydrodynamic simulation (Nelson et al. 2019), even with such low concentrations, are still much denser than these UDGs.

Baryonic feedback can steepen (sub)halo density profiles, via adiabatic contraction (Blumenthal et al. 1986; Gnedin et al. 2004), or core them via repeated cycles of supernova feedback (Navarro et al. 1996; Governato et al. 2010; Pontzen & Governato 2012). However, it is challenging for baryons to simultaneously explain the diverse observations described above. In particular, the SDSSJ0946+1006 perturber does not have a detectable luminous component; thus, its stellar mass is negligible and that adiabatic contraction is unlikely to explain its extremely steep density profile (Vegetti et al. 2010; M21). If baryonic feedback forms a core in the progenitor halo, the tension would be exacerbated.

At the opposite extreme, gas-rich UDGs have very shallow potentials; thus, it takes longer to replenish the gas necessary for bursty star formation. As a result, supernova feedback is expected to be inefficient in UDGs (Mancera Piña et al. 2022; K22), although we note that some hydrodynamic simulations with episodic feedback produce diffuse galaxies in median-concentration halos (e.g., Di Cintio et al. 2017; Chan et al. 2018). In addition, environmental effects such as tidal stripping—while important for UDGs in groups and clusters (Ogiya 2018; Carleton et al. 2019; Jiang et al. 2019; Sales et al. 2020; Tremmel et al. 2020; Yang et al. 2020; Benavides et al. 2021; Moreno et al. 2022)—are largely absent for the UDGs analyzed in K22.

We will show that both tensions may be resolved if DM has strong self-interactions that thermalize halos' inner regions (see Tulin & Yu 2018; Adhikari et al. 2022 for reviews). The key physics that enables this is the gravothermal evolution of self-interacting DM (SIDM) halos, which have two sequential phases: core expansion and collapse. In the first phase, interactions transport heat from outer to inner halo regions, lowering the central density; in the second, the direction of heat flow reverses and the inner halo becomes denser than its CDM counterpart (Balberg et al. 2002; Koda & Shapiro 2011). The timescale for the onset of core collapse scales as ${t}_{c}\propto {({\sigma }_{\mathrm{eff}}/m)}^{-1}{c}_{200}^{-7/2}{M}_{200}^{-1/3}$ (Essig et al. 2019; Kaplinghat et al. 2019), where σ_eff/m is an effective cross section (Yang & Yu 2022; Outmezguine et al. 2023; Yang et al. 2023c) and M₂₀₀ (c₂₀₀) are mass (concentration), respectively. If ${t}_{c}\sim { \mathcal O }(10)\,\mathrm{Gyr}$, the diversity of inner halo densities encoded by concentration (Kaplinghat et al. 2016; Kamada et al. 2017) is amplified at z ∼ 0 (Yang et al. 2023a).

In this Letter, we introduce a new velocity-dependent SIDM model, perform the first high-resolution cosmological N-body simulation with strong DM self-interactions on galaxy group scales, and identify simulated halos and subhalos that can self-consistently explain the observed lensing perturber and UDGs.

We consider velocity-dependent DM self-interactions with a Rutherford-like differential cross section (Feng et al. 2010; Ibe & Yu 2010; Tulin et al. 2013; Yang & Yu 2022)

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{d\cos \theta }=\displaystyle \frac{{\sigma }_{0}{w}^{4}}{2{\left[{w}^{2}+{v}^{2}{\sin }^{2}(\theta /2)\right]}^{2}},\end{eqnarray} \tag{ 1 }$$

where θ is the scattering angle in the center of mass frame, v is the relative velocity of DM particles, the cross section drops as v⁻⁴ for v > w, and σ₀ controls the normalization. We set the cross section per mass of DM particles σ₀/m = 147.1 cm² g⁻¹, w = 120 km s⁻¹ and refer to this as the "Group SIDM" model. A cross section with the same amplitude but a lower turnover scale of w = 24.33 km s⁻¹ was simulated in a Milky Way (MW) setting (Yang et al. 2023a; also see Zavala et al. 2019; Correa 2021; Turner et al. 2021; Correa et al. 2022).

In general, DM scattering is both velocity- and angle-dependent, as indicated in Equation (1). Yang & Yu (2022) demonstrated that the angular dependence can be absorbed by introducing a viscosity cross section, which weights the differential cross section with a factor of ${\sin }^{2}\theta$. Figure 1 (inset) shows the viscosity cross section σ_V for the SIDM model we consider; see Appendix A for details. Since σ_V does not explicitly depend on the scattering angle, it greatly simplifies the implementation of DM self-interactions in N-body simulations. Our SIDM simulation will be based on the velocity-dependent viscosity cross section.

Zoom In Zoom Out Reset image size

For a given halo, its gravothermal evolution can be captured by the constant effective cross section (Yang & Yu 2022)

$$\begin{eqnarray}&&{\sigma }_{\mathrm{eff}}=\displaystyle \frac{2\int {v}^{2}{dvd}\cos \theta \tfrac{d\sigma }{d\cos \theta }{\sin }^{2}\theta {v}^{5}\exp \left[-\tfrac{{v}^{2}}{4{\nu }_{\mathrm{eff}}^{2}}\right]}{\int {v}^{2}{dvd}\cos \theta {\sin }^{2}\theta {v}^{5}\exp \left[-\tfrac{{v}^{2}}{4{\nu }_{\mathrm{eff}}^{2}}\right]},\end{eqnarray} \tag{ 2 }$$

where ${\nu }_{\mathrm{eff}}\approx 0.64{V}_{\max }$ is a characteristic velocity dispersion for a CDM halo with maximum circular velocity ${V}_{\max }$. Figure 1 shows σ_eff/m versus ${V}_{\max }$ for our Group SIDM model (red), the model in Yang et al. (2023a, purple), which is similar to that in Turner et al. (2021), and the best-fit model from Kaplinghat et al. (2016, black) for comparison. The effective cross section of our SIDM model can be fitted with a functional form of $147.1\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}/{[1+{({V}_{\max }/80\,\mathrm{km}\,{{\rm{s}}}^{-1})}^{1.72}]}^{1.90}$. For a Navarro–Frenk–White (NFW) halo (Navarro et al. 1996), ${V}_{\max }\,=1.64{r}_{s}\sqrt{G{\rho }_{s}}$, where ρ_s is its scale density and r_s its scale radius, and G is Newton's constant. In this case, the effective cross section is ultimately related to the halo's mass M₂₀₀ and concentration c₂₀₀. For a halo with a given ${V}_{\max }$, the effective cross section provides a good proxy to evaluate the impact of the self-interactions on halo properties, although the scattering probability of individual DM particles is velocity- and angle-dependent; see Yang & Yu (2022) for further details.

Our Group SIDM model is motivated as follows. First, a cross section significantly larger than ∼1 cm² g⁻¹ is needed to produce a shallow density profile in low-concentration halos (e.g., see K22). For a typical UDG halo with ${V}_{\max }=40\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and c₂₀₀ 5σ below the CDM median, our Group SIDM model predicts σ_eff/m ≈ 86 km s⁻¹, and the collision rate at the scale radius is 1.3 per 10 Gyr, high enough to produce a density core. Second, for (sub)halos with ${V}_{\max }\approx 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$, relevant for the lensing perturber, σ_eff/m ≈ 25 cm² g⁻¹ and t_c ≲ 10 Gyr if c₂₀₀ is 1.5σ higher than the CDM median, where we have used the following formula to estimate the collapse time (Essig et al. 2019)

$$\begin{eqnarray}&&{t}_{c}=\displaystyle \frac{150}{C}\displaystyle \frac{1}{{r}_{s}{\rho }_{s}({\sigma }_{\mathrm{eff}}/m)}\displaystyle \frac{1}{\sqrt{4\pi G{\rho }_{s}}},\end{eqnarray} \tag{ 3 }$$

with C = 0.75 from calibrating against N-body simulations. Further using the relations

$$\begin{eqnarray}&&{r}_{s}={\left[\displaystyle \frac{3{M}_{200}}{(200{\rho }_{\mathrm{crit}})4\pi {c}_{200}^{3}}\right]}^{\tfrac{1}{3}},{\rho }_{s}=\displaystyle \frac{200{\rho }_{\mathrm{crit}}{c}_{200}^{3}}{3f({c}_{200})},\end{eqnarray} \tag{ 4 }$$

where ρ_crit is the critical density and $f({c}_{200})\equiv \mathrm{ln}({c}_{200}\,+1)-{c}_{200}/({c}_{200}+1)$, we have ${t}_{c}\propto {({\sigma }_{\mathrm{eff}}/m)}^{-1}{c}_{200}^{-7/2}{M}_{200}^{-1/3}$ as f has a very mild dependence on c₂₀₀ (Essig et al. 2019). Thus, strong DM self-interactions can further diversify inner halo densities encoded by concentration; for subhalos, tidal stripping can accelerate core collapse (Kahlhoefer et al. 2019; Nishikawa et al. 2020; Sameie et al. 2020; Correa 2021; Zeng et al. 2022).

Finally, our model evades constraints of σ_eff/m ≲ 1 (Sagunski et al. 2021) and 0.1 cm² g⁻¹ (Kaplinghat et al. 2016; Andrade et al. 2022) from groups and clusters (${V}_{\max }\approx 800$ and 1000 km s⁻¹, respectively). Cross section constraints from cluster mergers are typically weaker (Harvey et al. 2015; Jauzac et al. 2016; Wittman et al. 2018).

We also comment on constraints from MW satellites and spiral galaxies in the field. Previous studies have shown that MW satellites' DM distributions are diverse, and that it is difficult to accommodate the full diversity of these systems in SIDM if the subhalos that host these galaxies, with $10\,\mathrm{km}\,{{\rm{s}}}^{-1}\lesssim {V}_{\max }\lesssim 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, are in the core-expansion phase for ${\sigma }_{\mathrm{eff}}/m\lesssim { \mathcal O }(10)\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ (Read et al. 2018; Valli & Yu 2018; Kaplinghat et al. 2019; Silverman et al. 2023). This has motivated the exploration of SIDM scenarios with stronger DM self-interactions (Zavala et al. 2019; Correa 2021; Turner et al. 2021; Yang et al. 2023a), such that some subhalos could be in the core collapse phase, resulting in a high density as observed in, e.g., Draco (Nishikawa et al. 2020; Sameie et al. 2020). For our SIDM model, σ_eff/m ∼ 140 cm² g⁻¹ for the subhalos of the MW, such that many of them core collapse within the Hubble time.

For isolated spiral galaxies, SIDM models with σ_eff/m ∼ 3 cm² g⁻¹ can largely explain the diversity of galactic rotation curves over a wide mass range of $30\,\mathrm{km}\,{{\rm{s}}}^{-1}\,\lesssim {V}_{\max }\lesssim 300\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Kamada et al. 2017; Ren et al. 2019). This cross section amplitude should be regarded as a lower limit for most of these systems, as the fit will not change significantly when the cross section increases further (Kaplinghat et al. 2020). Our SIDM model has 2 cm² g⁻¹ ≲ σ_eff/m ≲ 100 cm² g⁻¹ over this ${V}_{\max }$ range, so we expect that the success of previous SIDM fits will largely remain, though a detailed reanalysis is needed. We will discuss the implications of our SIDM model for MW satellites and spiral galaxies further after presenting our results.

We present a pair of CDM and SIDM cosmological DM-only zoom-in simulations centered on a host halo with M₂₀₀ = 1.1 × 10¹³ M_⊙ using high-resolution initial conditions from the "Group" suite of Symphony simulations (Nadler et al. 2023); throughout, we measure masses that enclose 200 times the critical density of the universe at z = 0. These simulations are run with a high-resolution particle mass of 3 × 10⁵ h⁻¹ M_⊙ and a Plummer-equivalent gravitational softening length of 170 h⁻¹ pc using GADGET-2 (Springel 2005). Halos are identified with ROCKSTAR (Behroozi et al. 2013a) and merger trees are built using consistent-trees (Behroozi et al. 2013b). The cosmological parameters are Ω_M = 0.286, Ω_Λ = 0.714, n_s = 0.96, h = 0.7, and σ₈ = 0.82 (Hinshaw et al. 2013). DM self-interactions are implemented using the code developed and validated in Yang & Yu (2022), based on a velocity-dependent viscosity cross section; see Figure 1 (inset) and Appendix A for details.

We study (sub)halos within 6 Mpc of the host halo center or ≈10 times the virial radius of the group-mass host. High-resolution particles comprise >90% of the mass in this (6 Mpc)³ region, which is the effective volume of our zoom-in simulation. Furthermore, we only analyze (sub)halos that contain >2000 particles at z = 0, corresponding to a present-day mass threshold of 8.6 × 10⁸ M_⊙; there are ≈10³ such isolated halos throughout the high-resolution region, several hundred of which are in the mass range relevant for the UDGs studied in K22. Several large isolated halos are resolved in addition to the group-mass host, including a ∼10¹² M_⊙ (MW-mass) halo and eight ∼10¹¹ M_⊙ (Large Magellanic Cloud-mass) halos. Here, we focus on isolated halos and subhalos of the group-mass host, where subhalos are defined as halos whose centers are within the virial radius of the group host. Convergence tests and subhalo population analyses will be presented in a companion study (E. O. Nadler et al. 2024, in preparation).

At z = 0, the group-mass host in our SIDM simulation has a ≈25 kpc core, which is consistent with the core size expected semianalytically given our SIDM cross section (Kaplinghat et al. 2014, 2016; Robertson et al. 2021). Isolated halo mass functions are nearly identical in our CDM and SIDM simulations. However, as we will demonstrate, our Group SIDM model significantly diversifies isolated halos' and subhalos' inner density profiles, simultaneously yielding systems that resemble low-concentration UDG halos and the dense strong-lens perturber.

The left panel of Figure 2 shows the inner projected mass and density profile slope for all subhalos of our group-mass host with ${V}_{\max }\gt 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, in our CDM (black) and SIDM (blue) simulations. We measure projected logarithmic density slopes averaged between 0.75 and 1.25 kpc following M21; within this range, our results are insensitive to the exact radii chosen for this averaging procedure. For comparison, we show the parameter space inferred by M21 for the strong-lens perturber, assuming a truncated-NFW ("tNFW") profile for the perturber and including higher-order multipoles in the lens model ("tNFWmult").

Zoom In Zoom Out Reset image size

Our CDM simulation yields subhalos with inner projected density profile slopes clustered around −0.75, consistent with the findings in M21 based on IllustrisTNG (Nelson et al. 2019). Meanwhile, most surviving subhalos in our SIDM simulation have steeper inner density profiles than their CDM counterparts. Thus, the tension between the steepness of the strong-lens perturber's inner density profile and CDM predictions is significantly eased in our SIDM model. Only one of our SIDM subhalos approaches the required inner mass of 2 × 10⁹ M_⊙; however, our group host is on the low end of the halo mass distribution for SDSSJ0946+1006, between 1 and 6 × 10¹³ M_⊙ (Auger et al. 2010; M21). A modest increase in host halo mass is therefore likely to produce even more subhalos with enclosed masses similar to the SDSSJ0946+1006 perturber in our SIDM model. We highlight two SIDM subhalos that match the properties of the strong-lens perturber fairly well; both of these systems have high concentrations and early formation times relative to the rest of the SIDM subhalo population, and one of them undergoes a pericentric passage ≈4 Gyr ago, which may accelerate its gravothermal evolution (Nishikawa et al. 2020; Sameie et al. 2020).

We also find a population of isolated SIDM halos with properties similar to the SDSSJ0946+1006 perturber, shown as unfilled points in the left panel of Figure 2. M21 found that essentially no isolated halos in IllustrisTNG have sufficiently steep inner density slopes to explain the perturber's properties. Thus, our SIDM model also eases the ∼5σ tension between CDM predictions and the SDSSJ0946+1006 perturber when interpreted as a line-of-sight halo. This is intriguing as Şengül et al. (2022) and Şengül & Dvorkin (2022) reported that another lens, JVAS B1938+666 (King et al. 1998; Vegetti et al. 2012), is perturbed by a dense line-of-sight structure.

The right panel of Figure 2 compares subhalo density profiles from our CDM and SIDM simulations to that inferred for the observed perturber assuming a tNFW (magenta) or tNFWmult (cyan) profile. At the 1 kpc radius where the perturber profile is best constrained (corresponding to r/R_vir = 0.025 in Figure 2), our SIDM predictions span the inferred density profiles, which are most robustly constrained within ≈1 kpc, indicating a remarkable level of agreement that is difficult to achieve in CDM.

Most surviving subhalos in our SIDM simulation have higher-amplitude and steeper density profiles than their CDM counterparts in the same ${V}_{\max }$ range. By studying the evolution of matched CDM and SIDM subhalos, we find that initially high-mass subhalos with large cores at infall are quickly disrupted in our SIDM simulation, mainly due to tidal stripping (Vogelsberger et al. 2012; Dooley et al. 2016; Nadler et al. 2020; Yang & Yu 2021; Zeng et al. 2022; Slone et al. 2023). Meanwhile, subhalos with relatively low masses and high concentrations at infall often show signs of core collapse after their first pericenter, consistent with findings that tidal stripping accelerates subhalos' gravothermal evolution (Kahlhoefer et al. 2019; Nishikawa et al. 2020; Sameie et al. 2020; Correa 2021; Zeng et al. 2022). These core-collapsed subhalos are more resilient to disruption than cored systems; thus, most surviving subhalos in our Group SIDM model have steep inner density profiles.

The left panel of Figure 3 shows halo circular velocity, evaluated at a radius of 2 kpc, versus maximum circular velocity for the seven UDGs from K22 based on the kinematic data presented in Mancera Piña et al. (2019, 2020, 2022) and Shi et al. (2021) for AGC 114905 (red), 122966 (green), 219523 (blue), 248945 (orange), 334315 (magenta), 749290 (brown), and 242019 (black). The same quantities for isolated halos in our CDM (SIDM) simulation are shown by black (blue) points. We find that isolated halos with V_circ(2 kpc) ≲ 20 km s⁻¹ and $30\,\mathrm{km}\,{{\rm{s}}}^{-1}\lt {V}_{\max }\lt 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ are extremely rare in CDM, for which no halos in our sample fall within the 1σ uncertainties of the observed UDGs' V_circ(2 kpc) values in the relevant ${V}_{\max }$ range. Given the small volume of our high-resolution region, which only yields ≈10³ isolated halos in the UDG mass range, our result is consistent with the findings in K22 that (i) UDG-like halos are ≳3σ to 5σ away from the median CDM expectation in IllustrisTNG, and (ii) even extremely low-concentration CDM halos inner profiles are inconsistent with the circular velocities of the observed UDGs.

Zoom In Zoom Out Reset image size

Meanwhile, our Group SIDM model diversifies the V_circ(2 kpc)–${V}_{\max }$ relation and produces a population of halos that is consistent with UDGs' inner circular velocities. Consistent with K22, we find that these UDG analog halos preferentially form later and have higher spin parameters than the rest of the isolated SIDM halo population; thus, a detailed study of our UDG analogs' secondary halo properties is a promising area for future work.

Next, the right panel of Figure 3 shows the density profiles of all isolated halos with $30\,\mathrm{km}\,{{\rm{s}}}^{-1}\lt {V}_{\max }\lt 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and M₂₀₀ > 5.5 × 10⁹ M_⊙ in our CDM and SIDM simulations. This ${V}_{\max }$ range encompasses the observed UDGs, and the mass cut ensures that all UDGs' halo masses are larger than the universal DM-to-baryon ratio using the 1σ lower limits from K22. Roughly half of the resulting halos in our SIDM simulation have kiloparsec-scale cores. Dashed lines show best-fit Read et al. (2016) density profiles of the observed UDGs from K22; our SIDM predictions span the UDGs' inferred profiles, and this level of agreement is again difficult to achieve in CDM.

We note that there are detailed differences between our simulated SIDM profiles and the UDG profiles from K22. Among the UDGs, AGC 114905 (Mancera Piña et al. 2022) has the lowest inner halo density and the largest core size; our SIDM halos are still too dense in the inner regions to fit this system. It is possible that AGC 114905's halo concentration is even lower than any halos in our simulation. To test this, we used a parametric model to predict SIDM profiles as a function of halo concentration (Yang et al. 2023b); we found that our SIDM model can fit AGC 114905 for c₂₀₀ ≈ 2.2 with M₂₀₀ ≈ 1.4 × 10¹⁰ M_⊙ (${V}_{\max }\approx 34\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and σ_eff/m ≈ 99 cm² g⁻¹). This concentration would be "7σ" below the median according to the standard c₂₀₀–M₂₀₀ relation (Dutton & Macciò 2014). However, the low-concentration tail of this relation is a power law, rather than a Gaussian. Thus, the probability of such a low-concentration halo is small but not negligible as found in the IllustrisTNG simulation K22, which has a much larger volume than our zoom-in.

For the remaining UDGs, our simulated SIDM halos are slightly more cored than the observed systems K22. In fact, these UDGs have rather sharply rising rotation curves and they do not strongly prefer a cored density profile, although their halos must have a low concentration. It is important to investigate whether the trend holds with more kinematic data (the fits in K22 were based on two data points for each rotation curve, except for AGC 114905 and 242019). Furthermore, the presence of gas, which is dynamically important, may also bring our predicted SIDM density profiles into better agreement with observations.

A more detailed comparison between simulated and observed systems will require modeling the baryonic contribution and potentially expanding the Read et al. (2016) density profile to encompass our diverse SIDM halos. It would also be interesting to explore the formation of gas-rich galaxies in our Group SIDM model using hydrodynamic simulations, since the gas distribution could be extended due to the presence of large DM cores, potentially providing a novel explanation of UDGs' observed sizes. Note that the majority of cored field halos in our SIDM simulation are isolated throughout their histories, differentiating our model from scenarios in which environmental effects like tidal stripping produce shallow CDM density profiles (e.g., see Jiang et al. 2019; Benavides et al. 2021, 2023).

In addition to core-forming halos, we also find a sizable population of isolated SIDM halos with higher-amplitude inner density profiles and smaller values of ${R}_{\max }$ than any halos in our CDM simulation (where ${R}_{\max }$ is the radius at which ${V}_{\max }$ is achieved), indicating efficient core collapse in our SIDM model. These core-collapsed halos' inner logarithmic density profile slopes are often as steep as −3; they appear near the one-to-one relation in Figure 3 (Yang et al. 2023a). ⁴ The fraction of isolated core-collapsed halos in our SIDM simulation is consistent with analytic estimates of t_c in our SIDM model, which predict that higher-concentration halos collapse more quickly (Essig et al. 2019).

Because our SIDM scenario diversifies isolated halos in both directions relative to CDM, it predicts that extremely dense, isolated galaxies with rapidly rising rotation curves should be observed in addition to the UDGs. Rapidly rising rotation curves are particularly informative if the baryonic contribution to the rotation curve is small; existing data sets may hint at the existence of such galaxies (Kuzio de Naray et al. 2010; Santos-Santos et al. 2020), motivating detailed comparisons to simulations like ours.

We have focused on two opposite extremes of low-mass systems, and their corresponding halo properties—a dense strong-lens perturber and DM-deficient gas-rich UDGs. Using the first group-scale simulation with large-amplitude, velocity-dependent DM self-interactions, we showed that SIDM can simultaneously produce halos at both extremes.

Our study reveals exciting directions for the study of strong lensing and UDG (sub)halos. It is timely to revisit predictions for the diversity of galactic rotation curves, which have mainly been studied for SIDM halos in the core-expansion phase (Creasey et al. 2017; Kamada et al. 2017; Ren et al. 2019; Zentner et al. 2022), in the context of strong DM self-interactions. In our SIDM scenario, isolated halos are found in both core-expansion and core-collapse phases, yielding an even larger diversity of inner halo profiles than less extreme SIDM models. In parallel, this physics should be detectable in future strong lensing studies using both gravitational imaging and flux-ratio statistics (e.g., see Minor et al. 2017; Meneghetti et al. 2020, 2023; Gilman et al. 2021, 2023; Yang & Yu 2021; Loudas et al. 2022; Şengül & Dvorkin 2022; Dhanasingham et al. 2023; Zhang et al. 2023); observations of MW satellites will also probe our SIDM model, which predicts that a substantial fraction of the low-mass subhalos that host faint dwarf galaxies are core collapsed (e.g., Kaplinghat et al. 2019; Correa 2021; Turner et al. 2021; Silverman et al. 2023; Slone et al. 2023; Yang et al. 2023a).

Exploring a range of strong, velocity-dependent cross sections will be necessary to robustly constrain our SIDM scenario. We expect that models with σ₀/m smaller by an ${ \mathcal O }(1)$ factor, keeping w ∼ 100 km s⁻¹, can also produce halos at both density extremes. Thus, observational constraints on the fraction of core-collapsed halos will narrow the favored SIDM parameter space. Such constraints will be enabled by combining (semi)analytic predictions for SIDM halo populations (Jiang et al. 2023; Shah & Adhikari 2023; Yang et al. 2023b, 2023d), galaxy–halo connection modeling techniques (Nadler et al. 2019), and observations of DM structure and galaxies over a wide mass range.

We thank Manoj Kaplinghat for helpful discussions, Pavel Mancera Piña for comments on the manuscript, and Demao Kong for providing data from K22. This work was supported by the John Templeton Foundation under grant ID #61884 and the U.S. Department of Energy under grant No. de-sc0008541 (D.Y. and H.-B.Y.). The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. The computations presented here were conducted through Carnegie's partnership in the Resnick HPC Center, a facility supported by Resnick Sustainability Institute at the California Institute of Technology.

We assume Rutherford-like scattering in the Born limit (Feng et al. 2010; Ibe & Yu 2010; Tulin et al. 2013) and perform the SIDM simulation using the viscosity cross section defined in Yang & Yu (2022), i.e.,

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{V}(v) & = & \displaystyle \frac{3}{2}\displaystyle \int d\cos \theta {\sin }^{2}\theta \displaystyle \frac{d\sigma }{d\cos \theta }\\ & = & \displaystyle \frac{6{\sigma }_{0}{w}^{6}}{{v}^{6}}\left[\left(2+\displaystyle \frac{{v}^{2}}{{w}^{2}}\right)\mathrm{ln}\left(1+\displaystyle \frac{{v}^{2}}{{w}^{2}}\right)-\displaystyle \frac{2{v}^{2}}{{w}^{2}}\right],\end{array}\end{eqnarray} \tag{ A1 }$$

where we take σ₀/m = 147.1 cm² g⁻¹ to be the same as that used in the MW analog simulation from Yang et al. (2023a), while setting w = 120 km s⁻¹. Considering a scenario where DM particles interact with a dark photon, we have ${\sigma }_{0}=4\pi {\alpha }_{\chi }^{2}/({m}^{2}{w}^{4})$ and w = m_ϕ c/m, where α_χ is the fine structure constant in the dark sector, m_ϕ is the mass of the dark photon, and c is the speed of light. Yang & Yu (2022) showed that the viscosity cross section provides an excellent approximation for modeling angular- and velocity-dependent DM self-scattering. If DM is made of one particle species, the viscosity cross section for Møller scattering would be more appropriate (Girmohanta & Shrock 2022; Yang & Yu 2022), as it includes both t- and u-channel contributions, as well as their interference term. However, the difference is small and we can reinterpret the cross section of Rutherford scattering to that of Møller scattering by rescaling particle parameters within less than 10% (Girmohanta & Shrock 2023). In addition, our simulation assumes all particles participate in scattering; see Yang & Yu (2022) for a self-consistent interpretation in terms of quantum statistics.

For simplicity, we fit the exact viscosity cross section Equation (A1) with the function

$$\begin{eqnarray}&&{\sigma }_{V,\mathrm{fit}}(v)=\displaystyle \frac{{\sigma }_{0}}{{\left[1+{\left(v/{v}_{0}\right)}^{\delta }\right]}^{\beta }},\end{eqnarray} \tag{ A2 }$$

where σ₀/m = 147.1 cm² g⁻¹, v₀ = 190 km s⁻¹, δ = 1.72, and β = 2. Figure 4 shows that this fit matches the viscosity cross section for our Group SIDM model at the percent level over the relative velocity range of interest. Thus, our implementation accurately captures the underlying SIDM cross section.

Zoom In Zoom Out Reset image size

We compute surface mass density by integrating the 3D density profile along an arbitrary z-axis as

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={\int }_{-{R}_{\mathrm{vir}}}^{{R}_{\mathrm{vir}}}\rho (\sqrt{{R}^{2}+{z}^{2}})\,{dz},\end{eqnarray} \tag{ B1 }$$

assuming that the halo is spherically symmetric. Accordingly, the 2D mass enclosed within 1 kpc is computed as

$$\begin{eqnarray}&&{M}_{2{\rm{D}}}(1\,\mathrm{kpc})=2\pi {\int }_{0}^{1\,\mathrm{kpc}}{\rm{\Sigma }}(R)R\,{dR}.\end{eqnarray} \tag{ B2 }$$

We construct an interpolation function to smoothly model surface density versus projected radius by evaluating this relation on a log–log grid. The slope ${\gamma }_{2{\rm{D}}}\equiv d\mathrm{ln}{\rm{\Sigma }}/d\mathrm{ln}R$ is then numerically extracted from the interpolation function as a function of $\mathrm{ln}R$. To suppress numerical uncertainties, we present averaged γ_2D between 0.75 and 1.25 kpc, following M21. Figure 5 shows the 3D density profile (left), surface density profile (middle), and projected density profile slope for the tNFW strong-lensing perturber analog from our SIDM simulation (corresponding to the magenta density profile in the right panel of Figure 2). All three profiles are measured smoothly down to scales comparable to our simulation's softening length, and the results in the inner regions are robust to changes in the limits of integration in Equation (B1).

Zoom In Zoom Out Reset image size