Diversifying halo structures in two-component self-interacting dark matter models via mass segregation

Abstract

Self-interacting dark matter (SIDM), through gravothermal evolution driven by elastic self-scatterings, offers a compelling explanation for the observed diversity of inner halo densities. In this work, we investigate SIDM dynamics in a two-component dark matter model with mass ratios of order unity, motivated by an asymmetric dark matter framework that naturally evades constraints from relic abundance and mediator decay, while enabling strong, velocity-dependent self-interactions. We show that cross-component scatterings significantly enhance mass segregation, driving the formation of dense, core collapsed-like halos. This effect couples naturally to SIDM-induced diversity, introducing a new mechanism for generating structural variations beyond those arising from gravothermal evolution alone. Our results reveal a novel mechanism for reconciling SIDM with small-scale observational tensions by enabling shifts in central densities while preserving the flexibility to generate diverse halo structures. We further highlight that halo structural diversity may serve as a diagnostic of dark sector composition, opening a new observational window into the particle nature of SIDM.

Introduction. While the standard cold dark matter (CDM) paradigm has been well-tested at large scales, it faces challenges at small scales, particularly in accounting for the diversity of observed galactic rotation curves and inner halo densities Bullock and Boylan-Kolchin (2017); Salucci (2019); Gentile et al. (2004); van Dokkum et al. (2022); Bonaca et al. (2018). One promising solution is self-interacting dark matter (SIDM), where dark matter particles interact through elastic scatterings, leading to gravothermal evolution that results in both cored or cuspy density profiles in halos Tulin and Yu (2018); Adhikari et al. (2022); Spergel and Steinhardt (2000); Kochanek and White (2000); Kamada et al. (2017); Ren et al. (2019); Santos-Santos et al. (2020); Correa et al. (2022); Yang et al. (2023); Zavala et al. (2019); Kaplinghat et al. (2019); Turner et al. (2021); Slone et al. (2023); Silverman et al. (2022); Yang and Yu (2022); Yang et al. (2024a); Yang (2024); Hou et al. (2025). Recent studies Nadler et al. (2023); Kong et al. (2022); Zhang et al. (2024, 2025a, 2025b); Sameie et al. (2020); Yang et al. (2020); Kong and Yu (2025); Balberg and Shapiro (2002); Pollack et al. (2015); Choquette et al. (2019); Feng et al. (2022); Meshveliani et al. (2023); Alonso-Álvarez et al. (2024); Jiang et al. (2025) further show that it may explain a variety of extreme observations, including dark matter deficient galaxies, supermassive black holes at high redshifts, strong lensing perturbers, excessive small-scale lenses, and even the final parsec problem of binary black hole mergers. However, in minimal SIDM scenarios involving a single dark matter species coupled to a light mediator, efficient annihilation channels typically lead to signals in conflict with Cosmic Microwave Background (CMB) and other indirect detection constraints Bringmann et al. (2017). This tension motivates the exploration of more complex scenarios, such as asymmetric dark matter models Kaplan et al. (2009); Petraki et al. (2012). These models propose a dark sector with two particle species exhibiting an asymmetry, analogous to the matter-antimatter asymmetry in the visible sector Kaplan et al. (2009); Petraki et al. (2012); Cyr-Racine and Sigurdson (2013); Petraki and Volkas (2013); Petraki et al. (2014); Zurek (2014); Han et al. (2024). In the early universe, the symmetric components efficiently annihilate, leaving behind the asymmetric population that sets the relic dark matter abundance. Since the relic density is no longer tied to the annihilation cross section, these models allow greater flexibility to realize strong self-interactions. Moreover, just as the proton and electron have different masses but interact through the Coulomb force, the two dark matter components can have different masses and interact among themselves Petraki et al. (2014); Cyr-Racine and Sigurdson (2013).

The introduction of a light mediator naturally gives rise to a velocity-dependent SIDM cross section, a framework that has been extensively studied in recent years and has shown promise in addressing a range of anomalous observations Valli and Yu (2018); Zavala et al. (2019); Correa (2021); Ebisu et al. (2022); Silverman et al. (2022); Feng et al. (2009, 2010); Buckley and Fox (2010); Loeb and Weiner (2011); Tulin et al. (2013); Cyr-Racine et al. (2016); Kaplinghat et al. (2016); Agrawal et al. (2017); Vogelsberger et al. (2012, 2016); Robertson et al. (2019); Nadler et al. (2020); Turner et al. (2021); Bhattacharyya et al. (2022); Zeng et al. (2022); Robertson et al. (2017); Kahlhoefer et al. (2014); Fischer et al. (2021, 2022); Nadler et al. (2025). These studies typically assume either a single dark matter component that is its own antiparticle (i.e., a Majorana particle) or a two-component scenario where both components have equal or nearly equal masses. In this work, we present cosmological $N$-body simulations of a two-component SIDM model with a $3:1$ mass ratio. We show that interactions between the two species lead to pronounced mass segregation, modifying the range of inner density profiles compared to single-component SIDM. Previous simulations involving two-component dark matter have focused on either large mass ratios, as in atomic dark matter, or near-degenerate states Low et al. (2025); Roy et al. (2023); Schutz and Slatyer (2015), but neither has considered the effects of gravothermal evolution. Our framework, though simple, is broadly applicable and introduces a more versatile mechanism that generates novel structural features. These may help resolve existing tensions and ongoing debates stemming from recent precision constraints on SIDM from isolated galaxies, bullet clusters, halo shapes, and the Tully-Fisher relation Correa et al. (2025); Kong et al. (2024); Yang et al. (2024b); Ebisu et al. (2022); Vargya et al. (2022); Correa (2021).

Cosmological simulation of a two component SIDM model. We simulate a dark QED-like model consisting of two dark matter components, $\chi_{H}$ and $\chi_{L}$, with a mass ratio $m_{H}/m_{L}=3$ and equal number densities. For comparison, we also simulate a model with a single dark matter component, $\chi_{0}$, using the same total number of particles and assigning it a mass $m_{0}=2m_{H}/3$ to match the total mass of the two-component system. The scattering cross sections are modeled using the Møller and Rutherford equations parametrized by $\sigma_{0}/m$ and $w$, following Ref. Yang and Yu (2022). We assume a fixed mediator mass and a dark fine-structure constant in the perturbative regime, which results in different $\sigma_{0}/m$ and $w$ values for each scattering channel. For the $\chi_{0}\texttt{-}\chi_{0}$ case, we set $\sigma_{0}/m=147.1~{}\rm cm^{2}/g$ and $w=24.33~{}\rm km/s$, following the VD100 model in Ref. Yang et al. (2023), which yields a velocity-dependent cross section analogous to Ref. Turner et al. (2021). See Supplemental Material for further details.

Figure 1: Projected dark matter densities of MW analogs in the two-component SIDM2c simulation. Densities of the heavier (blue-to-white) and lighter (red-to-white) species are shown, with brighter regions indicating higher densities. The inset shows viscosity cross sections Tulin et al. (2013) for the $\chi_{H}-\chi_{H}$ (red), $\chi_{L}-\chi_{L}$ (orange), and $\chi_{H}-\chi_{L}$ (magenta) interactions in the SIDM2c simulation. For comparison, the $\chi_{0}-\chi_{0}$ (blue) interaction in the one-component case is shown as a dashed curve. All channels assume the same mediator mass and coupling. Model parameters for these interactions can be derived by fixing $\sigma_{0}/m=147.1~{}\rm cm^{2}/g$ and $w=24.33~{}\rm km/s$ in the $\chi_{0}\text{--}\chi_{0}$ case; see Supplemental Material for details.

We perform cosmological zoom-in simulations of the AGORA Milky Way (MW) analog system Kim et al. (2013); Roca-Fabrega et al. (2021) using the Gadget2 program Springel (2005); Springel et al. (2001), adopting $\Omega_{m}=0.272$, $\Omega_{\Lambda}=0.728$, $h=0.702$, and a box size of $L=60~{}\rm Mpc/h$. The initial condition at $z=100$ is generated using the MUSIC program Hahn and Abel (2011), with the transfer function computed by the CAMB program Lewis and Challinor (2011) based on these parameters. We randomly divide the initial particles into two components and conduct five cosmological simulations in both CDM and SIDM. These simulations include both single-component (1c) and two-component (2c) models, incorporating self- and cross-interaction (SIDMx) channels, as summarized in Table 1. The SIDM simulations are based on the module implemented in Refs. Yang and Yu (2022); Yang et al. (2023), which supports both types of differential cross sections considered in this work. In the high-resolution region, the particle mass for the lighter component is $12(1.5)\times 10^{4}~{}{\rm M_{\odot}}/h$ in SIDM (CDM) simulations, enabling us to probe structures down to a few times the softening length, $\epsilon=0.16(0.06)~{}{\rm kpc}/h$.

Simulation	Components	Interaction Type	Color
CDM1c	$\chi_{0}$	CDM	Green
SIDM1c	$\chi_{0}$	$\chi_{0}-\chi_{0}$	Orange
CDM2c	$\chi_{H}$, $\chi_{L}$	CDM	Blue
SIDM2c	$\chi_{H}$, $\chi_{L}$	$\chi_{H,L}-\chi_{H,L}$	Red
SIDMx	$\chi_{H}$, $\chi_{L}$	$\chi_{H}-\chi_{L}$	Magenta

Table 1: Summary of the cosmological simulations performed in this work. The one-component (1c) simulations use a particle mass of $m_{0}=(2/3)m_{H}$. The two-component SIDM2c simulations set $m_{H}=3m_{L}$ and include both self- and cross-component interactions. The SIDMx simulation considers only cross-component scattering between $\chi_{H}$ and $\chi_{L}$.

Figure 1 presents the projected dark matter density of the MW analog in the SIDM2c simulation at $z=0$. The heavier and lighter species are shown in blue-to-white and red-to-white, respectively, with brighter (darker) regions indicating higher (lower) densities. The inset panel illustrates the velocity-dependent interactions in the two-component (solid) and one-component (dashed) models through the viscosity cross sections $\sigma_{V}/m$ Tulin et al. (2013); Cline et al. (2014); Boddy et al. (2016); Blennow et al. (2017); Alvarez and Yu (2020); Colquhoun et al. (2021); Yang and Yu (2022). Compared to the one-component case, the two-component cross sections have lower normalization parameters $\sigma_{0}/m$ but larger velocity transition scales $w$. For quantitative analysis, we modify the Rockstar halo finder to accommodate two-component dark matter, which enables identification and characterization of halos in simulation snapshots. To facilitate density profile measurements, we restrict our analysis to MW subhalos with virial masses $M_{\rm vir}>5\times 10^{8}~{}\rm M_{\odot}/h$, ensuring a minimum of approximately $2000$ particles per halo. Convergence at this resolution has been demonstrated in Appendix B of Ref. Yang et al. (2023). We adopt $214\ \rm kpc/h$ as an estimate of the virial radius for MW analogs across all simulations.

Figure 2: Halo density profiles of MW subhalos of $M_{\rm vir}>5\times 10^{8}~{}\rm M_{\odot}/h$ in simulated models. The main panels show CDM1c (green) and SIDM1c (orange) on the left panel, SIDM2c (red) and SIDMx (magenta) on the middle and right panels, respectively. The CDM2c (blue) is shown in all panels for comparison. The smaller sub-panels display the fractional number density of the lighter species, $n_{L}/(n_{H}+n_{L})$. While the CDM2c case closely resembles the CDM1c scenario, with only slightly shallower inner regions, SIDM significantly redistributes the lighter species to larger radii, leading to denser inner regions dominated by the heavier species. In the SIDM2c case, self-interactions among the heavier species introduce both core and cusp profiles.

Mass segregation and diversity. In the two-component scenario, interactions between species of different masses tend to equalize their kinetic energies, causing the more massive species to sink into the inner halo regions—an effect known as mass segregation. In CDM, this effect is governed by the relaxation time due to gravitational scatterings, which scales as approximately $N/(10\ln N)$ times the particle crossing time. Since the particle number $N$ is extremely large, the relaxation time can be exceedingly long. In contrast, mass segregation by SIDM proceeds via collisional relaxation, with a timescale independent of $N$ and given by $1/(\rho\sigma v)$, where $\rho$ is the local density, $v$ the relative velocity, and $\sigma$ the scattering cross section. A larger $\sigma$ thus accelerates relaxation and enhances mass segregation.

Based on the particle positions and masses from the modified Rockstar program Behroozi et al. (2013), we compute the density profiles of subhalos in MW analogs and show them in Fig. 2 with radii normalized by the halos’ virial radii. The subpanels display the radial distribution of the lighter species, quantified by $n_{L}/(n_{H}+n_{L})$, as a measure of mass segregation. Compared to CDM2c (left), both SIDM2c (middle) and SIDMx (right) show a clear expulsion of the lighter species to larger radii, with its central fraction dropping to nearly zero. The nearly identical segregation patterns in SIDM2c and SIDMx confirm that cross-component interactions primarily drive this effect. Meanwhile, the agreement between CDM1c and CDM2c density profiles demonstrates the negligible impact of mass segregation in CDM.

To facilitate comparison with the one-component case, we include SIDM1c density profiles as dashed curves in the left panel. For the chosen cross section, SIDM1c develops prominent cores with central densities falling below observational estimates (e.g., Ref. Kaplinghat et al. (2019)). In contrast, mass segregation in the two-component models promotes the formation of denser central regions in many subhalos, resembling the core-collapse phase seen in single-component SIDM. As a result, both SIDM2c and SIDMx retain cuspy profiles. In SIDM2c, self-interactions among heavy particles further enhance structural diversity, producing both cores and cusps. As we will show, this variation in inner densities enables SIDM2c to better reproduce observational data.

The interplay between mass segregation and self-interactions thus broadens the phenomenology of SIDM. In particular, segregation allows some halos to exhibit core-collapsing features even at modest cross sections where one-component models would form overly large cores. The extent of this effect depends on multiple factors: the mass difference and cross-component scattering rates determine the degree of mass segregation, while the self-interaction cross section of the heavier component regulates the diversity in density profiles. These effects intertwine, collectively reshaping small-scale structures in the inner halo regions. In the limit where the two particle species have equal masses, mass segregation disappears, and the two-component model effectively reduces to the one-component SIDM scenario, although cross-species interactions remain governed by a Rutherford-like scattering cross section.

Inner halo densities and the fraction of lighter dark matter. In SIDM, gravothermal evolution has been identified as a key mechanism driving diversity in inner halo structures. Here, we demonstrate that cross-component scattering introduces an additional source of diversity via mass segregation, which can act in conjunction with the intra-component gravothermal evolution. To quantify this effect, we evaluate the inner halo density, $\rho_{\rm in}\equiv\rho(r=150~{}\rm pc)$, and the fractional number density of the lighter component, presenting the distributions for MW subhalos in Fig. 3. As the radius of $150$ pc approaches the softening length in our simulations, we fit the density profiles using a parametric model to obtain more robust results Yang et al. (2024a); see the Supplemental Material for details. In the literature, $\rho_{\rm in}$ has been reconstructed under various density profile models to test the diversity of observed MW satellites Andrade et al. (2024); Kaplinghat et al. (2019); Read et al. (2019, 2018); Hayashi et al. (2020). The fractional number density of the lighter component, defined as $f_{L}(<r)\equiv n_{L}(<r)/[n_{H}(<r)+n_{L}(<r)]$, and evaluated within $0.2R_{\rm vir}$, enables us to further probe the correlation between mass segregation and halo diversity. This provides a means to quantify cross-component SIDM as a new source of structural diversity.

Figure 3: Extrapolated inner halo density, $\rho_{\rm in}=\rho(r=150{\ \rm pc})$, vs fractional number density, $f_{L}(<0.2R_{\rm vir})$, for MW subhalos of $M_{\rm vir}>5\times 10^{8}~{}\rm M_{\odot}/h$ in two-component simulations SIDM2c (red), SIDMx (magenta), and CDM2c (blue). For comparison, the $\rho_{\rm in}$ of observed MW satellites, derived from Ref. Kaplinghat et al. (2019) under the assumptions of NFW (gray) and isothermal (brown) profiles, are shown as tick marks along the y-axis, with their Gaussian fits overlaid to highlight the probability distributions. The size of the points is proportional to the virial radius of the corresponding halos.

To assess how the diversity in our simulations compares with observations, we use data compiled in Ref. Kaplinghat et al. (2019), which analyzes MW satellite galaxies, including both classical dwarfs and ultrafaints, and extracts inner halo densities at $150$ pc using two extrapolation methods: one assuming Navarro-Frenk-White (NFW) profiles Navarro et al. (1997) and the other assuming cored isothermal (cISO) profiles Kaplinghat et al. (2016), the latter intended to model SIDM halos. We represent these extrapolated densities as horizontal ticks along the $y$-axis, with NFW-based values in gray and cISO-based values in brown. To compare data with our simulation results, we project the $\rho_{\rm in}$ of the selected subhalos onto the right panel, fitting their values in log space and display the resulting probability distributions as shaded Gaussians along the $y$-axis, using the same color scheme as the simulation benchmarks.

In the main panel of Fig. 3, we show the distribution of $\rho_{\rm in}$ vs $f_{L}(<0.2\,R_{\rm vir})$ for the three two-component simulations. In the CDM2c (blue) case, all subhalos have $f_{L}$ values clustered near $0.5$, whereas in SIDM2c and SIDMx, $f_{L}$ spans nearly the full range from $0$ to $0.5$, reflecting substantial variation in mass segregation across halos.

The SIDMx case includes only cross-component scatterings, and thus lacks gravothermal evolution in the conventional sense. Nonetheless, we observe a broader $\rho_{\rm in}$ distribution with a median shifted upward relative to CDM2c, indicating that cross-component interactions alone can generate significant diversity in inner halo structure.

When self-interactions are introduced among the same components (SIDM2c), the scatter further increases, yielding a larger population of core-like halos and a slightly lower median $\rho_{\rm in}$. Interestingly, the resulting diversity is consistent with CDM2c when inner densities are extrapolated using NFW profiles, but aligns more closely with the two SIDM cases under the assumption of cISO profiles. This highlights that current observational uncertainties limit our ability to discriminate between CDM and SIDM scenarios. Improved measurements and larger satellite samples will be essential to resolve the inner structure diversity in MW satellites.

Apart from the overall scatter, we find a moderate anti-correlation between $\rho_{\rm in}$ and $f_{L}(<0.2\,R_{\rm vir})$. In SIDMx, the Pearson correlation coefficient is $r=-0.36$ with a p-value of $0.14$, indicating a weak but coherent trend. The anti-correlation becomes stronger in SIDM2c, with $r=-0.60$ and $p=0.008$, suggesting that increased mass segregation is associated with enhanced central densities.

Discussion and conclusion. We have investigated the impact of mass segregation in two-component SIDM models and its implications for the diversity of inner halo structures. Mass segregation, a well-known phenomenon in astrophysical systems involving globular clusters and black holes, can also arise in multicomponent dark matter scenarios. In CDM, the segregation timescale is too long to yield observable consequences. By contrast, cross-component scatterings in SIDM can significantly accelerate this process, leading to substantial reshaping of halo density profiles.

Our results show that SIDM naturally connects mass segregation to structural diversity. In particular, we find that segregation can drive a subset of halos toward higher central densities, resembling core-collapse-like configurations. This behavior complements standard gravothermal evolution and introduces an additional mechanism for generating diversity in inner halo densities. While diversity in single-component SIDM is typically interpreted as a signature of gravothermal evolution, our findings highlight its potential as a probe of dark sector composition. Precise characterization of inner halo structures may offer clues to the presence of multiple self-interacting components, opening a new observational window into the particle nature of SIDM.

Although our analysis focuses on MW–like halos, the underlying mechanism is general and may operate in more massive systems. In such cases, dense, segregated substructures may leave observable imprints, such as enhanced gravitational lensing. More broadly, our results suggest a novel approach for resolving potential tensions in SIDM by enabling shifts in central densities without reducing the overall diversity.

Appendix A Implementation of differential scatterings

We consider a simple dark QED-like model with two dark matter species, $\chi_{H}$, $\chi_{L}$ of mass $m_{H}$ and $m_{L}$, interacting via a massive vector mediator of mass $m_{V}$ and coupling $e_{D}$ Kaplan et al. (2009); Petraki et al. (2012); Petraki and Volkas (2013); Petraki et al. (2014); Zurek (2014). The corresponding dark fine-structure constant is defined as $\alpha_{D}=e_{D}^{2}/(4\pi)$. We work in a weakly coupled perturbative regime, where differential cross sections can be calculated at tree level. To describe the interactions among $\chi_{L}$ or $\chi_{H}$ (self-interactions) and between $\chi_{L}$ and $\chi_{H}$ (cross-component interactions), we adopt the parameterizations of the Møller and Rutherford cross sections from Ref. Yang and Yu (2022).

Rutherford scattering applies to the interaction of two distinct species. When the masses of the two particles differ, the center-of-mass frame differential cross section can be generalized as Feng et al. (2010); Ibe and Yu (2010)

$$\frac{d\sigma}{d\cos\theta}=\frac{\sigma_{0}w^{4}}{2\left[w^{2}+{v^{2}}\sin^{2}(\theta/2)\right]^{2}},$$ (1)

where $\sigma_{0}\equiv\pi\alpha^{2}_{D}/(\mu^{2}w^{4})$, $w=m_{V}/(2\mu)$, and $\mu=m_{L}m_{H}/(m_{L}+m_{H})$ is the reduced mass of the two-particle system. Here, $v$ is the relative velocity and $\theta$ is the scattering angle. Integrating over $\theta\in(0,\pi)$ gives the total cross section

$$\sigma_{\rm tot}={\sigma_{0}}/(1+{v^{2}}/{w^{2}}),$$ (2)

which we use to determine scattering probabilities in our simulations. For further details, see Ref. Yang and Yu (2022).

Møller scattering describes the scattering of identical particles ($\chi\chi\rightarrow\chi\chi$). The differential cross section is Yang and Yu (2022); Girmohanta and Shrock (2022)

$$\dfrac{d\sigma}{d\cos\theta}=\frac{2\sigma_{0}w^{4}\left[\left(3\cos^{2}\theta+1\right)v^{4}+4v^{2}w^{2}+4w^{4}\right]}{\left(\sin^{2}\theta v^{4}+4v^{2}w^{2}+4w^{4}\right)^{2}},$$ (3)

where $\theta$ is restricted to $(0,\pi/2)$. The factor of $2$ ensures the total cross section matches that in Ref. Yang and Yu (2022). Upon integration, one obtains

$$\sigma_{\rm tot}=\sigma_{0}w^{4}\left[\frac{1}{v^{2}w^{2}+w^{4}}+\frac{1}{v^{4}+2v^{2}w^{2}}\ln\left(\frac{w^{2}}{v^{2}+w^{2}}\right)\right].$$ (4)

In our simulations, we include an additional $1/2$ phase space factor to account for the identical nature of the two initial-state particles.

In this work, we fix the mediator mass $m_{V}$, the coupling constant $\alpha_{D}$, and the total particle number across different simulation types, which result in varying $\sigma_{0}/m$ and $w$ parameters. In the single-component case, we denote the dark matter particles by $\chi_{0}$ and set their mass to $m_{0}=(2/3)m_{H}$. Table 2 summarizes the corresponding SIDM model parameters for each type of scattering.

Scatterings	Mass relation	$\frac{\sigma_{0}}{m}$ relation	$w\ \rm(km/s)$ relation
$\chi_{0}{\texttt{-}}\chi_{0}$	$m_{0}=\frac{2}{3}m_{H}$	$\frac{\sigma_{M}}{m_{0}}=147.1\ \rm\frac{cm^{2}}{g}$	$w_{M}=24.33\ \rm\frac{km}{s}$
$\chi_{H}-\chi_{H}$	$m_{H}=\frac{3}{2}m_{0}$	$\frac{\sigma_{H}}{m_{H}}=\frac{3}{2}\frac{\sigma_{M}}{m_{0}}$	$w_{1}=\frac{3}{2}w_{M}$
$\chi_{L}-\chi_{L}$	$m_{L}=\frac{1}{3}m_{H}$	$\frac{\sigma_{L}}{m_{L}}=\frac{1}{3}\frac{\sigma_{H}}{m_{H}}$	$w_{2}=3w_{1}$
$\chi_{H}-\chi_{H}$	$m_{H}=3m_{L}$	$\frac{\sigma_{x}}{m_{H}}=\frac{1}{4}\frac{\sigma_{H}}{m_{H}}$	$w_{x}=2w_{1}$

Table 2: SIDM model parameters for different scattering processes in a dark QED theory with fixed $m_{V}$ and $\alpha_{D}$. The mass assignments ensure the total particle number remains the same in all simulations.

Figure 4: The parametric model fitting results (solid black) for the density profiles of MW subhalos with $M_{\rm vir}>5\times 10^{8}~{}\rm M_{\odot}/h$ in CDM2c (left), SIDM2c (middle), and SIDMx (right) simulations. The relative differences between the simulated (Sim) and fitted (Fit) curves are measured as $2\rm(Fit-Sim)/(Fit+Sim)$ and shown in the sub-panels. Compared with the profiles in the main article, more bins are used here to improve the fit. The fitted curves are used to obtain the $\rho_{\rm in}=\rho(r=150{\ \rm pc})$, as referenced in the main text.

The viscosity cross section $\sigma_{V}$ is calculated by weighting the differential cross section with a kernel appropriate for kinetic-theory studies of viscosity and heat conductivity Tulin et al. (2013):

$$\sigma_{V}=\frac{3}{2}\int d\cos\theta\sin^{2}\theta\frac{d\sigma}{d\cos\theta}.$$ (5)

This quantity is of particular interest for studying the transport properties of SIDM halos and is used to generate the inset of Fig. 1 in the main text. See Ref. Yang and Yu (2022) for analytic results for the Møller and Rutherford scattering cases.

Appendix B Fitting the density profiles

The radius of $150$ pc, at which we evaluate the inner halo densities $\rho_{\rm in}$, lies below the softening length of the simulations. To obtain more reliable estimates of $\rho_{\rm in}$, we fit the simulated density profiles and extract values from the resulting fits.

We find that the parametric profile commonly used to model one-component SIDM halos Yang et al. (2024a); Yang (2024); Yang et al. (2025); Hou et al. (2025) can still be used to describe halos in two-component dark matter. Specifically, we fit the density profiles using the following $\beta 4$ profile:

$$\rho_{\rm\beta 4}(r)=\frac{\rho_{s}}{\frac{\left(r^{4}+r_{c}^{4}\right)^{1/{4}}}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}},$$ (6)

where $\rho_{s}$, $r_{s}$, and $r_{c}$ are free parameters that evolve with the normalized time variable $\tau\equiv t/t_{c}$. We treat $\tau$, $V_{\rm max}$, and $R_{\rm max}$ as fitting parameters and convert them to $\rho_{s}$, $r_{s}$, and $r_{c}$ using the relations provided by the parametric model for SIDM halos in Ref. Yang et al. (2024a). Specifically, given $\tau$, $V_{\rm max}$, and $R_{\rm max}$, we first evaluate

$$\hat{V}_{\rm max} = 1+0.1777\tau-4.399\tau^{3}+16.66\tau^{4}-18.87\tau^{5}$$

$$+9.077\tau^{7}-2.436\tau^{9}$$

$$\hat{R}_{\rm max} = 1+0.007623\tau-0.7200\tau^{2}+0.3376\tau^{3}$$

$$-0.1375\tau^{4},$$

obtaining NFW scale parameters as

$$r_{s,0} = \frac{R_{\rm max}}{\hat{R}_{\rm max}2.1626},$$ (9)

$$\rho_{s,0} = \frac{1}{G}\left(\frac{V_{\rm max}}{1.648r_{s,0}\hat{V}_{\rm max}}\right)^{2}.$$

Based on the $\rho_{s,0}$, $r_{s,0}$, and $\tau$, the three parameters in the $\beta 4$ profile is evaluated as

$$\frac{\rho_{s}}{\rho_{s,0}} = 2.033+0.7381\tau+7.264\tau^{5}$$

$$-12.73\tau^{7}+9.915\tau^{9}$$

$$+(1-2.033)(\ln 0.001)^{-1}\ln\left(\tau+0.001\right),$$

$$\frac{r_{s}}{r_{s,0}} = 0.7178-0.1026\tau+0.2474\tau^{2}-0.4079\tau^{3}$$

$$+(1-0.7178)(\ln 0.001)^{-1}\ln\left(\tau+0.001\right),$$

$$\frac{r_{c}}{r_{s,0}} = 2.555\sqrt{\tau}-3.632\tau+2.131\tau^{2}-1.415\tau^{3}$$

$$+0.4683\tau^{4}.$$

Figure 4 presents the fitted (Fit) and simulated (Sim) density profiles for the CDM2c, SIDM2c, and SIDMx simulations, from left to right. The relative differences, measured as $2{\rm(Fit-Sim)/(Fit+Sim)}$, are shown in the subpanels and indicate good agreement across all cases. While the relative differences occasionally reach $\sim 50\%$, these fluctuations are driven by statistical noise in the SIDM simulations and do not exhibit any systematic bias.