# A Parametric Model for Self-Interacting Dark Matter Halos

Daneng Yang,111Corresponding author. ID Ethan O. Nadler, ID Hai-Bo Yu, ID Yi-
Ming Zhong ID

###### Abstract

We propose a parametric model for studying self-interacting dark matter (SIDM)
halos. The model uses an analytical density profile, calibrated using a
controlled N-body SIDM simulation that covers the entire gravothermal
evolution, including core-forming and -collapsing phases. By normalizing the
calibrated density profile, we obtain a universal description for SIDM halos
at any evolution phase. The model allows us to infer properties of SIDM halos
based on their cold dark matter (CDM) counterparts. As a basic application, we
only require two characteristic parameters of an isolated CDM halo at $z=0$.
We then extend the model to incorporate effects induced by halo mass changes,
such as major mergers or tidal stripping, making it applicable to both
isolated halos and subhalos. The parametric model is tested and validated
using cosmological zoom-in SIDM simulations available in the literature.

###### Contents

  1. 1 Simulation data
  2. 2 A universal density profile
  3. 3 Transforming CDM halos into SIDM halos
    1. 3.1 The basic approach
    2. 3.2 The integral approach
  4. 4 Modeling the evolution of SIDM isolated halos
  5. 5 Modeling the evolution of SIDM subhalos
  6. 6 Exploring SIDM parameter space using the parametric model
  7. 7 Conclusions

In self-interacting dark matter (SIDM), dark matter particles can interact
with each other, in addition to their interactions through gravity; see refs.
[1, 2] for reviews. The self-interactions can thermalize the inner region of
dark matter halos and modify its structure [3]. An SIDM halo could experience
gravothermal core-forming and -collapsing phases sequentially, and its inner
density profile varies with time dynamically. Compared to collisionless cold
dark matter (CDM), SIDM predicts more diverse dark matter distributions in
both isolated halos and subhalos, see, e.g., [4, 5]. In particular, the
phenomenon of gravothermal collapse [6, 7, 8] could be probed using
observations of the Milky Way satellite galaxies [9, 10, 11, 12, 13, 14, 15,
16, 17, 4, 5] and strong gravitational lensing events [18, 19, 20, 21, 22].
The collapsed central halo could also provide a seed for supermassive black
holes [23, 24, 25, 26, 27, 28].

Cosmological N-body simulations are essential for studying the formation and
evolution of SIDM halos in realistic environments, see, e.g., [29, 30, 31, 32,
33, 34, 35, 36, 37, 38, 39]. These simulations are often computationally
expensive as they need a high resolution in order to resolve the central halo,
where the dark matter collision rate is the highest. In some cases, controlled
high-resolution N-body simulations can be used to study SIDM predictions, see,
e.g., [40, 41, 42, 43, 44, 45, 46]. A conducting fluid model is also broadly
used [7, 8, 47, 9, 25, 48, 49, 50], as it can resolve the central region of a
collapsed halo, but the model is based on an idealized setup. Another
complication is that dark matter self-scattering is generally angular- and
velocity-dependent [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61], and it is
challenging to directly use these simulation tools in the search for entire
particle parameter space.

Ref. [62] introduced a semi-analytical model for inferring SIDM density
profiles from their CDM counterparts. It assumes that the dark matter
distribution in the inner halo follows an isothermal distribution, and
determines the distribution parameters through a matching procedure. The model
has been tested and validated in both controlled and cosmological SIDM
simulations with and without baryons [40, 42, 63, 64]. It has been used to
analyze the rotation curves of a large sample of spiral galaxies [65, 66, 67,
68]. However, the basic version of this model is only valid for isolated halos
in the core-forming phase.

In this work, we propose a parametric model that can “transfer” CDM halos into
their SIDM counterparts for both core-forming and -collapsing phases, and it
works for isolated halos and subhalos. Our model is based on two main
observations. For velocity- and angular-dependent dark matter self-scattering,
there is a constant effective cross section [60, 48, 69] that can dynamically
incorporate the variation of the halo’s mass and concentration over its
evolution history. Furthermore, gravothermal evolution of an SIDM halo is
largely self-similar [7, 6, 8, 23, 47, 69, 50]. We will normalize halo
evolution time with its collapse timescale; see also [69, 50]. The rescaled
evolution trajectories of characteristic halo parameters become universal, as
they do not have an explicit dependence on the cross section.

We will introduce an analytical density profile and calibrate it using a
controlled N-body SIDM simulation from ref. [60] for the entire evolution
history. Then, we normalize the calibrated density profile and obtain a
universal description for SIDM halos. We consider two approaches in applying
the parametric model. In the basic application, we take density parameters of
an isolated CDM halo at $z=0$ as input and reconstruct the evolution history
of its SIDM counterpart. We extend the model to incorporate effects induced by
halo mass changes, such as major mergers or tidal stripping, making it
applicable to both isolated halos and subhalos. We will validate the
parametric model using the cosmological zoom-in SIDM simulation of a Milky Way
analog in ref. [5]. Since our model is accurate and computationally
inexpensive, it can be used to efficiently scan large parameter space of
particle physics models of SIDM. We will provide an example for such an
exercise.

The rest of the paper is organized as follows. In Sec. 1, we discuss the
simulation data used in this work. In Sec. 2, we introduce the universal
density profile and calibrate it with the controlled N-body SIDM simulation 2.
In Sec. 3, we discuss applications of the parametric model based on basic and
integral approaches. We validate the parametric model with isolated halos and
subhalos from the cosmological zoom-in simulation, in Secs. 4 and 5,
respectively. We present an example of using the parametric model to constrain
particle physics parameters of SIDM in Sec. 6 and conclude in Sec. 7.

##  1 Simulation data

We use two different simulation datasets in this study. We first take the
results from an isolated N-body SIDM simulation in ref. [60] to calibrate a
universal density profile over the course of gravothermal evolution, an
essential ingredient of our parametric model. Then, we apply the model to a
large population of isolated halos and subhalos from a cosmological zoom-in
SIDM simulation of a Milky Way analog from ref. [5], and show the parametric
model can successfully predict their properties. Both simulations were
performed using the Gadget2 program [70], with an SIDM module that adopts
similar techniques described in refs. [71, 30].

We consider the “BM2” halo in ref. [60]. The halo has an NFW initial condition
with the scale density and radius $\rho_{s}=2.74\times 10^{8}~{}\rm
M_{\odot}/kpc^{3}$ and $r_{s}=0.141~{}\rm kpc$, respectively. The halo
contains four million particles and the simulation assumes various SIDM
scenarios, including Rutherford- and Møller-like scatterings, both having
novel angular- and velocity-dependence, as well as isotropic scatterings with
a constant cross section. Ref. [60] proposed a constant effective cross
section to model the halo evolution; see also refs. [48, 69] for a similar
approach. For a differential cross section of $d\sigma/d\cos\theta$, the
effective cross section is evaluated as

$\displaystyle\sigma_{\rm eff}$ $\displaystyle=$ $\displaystyle\frac{2\int
v^{2}dvd\cos\theta\frac{d\sigma}{d\cos\theta}\sin^{2}\theta
v^{5}\exp\left[-\frac{v^{2}}{4\nu_{\rm eff}^{2}}\right]}{\int
v^{2}dvd\cos\theta\sin^{2}\theta v^{5}\exp\left[-\frac{v^{2}}{4\nu_{\rm
eff}^{2}}\right]}$ $\displaystyle=$ $\displaystyle\frac{1}{512\nu_{\rm
eff}^{8}}\int
v^{2}dvd\cos\theta\frac{d\sigma}{d\cos\theta}v^{5}\sin^{2}\theta\exp\left[-\frac{v^{2}}{4\nu_{\rm
eff}^{2}}\right],$

where $\nu_{\rm eff}=0.64V_{\rm max,NFW}\approx 1.05r_{\rm
eff}\sqrt{G\rho_{\rm eff}}$ is a characteristic velocity dispersion of dark
matter particles in the halo. We have used the relations $\rho_{\rm
eff}=(V_{\rm max}/(1.648r_{\rm eff}))^{2}/G$ and $r_{\rm eff}=R_{\rm
max}/2.1626$, where $V_{\rm max}$ is the maximal circular velocity of a halo
and $R_{\rm max}$ its corresponding radius. For an NFW halo, $r_{\rm eff}$ and
$\rho_{\rm eff}$ reduce to $r_{s}$ and $\rho_{s}$, respectively. In
cosmological simulations, $V_{\rm max}$ and $R_{\rm max}$ are more useful for
specifying a halo, as they do not replay on a specific form of density
profiles.

To be concrete, we will consider the Rutherford scattering benchmark in ref.
[60] with the differential cross section

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

where $\sigma_{0}/m=2.4\times 10^{4}~{}{\rm cm^{2}/g}$ and $w=1~{}{\rm km/s}$.
The corresponding effective cross section is $\sigma_{\rm eff}/m=7.1~{}{\rm
cm^{2}/g}$. In this work, we will calibrate our model to the simulation with
$\sigma_{\rm eff}/m=7.1~{}{\rm cm^{2}/g}$ [60], which is rescaled from the one
with $\sigma_{\rm eff}/m=10~{}{\rm cm^{2}/g}$. 222The simulation assumes all
particles participate in scattering. See the relevant discussion in ref. [60]
for a self-consistent interpretation in terms of quantum statistics.

The cosmological zoom-in simulation in ref. [5] contains a Milky Way analog
and the differential cross section is given in eq. (1.2), with
$\sigma_{0}/m=147.1~{}\rm cm^{2}/g$ and $w=24.33~{}\rm km/s$. We consider
simulated halos above a mass threshold of $10^{8}~{}{\rm M_{\odot}}/h$ and
each halo contains more than $3500$ simulation particles. Applying this
selection condition, we obtain $125$ and $115$ subhalos for CDM and SIDM,
respectively. For isolated halos at $z=0$, we further require them to reside
between $0.3\textup{--}3~{}{\rm Mpc}$ the center of the main halo so that
these candidates do not suffer from numerical contamination. There are $647$
CDM and $626$ SIDM isolated halos. Among them, some are slingshot and
splashback halos, which have experienced close encounters with other halos at
low redshifts. In addition, $32$ subhalos and $58$ isolated halos have their
collapse timescale less than $10~{}{\rm Gyr}$ [5].

##  2 A universal density profile

In SIDM, a halo first develops a shallow density core, and then the core size
shrinks and central density increases. We introduce the following function to
model the density profile of an isolated halo over the course of gravothermal
evolution,

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

where $\rho_{s}$ and $r_{s}$ are scale density and radius, respectively, and
$r_{c}$ is the core size. The transition between inner core and outer NFW
profile is controlled by $\beta$. We found that increasing $\beta$ from $1$ to
$4$, the fit to the simulated halo can be improved significantly, and its
remains almost unchanged for $\beta\geq 4$. We fix $\beta=4$, as in ref. [56],
where a similar trick was applied for halos with a Hernquist density profile.
In the limit $r_{c}\to 0$, eq. (2.1) smoothly reduces to an NFW profile. The
density profile in eq. (2.1) only applies to isolated halos. We will discuss
subhalos in Sec. 5.

Figure 1:  Density profiles of the simulated BM2 halo at different times for
an effective constant cross section of $\sigma_{\rm eff}/m=7.1~{}{\rm
cm^{2}/g}$ (dashed-red); data from ref. [60]. Each density profile is fitted
using the density profile in eq. (2.1) (solid-blue).

To develop a universal solution for parameterizing SIDM halo evolution, we
need to use a dimensionless timescale and eliminate explicit dependence on the
cross section. This can be achieved by normalizing the evolution time with the
collapse timescale $t_{c}$, which is approximately a fixed multiple of the
relaxation time. We estimate $t_{c}$ as [7, 23, 47]

$\displaystyle t_{\rm c}$ $\displaystyle=$
$\displaystyle\frac{150}{C}\frac{1}{(\sigma_{\rm eff}/m)\rho_{\rm eff}r_{\rm
eff}}\frac{1}{\sqrt{4\pi G\rho_{\rm eff}}},$ (2.2)

where $C$ is constant that can be calibrated with N-body simulations [8, 47,
9, 48], and we fix $C=0.75$. For the BM2 halo we consider, $t_{c}\approx
28.7~{}{\rm Gyr}$ for $\sigma_{\rm eff}=7.1~{}\rm cm^{2}/g$ [60].

Figure 1 shows the simulated density profile of the BM2 halo at
$t/t_{c}=0,\,0.09,\,0.18,\,0.25$ $0.50$ and $1$ (dashed-red), as well as the
one fitted using the density profile in eq. (2.1) (solid-blue). Our fit is
performed over the entire range, upper to the virial radius, and the agreement
is very good. Even at $t/t_{c}=1$, the agreement is within $10\%$. While
introducing additional parameters could further improve the fit, we prioritize
simplicity over sub-$10\%$ accuracy in our approach.

We further fit the density profile of the simulated BM2 halo at successive
time intervals of $0.2~{}{\rm Gyr}$, and obtain the evolution trajectories for
$r_{s}$, $\rho_{s}$ and $r_{c}$ in eq. (2.1) over the dimensionless time
variable $\tilde{t}=t/t_{c}$. They can be expressed as the following functions

$\displaystyle\frac{\rho_{s}}{\rho_{s,0}}$ $\displaystyle=$ $\displaystyle
2.033+0.7381\tilde{t}+7.264\tilde{t}^{5}-12.73\tilde{t}^{7}+9.915\tilde{t}^{9}-(1-2.033)(\ln
0.001)^{-1}\ln\left(\tilde{t}+0.001\right),$
$\displaystyle\frac{r_{s}}{r_{s,0}}$ $\displaystyle=$ $\displaystyle
0.7178-0.1026\tilde{t}+0.2474\tilde{t}^{2}-0.4079\tilde{t}^{3}-(1-0.7178)(\ln
0.001)^{-1}\ln\left(\tilde{t}+0.001\right),$
$\displaystyle\frac{r_{c}}{r_{s,0}}$ $\displaystyle=$ $\displaystyle
2.555\sqrt{\tilde{t}}-3.632\tilde{t}+2.131\tilde{t}^{2}-1.415\tilde{t}^{3}+0.4683\tilde{t}^{4},$
(2.3)

where the subscript $``0"$ denotes the corresponding value of the initial NFW
profile. To enforce $\rho_{s}/\rho_{s,0}=1$ and $r_{s}/r_{s,0}=1$ at
$\tilde{t}=0$, we have chosen a term of $\alpha-(1-\alpha)(\ln
0.001)^{-1}\ln\left(\tilde{t}+0.001\right)$ and $\alpha$ is the fitting
parameter. The numerical factor $0.001$ inside the logarithm is introduced to
avoid singularity at $\tilde{t}=0$. Although $t_{c}$ may vary slightly with a
different calibration of $C$, the functional forms in eq. (2), including the
coefficients, are robust to this variation. Thus, our SIDM density described
in eqs. (2.1) and (2) is highly universal.

In figure 2, we present the evolution of $\rho_{s}/\rho_{s,0}$ (left),
$r_{s}/r_{s,0}$ (middle), and $r_{c}/r_{c,0}$ (right). In all cases, our
density model in eq. (2) (solid-blue) accurately describes the simulation
results ($\sigma_{\rm eff}/m=7.1~{}\rm cm^{2}/g$, solid-red). We further show
the results from the N-body simulation [60] with the differential cross
section ($d\sigma/d\cos\theta$, dotted-red), which agree with those from the
simulation with $\sigma_{\rm eff}/m=7.1~{}\rm cm^{2}/g$ (solid-red). This
comparison serves as a systematic test, showing that a differential cross
section can indeed be approximated by an effective constant cross section for
an isolated halo as proposed in ref. [60].

In figure 3, we show the density (left) and circular velocity (middle)
profiles obtained from the N-body simulation (dashed) and our model in eq. (2)
(solid) at three representative times, $\tilde{t}=0$, $0.18$, and $1$. The
overall agreement is good, with only a minor discrepancy at $\tilde{t}\approx
1$, which is expected from the fit shown in figure 1. Figure 3 (right) further
shows the evolution of central densities from the N-body simulation with the
effective cross section ($\sigma/m=7.1~{}\rm cm^{2}/g$; solid-red) and the one
with $d\sigma/d\cos\theta$ (dotted-red), as well as that obtained using eq.
(2) (solid-blue). We see that they agree with each other very well up to
$\tilde{t}=25/28.7\approx 0.87$, beyond which eq. (2) underestimates the
central density by $\sim 15\%$ at $\tilde{t}=1$.

We can use eq. (2) and obtain $V_{\rm max}$ and $R_{\rm max}$ from the
reconstructed halo rotation curve at each snapshot. For the BM2 halo with
$\sigma/m=7.1~{}{\rm cm^{2}/g}$, we find that the evolution trajectories of
$V_{\rm max}$ and $R_{\rm max}$ can be well fitted by

$\displaystyle\frac{V_{\rm max}}{V_{\rm max,0}}$ $\displaystyle=$
$\displaystyle
1+0.1777\tilde{t}-4.399\tilde{t}^{3}+16.66\tilde{t}^{4}-18.87\tilde{t}^{5}+9.077\tilde{t}^{7}-2.436\tilde{t}^{9}$
$\displaystyle\frac{R_{\rm max}}{R_{\rm max,0}}$ $\displaystyle=$
$\displaystyle
1+0.007623\tilde{t}-0.7200\tilde{t}^{2}+0.3376\tilde{t}^{3}-0.1375\tilde{t}^{4},$
(2.4)

where the subscript $``0"$ denotes the corresponding value of the initial NFW
profile.

Figure 2:  Evolution of normalized density parameters $\rho_{s}/\rho_{s,0}$
(left), $r_{s}/r_{s,0}$ (middle), and $r_{c}/r_{s,0}$ (right) from the
calibrated functions in eq. (2) (solid-blue) and the N-body SIDM simulation
with the effective constant cross section $\sigma/m=7.1~{}{\rm cm^{2}/g}$ [60]
(solid-red). For comparison, the SIDM simulation [60] based on the
differential scattering cross section $d\sigma/d\cos\theta$ is also shown
(dotted-red).

Figure 3:  Density (left) and circular velocity (middle) profiles from eq. (2)
(solid) and the SIDM simulation with the effective cross section
$\sigma/m=7.1~{}\rm{cm^{2}/g}$ [60] (dashed) at $\tilde{t}=0,~{}0.18$ and $1$.
The evolution curves of the central dark matter density (right) from eq. (2)
(solid-blue) and the simulations [60] with $\sigma/m=7.1~{}\rm cm^{2}/g$
(solid-red) and $d\sigma/d\cos\theta$ (dashed-red) are also shown.

##  3 Transforming CDM halos into SIDM halos

To apply the functions in eq. (2) for halos in cosmological simulations, we
need to determine $\rho_{s,0}$, $r_{s,0}$, $t_{c}$, and evolution time after
halo formation. In this section, we first introduce a basic approach for the
parametric model that simply utilizes the corresponding CDM halo parameters at
$z=0$ as input to make predictions for a given cross section specified in eq.
(1.2). We then present a more synthesized approach that makes use of the CDM
halo properties over the evolution history, yielding more detailed and
reliable SIDM predictions. To demonstrate the accuracy and effectiveness, we
will take the CDM halos from the cosmological simulation in ref. [5], use the
parametric model to obtain their corresponding SIDM halos, which will be
further compared to the SIDM halos directly from the N-body simulation.

To facilitate discussion, we rank the simulated halos in ref. [5] by their
masses in descending order. We identify the corresponding CDM counterpart for
each simulated SIDM halo by examining the evolution history. We label these
pairs as “Cosmo-#.” For example, Cosmo-501 represents the halo ranked 501st
from the top in terms of mass.

###  3.1 The basic approach

We propose to use the NFW halo parameters $r_{s}=R_{\rm max}/2.1626$ and
$\rho_{s}=(V_{\rm max}/(1.648r_{s}))^{2}/G$ evaluated for simulated CDM halos
at $z=0$ as input for the functions in eq. (2). To justify this approach, we
show in figure 4 the evolution history of the density profile (left), scale
density (middle) and radius (right) for the Cosmo-501 CDM halo in physical
coordinates. For this halo, its inner NFW profile was established $10~{}{\rm
Gyr}$ ago and has remained largely unchanged since then. In particular,
$\rho_{s}$ and $r_{s}$ are approximately constants after first $3~{}{\rm
Gyr}$.

In general, isolated halos undergo exponential mass growth at high redshifts,
during which their inner density profiles are established. At later stages,
minor mergers primarily contribute mass to the outer regions, resulting in
moderate mass increases, and the inner density profile remains largely
unchanged. Thus we can use the NFW halo parameters at $z=0$ for the functions
in eq. (2) to make SIDM predictions, if we know the halo formation time
$t_{f}$ at which the $\rho_{s}$ and $r_{s}$ become stable.

Figure 4:  Density profiles (left) of an isolated CDM halo at different
lookback times, from the cosmological zoom-in simulation [5]. For
$t_{L}\lesssim 0.2~{}{\rm Gyr}$, the halo density can be described by a single
NFW profile. Evolution of the NFW parameters $\rho_{s}$ (middle) and $r_{s}$
(right) constructed from the simulation [5]. The subscript “phys” denotes that
the relevant quantities are evaluated in physical coordinates.

To estimate $t_{f}$, we use the following relation [72],

$\displaystyle z_{f}=-0.0064(\log_{10}(M_{\rm vir,0}/{\rm
M_{\odot}}))^{2}+0.0237\log_{10}(M_{\rm vir,0}/{\rm M_{\odot}})+1.8837,$ (3.1)

where $z_{f}$ is the redshift of halo formation and $M_{\rm vir,0}$ is the
halo mass at present. The formation redshift is determined through the
condition that at $z_{f}$ the halo mass is $M_{\rm vir,0}/q$, where $q\approx
4.137z^{-0.9476}_{f}$ [72]. The corresponding lookback time is computed as
$t_{L}(z_{f})=\int_{0}^{z_{f}}dz/(H\times(1+z))$, where $H$ is the Hubble rate
as a function of $z$. For the cosmology with $h=0.7$, $\Omega_{m,0}=0.286$,
and $\Omega_{\Lambda_{0}}=0.714$, as in ref. [73], the lookback time can be
evaluated as

$\displaystyle
t_{L}(z)=13.647-11.020\ln\left(\frac{1.5800}{(1+z)^{1.5}}+\sqrt{1+\frac{2.4965}{(1+z)^{3}}}\right)\text{
Gyr}.$ (3.2)

For convenience, we define $t_{f}=13.647~{}{\rm Gyr}-t_{L}(z_{f})$. For the
Cosmo-501 isolated halo, $M_{\rm vir,0}=1.81\times 10^{9}\rm~{}M_{\odot}$, we
have $z_{f}\approx 1.55$ and $t_{L}(z_{f})\approx 9.48~{}{\rm Gyr}$ estimated
using eqs. (3.1) and (3.2), respectively. This is well consistent with the
evolution history shown in figure 4.

In short, we take the following steps for applying the parametric model with
the basic approach. For an isolated CDM halo, we take its NFW scale parameters
at $z=0$, i.e., $\rho_{s,0}$ and $r_{s,0}$ in eq. (2), and calculate the
virial mass and sequentially the formation time $t_{L}(z_{f})$ using eqs.
(3.1) and (3.2), which could also be obtained from halo catalogs if available.
Then, we evaluate the effective cross section eq. (1) and collapse time eq.
(2.2) for a given particle physics model of SIDM, see eq. (1.2). Lastly, we
use the functions in eq. (2) and obtain $\rho_{s}$, $r_{s}$ and $r_{c}$ as a
function of time $t$, and reconstruct the density profile at $t$ based on eq.
(2.1). In practice, we also incorporate a Gaussian scatter of $0.16$ dex in
$\log_{10}(t_{L}(z_{f})/{\rm Gyr})$ [74], when applying to a population of
halos. If $t_{L}(z_{f})/t_{c}>1$, we truncate the halo’s evolution at $t_{c}$
to avoid extrapolation. This truncation only affects the value of the central
density of halos that would be deeply collapsed.

###  3.2 The integral approach

The basic approach is straightforward, but it has limitations when applied to
halos that have undergone significant mass accretion or loss at relatively
late times. Isolated halos, for instance, can experience mass increases due to
minor mergers, or mass loss due to slingshot and splashback events, which
cause spikes and valleys in $V_{\rm max}$ and $R_{\rm max}$. Furthermore,
subhalos evolve in the tidal field and lose mass. Here, we propose an extended
approach that can effectively “integrate” the SIDM effects over the period of
mass growth or reduction.

We assume the mass accretion history in SIDM is similar to that in CDM. This
is a reasonable assumption for most cases, as the self-interactions mainly
modify the inner halo structure; we will discuss later about subhalos where
the assumption is violated. Given the evolution of $V_{\rm max,CDM}$ and
$R_{\rm max,CDM}$, we propose the following integral functions

$\displaystyle V_{\rm max}(t)$ $\displaystyle=$ $\displaystyle V_{\rm
max,CDM}(t_{f})+\int_{t_{f}}^{t}dt^{\prime}\frac{dV_{\rm
max,CDM}(t^{\prime})}{dt^{\prime}}+\int_{t_{f}}^{t}\frac{dt^{\prime}}{t_{c}(t^{\prime})}\frac{dV_{\rm
max,Model}(\tilde{t^{\prime}})}{d\tilde{t^{\prime}}}$ $\displaystyle R_{\rm
max}(t)$ $\displaystyle=$ $\displaystyle R_{\rm
max,CDM}(t_{f})+\int_{t_{f}}^{t}dt^{\prime}\frac{dR_{\rm
max,CDM}(t^{\prime})}{dt^{\prime}}+\int_{t_{f}}^{t}\frac{dt^{\prime}}{t_{c}(t^{\prime})}\frac{dR_{\rm
max,Model}(\tilde{t^{\prime}})}{d\tilde{t^{\prime}}},$ (3.3)

where the first integral on the right-hand side accommodates for the change
due to the evolution events, such as mergers and tidal mass loss, while the
second integral captures the SIDM effects by accounting for variations of the
effective cross section and the collapse timescale, resulting from the change
in the halo properties.

At each moment, the collapse time $t_{c}(t^{\prime})$ is computed using
$V_{\rm max,CDM}(t^{\prime})$ and $R_{\rm max,CDM}(t^{\prime})$. The $dV_{\rm
max,Model}(\tilde{t^{\prime}})/d\tilde{t^{\prime}}$ and $dR_{\rm
max,Model}(\tilde{t^{\prime}})/d\tilde{t^{\prime}}$ terms are derivatives of
the corresponding functions in eq. (2). Once obtaining $V_{\rm max}(t)$,
$R_{\rm max}(t)$ and $\tilde{t}=t/t_{c}(t)$ for time $t$, we further use eq.
(2) to solve for $V_{\rm max,0}$ and $R_{\rm max,0}$, which in turn give rise
to $\rho_{s,0}$ and $r_{s,0}$ assuming an NFW density profile. In this way, we
“construct” an isolated CDM halo at a given moment $t$, to which we can apply
the functions in eq. (2) for obtaining $\rho_{s}$, $r_{s}$ and $r_{c}$. Note
that in the integral approach, $\rho_{s,0}$ and $r_{s,0}$ in eq. (2) are not
necessary their corresponding values found in the simulation at $z=0$. With
the effects of mass changes incorporated, we can apply the model to the halo
earlier than the formation time used in the basic approach. For the integral
approach, we take it to be half of $13.647~{}{\rm Gyr}-t_{L}(z_{f})$, see eq.
(3.2), although the difference is small. In practice, we discretize the
lookback time for a simulated halo into $1000$ intervals with equal spacing
and obtain its $V_{\rm max}$ and $R_{\rm max}$ at all times incrementally.
Then, we replace the integrals in eq. (3.2) with discrete summation.

The proposed integrals in eq. (3.2) can be understood as follows. For a CDM
halo, at each small time interval, the change in the mass due to accretion or
tidal stripping, is reflected in the change of $\rho_{s}$ and $r_{s}$. The
effective cross section eq. (1) with the velocity dispersion $\nu_{\rm eff}$,
as well as the collapse time $t_{c}$, being computed using the obtained
$\rho_{s}$ and $r_{s}$ values, effectively incorporates the SIDM effects
during this time interval. Thus, the integral approach in eq. (3.2) gives rise
to the net change in $V_{\rm max}$ and $R_{\rm max}$ for the corresponding
SIDM during the interval. If the evolution history of CDM halos is known, the
integral approach is easy to implement and it has broader applications than
the basic approach.

We also note that the integral approach becomes the basic one in the limit
that $V_{\rm max,CDM}$ and $R_{\rm max,CDM}$ do not change after the halo
formation time $t_{f}$. From eq. (3.2), we see that in the limit the first
integral vanishes, and the second integral simply reduces to two boundary
terms at $t$ and $t_{f}$ as $t_{c}$ is a constant. The latter term carries a
minus sign and it cancels the first term in the right-hand side of eq. (3.2).

##  4 Modeling the evolution of SIDM isolated halos

We first consider a few representative halos from the simulation [5] that have
evolved into late stages; see table 1 for their characteristic properties.
Among them, Cosmo-26 and Cosmo-32 are core-forming, and both have
$t_{L}(z_{f})/t_{c}\ll 0.2$. All others are deeply collapsed at $z=0$, with
$t_{L}(z_{f})/t_{c}\gtrsim 1$, and they will be truncated at
$t_{L}(z_{f})/t_{c}=1$. Cosmo-796 and Cosmo-1128 have relatively complex
accretion history. The former had a late merger and the latter experienced a
splashback event (i.e., it passed through a larger halo before $z=0$).

Cosmo-ID | $M_{\rm vir,0}$ | $t_{L}(z_{f})$ | $t_{L}(z_{f})/t_{c}$ | $V_{\rm max,SIDM,0}$ | $R_{\rm max,SIDM,0}$ | $V_{\rm max,Model,0}$ | $R_{\rm max,Model,0}$
---|---|---|---|---|---|---|---
Isolated | $\rm(10^{8}\ M_{\odot})$ | $\rm(Gyr)$ | - | $\rm(km/s)$ | $\rm(kpc)$ | $\rm(km/s)$ | $\rm(kpc)$
26 | $186$ | 9.23 | 0.035 | 45.4 | 9.11 | 45.2 | 9.27
32 | $146$ | 9.26 | 0.018 | 40.7 | 10.4 | 40.4 | 12.6
501 | $6.47$ | 9.58 | 1.06 | 21.9 | 0.63 | 22.8 | 0.61
796 | $3.92$ | 9.62 | 1.00 | 21.1 | 0.56 | 21.7 | 0.61
800 | $4.03$ | 9.62 | 1.17 | 21.7 | 0.40 | 21.7 | 0.57
1128 | $2.64$ | 9.65 | 1.27 | 19.6 | 0.60 | 20.9 | 0.53
Table 1:  Characteristic properties of isolated halos at $z=0$, selected from
the cosmological zoom-in SIDM simulation [5] for testing the parametric model.
The virial mass, formation time, $t_{L}(z_{f})/t_{c}$, $V_{\rm max}$ and
$R_{\rm max}$ are reported for the halos. For comparison, their $V_{\rm max}$
and $R_{\rm max}$ values predicted in the parametric model eq. (2) are also
listed in the last two columns.

Figure 5:  Evolution of $V_{\rm max}$ (top) and $R_{\rm max}$ (bottom) for
three isolated halos from the parametric model with the basic (solid-blue) and
integral (solid-green) approaches, as well as the N-body simulation [5]
(solid-magenta). Cosmo-32 is in the core-forming phase at $z=0$, while the
other two are deeply collapsed. Cosmo-796 has a late major merger at
$t_{L}\sim 7~{}{\rm Gyr}$. The subscript “phys” denotes that the relevant
quantities are evaluated in physical coordinates.

Figure 6:  Density profiles of the isolated halos in table 1 from the
parametric model with the basic (solid-blue) and integral (solid-green)
approaches, and the SIDM simulation [5] (solid-magenta), as well as the CDM
simulation (dashed-black) for comparison.

In figure 5, we show the $V_{\rm max}$ (top) and $R_{\rm max}$ (bottom)
evolution of three halos: Cosmo-32, Cosmo-501, and Cosmo-796. Both basic
(solid-blue) and integral (solid-green) approaches agree with the SIDM
simulation (solid-magenta) within $\sim 10\%$. Compared to their CDM
counterparts (dashed-black), Cosmo-501 and Cosmo-796 have larger $V_{\rm max}$
and smaller $R_{\rm max}$ at later stages due to core collapse. Although
Cosmo-796 has a significant merger at $t_{L}\sim 7~{}{\rm Gyr}$, both
approaches make similar predictions, and they agree with the simulation. For
Cosmo-32, the evolution is almost identical in SIDM and CDM. This halo is in
the core-forming phase over the entire evolution history, as it has a higher
mass and a smaller effective cross section accordingly. As a good
approximation, dark matter self-interactions do not change $V_{\rm max}$ and
$R_{\rm max}$ ($z=0$) for core-forming isolated SIDM halos [5].

In figure 6, we show density profiles of all isolated halos at $z=0$ listed in
table 1. We again see that both basic (solid-blue) and integral (solid-green)
approaches agree well with the SIDM simulation (solid-magenta). Cosmo-26 and
Cosmo-32 are in the core-forming phase at $z=0$ and they have a cored density
profile. The four others are deeply collapsed, and their central densities are
higher compared to CDM (dashed-black).

We have tested both approaches with a sample of $647$ isolated CDM halos
resolved in the simulation [5]. As discussed previously, the integral approach
requires well-resolved evolution history of a CDM halo, but numerical noises
could be large for low-mass halos in the sample that are close to the
resolution limit. In addition, the basic and integral approaches agree well
for the representative examples shown in figures 5 and 6. Thus, we show the
comparison between N-body simulation and basic approach for the whole sample
of simulated halos.

Figure 7 (top) shows $R_{\rm max}\textup{--}V_{\rm max}$ and $V_{\rm
circ}(r_{\rm fid})\textup{--}V_{\rm max}$ distributions predicted using the
basic approach, where $r_{\rm fid}=2V_{\rm max}/(70~{}{\rm km/s})\ {\rm kpc}$
and each halo is color coded according to $t_{L}(z_{f})/t_{c}$. We see that
the predicted distributions agree well with those directly from the N-body
SIDM simulation [5] as shown in figure 7 (bottom), where the simulated SIDM
halo is colored based on its effective concentration as in [5]. In particular,
our parametric model successfully reproduces halos that are deeply collapsed.
These halos have the lowest $R_{\rm max}$ and highest $V_{\rm circ}(r_{\rm
fid})$ values for given $V_{\rm max}$.

Figure 7:  The $R_{\rm max}\textup{--}V_{\rm max}$ (left) and $V_{\rm
circ}(r_{\rm fid})\textup{--}V_{\rm max}$ (right) distributions of isolated
SIDM halos at $z=0$ predicted using the parametric model with the basic
approach (top), where each halo is color coded according to
$t_{L}(z_{f})/t_{c}$. For comparison, the corresponding distributions from the
N-body SIDM simulation are also shown (bottom), and the color coding is based
on the effective concentration $c_{\rm eff}=R_{\rm vir}/(R_{\rm max}/2.1626)$;
taken from [5]. As in [5], the median CDM isolated halo $V_{\rm
max}\textup{--}R_{\rm max}$ relation (solid-black) together with a $\pm 0.6$
dex band (dashed-black), and resolution limit (gray-shaded) are shown in the
left panel; the $1\textup{:}1$ relation (dashed-black) in the right panel.

Figure 8:  The $\log_{10}({R_{\rm max}}/{R_{\rm max,median}})$ (left), $V_{\rm max}$ (middle), and $(V_{\rm max}-V_{\rm circ}(r_{\rm fid}))$ (right) distributions for the sample of isolated halos from the parametric model (solid-blue), SIDM (solid-magenta) and CDM (dashed-black) simulations [5]. Properties | CDM simulation | SIDM simulation | Model prediction
---|---|---|---
$\mu\left(\log_{10}\frac{R_{\rm max}}{R_{\rm max,median}}\right)$ | 0.0144 | -0.0763 | -0.0918
$\sigma\left(\log_{10}\frac{R_{\rm max}}{R_{\rm max,median}}\right)$ | 0.173 | 0.224 | 0.272
$\mu\left(V_{\rm max}\right)$ $\rm(km/s)$ | 15.9 | 16.9 | 16.7
$\sigma\left(V_{\rm max}\right)$ $\rm(km/s)$ | 7.40 | 7.75 | 7.46
$\mu\left(V_{\rm max}-V_{\rm circ}(r_{\rm fid})\right)$ $\rm(km/s)$ | 3.00 | 5.40 | 4.16
$\sigma\left(V_{\rm max}-V_{\rm circ}(r_{\rm fid})\right)$ $\rm(km/s)$ | 2.32 | 3.74 | 3.90
Table 2:  The mean and standard deviation of the $\log_{10}({R_{\rm
max}}/{R_{\rm max,median}})$, $V_{\rm max}$, and $(V_{\rm max}-V_{\rm
circ}(r_{\rm fid}))$ distributions.

Figure 8 further shows the distribution of $\log_{10}(R_{\rm max}/R_{\rm
max,median})$ (left), $V_{\rm max}$ (middle), and $(V_{\rm max}-V_{\rm
circ}(r_{\rm fid}))$ (right) for the sample, from the parametric model (solid-
blue), as well as the N-body SIDM (solid-magenta) and CDM (dashed-black)
simulations [5]. We see good agreement between parametric model and
simulation. In addition, compared to CDM, SIDM predicts a higher population
for $\log_{10}(R_{\rm max}/R_{\rm max,median})\lesssim-0.5$ due to core
collapse, as well as a higher population for $(V_{\rm max}-V_{\rm circ}(r_{\rm
fid}))>5~{}{\rm km/s}$ because of core formation. For a more quantitative
comparison, we show the mean and standard deviation of the three distributions
in table 2.

##  5 Modeling the evolution of SIDM subhalos

Subhalos can be significantly affected by the host halo due to tidal
interactions. As a result, their density profiles are different from those of
isolated halos. To account for the tidal effects, we modify the density
profile in eq. (2.1) as

$\displaystyle\rho_{\rm
tSIDM}(r)=\frac{\rho_{s}}{\frac{(r^{\beta}+r_{c}^{\beta})^{1/{\beta}}}{r_{s}}\left(1+\frac{r}{r_{s}}\right)^{2}\left(1+\left(\frac{r}{r_{t}}\right)^{2-u}\right)^{1+3u}},$
(5.1)

where the additional parameters are the tidal radius $r_{t}$ and the numerical
factor $u$ controlling the logarithmic slope of the truncated density profile
in the outer region. We evaluate the tidal radius $r_{t}$ as [75, 76]

$\displaystyle r_{t}=d\left(\frac{M_{\rm sub}}{M_{\rm
host}(r<d)}\right)^{1/3},$ (5.2)

where $d$ is the distance of the subhalo from its host halo, $M_{\rm sub}$ is
the subhalo virial mass and $M_{\rm host}(r<d)$ is the host halo mass enclosed
within $d$. We further take into account the evolution of the tidal radius in
our study. The smooth truncation is controlled by the term
$1/(1+(r/r_{t})^{2-u})^{1+3u}$, which falls off at large radii as $r^{-2}$ at
$u=0$ and as $r^{-4}$ at $u=1$, as in the case of a truncated NFW profile for
CDM subhalos, see, e.g., [77, 18]. In this work, we take $u(c_{\rm eff})={\rm
Min}(1,0.0004c_{\rm eff}^{2.2})$, where $c_{\rm eff}=R_{\rm vir}/(R_{\rm
max}/2.1626)$ is the effective concentration [5], and the corresponding
truncated profile in eq. (5.1) can fit both core-forming and -collapsing
subhalos. We have also checked that the SIDM evaporation effect [78] is
negligible for the subhalos in the simulation [5], although such an effect
could be incorporated using a semi-analytical method [79].

Since subhalos undergo tidal mass loss after falling into the host halo, the
basic approach is insufficient for modeling their evolution. However, the
integral approach based on eq. (3.2) still works. As a demonstration, we
select $6$ representative cases, ranging from core-forming subhalos to deeply
collapsing ones, as summarized in table 3. The evolution stage is quantified
by $\tilde{t}_{0}=\int_{t_{f}}^{t_{0}}dt/t_{c}(t)$, where $t_{0}=13.647~{}{\rm
Gyr}$ and we take the formation time $t_{f}$ to be half of $13.647~{}{\rm
Gyr}-t_{L}(z_{f})$, see eq. (3.2). In addition, the collapse time $t_{c}$
changes along with tidal mass loss, and we use eq. (2.2) to evaluate it with
$\rho_{\rm eff}=(V_{\rm max,CDM}/1.648r_{\rm eff})^{2}/G$ and $r_{\rm
eff}=(R_{\rm max,CDM}/2.1626)$ for a given moment.

Cosmo-ID | $M_{\rm vir,0}$ | $\tilde{t}_{0}$ | $V_{\rm max,SIDM,0}$ | $R_{\rm max,SIDM,0}$ | $V_{\rm max,Model,0}$ | $R_{\rm max,Model,0}$
---|---|---|---|---|---|---
Subhalos | $\rm(10^{8}\ M_{\odot})$ | - | $\rm(km/s)$ | $\rm(kpc)$ | $\rm(km/s)$ | $\rm(kpc)$
71 | 57.1 | 0.092 | 44.6 | 4.93 | 43.0 | 5.60
186 | 20.7 | 0.27 | 34.4 | 2.46 | 31.4 | 2.00
567 | 6.04 | 0.85 | 31.7 | 0.58 | 30.3 | 0.82
665 | 4.61 | 0.12 | 13.7 | 3.72 | 13.6 | 3.43
945 | 3.92 | 1.10 | 38.8 | 0.24 | 38.3 | 0.045
1357 | 2.38 | 1.18 | 26.0 | 0.37 | 24.6 | 0.42
Table 3:  Characteristic properties of subhalos at $z=0$, selected from the
cosmological zoom-in SIDM simulation [5] for testing the parametric model. The
virial mass, $\tilde{t}_{0}=\int_{t_{f}}^{t_{0}}dt/t_{c}(t)$, $V_{\rm max}$
and $R_{\rm max}$ are reported for the subhalos. For comparison, their $V_{\rm
max}$ and $R_{\rm max}$ values predicted in the parametric model with the
integral approach are also listed in the last two columns; see the relevant
discussion after eq. (3.2).

Figure 9:  The $V_{\rm max}$ (top) and $R_{\rm max}$ (bottom) evolution of
three SIDM subhalos that are in the deeply collapse phase at $z=0$ from the
parametric model with the integral approach (solid-green) and the SIDM
simulation [5] (solid-magenta), as well as their CDM counterparts (dashed-
black).

Figure 9 shows the $V_{\rm max}$ and $R_{\rm max}$ evolution of three SIDM
subhalos that are in the collapse phase at $z=0$ from the parametric model
(solid-green) and the N-body simulation [5] (solid-magenta), as well as their
CDM counterparts (dashed-black). We see that the model well reproduces the
evolution history of $V_{\rm max}$ and $R_{\rm max}$ of the simulated SIDM
subhalos, and the agreement is within $10\%$. For the SIDM subhalos, $V_{\rm
max}$ increases, while $R_{\rm max}$ decreases continuously for $t_{L}\gtrsim
6~{}{\rm Gyr}$, a significant deviation from CDM. In figure 11, we further
show the density profiles for all subhalos listed in table 3 at $z=0$. We
again see the agreement between parametric model (solid-green) and simulation
(solid-magenta) predictions for both core-forming and -collapsing SIDM
subhalos. Their density profiles different from the corresponding CDM ones
(dashed-black).

Figure 10:  The $\sigma_{\rm eff}/m$, $t_{c}(t)$, and $t/t_{c}(t)$ evolution
of the SIDM subhalos in table 3, where $t$ is the time duration the halo has
evolved since its formation.

As discussed in Sec. 3.2, $\sigma_{\rm eff}/m$ and $t_{c}$ evolve with the
mass change, which is captured in the integral approach. Figure 10 shows the
$\sigma_{\rm eff}/m$, $t_{c}$ and $t/t_{c}$ evolution of the $6$ subhalos
listed in table 3. Apart from fluctuations, both $\sigma_{\rm eff}/m$ and
$t_{c}(t)$ changed moderately on average, within an order of magnitude. The
ratio $t/t_{c}$ has large fluctuations but gradually increases since the halo
formation.

Figure 11:  Density profiles of core-collapsing (top) and -forming (bottom)
SIDM subhalos at $z=0$ from the parametric model with the integral approach
(solid-green) and the SIDM simulation [5] (solid-magenta), as well as their
CDM counterparts [5] (dashed-black).

Figure 12:  The $V_{\rm max}\textup{--}R_{\rm max}$ distribution of subhalos
predicted using the parametric model with the integral approach (left), where
each subhalo is colored according to its $t_{L}(z_{f})/t_{c}$ value. For
comparison, the distributions from the SIDM (middle) and CDM (right)
simulations are shown; taken from ref. [5]. As in [5], the median CDM subhalo
$V_{\rm max}\textup{--}R_{\rm max}$ relation (solid-black) together with a
$\pm 0.6$ dex band (dashed-black), and resolution limit (gray-shaded) are
shown.

We have applied the parametric model to all $125$ resolved CDM subhalos in the
cosmological simulation of a Milky Way analog in ref. [5] and obtain their
SIDM predictions in the $R_{\rm max}\textup{--}V_{\rm max}$ plane. Since some
of the low-mass simulated subhalos that undergo strong stripping may suffer
from numerical artifacts, we set $V_{\rm max}=2~{}{\rm km/s}$ and $R_{\rm
max}=0.1~{}{\rm kpc}$, if their values drop below the corresponding ones when
performing the discrete summation in eq. (3.2). With this condition, we avoid
producing unphysical results while preserving the numerical error in the same
order as the numerical noise induced by the resolution limit.

Figure 12 shows the $R_{\rm max}\textup{--}V_{\rm max}$ distributions of the
SIDM subhalos predicted using the parametric model (left) and the N-body
simulation [5] (middle), compared to the simulated CDM subhalos [5] (right).
We see that our model successfully reproduces the main trends of the simulated
SIDM subhalos, especially the spread towards the lower-right region in the
$R_{\rm max}\textup{--}V_{\rm max}$ plane, corresponding to a population of
collapsed subhalos. In the upper-left region, the model prediction misses a
few subhalos with large $R_{\rm max}$ seen in the SIDM simulation. For these
subhalos, the tidal mass loss is amplified due to SIDM core formation and they
become more diffuse, compared to CDM. Our current model does not capture this
effect, as we have assumed that the tidal mass loss rate of an SIDM subhalo is
similar to that of its CDM counterpart. The assumption could be violated in
some extreme cases, see, e.g., [10, 44, 80, 5, 81]. A dedicated study on this
aspect is beyond the scope of this work, and we leave it for future work.

##  6 Exploring SIDM parameter space using the parametric model

In Sec. 4, we applied the parametric model to the isolated CDM halos in ref.
[5] and obtained their SIDM counterparts, whose properties well agree with
those from the cosmological N-body simulation. The simulation in ref. [5]
assumes $\sigma_{0}/m=147.1~{}{\rm cm^{2}/g}$ and $w=24.33~{}{\rm km/s}$.
Since our model is flexible, it becomes possible to efficiently sample large
parameter space of $\sigma_{0}/m$ and $w$, which can be in turn mapped to
fundamental particle physics parameters, such as mediator mass and coupling
constant; see ref. [60] for details. As a demonstration, we take the $647$
isolated CDM halos in [5] as input and apply the parametric model, while
varying $\sigma_{0}/m$ and $w$ over a wide range.

We first consider four representative constant cross sections to gain
insights: $\sigma/m=3~{}\rm cm^{2}/g$, $10~{}\rm cm^{2}/g$, $50~{}\rm
cm^{2}/g$, and $100~{}\rm cm^{2}/g$. In figure 13, we show the $V_{\rm
circ}(r_{\rm fid})\textup{--}V_{\rm max}$ distribution predicted using the
parametric model with the basic approach (circle), where each halo is colored
according to its $t_{L}(z_{f})/t_{c}$ value. For comparison, we also include
observed galaxies with data compiled in [67] (cross). For $\sigma/m=3~{}\rm
cm^{2}/g$, all halos are in the core-forming phase, $V_{\rm circ}(r_{\rm
fid})$ is systematically lower and its spread for fixed $V_{\rm max}$ becomes
larger, compared to CDM [5]. As the cross section increases, the number of
core-collapsing halos increases accordingly. When
$\sigma/m=50\textup{--}100~{}{\rm cm^{2}/g}$, a significant number of halos
are in the collapse phase and they have a high inner density, close to the
$1\textup{:}1$ line in the $V_{\rm circ}(r_{\rm fid})\textup{--}V_{\rm max}$
plane. Meanwhile, the most diffused halos are still in the core-forming phase
and remain far off the $1\textup{:}1$ line, even for $\sigma/m=100~{}\rm
cm^{2}/g$. These halos could potentially host ultra-diffuse galaxies in the
field [82].

The deviation of the $V_{\rm circ}(r_{\rm fid})\textup{--}V_{\rm max}$
distribution from the $1\textup{:}1$ line provides a useful measure to
characterize the SIDM effects. We fit the halos with a relation of $V_{\rm
circ}(r_{\rm fid})=c~{}V_{\rm max}$ for the range of $15~{}{\rm km/s}<V_{\rm
max}<50~{}{\rm km/s}$, consistent with that for observed galaxies shown in
figure 13, and minimize the absolute mean distance from the $1\textup{:}1$
line to determine the parameter $c$; see figure 13 (solid-black line). Our fit
tends to split the selected halos into two equal halves on two sides of the
fitted line. For the CDM counterparts, $c=0.82$. As the cross section
increases, the $c$ value decreases first as the core size increases for most
halos, but it increases again as more and more halos are collapsed. For
$\sigma/m=100~{}{\rm cm^{2}/g}$, $c=0.96$. As a reference, $c=0.60$ for the
observed galaxies shown in figure 13. We expect that there is significant
coverage bias for the sample of observed galaxies, and hence it is challenging
to make a quantitative comparison between observations and SIDM predictions
based on the $c$ value. Nevertheless, it is a useful indicator for the overall
significance of dark matter self-interactions over a population of halos.

Figure 13:  The $V_{\rm circ}(r_{\rm fid})\textup{--}V_{\rm max}$ distribution
of isolated halos predicted using the parametric model with the basic approach
for the constant cross sections $\sigma/m=3~{}{\rm cm^{2}/g}$, $10~{}{\rm
cm^{2}/g}$, $50~{}{\rm cm^{2}/g}$ and $100~{}{\rm cm^{2}/g}$ (circle). Each
halo is colored according to $t_{L}(z_{f})/t_{c}$. The $V_{\rm circ}(r_{\rm
fid})\textup{--}V_{\rm max}$ relation from fitting to the halos in the range
$15~{}{\rm km/s}<V_{\rm max}<50~{}\rm km/s$ (solid-black), the $1\textup{:}1$
relation (dotted-black), and the data points from observed galaxies compiled
in [67] (red cross) are shown for comparison.  Figure 14:  Contours of the
fitted $c$ value in the $w\textup{--}\sigma_{0}/m$ plane, where $\sigma_{0}/m$
is the overall normalization factor of the self-scattering cross section and
$w$ controls the velocity dependence; see eq. (1.2).

In figure 14, we show contours of the fitted $c$ value based on a scan with a
grid of $100\times 100$ points in the $w\textup{--}\sigma_{0}/m$ plane, where
$\sigma_{0}/m$ is the overall normalization factor of the cross section and
$w$ controls the velocity dependence. We see the distribution of $c$ has
nontrivial features, including a valley (darker green) and a hill (darker
pink). In the lower-left region, the SIDM effect is small and the $c$ value
reduces to the CDM case around $0.82$. For $w\sim{\cal O}(10)~{}\rm km/s$, the
increase of $\sigma_{0}/m$ first decreases $c$, enlarging the spread of halos
in the $V_{\rm circ}(r_{\rm fid})\textup{--}V_{\rm max}$ from the
$1\textup{:}1$ line. When $\sigma_{0}/m$ further increases, $c$ increases as
well because more and more halos enter the collapse phase. These nontrivial
features could help us distinguish different particle physics models of SIDM.

##  7 Conclusions

In this work, we have proposed a parametric model to transfer CDM (sub)halos
into their SIDM counterparts. The model is based on an analytical universal
density profile, which can accurately describe the dark matter distribution of
halos over the course of the evolution. We calibrated the density profile
using controlled N-body simulations of an isolated halo and obtained a set of
equations capturing the time evolution of the parameters of the density
profile. We then introduced two ways of applying the model. The basic approach
uses the scale radius and density of an isolated CDM halo ($z=0$) as input and
reconstructs the halo’s evolution history in the presence of dark matter self-
interactions. The integral approach further extends the basic one to
incorporate the mass change, and integrates the SIDM effects along the
evolution history. Thus it can be used to model gravothermal evolution of
subhalos, which suffer from tidal mass loss.

We have tested and validated the parametric model using zoom-in cosmological
SIDM simulations of a Milky Way analog for both isolated halos and subhalos.
We further applied it to map out the relevant parameter space for a particle
physics model of SIDM. In the future, we could apply the parametric model to
state-of-the-art cosmological CDM simulations, e.g., FIRE2 [83, 84] and
IllustrisTNG [85, 86], and transform existing CDM (sub)halo properties into
SIDM predictions. We could integrate the parametric model together with
existing semi-analytical models with halo merger trees [87, 88, 89] and
subhalo evolution trajectories [90]. The parametric model is flexible and it
provides an efficient tool for identifying favored SIDM parameter regions in
light of latest observations of galactic systems, such as satellite galaxies
of the Milky Way and ultra-diffuse galaxies.

We could further improve the model. While this work focuses on the dark
matter-only case, it is crucial to include baryonic effects on gravothermal
evolution in order to make better connections with observations of galaxies
whose baryon concentration is high. For the integral approach, extra terms
could be included to incorporate SIDM effects that couple to mass changes. In
addition, we use the collapse timescale as a normalization factor to derive
the universal density profile. It would be interesting to explore if there is
a more accurate way to estimate the timescale for halos in cosmological
environments. We will leave these interesting topics for future investigation.

## Acknowledgments

We thank Manoj Kaplinghat for helpful discussion. D.Y. and H.-B. Y were
supported by the John Templeton Foundation under grant ID #61884 and the U.S.
Department of Energy under grant No. de-sc0008541. Y.Z. was supported by the
Kavli Institute for Cosmological Physics at the University of Chicago through
an endowment from the Kavli Foundation and its founder Fred Kavli. The
opinions expressed in this publication are those of the authors and do not
necessarily reflect the views of the John Templeton Foundation.

## References

  * Tulin and Yu [2018] S. Tulin and H.-B. Yu, Physics Reports 730, 1 (2018), 1705.02358.
  * Adhikari et al. [2022] S. Adhikari, A. Banerjee, K. K. Boddy, F.-Y. Cyr-Racine, H. Desmond, C. Dvorkin, B. Jain, F. Kahlhoefer, M. Kaplinghat, A. Nierenberg, et al., arXiv e-prints arXiv:2207.10638 (2022), 2207.10638.
  * Spergel and Steinhardt [2000] D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84, 3760 (2000), astro-ph/9909386.
  * Correa et al. [2022] C. A. Correa, M. Schaller, S. Ploeckinger, N. Anau Montel, C. Weniger, and S. Ando (2022), 2206.11298.
  * Yang et al. [2022] D. Yang, E. O. Nadler, and H.-B. Yu (2022), 2211.13768.
  * Balberg et al. [2002] S. Balberg, S. L. Shapiro, and S. Inagaki, APJ 568, 475 (2002), astro-ph/0110561.
  * Balberg et al. [2002] S. Balberg, S. L. Shapiro, and S. Inagaki, Astrophys. J. 568, 475 (2002), astro-ph/0110561.
  * Koda and Shapiro [2011] J. Koda and P. R. Shapiro, MNRAS 415, 1125 (2011), 1101.3097.
  * Nishikawa et al. [2020] H. Nishikawa, K. K. Boddy, and M. Kaplinghat, Phys. Rev. D 101, 063009 (2020), 1901.00499.
  * Sameie et al. [2019] O. Sameie, H.-B. Yu, L. V. Sales, M. Vogelsberger, and J. Zavala (2019), 1904.07872.
  * Zavala et al. [2019] J. Zavala, M. R. Lovell, M. Vogelsberger, and J. D. Burger, Phys. Rev. D 100, 063007 (2019), 1904.09998.
  * Kaplinghat et al. [2019] M. Kaplinghat, M. Valli, and H.-B. Yu, Mon. Not. Roy. Astron. Soc. 490, 231 (2019), 1904.04939.
  * Kahlhoefer et al. [2019] F. Kahlhoefer, M. Kaplinghat, T. R. Slatyer, and C.-L. Wu, JCAP 12, 010 (2019), 1904.10539.
  * Correa [2021] C. A. Correa, Mon. Not. Roy. Astron. Soc. 503, 920 (2021), 2007.02958.
  * Turner et al. [2021] H. C. Turner, M. R. Lovell, J. Zavala, and M. Vogelsberger, Mon. Not. Roy. Astron. Soc. 505, 5327 (2021), 2010.02924.
  * Slone et al. [2023] O. Slone, F. Jiang, M. Lisanti, and M. Kaplinghat, Phys. Rev. D 107, 043014 (2023), 2108.03243.
  * Silverman et al. [2022] M. Silverman, J. S. Bullock, M. Kaplinghat, V. H. Robles, and M. Valli, Mon. Not. Roy. Astron. Soc. 518, 2418 (2022), 2203.10104.
  * Minor et al. [2021] Q. E. Minor, S. Gad-Nasr, M. Kaplinghat, and S. Vegetti, Mon. Not. Roy. Astron. Soc. 507, 1662 (2021), 2011.10627.
  * Yang and Yu [2021] D. Yang and H.-B. Yu, Phys. Rev. D 104, 103031 (2021), 2102.02375.
  * Gilman et al. [2021] D. Gilman, J. Bovy, T. Treu, A. Nierenberg, S. Birrer, A. Benson, and O. Sameie, Mon. Not. Roy. Astron. Soc. 507, 2432 (2021), 2105.05259.
  * Gilman et al. [2023] D. Gilman, Y.-M. Zhong, and J. Bovy, Phys. Rev. D 107, 103008 (2023), 2207.13111.
  * [22] E. O. Nadler, D. Yang, and H.-B. Yu.
  * Pollack et al. [2015] J. Pollack, D. N. Spergel, and P. J. Steinhardt, Astrophys. J. 804, 131 (2015), 1501.00017.
  * Choquette et al. [2019] J. Choquette, J. M. Cline, and J. M. Cornell, JCAP 07, 036 (2019), 1812.05088.
  * Feng et al. [2021] W.-X. Feng, H.-B. Yu, and Y.-M. Zhong, Astrophys. J. Lett. 914, L26 (2021), 2010.15132.
  * Feng et al. [2022] W.-X. Feng, H.-B. Yu, and Y.-M. Zhong, JCAP 2022, 036 (2022), 2108.11967.
  * Xiao et al. [2021] H. Xiao, X. Shen, P. F. Hopkins, and K. M. Zurek, JCAP 07, 039 (2021), 2103.13407.
  * Meshveliani et al. [2023] T. Meshveliani, J. Zavala, and M. R. Lovell, Phys. Rev. D 107, 083010 (2023), 2210.01817.
  * Vogelsberger et al. [2012] M. Vogelsberger, J. Zavala, and A. Loeb, Mon. Not. Roy. Astron. Soc. 423, 3740 (2012), 1201.5892.
  * Rocha et al. [2013] M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Onorbe, and L. A. Moustakas, Mon. Not. Roy. Astron. Soc. 430, 81 (2013), 1208.3025.
  * Peter et al. [2013] A. H. G. Peter, M. Rocha, J. S. Bullock, and M. Kaplinghat, MNRAS 430, 105 (2013), 1208.3026.
  * Zavala et al. [2013] J. Zavala, M. Vogelsberger, and M. G. Walker, Mon. Not. Roy. Astron. Soc. 431, L20 (2013), 1211.6426.
  * Vogelsberger et al. [2016] M. Vogelsberger, J. Zavala, F.-Y. Cyr-Racine, C. Pfrommer, T. Bringmann, and K. Sigurdson, Mon. Not. Roy. Astron. Soc. 460, 1399 (2016), 1512.05349.
  * Fitts et al. [2019] A. Fitts et al., Mon. Not. Roy. Astron. Soc. 490, 962 (2019), 1811.11791.
  * Robles et al. [2019] V. H. Robles, T. Kelley, J. S. Bullock, and M. Kaplinghat, Mon. Not. Roy. Astron. Soc. 490, 2117 (2019), 1903.01469.
  * Banerjee et al. [2019] A. Banerjee, S. Adhikari, N. Dalal, S. More, and A. Kravtsov (2019), 1906.12026.
  * Sameie et al. [2021] O. Sameie, M. Boylan-Kolchin, R. Sanderson, D. Vargya, P. F. Hopkins, A. Wetzel, J. Bullock, A. Graus, and V. H. Robles, Mon. Not. Roy. Astron. Soc. 507, 720 (2021), 2102.12480.
  * Nadler et al. [2020] E. O. Nadler, A. Banerjee, S. Adhikari, Y.-Y. Mao, and R. H. Wechsler (2020), 2001.08754.
  * Ebisu et al. [2022] T. Ebisu, T. Ishiyama, and K. Hayashi, Phys. Rev. D 105, 023016 (2022), 2107.05967.
  * Elbert et al. [2018] O. D. Elbert, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, A. S. Graus, and M. Rocha, Astrophys. J. 853, 109 (2018), 1609.08626.
  * Creasey et al. [2017] P. Creasey, O. Sameie, L. V. Sales, H.-B. Yu, M. Vogelsberger, and J. Zavala, Mon. Not. Roy. Astron. Soc. 468, 2283 (2017), 1612.03903.
  * Sameie et al. [2018] O. Sameie, P. Creasey, H.-B. Yu, L. V. Sales, M. Vogelsberger, and J. Zavala, Mon. Not. Roy. Astron. Soc. 479, 359 (2018), 1801.09682.
  * Huo et al. [2020] R. Huo, H.-B. Yu, and Y.-M. Zhong, JCAP 06, 051 (2020), 1912.06757.
  * Yang et al. [2020] D. Yang, H.-B. Yu, and H. An, Phys. Rev. Lett. 125, 111105 (2020), 2002.02102.
  * Zeng et al. [2022] Z. C. Zeng, A. H. G. Peter, X. Du, A. Benson, S. Kim, F. Jiang, F.-Y. Cyr-Racine, and M. Vogelsberger, Mon. Not. Roy. Astron. Soc. 513, 4845 (2022), 2110.00259.
  * Burger et al. [2022] J. D. Burger, J. Zavala, L. V. Sales, M. Vogelsberger, F. Marinacci, and P. Torrey, Mon. Not. Roy. Astron. Soc. 513, 3458 (2022), 2108.07358.
  * Essig et al. [2019] R. Essig, S. D. Mcdermott, H.-B. Yu, and Y.-M. Zhong, Phys. Rev. Lett. 123, 121102 (2019), 1809.01144.
  * Yang et al. [2023a] S. Yang, X. Du, Z. C. Zeng, A. Benson, F. Jiang, E. O. Nadler, and A. H. G. Peter, Astrophys. J. 946, 47 (2023a), 2205.02957.
  * Yang et al. [2023b] S. Yang, F. Jiang, A. Benson, Y.-M. Zhong, C. Mace, X. Du, Z. C. Zeng, and A. H. G. Peter (2023b), 2305.05067.
  * [50] Y.-M. Zhong, D. Yang, and H.-B. Yu.
  * Feng et al. [2010] J. L. Feng, M. Kaplinghat, and H.-B. Yu, Phys. Rev. Lett. 104, 151301 (2010), 0911.0422.
  * Feng et al. [2009] J. L. Feng, M. Kaplinghat, H. Tu, and H.-B. Yu, JCAP 2009, 004 (2009), 0905.3039.
  * Tulin et al. [2013] S. Tulin, H.-B. Yu, and K. M. Zurek, Phys. Rev. D 87, 115007 (2013), 1302.3898.
  * Kahlhoefer et al. [2014] F. Kahlhoefer, K. Schmidt-Hoberg, M. T. Frandsen, and S. Sarkar, Mon. Not. Roy. Astron. Soc. 437, 2865 (2014), 1308.3419.
  * Agrawal et al. [2017] P. Agrawal, F.-Y. Cyr-Racine, L. Randall, and J. Scholtz, JCAP 2017, 022 (2017), 1610.04611.
  * Robertson et al. [2017a] A. Robertson, R. Massey, and V. Eke, Mon. Not. Roy. Astron. Soc. 467, 4719 (2017a), 1612.03906.
  * Kahlhoefer et al. [2017] F. Kahlhoefer, K. Schmidt-Hoberg, and S. Wild, JCAP 08, 003 (2017), 1704.02149.
  * Kummer et al. [2018] J. Kummer, F. Kahlhoefer, and K. Schmidt-Hoberg, Mon. Not. Roy. Astron. Soc. 474, 388 (2018), 1706.04794.
  * Colquhoun et al. [2021] B. Colquhoun, S. Heeba, F. Kahlhoefer, L. Sagunski, and S. Tulin, Phys. Rev. D 103, 035006 (2021), 2011.04679.
  * Yang and Yu [2022] D. Yang and H.-B. Yu, JCAP 2022, 077 (2022), 2205.03392.
  * Girmohanta and Shrock [2022] S. Girmohanta and R. Shrock, Phys. Rev. D 106, 063013 (2022), 2206.14395.
  * Kaplinghat et al. [2016] M. Kaplinghat, S. Tulin, and H.-B. Yu, Phys. Rev. Lett. 116, 041302 (2016), 1508.03339.
  * Robertson et al. [2021] A. Robertson, R. Massey, V. Eke, J. Schaye, and T. Theuns, Mon. Not. Roy. Astron. Soc. 501, 4610 (2021), 2009.07844.
  * Jiang et al. [2023] F. Jiang et al., Mon. Not. Roy. Astron. Soc. 521, 4630 (2023), 2206.12425.
  * Kamada et al. [2017] A. Kamada, M. Kaplinghat, A. B. Pace, and H.-B. Yu, Phys. Rev. Lett. 119, 111102 (2017), 1611.02716.
  * Ren et al. [2019] T. Ren, A. Kwa, M. Kaplinghat, and H.-B. Yu, Physical Review X 9, 031020 (2019), 1808.05695.
  * Santos-Santos et al. [2020] I. M. E. Santos-Santos, J. F. Navarro, A. Robertson, A. Benítez-Llambay, K. A. Oman, M. R. Lovell, C. S. Frenk, A. D. Ludlow, A. Fattahi, and A. Ritz, MNRAS 495, 58 (2020), 1911.09116.
  * Zentner et al. [2022] A. Zentner, S. Dandavate, O. Slone, and M. Lisanti, JCAP 07, 031 (2022), 2202.00012.
  * Outmezguine et al. [2022] N. J. Outmezguine, K. K. Boddy, S. Gad-Nasr, M. Kaplinghat, and L. Sagunski (2022), 2204.06568.
  * Springel [2005] V. Springel, MNRAS 364, 1105 (2005), astro-ph/0505010.
  * Robertson et al. [2017b] A. Robertson, R. Massey, and V. Eke, Mon. Not. Roy. Astron. Soc. 465, 569 (2017b), 1605.04307.
  * Correa et al. [2015] C. A. Correa, J. S. B. Wyithe, J. Schaye, and A. R. Duffy, Mon. Not. Roy. Astron. Soc. 450, 1514 (2015), 1409.5228.
  * Hinshaw et al. [2013] G. Hinshaw, D. Larson, E. Komatsu, D. N. Spergel, C. L. Bennett, J. Dunkley, M. R. Nolta, M. Halpern, R. S. Hill, N. Odegard, et al., APJS 208, 19 (2013), 1212.5226.
  * Giocoli et al. [2012] C. Giocoli, G. Tormen, and R. K. Sheth, Mon. Not. Roy. Astron. Soc. 422, 185 (2012), 1111.6977.
  * van den Bosch et al. [2018] F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, Mon. Not. Roy. Astron. Soc. 474, 3043 (2018), 1711.05276.
  * Klypin et al. [1999] A. A. Klypin, S. Gottlober, and A. V. Kravtsov, Astrophys. J. 516, 530 (1999), astro-ph/9708191.
  * Baltz et al. [2009] E. A. Baltz, P. Marshall, and M. Oguri, JCAP 01, 015 (2009), 0705.0682.
  * Nadler et al. [2020] E. O. Nadler, A. Banerjee, S. Adhikari, Y.-Y. Mao, and R. H. Wechsler, APJ 896, 112 (2020), 2001.08754.
  * Shirasaki et al. [2022] M. Shirasaki, T. Okamoto, and S. Ando, MNRAS 516, 4594 (2022), 2205.09920.
  * Moreno et al. [2022] J. Moreno, S. Danieli, J. S. Bullock, R. Feldmann, P. F. Hopkins, O. çatmabacak, A. Gurvich, A. Lazar, C. Klein, C. B. Hummels, et al., Nature Astronomy 6, 496 (2022), 2202.05836.
  * Carleton et al. [2019] T. Carleton, R. Errani, M. Cooper, M. Kaplinghat, J. Peñarrubia, and Y. Guo, MNRAS 485, 382 (2019), 1805.06896.
  * Kong et al. [2022] D. Kong, M. Kaplinghat, H.-B. Yu, F. Fraternali, and P. E. Mancera Piña, APJ 936, 166 (2022), 2204.05981.
  * Hopkins et al. [2018] P. F. Hopkins, A. Wetzel, D. Kereš, C.-A. Faucher-Giguère, E. Quataert, M. Boylan-Kolchin, N. Murray, C. C. Hayward, S. Garrison-Kimmel, C. Hummels, et al., MNRAS 480, 800 (2018), 1702.06148.
  * Wetzel et al. [2022] A. Wetzel et al. (2022), 2202.06969.
  * Pillepich et al. [2018] A. Pillepich et al., Mon. Not. Roy. Astron. Soc. 473, 4077 (2018), 1703.02970.
  * Nelson et al. [2018] D. Nelson et al. (2018), 1812.05609.
  * Behroozi et al. [2013] P. S. Behroozi, R. H. Wechsler, H.-Y. Wu, M. T. Busha, A. A. Klypin, and J. R. Primack, APJ 763, 18 (2013), 1110.4370.
  * Somerville et al. [2008] R. S. Somerville, P. F. Hopkins, T. J. Cox, B. E. Robertson, and L. Hernquist, Mon. Not. Roy. Astron. Soc. 391, 481 (2008), %****␣main.bbl␣Line␣850␣****0808.1227.
  * Benson [2012] A. J. Benson, New Astron. 17, 175 (2012), 1008.1786.
  * Ogiya et al. [2019] G. Ogiya, F. C. van den Bosch, O. Hahn, S. B. Green, T. B. Miller, and A. Burkert, Mon. Not. Roy. Astron. Soc. 485, 189 (2019), 1901.08601.