# Exploring Self-Interacting Dark Matter Halos with Diverse Baryonic
Distributions:
A Parametric Approach

Daneng Yang<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of California, Riverside, California 92521, USA

###### Abstract

Galaxies residing in dark matter halos exert significant gravitational effects
that alter halo structure and dynamics. The complexity of these interactions
escalates with the diversity of galactic structures and the variability in
dark matter halo profiles under self-interacting dark matter (SIDM) models.
This work extends the parametric model for dark matter-only halos presented in
Ref. Yang  _et al._ (2024a) to incorporate baryons. We adapt this model to
consistently represent the SIDM halo density profile over time, highlighting
the role of a gravothermal phase in characterizing the state of an SIDM halo.
Given this phase, the density profile in SIDM is determined by a fictitious
progenitor — consisting of an NFW halo influenced by a baryonic potential —
that has evolved to its present state. In the temporal dimension, the model
incorporates a form factor that rescales the evolution time in the dark
matter-only case, thereby enabling the introduction of a universal phase. In
the radial dimension, the halo density profile is parametrized to reflect the
influences of baryons. We calibrate the model through N-body simulations with
baryon potentials to fit various stellar-to-halo mass ratios and size-mass
relationships. Our parametric approach is numerically efficient, enabling the
exploration of SIDM effects across a diverse set of halos, as exemplified by a
case study using an illustrative sample that spans five orders of magnitude in
the mass range. We also demonstrate that the effects of evolution history and
the specific SIDM model can be separated from the current states of galaxies
and halos, leaving the task of identifying consistent SIDM models to dedicated
post-processing analyses.

###### Contents

  1. I Introduction
  2. II Simulation data
  3. III Normalized gravothermal evolution with baryons
    1. III.1 The universal gravothermal evolution
    2. III.2 A universal form factor for normalizing the evolution time with baryons
    3. III.3 The gravothermal phase of an SIDM halo with baryons
  4. IV Modeling baryon effects on dark matter density profiles
    1. IV.1 The Contracted $\beta 4$ model
    2. IV.2 The Cored-DZ model
  5. V Model validation
  6. VI Effect of SIDM on the baryon profiles
  7. VII Applications
    1. VII.1 Probing SIDM through sampled halos
    2. VII.2 Deciphering SIDM models from gravothermal phases and accretion histories
  8. VIII Conclusion and discussion
  9. A The DM-only limit of the Cored-DZ model
  10. B Testing the hybrid approach
  11. C Adiabatic contraction in different models

##  I Introduction

Cold dark matter (CDM) aligns well with observations at scales larger than
megaparsecs. However, discrepancies arise at smaller scales in the inner
structures of dark matter (DM) halos, where the density profiles can be either
cored or cuspy Flores and Primack (1994); Spergel and Steinhardt (2000); de
Blok  _et al._ (2001); Kuzio de Naray  _et al._ (2010); Oh  _et al._ (2011);
Oman  _et al._ (2015). The cored cases may exhibit fainter or more diffuse
profiles than most models predicted, while the cuspy cases can display steeper
inner-density slopes. To explain these variations, researchers have elaborated
on both baryonic feedback models and self-interacting dark matter (SIDM)
models. Baryonic effects diversify the inner halo structures because it alters
the structure by leading to compact stellar distributions through dissipation
while simultaneously dispersing these distributions through feedback. The
effectiveness of these processes varies with scale and depends on complex
mechanisms including star formation and supernova feedback Navarro  _et al._
(1996); Gnedin and Zhao (2002); Read and Gilmore (2005); Governato  _et al._
(2010); Di Cintio  _et al._ (2014); Chan  _et al._ (2015); Tollet  _et al._
(2016); Lazar  _et al._ (2020); Macciò  _et al._ (2020); Santos-Santos  _et
al._ (2020a). Similarly, SIDM contributes to diversity in inner halo
structures through elastic scatterings that facilitate energy transfer within
halos. This results in core formation and core collapse stages in the
thermodynamic evolution of self-gravitating halos, known as gravothermal
evolution 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._ (2020b); Spergel and Steinhardt
(2000); Correa  _et al._ (2022); Yang  _et al._ (2023a); Balberg  _et al._
(2002, 2002); Koda and Shapiro (2011); Nishikawa  _et al._ (2020); Sameie  _et
al._ (2020); Zavala  _et al._ (2019); Kaplinghat  _et al._ (2019); Kahlhoefer
_et al._ (2019); Correa (2021); Turner  _et al._ (2021); Slone  _et al._
(2023); Silverman  _et al._ (2022); Correa  _et al._ (2022); Yang  _et al._
(2023a); Minor  _et al._ (2021); Yang and Yu (2021); Gilman  _et al._ (2021,
2023); Nadler  _et al._ (2023); Pollack  _et al._ (2015); Choquette  _et al._
(2019); Feng  _et al._ (2021, 2022); Meshveliani  _et al._ (2023); Fischer
_et al._ (2024); Dhanasingham  _et al._ (2023); Zhang  _et al._ (2024). In
probing SIDM, the diversity from baryon physics naturally integrate through
gravitational interactions, resulting in a nonlinear and complex combined
evolution. Nonetheless, understanding the interplay between baryons and dark
matter is vital for making realistic SIDM predictions, as we depend on visible
baryonic matter to probe the elusive dark matter structures Vogelsberger  _et
al._ (2014); Robertson  _et al._ (2019, 2021); Yang and Yu (2021); Rose  _et
al._ (2023); Mastromarino  _et al._ (2023); Rahimi  _et al._ (2023); Zhong
_et al._ (2023); Creasey  _et al._ (2017); Fry  _et al._ (2015); Kong  _et
al._ (2024); Fischer  _et al._ (2023).

To date, there have been only a handful of hydrodynamical cosmological
simulations that incorporate SIDM physics Robertson  _et al._ (2018); Robles
_et al._ (2017); Elbert  _et al._ (2018); Despali  _et al._ (2019); Robertson
_et al._ (2019). These simulations have revealed a complex interplay dependent
on mass, galaxy type, and accretion history. Meanwhile, theoretical studies
have used conducting fluid models, semi-analytic models, and idealized N-body
simulations to model the combined effects of a baryonic potential and
gravothermal evolution Balberg  _et al._ (2002); Koda and Shapiro (2011);
Essig  _et al._ (2019); Nishikawa  _et al._ (2020); Feng  _et al._ (2021);
Yang  _et al._ (2023b); Jiang  _et al._ (2023); Yang  _et al._ (2023c); Zhong
_et al._ (2023). It has been illuminated that galaxy formation can
adiabatically contract dark matter, increase the DM scattering rate, and thus
accelerate the gravothermal evolution. Some empirical relations have also been
proposed to enhance our understanding of these complex processes Elbert  _et
al._ (2018); Sameie  _et al._ (2018); Feng  _et al._ (2021); Yang  _et al._
(2023c); Zhong  _et al._ (2023).

In this study, we explore the gravothermal evolution of DM halos induced by
SIDM in the presence of baryons, utilizing a parametric approach. Within this
framework, each halo is posited to reside in a gravothermal state
characterized by a dimensionless gravothermal phase, $\tau$, defined as the
ratio of the elapsed time since formation to the core collapse time
($t/t_{c}$). The halo’s density profile is then determined from parameters
that describe its characteristics under the CDM scenario, accounting for both
halo and baryonic properties. The model for the density profile is parametric,
with time-dependent parameters represented as functions of $\tau$. These
parameters are calibrated against high-resolution N-body simulations, as
detailed in Ref. Yang  _et al._ (2024a). This method provides an efficient way
to achieve accurate results. It is founded on the universal characteristics of
gravothermal evolution and integrates methods to manage differential SIDM
cross sections and halos with realistic accretion histories Outmezguine  _et
al._ (2023); Yang  _et al._ (2024a); Yang and Yu (2022); Yang  _et al._
(2023b).

Building upon the model of Ref. Yang  _et al._ (2024a) for dark matter-only
halos, we derive and calibrate a new formula for the core collapse time,
making it applicable for a broad range of stellar-to-halo mass and size
ratios. We also introduce procedures to modify the SIDM density profile in DM-
only scenarios to account for the presence of baryons. Unlike the DM-only
case, which has an initial Navarro-Frenk-White (NFW) profile Navarro  _et al._
(1997), our new profile is based on the Dekel-Zhao (DZ) profile Zhao (1996);
Dekel  _et al._ (2017); Freundlich  _et al._ (2020), which has been shown
applicable for halos in hydrodynamical simulations, with the effects of
stellar feedback and adiabatic contraction accounted for. We validate the
model through eight controlled N-body simulations, each employing distinct yet
representative sets of baryon and halo parameters.

Figure 1:  The diverse baryonic structures in isolated halos from the TNG-50-1
simulation, illustrated by the distribution of $\rho_{H}/\rho_{s}$ versus
$r_{H}/r_{s}$ Nelson  _et al._ (2019a, b); Pillepich  _et al._ (2019). Here,
$\rho_{s}$ and $r_{s}$ represent the NFW scale density and radius of the halo
component, respectively, while $\rho_{H}$ and $r_{H}$ denote the corresponding
Hernquist parameters for the stellar component. The color coding represents
the stellar-to-halo mass ratio on a logarithmic scale, indicating that higher
mass halos tend to have higher $\rho_{H}/\rho_{s}$ values on average. The red
dots correspond to the simulated benchmark halos listed in Table 1, which are
used to calibrate the analytic model.

In Fig. 1, we highlight the diverse baryonic structures across a wide range of
samples depicted on the plane of scale density ratio $\rho_{H}/\rho_{s}$
versus scale radius ratio $r_{H}/r_{s}$. Here, the halos and their stellar
components are modeled using NFW and Hernquist Hernquist (1990) profiles,
characterized by scale parameters $\rho_{s}$, $r_{s}$ for halos and
$\rho_{H}$, $r_{H}$ for stellar components. The color codes on the plot
correspond to the stellar-to-halo mass ratio, represented on a logarithmic
scale, revealing a trend where higher $M_{s}/M_{h}$ halos move up along the
$\ln(y)\propto\ln(y)$ direction, exhibiting larger $\rho_{H}/\rho_{s}$ and
$r_{H}/r_{s}$ ratios. Red dots indicate benchmark halos from Table 1, utilized
for validating the model. These halos cover a broad spectrum in the
$\rho_{H}/\rho_{s}$ and $r_{H}/r_{s}$ plane, demonstrating the model’s broad
applicability.

Benchmark | $M_{h}({\rm M_{\odot}})$ | $\rho_{s}({\rm M_{\odot}/kpc^{3}})$ | $r_{s}({\rm kpc})$ | $M_{s}/M_{h}$ | $\rho_{H}/\rho_{s}$ | $r_{H}/r_{s}$ |  $\sigma/m$ ($\rm cm^{2}/g$) | $t_{c,b}~{}\rm(Gyr)$ | $t_{c,b}/t_{c,0}$
---|---|---|---|---|---|---|---|---|---
DM11 | $10^{11}$ | $6.89\times 10^{6}$ | $9.1$ | 0 | 0 | 0 | 10 | 66 | 1.0
DM11+baryon | $10^{11}$ | $6.89\times 10^{6}$ | $9.1$ | 0.010 | 51.7 | 0.0844 | 10 | 25 | 0.38
DM11+baryon R2 | $10^{11}$ | $6.89\times 10^{6}$ | $9.1$ | 0.010 | 13.0 | 0.134 | 10 | 32 | 0.49
DM11+baryon M2 | $10^{11}$ | $6.89\times 10^{6}$ | $9.1$ | 0.020 | 77.5 | 0.0929 | 10 | 12 | 0.18
DM12+baryon | $10^{12}$ | $1.23\times 10^{7}$ | $15.5$ | 0.034 | 149 | 0.0926 | 1 | 9.0 | 0.055
DM13 compact | $10^{13}$ | $1.43\times 10^{7}$ | 42.1 | 0.013 | 2972 | 0.0187 | 1 | 9.0 | 0.19
DM13+baryon | $2.2\times 10^{13}$ | $1.17\times 10^{6}$ | 115 | 0.010 | 1544 | 0.0234 | 1 | 173 | 0.23
DM13+baryon2 | $3\times 10^{13}$ | $5.64\times 10^{6}$ | 67.1 | 0.0074 | 321 | 0.0402 | 1 | 42 | 0.35
DM13 extreme | $3\times 10^{13}$ | $7.80\times 10^{6}$ | 70.0 | 0.0043 | 5442 | 0.0112 | 1 | 27 | 0.38
Table 1:  Properties of N-body simulated benchmarks: The first four
simulations, introduced in Ref.Zhong  _et al._ (2023), correspond to the
initial block of cases in Table 1 therein. The final two benchmarks,
DM13+baryon2 and DM13 extreme are analogous to J1636+4707 and JWST-ER1, as
described in Ref. Kong  _et al._ (2024). Additionally, the DM12+baryon, DM13
compact, and DM13+baryon benchmarks represent further simulations conducted to
broaden the range of benchmark scenarios. The core collapse times are computed
considering $\sigma/m=10~{}\rm cm^{2}/g$ for the DM11 cases in the first block
and $\sigma/m=1~{}\rm cm^{2}/g$ for the DM12 and DM13 cases, to be consistent
with the curves in Fig. 4 (left panel). We use $t_{c,b}$ and $t_{c,0}$ to
represent the core collapse times in the DM+baryon and DM-only scenarios. The
halo masses are calculated as the mass enclosed within radii where the density
is 200 times the critical density ($R_{\rm 200}$), using the scale density
($\rho_{s}$) and scale radius ($r_{s}$) of the NFW profile. The stellar masses
are total masses computed using the scale density ($\rho_{H}$) and scale
radius ($r_{H}$) of the Hernquist profile.

Recognizing the mutual gravitational influence between dark matter and
baryons, we employ the parametric model to investigate the conditions under
which the effects of SIDM on baryons are pronounced. Our findings indicate
that SIDM can notably affect baryons in cases where the maximum core sizes
exceed the baryons’ scale radii. This core formation tends to expand the
baryon distribution, whereas core collapse generally compresses it.
Conversely, when the maximum core size is smaller than the baryon scale
radius, the baryonic potential strongly suppresses distinct core formation
within the halo. As a result, the baryon profile becomes more compact as the
halo contracts. Our study on this topic is simplified; therefore, we leave the
modeling of consistent coevolution of dark matter and baryons to future work.

To demonstrate the effectiveness of our model, we present two illustrative
applications. We first apply the model to a population of sampled isolated
halos, utilizing generic relationships between halos and galaxies, and
focusing on the distribution of maximum circular velocity ($V_{\rm max}$) and
its corresponding radius ($R_{\rm max}$). Next, we illustrate how the
gravothermal evolution induced by SIDM is primarily captured in the
gravothermal phases, which can, in principle, be reconstructed from data. The
CDM and SIDM density profiles are linked via the gravothermal phases, and by
reversing these phases to zero, the outcomes are expected to align with our
existing knowledge of $\Lambda$CDM. Notably, these gravothermal phases are
independent of the halos’ accretion history as well as the underlying SIDM
model. By post-processing these results, we can discern SIDM interpretations
based on stellar accretion histories, aiding in the identification of the
preferred SIDM model that coherently explains all observations. These
applications highlight the model’s adaptability and effectiveness in
addressing the complex interplay between SIDM halos and their baryonic
counterparts.

We organize the paper as follows. Section II provides a summary of the
simulation samples utilized in this study. In Section III, we first discuss
the theoretical basis for the parametric model, following that we introduce a
core collapse time formula that incorporates baryons. Section IV presents
parametric SIDM models based on a contracted $\beta 4$ and Cored-DZ profiles.
The latter is more intricate but consistently incorporates the baryon-induced
contraction for halos hosting massive galaxies. Section V is dedicated to the
validation of the proposed model. In Section VI, we investigate the effect of
SIDM on the baryon distributions using the proposed model. We provide two
example applications in Section VII: applying the parametric model to sampled
halos, and a post-processing of the gravothermal state for the SIDM
interpretations. Finally, Section VIII summarizes the paper and discusses
directions for future research.

##  II Simulation data

As illustrated in Fig. 1, there is a significant diversity in the baryonic
structures within halos. The range of the $\rho_{H}/\rho_{s}$ ratio, in
particular, spans approximately five orders of magnitude. For a better
presentation of physically relevant cases, we selected isolated halos with
masses greater than $10^{10}\rm M_{\odot}$. This criterion is based on the
consideration that the stellar-to-halo mass ratio reaches its peak at halos of
around $10^{12}\rm M_{\odot}$, and the overall halo population increases
towards the lower mass, inversely proportional to the square of the mass.

Given this diversity, it is crucial to test our model against a broad spectrum
of scenarios. To this end, we make use of eight N-body simulations where
baryons are represented using a Hernquist potential. See Table 1 for details
of their properties. Although our sample size is limited, these examples span
a substantial range in both $\rho_{H}/\rho_{s}$ and $r_{H}/r_{s}$ ratios,
ensuring a comprehensive validation of our model across different baryonic
structures and halo configurations. We depict these samples as red dots in
Fig. 1.

Our simulations are idealized in the treatment of baryon potentials, which are
fixed at all times and take the form of the Hernquist potential. While this
setup overlooked baryonic dynamics, it fits into the context of the integral
approach in Ref. Yang  _et al._ (2024a), where each SIDM halo is assumed to
arise from a fictitious CDM halo. In application, the effect of baryon
dynamics, including its interplay with the gravothermal evolution, is encoded
into the evolution of the potential and can be consistently incorporated in
the integral approach, given the parametric model with baryons.

Among the eight samples analyzed, four were previously utilized in Ref. Zhong
_et al._ (2023). In these simulations, each consists of $4\times 10^{6}$
particles, with a force softening length set to $h=0.13$ kpc. The Hernquist
baryonic potentials in these simulations are introduced instantaneously at
$t=0$ Gyr, and analyses have shown that this approach yields results
comparable to those obtained by gradually developing the potential over 4
Gyr—a period significantly shorter than the core collapse time for the
simulated halos.

The remaining samples are simulated first in CDM, with the baryon potential
increased linearly in $\rho_{H}$ over 4 Gyr. The snapshots at $5$ Gyr are
subsequently used as initial conditions for the SIDM simulations. The
simulation parameters vary in the samples. The DM13+baryon and DM13+baryon2
samples are simulated using $2\times 10^{6}$ particles with a force softening
length $h=2.8~{}$kpc, while the other samples are simulated using $10^{6}$
particles with $h=0.42$ kpc. The DM13+baryon2 and DM13 extreme samples
correspond to the J1636+4707 and JWST-ER1 systems, respectively, as detailed
in Ref. Kong  _et al._ (2024). These samples have been re-simulated in this
study to reach the deep core collapse phases.

##  III Normalized gravothermal evolution with baryons

In this section, we first review the theoretical rationale that supports the
construction of parametric models. Following that, we present a new core
collapse formula that incorporates the effects of baryons.

###  III.1 The universal gravothermal evolution

We start with CDM halos that are spherically symmetric and isolated. Their
density profiles are known to be well-described by the NFW profile, which can
be written in terms of two scale parameters $\rho_{s}$ and $r_{s}$ as
$\rho_{\rm NFW}=\rho_{s}/\left((r/r_{s})(1+r/r_{s})^{2}\right)$. This profile
is self-similar in nature: all halos exhibit the same density profile once
normalized by the $\rho_{s}$ and $r_{s}$. In SIDM, scatterings introduce an
“arrow of time” which breaks the self-similarity of the profile. Intriguingly,
such self-similarity is replaced by so-called “universality”, after
incorporating the time dimension into the modeling and with explicit SIDM
dependence absorbed into the evolution time Outmezguine  _et al._ (2023); Yang
_et al._ (2024a); Zhong  _et al._ (2023).

To see how this is achieved, we examine the three fluid equations describing
the halo evolution: the continuity equation, the momentum conservation
equation, and the transport equation. In the context of SIDM, the continuity
equation remains unaffected since the interactions are elastic. The momentum
and transport equations will both be affected by viscosity through the viscous
tensor, denoted as $\Pi_{ij}^{\rm vis}$, which is proportional to
$\eta(\nabla_{i}u_{j}-\nabla_{j}u_{i})$, with $\eta$ representing the
viscosity coefficient and $\mathbf{u}$ denoting the average particle velocity
in a small given volume. However, in spherical halos that maintain isotropic
velocity distributions, the average velocities in both polar and azimuthal
directions become negligible, eliminating viscosity’s impact. Under these
conditions, the transport equation emerges as particularly crucial and drives
the gravothermal evolution. This equation is typically represented as:

$\frac{3\rho}{2m}\frac{DT}{Dt}=-\mathbf{\nabla}\cdot\mathbf{q}-P\mathbf{\nabla}\cdot\mathbf{u},$
(1)

where $D/Dt=(\partial/\partial t+\mathbf{u\cdot\nabla})$ refers to the
material derivative, $T\equiv m\nu^{2}$ denotes the “temperature,”
$P\equiv\rho\nu^{2}$ indicates the “pressure,” and
$\mathbf{q}=-\kappa\mathbf{\nabla}T$ describes the heat flux, linked to
thermal conductivity $\kappa$. In the spherical coordinates, this equation
simplifies to

$\displaystyle\frac{\partial}{\partial r}\left(r^{2}\kappa
m\frac{\partial\nu^{2}}{\partial
r}\right)=r^{2}\rho\nu^{2}\frac{D}{Dt}\ln\frac{\nu^{3}}{\rho},$ (2)

where the left-hand-side depends on SIDM through the heat conductivity
$\kappa$, and the right-hand-side introduces evolution through the material
derivative $D/Dt$. Suppose that $\kappa$ depends on the scattering cross
section per mass through a function $f(\sigma/m)$, which has no $r$ or $t$
dependencies, one can absorb the SIDM dependence into the evolution time by
redefining $t$ as $tf(\sigma/m)$. Notably, the continuity and the momentum
equations primarily represent mass conservation and hydrostatic equilibrium,
both of which are inherently independent of time. Therefore, the evolution of
a dark matter halo’s density profile, starting from an initial NFW profile,
follows a unique trajectory. This uniqueness implies a universal solution for
the gravothermal evolution of dark matter halos in SIDM.

We are particularly interested in how the the conductivity in SIDM halos
depends on the cross section $\sigma/m$. While dark matter halos can be
modeled using fluid equations, their properties differ from ordinary fluid in
several crucial aspects, among which one of the most important differences is
the sparsity in the number of scatterings. In the bulk region of realistic
halos, the scattering rate is suppressed by the dark matter density, making
their mean-free-path longer than the scale height $H=\sqrt{\nu^{2}/(4\pi
G\rho)}$ of the halos. Consequently, a post-scattering particle can orbit
around the halo center multiple times before colliding again. It implies that
the energy transferred through the collision is diluted by the orbital
evolution, reducing the effect of heat conduction. In the literature, this
scenario is called the long-mean-free-path (LMFP) regime, where the heat
conductivity is proportional to the scattering cross section and can be
expressed as Pollack  _et al._ (2015); Choquette  _et al._ (2019); Feng  _et
al._ (2021, 2022)

$\displaystyle\kappa_{\rm LMFP}\propto\frac{\sigma}{m}\frac{\rho\nu^{3}}{mG}.$
(3)

This can be contrasted to the ordinary fluid, where the short-mean-free-path
(SMFP) inhibits heat conduction, resulting in a conductivity $\kappa_{\rm
SMFP}\propto\nu/\sigma$, being inversely proportional to the cross section. We
see that in both the LMFP and SMFP regimes, the conductivities have simple
dependences on the scattering cross section and universality can be achieved.

In practice, SIDM halos that conform to astrophysical constraints
predominantly exist within the LMFP regime across most of their volume and
throughout much of their evolutionary history. Gravothermal core collapse can
intensify their inner densities, potentially initiating the SMFP regime;
however, this occurs in the inner regions, which largely evolve independently
of the outer regions. These decoupled inner regions exhibit a self-similar
density profile described by $r^{-u}$ with $u$ ranges between $2$ and $2.5$.
Numerical simulations suggest a preferred value of $u=2.21$ Lynden-Bell and
Eggleton (1980).

###  III.2 A universal form factor for normalizing the evolution time with
baryons

We have demonstrated the universality in the evolution of DM-only halos. In
the presence of baryons, the contraction in response to the potential
necessarily distorts the halo density profiles, violating the universality in
the spatial dimension. Despite this, the underlying principle that the effects
of SIDM can be integrated into the evolution timeline remains valid.

We assume baryonic components, particularly stellar ones, interact
predominantly through gravity, blending with dark matter in phase space to
achieve a quasi-hydrostatic equilibrium. In this setup, the fluid equations
are still applicable, with the density and velocity dispersion including
contributions from both dark matter and baryons. However, because the baryons
do not participate in this self-scatterings, the heat conductivity is adjusted
to account only for the dark matter density. This adjustment does not affect
the form of Eq. (2) and the SIDM dependence of the conductivity can still be
absorbed into the evolution time through the same redefinition. In the LMFP
regime, we have

$\displaystyle\frac{\partial}{\partial r}\left(r^{2}\frac{\rho_{\rm
DM}\nu_{\rm tot}^{3}}{G}\frac{\partial\nu_{\rm tot}^{2}}{\partial
r}\right)\propto r^{2}\rho_{\rm tot}\nu_{\rm
tot}^{2}\frac{D}{D(t\frac{\sigma}{m})}\ln\frac{\nu_{\rm tot}^{3}}{\rho_{\rm
tot}},$ (4)

where $\rho_{\rm DM}$ refers to the dark matter density and $\rho_{\rm tot}$,
$\nu_{\rm tot}$ represent the total density and velocity dispersion of the
baryon and dark matter components.

Figure 2:  The unbroken universality in the evolution timeline with baryons,
depicted through the evolution of normalized density profiles. For clarity,
three benchmark cases—DM11 (blue), DM11+baryon (red), and DM13 compact
(gray)—from Table 1 are illustrated to demonstrate the predictions of the
parametric model. All the cases evolve into a fixed configuration as $t/t_{c}$
approaching $1.1$. For more details, refer to Section V and Fig. 4.

Before delving into the specifics, we illustrate in Fig. 2 the aspects just
discussed. In this figure, the DM-only case is depicted in blue. The other
cases include different amounts of baryons, and their inner density evolution
and density profiles are compared with the DM-only case in the left and right
panels, respectively. Both panels are normalized to depict the universality
and its breaking clearly. From this figure, it is evident that all scenarios
follow a unified evolution timeline, represented as $t/t_{c}$ and ranging from
$0$ to $1.1$. We define each point along this timeline as a “gravothermal
phase.” The density profiles, conversely, vary according to the amount of
baryonic matter and show different deviations from the DM-only scenario.
Subsequently, we will discuss methods to normalize the evolution time and
adjust the DM density profile for a parametric model.

Equation (4) suggests that the SIDM dependence can be absorbed into the
evolution time. Identifying an equation for the core collapse time would allow
us to normalize this timeline effectively. We introduce the core collapse time
equation here and will discuss the modeling of density profiles in the
following section.

In the literature, the core collapse time for dark matter-only halos with NFW
profiles has been estimated as follows Balberg  _et al._ (2002); Pollack  _et
al._ (2015); Essig  _et al._ (2019)

$\displaystyle
t_{c,0}=\frac{150}{C}\frac{1}{\frac{\sigma}{m}\rho_{s}}\left(\frac{1}{4\pi
G\rho_{s}r_{s}^{2}}\right)^{\frac{1}{2}},$ (5)

where $C\approx 0.9$ is a normalization constant chosen to match the core
collapse time of the simulated DM-only benchmark in Table 1 and it differs
from the $0.75$ used in Ref. Yang  _et al._ (2024a). This approximately $20$%
difference is at the same level of the uncertainty associated with the core
collapse time estimated from N-body simulations.

Equation (5) can be derived assuming the core collapse time to be proportional
to the collisional relaxation time $t_{r}\approx 1/(n\sigma v)$, where
$n=\rho/m$ and $v$ approximately scales as the square root of the inner halo
potential $\lim_{r\to 0}\Phi_{\rm NFW}(r)=4\pi
G\rho_{s}r_{s}^{2}\equiv\Phi_{0,\rm NFW}$. However, the scattering rate does
not always linearly represent the SIDM effect on the halo density profile, and
a satisfactory extension to incorporate baryons has been found challenging. To
compute the relaxation time more accurately, we consider the following
equation from the perspective of energy transport

$\displaystyle t_{R}(r)\propto\frac{M(r)|\Phi(r)|/2}{4\pi r^{2}\kappa|\nabla
T|},$ (6)

where the numerator represents the total energy enclosed in radius $r$ and the
denominator computes the absolute luminosity at $r$. In the dark matter only
case, one can estimate $t_{R}$ as

$\displaystyle t_{R}$ $\displaystyle\propto$
$\displaystyle\frac{4\pi\rho_{s}r_{s}^{3}|\Phi_{0,\rm NFW}|/2}{4\pi
r_{s}^{2}\rho_{s}\frac{\sigma}{m}\frac{|\Phi_{0,\rm NFW}|^{3/2}}{G}|\lim_{r\to
0}\nabla\Phi_{\rm NFW}(r)|},$ $\displaystyle\approx$
$\displaystyle\frac{1}{4\pi\rho_{s}\frac{\sigma}{m}}\left(\frac{1}{4\pi
G\rho_{s}r_{s}^{2}}\right)^{\frac{1}{2}}\propto t_{c,0}$

which is consistent with the equation derived based on the scattering rate in
Eq. (5).

In the presence of baryons, we observe that the $\nabla\Phi_{\rm
Hern}(r)=(2\pi G\rho_{H}r_{H}^{3})/(r+r_{H})^{2}$ for the Hernquist profile
has a large dynamical range between $r_{H}$ and $r_{s}$. Therefore, it is
crucial to find an effective radius where the relative contributions of the
two potentials are appropriately accounted for in estimating the core collapse
time. For this purpose, we introduce

$\displaystyle r_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{r_{s}\Phi_{0,\rm NFW}+\alpha r_{H}\Phi_{\rm
Hern}(0)}{\Phi_{0,\rm NFW}+\alpha\Phi_{\rm Hern}(0)}$ $\displaystyle=$
$\displaystyle\frac{\rho_{s}r_{s}^{3}+\alpha\rho_{H}r_{H}^{3}/2}{\rho_{s}r_{s}^{2}+\alpha\rho_{H}r_{H}^{2}/2}$
$\displaystyle=$ $\displaystyle
r_{s}\frac{1+\alpha\hat{\rho}_{H}\hat{r}_{H}^{3}/2}{1+\alpha\hat{\rho}_{H}\hat{r}_{H}^{2}/2}\equiv
r_{s}\hat{r}_{\rm eff},$

where $\alpha=1.6$ is a parameter tuned against simulations with baryons,
$\hat{\rho}_{H}\equiv\rho_{H}/\rho_{s}$ and $\hat{r}_{H}\equiv r_{H}/r_{s}$
are dimensionless scale density and radius introduced to parameterize the form
factor later. In our discussion, all quantities shown with hats will be
dimensionless.

Base on this $r_{\rm eff}$, we estimate the $t_{R}$ as

$\displaystyle t_{R}$ $\displaystyle\propto$ $\displaystyle\frac{r_{\rm
eff}G}{\frac{\sigma}{m}\sqrt{\Phi_{0,\rm NFW}+\alpha\Phi_{\rm
Hern}(0)}(\lim_{r\to 0}\nabla\Phi_{\rm NFW}+\gamma(\nabla\Phi_{\rm
Hern})(0))}$ $\displaystyle\propto$
$\displaystyle\frac{150}{C}\frac{1}{\frac{\sigma}{m}\left(\frac{\rho_{s}r_{s}}{r_{\rm
eff}}+\gamma\frac{\rho_{H}r_{H}^{3}}{r_{\rm eff}(r_{\rm
eff}+r_{H})^{2}}\right)}\left(\frac{1}{4\pi
G(\rho_{s}r_{s}^{2}+\alpha\rho_{H}r_{H}^{2}/2)}\right)^{\frac{1}{2}}\equiv
t_{c,b},$

where $\gamma=20$ is introduced to control the relative importance of baryons
and dark matter to the temperature gradient. In the last line, we posited that
the core collapse time in the presence of baryons, denoted as $t_{c,b}$, is
proportional to $t_{R}$.

A comparison between Eqs. (III.2) and (5) reveals their similar structures. In
particular, the removal of baryonic components from $t_{c,b}$ automatically
simplifies the equation to the DM-only case. This relation is made explicit by
writting $t_{c,b}=t_{c,0}{\cal F}_{t}(\hat{\rho}_{H},\hat{r}_{H})$, where

$\displaystyle{\cal F}_{t}$ $\displaystyle=$
$\displaystyle\left(\frac{1}{\hat{r}_{\rm
eff}}+\frac{\gamma\hat{\rho}_{H}\hat{r}_{H}^{3}}{\hat{r}_{\rm
eff}(\hat{r}_{\rm
eff}+\hat{r}_{H})^{2}}\right)^{-1}\left(1+\alpha\frac{\hat{\rho}_{H}\hat{r}_{H}^{2}}{2}\right)^{-\frac{1}{2}}$
(10)

is a form factor that only depends on the hatted quantities and contains all
the baryon dependencies.

Figure 3:  Contours of ${\cal F}_{t}\equiv t_{c,b}/t_{c,0}$ on the
$\rho_{H}/\rho_{s}$ ($\hat{\rho}_{H}$) versus $r_{H}/r_{s}$ ($\hat{r}_{H}$)
plane, depicted together with the isolated halos in the TNG-50-1 simulation
(black dots) Nelson  _et al._ (2019a, b); Pillepich  _et al._ (2019) and the
ones we utilized for calibration (red dots, detailed in Table 1).

In Fig. 3, we plot contours of ${\cal F}_{t}$ in the
$\hat{\rho}_{H}-\hat{r}_{H}$ plane and overlay isolated TNG-50 halos and those
used to calibrate the $\alpha$ and $\gamma$ parameters to delineate the
parameter region of interest. Interestingly, the contours align well with the
median behavior of the simulated halos. Halos that scatter towards the upper-
right can have core collapse times several orders of magnitude shorter than
their dark matter-only counterparts. Conversely, towards the lower-left, a
large population of halos remains nearly unaffected by baryons.

The broad range of ${\cal F}_{t}$ values covered by the TNG-50 halos
demonstrates that SIDM can significantly amplify diversity in inner halo
structures. Since the SIDM dependence does not manifest in ${\cal F}_{t}$, our
discussion is generic; physical results should be interpreted by factoring in
the effective SIDM cross sections.

###  III.3 The gravothermal phase of an SIDM halo with baryons

The integral approach, as detailed in Ref. Yang  _et al._ (2024a), is built
upon the assumption that each halo resides in a specific gravothermal state,
which can arise from the gravothermal evolution of a fictitious progenitor of
the same halo and baryon parameters in CDM, regardless of its realistic
accretion history. In this approach, the SIDM effect in a halo is accumulated
through the following equations for the $V_{\rm max}$ and $R_{\rm max}$ of the
dark matter component

$\displaystyle V_{\rm max}(t)$ $\displaystyle=$ $\displaystyle V_{\rm
max,CDM}(t)+\int_{0}^{\tau(t)}d\tau^{\prime}\frac{dV_{\rm
max,Model}(\tau^{\prime})}{d\tau^{\prime}},$ $\displaystyle R_{\rm max}(t)$
$\displaystyle=$ $\displaystyle R_{\rm
max,CDM}(t)+\int_{0}^{\tau(t)}d\tau^{\prime}\frac{dR_{\rm
max,Model}(\tau^{\prime})}{d\tau^{\prime}},$

where the $dV_{\rm max,Model}(\tau^{\prime})/d\tau^{\prime}$ and $dR_{\rm
max,Model}(\tau^{\prime})/d\tau^{\prime}$ terms are analytic equations taken
from the parametric model that describe how fast the dark matter $V_{\rm max}$
and $R_{\rm max}$ evolve along the gravothermal phase $\tau(t)=t/t_{c}$. Based
on the obtained $V_{\rm max}(t)$ and $R_{\rm max}(t)$, one can identify the
fictitious CDM progenitor that gives rise to these values $\tau(t)$. This
fictitious CDM halo only moderately differs from the CDM halo at time $t$;
hence, replacing it with the CDM halo would only decrease a small amount of
accuracy. In Ref. Yang  _et al._ (2024b), a comparison of results with and
without this simplification is performed for dark matter-only halos. In
Appendix B, we verify this point through quantitative examples with baryons.
Applications in this work will take advantage of this observation and use the
CDM halo parameters to replace the fictitious ones.

In practice, the integral can be calculated by summing contributions over
small, incremental time steps. This method effectively yields the accretion
history in SIDM as a direct result of the evaluation process. Moreover, this
method offers heuristic insights: the SIDM effect with realistic accretion
history is obtained by summing over the incremental evolution of many isolated
halos. At any time $t$ in the accretion history, the halo density profile in
SIDM can be determined parametrically, taking the instantaneous CDM halo and
baryon parameters, along with the gravothermal phase, as input. An incremental
change in this evolution can be expressed as follows:

$\displaystyle\rho_{\rm SIDM}(r,\text{``CDM'' halo \& baryon params at
}t,\tau)$ $\displaystyle\Downarrow\ \ t\to t+\delta t$ $\displaystyle\rho_{\rm
SIDM}(r,\text{``CDM'' halo \& baryon params at }t+\delta t,\tau+\delta\tau)$

where the use of quotes around CDM underscores the use of a fictitious CDM
progenitor is theoretically most plausible and can enhance the accuracy. In
the presence of baryons, the fictitious progenitor should be solved under a
fixed baryon configuration at the same time. Therefore, we added baryon
parameters to the argument list.

The gravothermal phase in the integral approach can be computed independently.
Taking the evolution time to commence at the time of the halo formation, it
can be calculated as

$\tau(t)=\int_{0}^{t}\frac{dt}{t_{c,b}[\sigma_{\rm
eff}(t)/m,\rho_{s}(t),r_{s}(t),\rho_{H}(t),r_{H}(t)]},$

where $\rho_{s}(t)$ and $r_{s}(t)$ specifies an NFW halo, and $\rho_{H}(t)$,
$r_{H}(t)$ specifies the baryon profile, both at the time $t$. The SIDM cross
section can be time-dependent for a differential cross section. At each time,
we average out the angular and velocity dependencies by computing the
effective constant cross section as Yang and Yu (2022); Yang  _et al._ (2023b)

$\displaystyle\sigma_{\rm eff}$ $\displaystyle=$ $\displaystyle\frac{2\int
dvd\cos\theta\frac{d\sigma}{d\cos\theta}\sin^{2}\theta v^{5}f_{\rm
MB}(v,\nu_{\rm eff})}{\int dvd\cos\theta\sin^{2}\theta v^{5}f_{\rm
MB}(v,\nu_{\rm eff})},$ (12)

where $v$ denotes the relative velocity between the two incoming particles,
$\theta$ refers to the polar angle that takes values in $[0,\pi]$, and

$\displaystyle f_{\rm MB}(v,\nu_{\rm eff})\propto
v^{2}\exp\left[-\frac{v^{2}}{4\nu_{\rm eff}^{2}}\right],$ (13)

is a Maxwell-Boltzmann velocity distribution that approximates the dark matter
velocity distribution. In the presence of baryons, the effective velocity
dispersion can be obtained by rescaling the dark matter-only case result
($\approx 0.64V_{\rm max,NFW}\approx 1.05r_{s}\sqrt{G\rho_{s}}$) by
$\sqrt{\Phi_{\rm 0,NFW}+\Phi_{\rm Hern}(0)}/\sqrt{\Phi_{\rm 0,NFW}}$.

For the integral approach to yield reliable results, the gravothermal
evolution must proceed slowly during a period of a merger event. This
condition implies that the CDM halo profile is predominantly described by the
NFW profile during a notable change in gravothermal evolution. To justify this
condition, one can compare the gravitational relaxation time with a small
increment in the gravothermal phase $\Delta\tau$. Given that these parameters
vary depending on the halo and the specific SIDM model, we perform a
quantitative assessment for a particularly challenging case. We consider the
BM2 halo benchmark detailed in Ref. Yang  _et al._ (2024a), which is
characterized by a low mass of $2\times 10^{7}~{}\rm M_{\odot}$, a high
concentration of $c_{200}=19.7$ assuming $z=2$, and has been used to calibrate
the parametric model. By adopting an effective constant cross section of
$68~{}\rm cm^{2}/g$, we anticipate this halo to undergo core collapse within
$3$ Gyr. The corresponding gravitational relaxation time is estimated to be
$\sqrt{3\pi/(16G\rho_{s})}\approx 0.02~{}\rm Gyr$, which is about $150$ times
shorter than the core collapse duration, allowing for $\Delta\tau$ to be
smaller than one percent in one gravitational relaxation time.

In Ref. Yang  _et al._ (2024a), the integral approach is validated using halos
from a cosmological SIDM simulation. The consistent agreement across halos of
various mass scales and differing effective SIDM cross sections provides
robust support for employing the idea of the gravothermal state. Notably, most
merger events do not compromise the performance of the model because the
changes in the gravothermal phases are negligible during these events. Even if
a halo is in the collapsing phase during a merger, the effectiveness of the
parametric model does not seem to decrease significantly. Such violation
effect is explored in Ref. Yang  _et al._ (2024b) using an extreme SIDM cross
section simulation.

##  IV Modeling baryon effects on dark matter density profiles

When baryons are present, the contracted halo profile is more accurately
described by the DZ profile than by the NFW profile. However, as SIDM
thermalizes the inner regions of a halo, this initial discrepancy diminishes,
allowing for the calibration of a parametric model using the cored NFW
profiles from the DM-only scenario. Specifically, Ref. Zhong  _et al._ (2023)
demonstrated that the differences between an instantaneously inserted and a
gradually grown baryon potential become negligible in halos undergoing core
formation.

Therefore, we introduce two versions of the parametric model. The first is
based on the SIDM profile from Ref. Yang  _et al._ (2024a), hereafter $\beta
4$ profile, which we adapt to include baryon effects. We name it as the
Contracted $\beta 4$ model. The second, based on the DZ profile and adjusted
for SIDM effects, will be termed as the Cored-DZ model. This latter model is
more complex but effectively accounts for baryon-induced contraction in both
early gravothermal phases and halos containing massive galaxies.

###  IV.1 The Contracted $\beta 4$ model

When the effect of adiabatic contraction is not significant, the
thermalization caused by SIDM would be most important for determining the
shape of the density profile. In this case, we found the $\beta 4$ profile
still applies, which reads

$\displaystyle\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}},$
(14)

and the SIDM halo profile at $\tau$ is specified by evaluating the
$\rho_{s}(\tau)$, $r_{s}(\tau)$, and $r_{c}(\tau)$. For completeness, we
rewrite these three trajectories in the DM-only case as
$\rho_{s}(\tau)=\rho_{s,0}g_{\rho}(\tau)$, $r_{s}(\tau)=r_{s,0}g_{r}(\tau)$,
and $r_{c}(\tau)=r_{s,0}g_{c}(\tau)$, where $\rho_{s,0}$ and $r_{s,0}$ are
initial NFW parameters and

$\displaystyle g_{\rho}(\tau)$ $\displaystyle=$ $\displaystyle
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),$ $\displaystyle g_{r}(\tau)$
$\displaystyle=$ $\displaystyle
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),$ $\displaystyle g_{c}(\tau)$
$\displaystyle=$ $\displaystyle
2.555\sqrt{\tau}-3.632\tau+2.131\tau^{2}-1.415\tau^{3}+0.4683\tau^{4}.$ (15)

For core collapsed halos having $\tau>1.1$, we freeze the density profile
taking $\tau=1.1$, assuming the post-collapse evolution to proceed in the
inner halo region, which separates from the outer region evolution.

In the presence of baryons, one only needs two simple adjustments to
incorporate the baryon effect and the SIDM profile is obtained as

$\displaystyle\rho_{s}(\tau)$ $\displaystyle=$
$\displaystyle\rho_{s,0}g_{\rho}(\tau),$ (16) $\displaystyle r_{s}$
$\displaystyle=$ $\displaystyle r_{s,0}g_{r}(\tau),$ $\displaystyle r_{c}$
$\displaystyle=$ $\displaystyle r_{s,0}g_{c}(\tau)({\cal F}_{t})^{2},$
$\displaystyle k$ $\displaystyle=$ $\displaystyle 4({\cal F}_{t})^{1/2},$

where $\tau=t/t_{c,b}$ is the gravothermal phase with baryon effect
incorporated. We have adjusted the core size and its transition by the form
factor ${\cal F}_{t}=t_{c,b}/t_{c,0}$, with the indices calibrated based on
the simulations.

In Appendix A, we demonstrate the quality of this simplified scenario in Fig.
15 and discuss the necessity of incorporating adiabatic contraction for more
massive halos.

###  IV.2 The Cored-DZ model

As depicted in Fig. 2, the dark matter density profile changes in response to
the baryon potential, exhibiting nonuniversal behavior. To model this effect
in a manner that depends on the baryon content, we utilize the Dekel-Zhao (DZ)
profile, which has been demonstrated to accurately fit the density profiles of
dark matter halos in cosmological hydrodynamical simulations Freundlich  _et
al._ (2020). We introduce a core modification to the DZ profile, similar to
our approach in the DM-only scenario, and incorporate additional functions to
ensure it aligns with the key characteristics of the DZ profile during
gravothermal evolution.

The DZ profile takes the following form

$\displaystyle\rho_{\rm
DZ}(r)=\frac{\rho_{x}}{\left(\frac{r}{r_{x}}\right)^{a}\left(1+\left(\frac{r}{r_{x}}\right)^{1/2}\right)^{2(3.5-a)}}$
(17)

where $\rho_{x}$ and $r_{x}$ are scale density and radius that can be obtained
through a concentration parameter $c=R_{\rm cut}/r_{x}$, with $R_{\rm cut}$
being the virial radius of the halo, and a density slope parameter
$a=-\lim_{{r\to 0}}d\ln\rho/d\ln r$. In terms of $a$ and $c$,

$\displaystyle\rho_{x}$ $\displaystyle=$
$\displaystyle\left(1-a/3\right)c^{3}\mu\bar{\rho}_{\rm h},$ (18)
$\displaystyle r_{x}$ $\displaystyle=$ $\displaystyle R_{\rm cut}/c,$

and $\mu=c^{a-3}\left(1+c^{1/2}\right)^{2(3-a)}$, $\bar{\rho}_{\rm
h}=M_{h}/\left(\frac{4\pi}{3}R_{\rm cut}^{3}\right)$.

In the presence of baryons, the $a$ and $c$ parameters can be determined by
fitting the contracted profile, which can be obtained following the method of
Ref. Gnedin  _et al._ (2004). For simplified studies, an analytic method in
Ref. Freundlich  _et al._ (2020) could be used, where these quantities are
computed as a function of stellar to halo mass ratio $M_{s}/M_{h}$. However,
as we illustrated in Fig. 1, there is a significant spread in the
$\rho_{H}/\rho_{s}$ and $r_{H}/r_{s}$ plane with increasing $M_{s}/M_{h}$
ratios, indicating the need for enhanced theoretical modeling to achieve
greater accuracy. In this study, we employ the method detailed in Ref. Gnedin
_et al._ (2004) to model the contracted profiles accurately. The $a$ and $c$
parameters are extracted by fitting the contracted profiles.

The CoredDZ profile is parameterized as

$\displaystyle\rho_{\rm CoredDZ}(r)$ $\displaystyle=$
$\displaystyle\frac{f_{\rm in}(r)\rho_{x}f_{\rm
out}(r)}{\frac{\left(r^{k}+r_{c}^{k}\right)^{1/{k}}}{r_{x}}\left(1+\left(\frac{r}{r_{x}}\right)^{1/2}\right)^{2(3.5-a)}},$
(19)

where we have introduced two functions to reshape the inner and outer profiles

$\displaystyle f_{\rm
in}(r)=\left(\frac{r}{r_{x}}+\frac{r_{c}}{0.4r_{s}}\left(\frac{\rho_{x,0}r_{x,0}}{\rho_{s}r_{s}+0.4\rho_{H}r_{H}}\right)^{1/(a-1)}\right)^{1-a}$
(20) $\displaystyle f_{\rm out}(r)=\left(1+\frac{r}{R_{\rm
cut}}\left(\frac{r_{x}}{r_{x,0}}-1\right)\right)^{-1/2}.$

While one can formally describe the form factor that cores the DZ profile as
“$\rho_{\rm CoredDZ}(r)={\cal F}_{r}(r)\rho_{\rm DZ}(r)$,” in practice, we
have found it more convenient to use a hybrid density profile. This profile
introduces a core analogous to the contracted $\beta 4$ profile while
preserving the second factor in the denominator as in the DZ profile. The
function $f_{\rm in}(r)$ ensures that in the limit of a vanishing core, the
profile reverts to the DZ profile, with an additional constant term
proportional to the core size to maintain a cored configuration in regions
much smaller than $r_{c}$. Meanwhile, $f_{\rm out}(r)$ addresses the
distortion caused by differing power slopes at large radii between the NFW and
DZ profiles, ensuring minimal change in the total halo mass during
gravothermal evolution. A key aspect of “${\cal F}_{r}$” is that it depends
not only on the radius but also on the interplay between baryon and dark
matter properties. This feature reflects the loss of universality in the
radial direction: rather than a single profile that fits all, we need a family
of profiles parametrized by factors $a$ and $c$, which are coupled to the
baryon parameters.

To model gravothermal evolution with baryons based on the Cored-DZ profile, we
found it sufficient to start with the trajectories of the $\beta 4$ profile in
Ref. Yang  _et al._ (2024a). The effect induced by the baryons can be
incorporated by adjusting the core size $r_{c}$ and the power index $k$ as a
function of the gravothermal phase $\tau$. The parameter trajectories for the
Cored-DZ profile are

$\displaystyle\rho_{x}$ $\displaystyle=$
$\displaystyle\rho_{x,0}g_{\rho}(\tau),$ $\displaystyle r_{x}$
$\displaystyle=$ $\displaystyle r_{x,0}g_{r}(\tau),$ $\displaystyle r_{c}$
$\displaystyle=$ $\displaystyle r_{s,0}g_{c}(\tau)({\cal F}_{t})^{2},$ (21)
$\displaystyle k$ $\displaystyle=$ $\displaystyle 4({\cal F}_{t})^{1/2},$

where the $\rho_{x,0}$ and $r_{x,0}$ correspond to the initial DZ profile
parameters and the $r_{s,0}$ is the initial NFW scale radius.

This modified profile facilitates the initiation of gravothermal evolution
from a profile that more closely mirrors a realistic CDM context,
incorporating the contraction induced by baryons. The evolution of this
response in the context of SIDM has been calibrated using a series of N-body
simulations listed in Table 1 with a fixed baryon potential.

##  V Model validation

Figure 4:  The averaged inner halo density evolution from N-body simulations
(colored bands) and model predictions (colored curves) for the benchmarks in
Table 1. For clearer visualization, inner densities are calculated within
$0.3$ kpc for the DM11 cases and $1$ kpc for the DM12 and DM13 cases. Core
collapse timings are determined based on an SIDM of
$\sigma/m=10~{}\mathrm{cm}^{2}/\mathrm{g}$ for DM11 cases and
$\sigma/m=1~{}\mathrm{cm}^{2}/\mathrm{g}$ for the DM12 and DM13 cases. In the
DM13 cases (dashed), the presence of heavy, compact baryons leads to the
compression of the innermost density profiles, whereas the outer regions are
less sensitive to the baryon potential and evolve into shallower densities in
N-body simulations.

Based on the samples in Table 1, we adjust the $\alpha$ and $\gamma$
parameters in Eq. III.2 to obtain an overall agreement between the core
collapse times predicted by the equation and derived from simulations. The
optimal values we found are $\alpha=1.6$ and $\gamma=20$ and will be used
throughout the paper. The additional factors, including the $({\cal
F}_{t})^{2}$ for the core size and $({\cal F}_{t})^{1/2}$ for the $k$, are
also introduced based on the quality of model predictions versus the simulated
density profiles.

In Fig. 4, we present the evolution of the averaged inner density for all
eight samples, alongside the DM-only scenario (DM11, shown in blue). The
results from N-body simulations are depicted with colored bands, while the
predictions from our model are represented by colored curves for comparison.
To enhance clarity, we have calculated inner densities within $0.3$ kpc for
the DM11 scenarios and within $1$ kpc for the DM12 and DM13 scenarios. The
core collapse times have been evaluated assuming an SIDM cross section of
$\sigma/m=10~{}\mathrm{cm}^{2}/\mathrm{g}$ for the DM11 cases and
$\sigma/m=1~{}\mathrm{cm}^{2}/\mathrm{g}$ for the DM12 and DM13 cases.
Additionally, to incorporate the effect of adiabatic contraction, which is
significant for the DM11+baryon M2, DM12 and DM13 cases, we fit the contracted
initial conditions in the simulation using the DZ profile and used them as
input for the initial condition of the Cored-DZ parametric model for the DM12
and DM13 cases. For the DM11+baryon M2 case where the simulated initial
condition is NFW, we use the equations for $s_{1}$ and $c_{2}$ in Ref.
Freundlich  _et al._ (2020) (Eq. 45 and 46, with coefficients in Table 1) to
incorporate this effect.

As previously discussed in terms of Fig. 2, baryon potentials alter dark
matter density profiles in ways that are not universally predictable. Our
model, which has been calibrated using controlled N-body simulations, is
effective in accounting for the impact of adiabatic contraction by assuming
that the contracted profiles can consistently be described by the Cored-DZ
profile at all times. Nonetheless, when the baryon content is extremely
compact—indicated by high $\rho_{H}/\rho_{s}$ ratios and low $r_{H}/r_{s}$
ratios—the accuracy of the Cored-DZ profile diminishes, leading to
discrepancies between the model predictions and the N-body simulations. In
Fig. 4, the DM13 scenarios have dense and compact baryons that exert a
stronger contraction on the dark matter within the inner regions of the halo,
whereas the outer regions remain relatively less affected by the baryon
potential. This discrepancy results in the formation of larger cores and
steeper density minima in N-body simulations compared to what our parametric
model predicts. In the cases of DM11+baryon M2 (orange) and DM12+baryon
(magenta), the early phases of core formation are accurately captured by the
Cored-DZ profile. However, as they evolve into the core collapse phase, the
relatively high $\rho_{H}/\rho_{s}$ and $r_{H}/r_{s}$ values lead to late-
stage distortions in the density profile during the gravothermal evolution.
While the large $\rho_{H}/\rho_{s}$ compacts the inner halo earlier on, a
larger $r_{H}$ reduces the degree of enhancement within this region.
Nevertheless, the core collapse times predicted by models align closely with
those from simulations across all cases.

In Fig. 5, we present the simulated (colored bands) and Cored-DZ model
(colored curves) predicted halo profiles at three representative gravothermal
phases: $t/t_{c}=0,0.2$, and $1$. In the left panel ($t/t_{c}\approx 0$), the
DM11 cases in the simulations are close to NFW. This is because they are
simulated with an instant insertion of baryon potentials. We show the results
of these cases at $t=0.25$ Gyr to allow some initial evolution away from the
original NFW profile. Nevertheless, the impact of adiabatic contraction in
these instances remains minimal. In comparison, the initial conditions of DM12
and DM13 cases are very well modeled by the DZ profile. In the middle panel,
the agreement is good for most cases, except that the simulated DM13+baryon2
and DM13 extreme cases demonstrate shallower profiles around the maximum
cores, being consistent with what we discussed for Fig. 4. At $t/t_{c}\approx
1$, the core collapse time, as calculated using Eq. (III.2), is found to be
10% (30%) shorter than the simulated DM13+baryon2 (DM13 extreme). To align the
profiles for equivalent gravothermal phases, we adjust the timing of the
simulated curves accordingly in these specific cases. While distortions in the
shape exist, the parametric model correctly captures most features even in
this collapsed phase.

##  VI Effect of SIDM on the baryon profiles

Figure 5:  The simulated (colored bands) and Core-DZ model predicted (colored
curves) halo density profiles at three representative gravothermal phases:
$t/t_{c}\approx 0,0.2$, and $1$. The DM12 and DM13 scenarios use a contracted
CDM profile as the initial condition, whereas the DM11 scenarios commence with
an instant insertion method. In the left panel ($t/t_{c}\approx 0$), the DM11
cases are depicted at $t=0.25$ Gyr to allow some initial evolution away from
the original NFW profile. At $t/t_{c}\approx 1$, the core collapse time, as
calculated using Eq. (III.2), is found to be 10% (30%) shorter than the
simulated DM13+baryon2 (DM13 extreme). To align the profiles for equivalent
gravothermal phases, we adjust the timing of the simulated curves accordingly
in these specific cases.

In this work, we focus on the impact of baryons on SIDM halo density profiles.
However, the gravitational influence between baryons and dark matter is
mutual, and understanding the feedback of SIDM gravothermal evolution on
baryonic evolution is crucial. In this section, we explore the role of such
feedback in the $\rho_{H}/\rho_{s}-r_{H}/r_{s}$ plane.

Figure 6:  Contours of the maximum core size $r_{c,\rm max}$ relative to
$r_{H}$ on the $\rho_{H}/\rho_{s}$ ($\hat{\rho}_{H}$) versus $r_{H}/r_{s}$
($\hat{r}_{H}$) plane, obtained using the parametric model with baryons. The
isolated halos in the TNG-50-1 simulation (black dots) Nelson  _et al._
(2019a, b); Pillepich  _et al._ (2019) and the ones utilized for calibration
(red dots, detailed in Table 1), are shown to indicate the region of interest.
The red-dashed line highlights the contour at $r_{c,\rm max}=r_{H}$, which is
close to the median of samples.

Figure 7:  Evolution of the baryon scale radius (blue) and density (red) for
the DM11 halo, depicted across three panels. The left panel corresponds to the
median baryon distribution from the DM11+baryon benchmark (Table 1), where
$r_{c,\rm max}\approx r_{H}$, aligning with the dashed red curve in Fig. 6.
The middle and right panels show scenarios with the scale radius decreased and
increased by 50% concerning the case on the left while maintaining the
$\rho_{H}$. These adjustments position the middle panel below and the right
panel above the dashed red curve of $r_{c,\rm max}\approx r_{H}$.

We are particularly interested in comparing the maximum core size, $r_{c,\rm
max}$, with the radial scale $r_{H}$ of the baryon profile, as we anticipate
that core expansion will make the baryon distribution more diffuse, while core
collapse will compact it. In Fig. 6, we present the contours of $r_{c,\rm
max}=0.49r_{s,0}({\cal F}_{t})^{2}$ relative to $r_{H}$ on the
$\rho_{H}/\rho_{s}-r_{H}/r_{s}$ plane. Similar to the findings in Fig. 3 for
${\cal F}_{t}$, the ratio $r_{c,\rm max}/r_{H}$ varies by several orders of
magnitude, with values as low as $10^{-5}$ towards the upper-right points and
exceeding $100$ towards the lower-left. These results also indicate a
diversity in the effect of SIDM feedback on baryon distributions.

We also observe that the expression $r_{c,\rm max}=0.49r_{s,0}({\cal
F}_{t})^{2}$ introduces dependencies on $\hat{\rho}_{H}$ and $\hat{r}_{H}$
that align well with the median data points. This alignment, which contrasts
with the dark matter-only scenario’s lack of scale density dependence, further
justifies our choice of the $({\cal F}_{t})^{2}$ factor.

To better understand the influence of SIDM on baryons for varying ratios of
maximum core radius to baryon scale radius $r_{c,\rm max}/r_{H}$, we employ
the approximate conservation principle from the simple adiabatic contraction
model described in Ref. Ryden and Gunn (1987); Gnedin  _et al._ (2004). We
refine the integral approach to allow variation in the baryon profile while
maintaining a Hernquist profile. Given a small timestep forward, the change in
the baryon profile is solved using

$\displaystyle[M_{h}(r_{H})+M_{s}(r_{H})]r_{H}=[M_{h}(r_{H}^{\prime})+M_{s}(r_{H}^{\prime})]r_{H}^{\prime},$
(22)

where $r_{H}^{\prime}$ denotes the updated Hernquist radius. The scale density
can be derived from $r_{H}^{\prime}$ and the total baryon mass
$M_{s}^{\prime}$ as $\rho_{H}^{\prime}=M_{s}^{\prime}/(2\pi r_{H}^{\prime
3})$.

For this analysis, we take the DM11+baryon system as an example and assume the
dark matter and baryon masses to both be conserved. We first take the initial
condition as given in Table 1. We modify the cross section to be $25.3~{}\rm
cm^{2}/g$ such that $t_{c,b}=10~{}$Gyr. Evolving this initial condition for 10
Gyr, we obtain the results in Fig. 7 (left panel). We observe the expected
features in this figure: the baryon scale radius (blue) increases during the
core formation phase $0-2~{}$Gyr, after which it decreases in response to the
collapse. The evolution of scale density (red) is solved from mass
conservation. It looks much more significant than the size evolution because
$\rho_{H}\propto r_{H}^{-3}$.

Notably, we choose the baryon content of the DM11 system to be close to the
median in both the stellar to halo mass ratio and the stellar size-mass
relation, which makes $r_{c,\rm max}\approx r_{H}$ and resulting in a point
residing on the dashed-red line of Fig. 6. To explore the effect of SIDM on
baryons in regions above and below the $r_{c,\rm max}=r_{H}$ contour (dashed
red in Fig. 6), we perform analogous studies for two cases: one with $r_{H}$
reduced by half and one with $r_{H}$ increased by half; the $\rho_{H}$ is
fixed and hence the total stellar masses and core collapse times are modified
accordingly.

The obtained results are illustrated in the middle and right panels of Fig. 7.
The middle panel case corresponds to the region below the dashed red curve. In
this region, $r_{c,\rm max}>r_{H}$ and hence the baryon content can be more
sensitive to the SIDM core evolution. In the simulated case, even if the
stellar mass is reduced by a factor of $2^{3}$, the amount of size growth is
comparable (slightly larger) to the case with the original $r_{H}$ value. The
right panel case, instead, represents a case above the dashed red curve, where
$r_{c,\rm max}<r_{H}$. This region is featured by a strong baryon effect on
SIDM halos that can be cross-checked with contours in Fig. 3. In the case that
we simulate, the total baryon mass is about $3.4$ times the original case and
the core collapse time is reduced to about 3 Gyr. Interestingly, we found the
baryon size evolution to have bypassed the expansion phase and go directly
into the collapsing phase. This result can be correlated to the behavior in
the dark matter distribution. For example, the DM12+baryon case has
$\rho_{H}/\rho_{s}=149$ and $r_{H}/r_{s}=0.09$, which is well above the dashed
red curve. Its inner density evolution corresponds to the magenta curve in
Fig. 4, where we see no signs of core formation. Similarly, the DM11+baryon M2
case also resides above the dashed red curve and demonstrates no core
formation.

##  VII Applications

Figure 8:  Distributions of maximum total circular velocity ($V_{\rm max}$)
and radius at which this velocity occurs ($R_{\rm max}$) for sampled halos.
The left panels depict results for the CDM model. The middle and right panels
illustrate outcomes from two SIDM models: the first adopts a constant cross
section per mass $\sigma/m=10~{}\rm cm^{2}/g$, and the second features a
velocity-dependent cross section per mass, parameterized through a Rutherford-
like cross section, with $\sigma_{0}/m=100~{}\rm cm^{2}/g$ and $w=100~{}\rm
km/s$ in Eq. (23). The upper row displays results in a dark matter-only
scenario, while the middle and lower rows show the impacts of baryonic
presence: the middle row details the $V_{\rm max,DM}$ and $R_{\rm max,DM}$ for
the dark matter component alone, and the lower row includes the contributions
from both dark matter and baryons. The SIDM predictions with baryons are
obtained through the integral approach. All halos are uniformly sampled using
a flat halo mass function over the mass range $M_{h}\in[10^{9},10^{14}]\rm\
M_{\odot}$.

As with the DM-only parametric model, the baryonic extension of this model can
be effectively applied to CDM simulations or semi-analytic models to forecast
SIDM outcomes. The dark matter-only version of the model has recently be
implemented into the Semi-Analytical Sub-Halo Inference ModelIng (SASHIMI)
program Ando  _et al._ (2024). Here, we demonstrate such an application using
the Cored-DZ model on a sample of isolated halos. This application simplifies
the modeling by considering generic accretion histories that do not
differentiate between late and early-type galaxies. Additionally, we disregard
the SIDM impacts on baryons discussed in Section VI, due to the absence of a
calibration and validation using SIDM simulation data. Despite this omission,
the influence on the maximum circular velocity ($V_{\rm max}$) and its
corresponding radius ($r_{\rm max}$) is minimal when considering combined dark
matter and baryon contributions. Thus, our results focus solely on these
parameters for clarity.

The parametric model also facilitates extracting the profile parameters from
observed galactic rotation curves. It is advantageous because it parametrizes
a halo’s gravothermal state, bypassing the complexities in the accretion
histories and the effect of SIDM feedback on baryon profiles. As an
illustration, we show that SIDM parameters can be effectively retrieved given
an assumed accretion history via post-processing, underscoring the model’s
utility in analyzing complex astrophysical phenomena.

###  VII.1 Probing SIDM through sampled halos

We demonstrate the effectiveness of the parametric model with baryons using a
sample of isolated halos. These halos are constructed based on established
relations in the literature for exploring the SIDM impact on the scatter
within the $V_{\rm max}-R_{\rm max}$ plane. To efficiently explore a wide
spectrum of halos, we sample the halo mass in $[10^{9},10^{14}]\rm\ M_{\odot}$
using a flat mass function instead of the approximately $M^{-2}$ mass
function. Based on the halo masses, we sample the halo concentrations based on
the median $z=0$ concentration-mass relation from Ref. Dutton and Macciò
(2014), incorporating a mass-independent lognormal scatter of 0.11 dex. All
other halo properties can be computed based on the sampled masses and
concentrations.

For baryonic components, we sample the stellar masses according to the
stellar-to-halo-mass relation from Ref. Moster  _et al._ (2013), with an
intrinsic scatter of 0.15 dex. The stellar sizes are determined using the
size-mass relation for red galaxies from Ref. Carleton  _et al._ (2019), with
a 0.1 dex scatter incorporated and the half-light radii ($r_{e}$) converted to
Hernquist scale radii using $r_{H}=r_{e}(\sqrt{2}-1)4/3$. The inner dark
matter scale densities and radii, represented through the NFW $\rho_{s}$ and
$r_{s}$, are assumed to have been established 10 Gyr ago, a simplification
supported by hierarchical structure formation theories (see section 3.1 of
Ref. Yang  _et al._ (2024a) for associated discussion). In contrast,
significant growth is assumed in stellar parameters within the last 10 Gyr,
modeled using the average stellar accretion histories from Ref. Behroozi  _et
al._ (2013). We simplify the stellar size evolution considering
$r_{H}=r_{H,z=0}(1+z)^{-0.75}$, which ignores dependencies on mass and galaxy
type; for a comprehensive treatment, equations from Ref. van der Wel  _et al._
(2014) can be employed.

To calculate $V_{\rm max,DM}$ and $R_{\rm max,DM}$ in the presence of an
accretion history, we apply the integral approach based on the discretized
version of Eq.(LABEL:eq:int). In this approach, each incremental step
calculates the core collapse time using the current configurations of the halo
and baryons. The Cored-DZ model is applied to ensure accuracy.

When incorporating baryons using the integral approach, subtle differences
emerge compared to the dark matter-only scenario. One is that the changes in
$V_{\rm max,DM}$ and $R_{\rm max,DM}$ must be calculated based on the
differences in the density profiles of two successive snapshots. The second is
that we minimize the discrepancies between the predicted and actual values of
$V_{\rm max,DM}$ and $R_{\rm max,DM}$ to derive the parameters of the
fictitious CDM halo. An efficient code for doing this is provided on the
program page.

These subtleties can significantly increase the numerical complexity of the
integral approach compared with the dark matter-only case. Here we introduce a
hybrid approach that yields results reasonably close to those obtained using
the integral approach but significantly reduces the computational requirement.
In this hybrid model, we disregard the differences between fictitious CDM
halos and actual simulated CDM halos at $z=0$. Instead, we incorporate all
accretion effects through the gravothermal phase at $z=0$, calculated as

$\tau(z=0)=\int_{0}^{10~{}\rm Gyr}\frac{dt}{t_{c,b}[\sigma_{\rm
eff}/m,\rho_{s},r_{s},\rho_{H}(t),r_{H}(t)]}.$

Given this gravothermal phase, obtaining the density profile becomes identical
to that of the basic approach. In Appendix B, we provide detailed comparisons
between the integral and hybrid approaches. Our analysis shows that the hybrid
approach effectively captures the overall trends and variations associated
with SIDM, with secondary effects from the integral approach being suppressed
once baryonic contributions to $V_{\rm max}$ and $R_{\rm max}$ are included.
When applying our model, we recommend starting with the hybrid approach and
then assessing whether the more complex integral approach is necessary based
on the problem’s required precision.

For the study here, we analyze the CDM and two variants of SIDM models. The
first SIDM model employs a constant cross section per mass, set at
$\sigma/m=10~{}\rm cm^{2}/g$. The second SIDM model uses a velocity-dependent
cross section, modeled assuming Rutherford-like scatterings, as described by

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

where $\sigma_{0}/m=100~{}\rm cm^{2}/g$ and $w=100~{}\rm km/s$.

Our findings are depicted in Fig. 8, organized from left to right, showcasing
results for the CDM, constant cross section SIDM, and velocity-dependent SIDM
models. Each model is examined across three rows of figures to elucidate the
effects of baryons: the top row presents the dark matter-only scenario, while
the middle and bottom rows explore the combined influence of dark matter and
baryons. Specifically, the middle row visualizes $V_{\rm max}$ and $R_{\rm
max}$ for the dark matter component alone, and the bottom row incorporates the
baryonic contributions.

In the CDM case, results for the dark matter-only scenario are color-coded by
halo concentration, revealing that higher concentrations are associated with
lower values of $V_{\rm max}$ and $R_{\rm max}$. When baryons are included,
results are color-coded by the stellar-to-halo mass ratio to highlight the
impact of stellar mass, which tends to compactify and increase the mass of
systems, especially in halos around the Milky Way scale ($\sim 10^{12}\rm
M_{\odot}$). The middle row indicates that baryons have a moderate effect on
the dark matter $V_{\rm max,DM}$ and $R_{\rm max,DM}$. This is because the NFW
halos are already cuspy.

For the SIDM models, results are color-coded based on the halos’ gravothermal
phases at $z=0$. In both models, baryonic effects are pronounced: data
distribution spreads wider, and halos hosting massive galaxies can exhibit
core-collapsing phases characterized by decreased $R_{\rm max,DM}$ and
slightly increased $V_{\rm max,DM}$. We also explore novel correlations
between SIDM model parameters and halo mass scales. In the constant cross
section model ($\sigma/m=10~{}\rm cm^{2}/g$), a significant number of halos
are in core-collapsing phases at $V_{\rm max}\gtrsim 300~{}\rm km/s$.
Conversely, the velocity-dependent model shows a higher incidence of core
collapse at lower velocities ($V_{\rm max}\lesssim 60~{}\rm km/s$).

We have provided scripts for this analysis on GitHub, accessible at
https://github.com/DanengYang/parametricSIDM. Readers can employ the codes to
investigate various SIDM models.

###  VII.2 Deciphering SIDM models from gravothermal phases and accretion
histories

Figure 9:  Distributions of maximum circular velocity ($V_{\rm max,DM}$) and
radius at which this velocity occurs ($R_{\rm max,DM}$) for CDM halos, color-
coded by the $\tau$ for the $\sigma_{0}/m=100~{}\rm cm^{2}/g$ and $w=100~{}\rm
km/s$ model obtained through the integral approach. In reality, these
quantities could be uncovered by analyzing rotation curves. The presented
results are actual values that correspond to the analysis in the middle right
panel of Fig. 8, where the halos are uniformly sampled using a flat halo mass
function over the mass range $M_{h}\in[10^{9},10^{14}]\rm\ M_{\odot}$.

In Fig. 8, we have presented predictions from two SIDM models for the
parameters $V_{\rm max}$, $R_{\rm max}$, and the gravothermal phase $\tau$. It
appears that different SIDM models influence these quantities in varied ways.
Assuming these parameters can be uncovered from observations, such as through
fitting rotation curves, identifying the most plausible SIDM model might be
feasible. However, extracting reliable data from observations and generating
precise theoretical predictions pose significant challenges. A notable issue
is that observations of individual galaxies capture only a single moment in
their complex evolutionary histories. Consequently, one must test theoretical
predictions across a broad range of possible accretion histories and models,
which requires substantial effort.

The concept of gravothermal states facilitates contending with this challenge
through a model-independent insight. By rewinding the current phase $\tau$ of
each halo’s profile—modeled using a parametric model with baryons—to zero, we
effectively connect it to a fictitious CDM halo. This framework presents
several advantages when interpreting observational data. First, the
gravothermal phases derived from observations are not influenced by the
galaxies’ evolutionary histories, allowing researchers to focus solely on
fitting the halos’ gravothermal states without concern for the underlying SIDM
model dynamics. Instead, the search for a plausible SIDM model is left to a
post-processing process, where one can focus on the subtle correlations
between SIDM models and accretion histories. Additionally, given our extensive
knowledge about CDM halos, we could leverage them to constrain the parametric
model’s CDM input (fictitious CDM halos), thereby aiding the extraction of
gravothermal phases from data.

In Fig.9, we illustrate a heuristic representation of “uncovered” results,
which should come from data analysis in reality, but are depicted using the
predictions of the velocity-dependent SIDM model with $\sigma_{0}/m=100\rm
cm^{2}/g$ and $w=100~{}\rm km/s$. The figure demonstrates $V_{\rm max,DM}$ and
$R_{\rm max,DM}$ of fictitious CDM halos, color-coded by the phase $\tau$.
These quantities can all be derived from the halos’ current gravothermal
states. The color coding of the phase already illustrates distinguished
features associated with the velocity-dependent cross section exemplified
here, and the reliability and representativity can be further explored by
post-processing the results for individual halos.

Figure 10:  The normalized evolution histories of the stellar mass (red) and
the gravothermal phase (green). An example post-processing is applied to
extract the SIDM model for the snapshot of DM11+baryon simulation at
$\tau=0.8$. The stellar mass and size evolution are modeled as
$M_{s}(z)=M_{s,0}e^{-2z}$ and $r_{H}=r_{H,0}(1+z)^{-0.75}$, respectively. The
evolution of $\tau$ is modeled in the integral approach to reach the value
$0.8$ from the Cored-DZ fit. The case without stellar accretion is shown with
dashed lines for comparison.

To illustrate this post-processing process, we take the snapshot of the
DM11+baryon simulation at $\tau_{0}=0.8$ for a concrete example. This halo has
$\rho_{s}=6.9\times 10^{6}~{}\rm M_{\odot}/kpc^{3}$ and $r_{s}=9.1~{}$kpc.
Assuming the halo formed $10$ Gyr ago, we determine its core collapse time to
be $t_{c}=10~{}{\rm Gyr}/\tau_{0}\approx 12.5~{}$Gyr. The SIDM models capable
of yielding this $t_{c}=12.5$ Gyr depend on the modeling of the baryon
contents. In scenarios without baryons, we use Eq. (5) to calculate the SIDM
cross section, yielding $\sigma_{\rm eff}/m=53~{}\rm cm^{2}/g$. Alternatively,
assuming the baryon potential remains fixed based on present-day observations,
we utilize Eq. (III.2) with $\rho_{H}=3.6\times 10^{8}\rm M_{\odot}/kpc^{3}$
and $r_{H}=0.77~{}$kpc. This results in a reduced effective cross section of
$20~{}\rm cm^{2}/g$.

To examine how accretion history influences the interpretation, we analyze the
influence of a toy accretion history. Here, we assume the halo to have
established the $\rho_{s}$ and $r_{s}$ since 10 Gyr ago and stellar mass grow
according to $M_{s}(z)=M_{s,0}e^{-2z}$, respectively, with subscript “0”
indicating present day values. The stellar size evolution is modeled as
$r_{H}=r_{H,0}(1+z)^{-0.75}$, as we choose for the sampled halos. To determine
the effective cross section, we solve for

$\tau_{0}=\int_{0}^{10~{}\rm Gyr}\frac{dt}{t_{c,b}[\sigma_{\rm
eff}/m,\rho_{s},r_{s},\rho_{H}(t),r_{H}(t)]},$

where $\rho_{H}$ is derived from $M_{s}$ and $r_{H}$ assuming a Hernquist
profile at each time. The resulting effective cross section is $32~{}\rm
cm^{2}/g$. In Fig. 10, we depict the normalized stellar mass, and gravothermal
phase evolution using red, and green curves, respectively. The case without
stellar accretion is shown with dashed lines for comparison.

Through this analysis, it becomes evident that the inclusion of a baryon
potential significantly influences the derivation of an SIDM model, even when
assuming a constant cross section. Neglecting the baryon effect results in an
SIDM cross section approximately $66$% larger than that from the integral
model. Conversely, when incorporating the current baryon profile without
considering accretion history, the derived SIDM cross section is approximately
38% smaller.

This post-processing approach integrates specific SIDM models with various
accretion histories to model gravothermal evolution. For a given SIDM model,
one can evaluate the predictions using a range of accretion history scenarios
to determine if the uncovered phases are consistently explained under certain
scenarios. If a model outperforms others when paired with well-understood
accretion histories, it becomes a viable SIDM scenario that could explain
observational data and warrants further investigation, especially using
cosmological simulations. Importantly, the post-processing of many halos under
vast possibilities of accretion histories and SIDM models is most effectively
conducted using a parametric approach due to its numerical efficiency.

##  VIII Conclusion and discussion

In this study, we expand upon the parametric model introduced in Ref. Yang
_et al._ (2024a) to encompass the influence of baryons. Our model is based on
two key elements: a formula for calculating core collapse time (Eq. (III.2))
and a density profile model (Cored-DZ) calibrated through controlled N-body
simulations. These elements successfully incorporated the boost of
gravothermal evolution and adiabatic contraction induced by the presence of
baryons.

Utilizing this model, we investigate the interplay between SIDM and baryonic
distributions within the $\rho_{H}/\rho_{s}$ versus $r_{H}/r_{s}$ plane. As
depicted in Fig. 3, the influence of baryons on gravothermal evolution can
vary significantly; halos in the upper-right region, characterized by larger
$\rho_{H}/\rho_{s}$ and $r_{H}/r_{s}$ ratios, experience substantial boost in
the evolution, while those in the lower-left region remain relatively
unaffected. In Fig. 6, we further examine the effect of SIDM on the baryon
distribution by analyzing the ratio of maximum core size to baryon scale
radius $r_{c,\rm max}/r_{H}$. We demonstrate that baryons corresponding to the
lower-left region are more susceptible to gravothermal evolution. Above the
$r_{c,\rm max}=r_{H}$ contour, halos may have suppressed or omitted core
formation behavior, accompanied by contracted baryon profiles. These
interactions could be systematically incorporated into the framework of the
integral approach, which we propose to explore in-depth in future studies.

In this work, we present a simplified application by considering sampled halos
in isolation. We assume that the NFW parameters $\rho_{s}$ and $r_{s}$ have
remained unchanged in the sampled halos over the past $10$ Gyr. Additionally,
we implement simplified accretion histories for the stellar components of all
candidates. This method allows us to concentrate on the influence of baryons
on gravothermal evolution. It also serves as a reasonable approximation of the
behavior observed in isolated halos within cosmological simulations. During
the hierarchical structure formation, the inner configurations of isolated
halos, characterized by the $\rho_{s}$ and $r_{s}$, are established early and
remain relatively stable in the context of CDM Ludlow  _et al._ (2013); Correa
_et al._ (2015). See Fig. 4 of Ref. Yang  _et al._ (2024a) for a visual
representation. Conversely, the baryonic component can exhibit more pronounced
growth, particularly in recent epochs. As depicted in Fig. 7 of Ref. Behroozi
_et al._ (2013), halos within the mass range $10^{11}-10^{12}~{}\rm M_{\odot}$
demonstrate an increasing stellar-to-halo mass ratio ($M_{s}/M_{h}$) from
redshift two until the present day. For more massive halos, the evolution of
$M_{s}/M_{h}$ becomes more moderate, reflecting varied impacts from baryonic
components. Our simplified model features the growth of baryons atop
stationary CDM halos, which is more suitable for halos in the
$10^{11}-10^{12}~{}\rm M_{\odot}$ range.

We also illustrate the effectiveness of the parametric model with baryons in
probing SIDM models. By assuming that each halo resides in a specific
gravothermal phase, retrievable through observational analysis, we illustrate
how the SIDM interpretations depend on the accretion histories and SIDM
models. Exploring the SIDM models across a population of galaxies could enable
the identification of SIDM models that consistently explain all the
observations.

Crucially, the assumption that each SIDM halo is in a gravothermal state
defined by $\tau$ operates independently of the parametric model. This
assumption allows for the use of controlled N-body simulations to accurately
model any halo’s state within a cosmological simulation, as long as their
$\tau$ phases match. For dark matter-only halos, the validity of this
assumption is corroborated by comparisons between the model and simulations,
as detailed in Refs. Yang  _et al._ (2024a) and Yang  _et al._ (2024b), the
latter of which examines the relative differences on a case-by-case basis. In
scenarios involving baryons, however, verification via N-body simulations
remains to be shown. Nonetheless, as shown in Appendix B, the reliability of
this assumption supports theoretical explorations. Conducting dedicated N-body
tests of this assumption would facilitate obtaining accurate theoretical
predictions in future studies.

Applications illustrated in this work are mostly proof-of-principle. Despite
their simplified nature, they provide valuable insights into the complex
interactions between SIDM halos and baryonic distributions. These preliminary
findings underscore the potential of further refining parametric models to
deepen our understanding of the self-interacting characteristics of dark
matter. Continued advancements in this direction can enhance our ability to
decipher the underlying dynamics and implications of SIDM in cosmological
contexts.

## Acknowledgments

I would like to thank Hai-Bo Yu, Ethan O. Nadler, Yi-Ming Zhong, Shin’ichiro
Ando, and Demao Kong, for their helpful discussions. Special appreciation goes
to Hai-Bo Yu for helpful comments on fitting rotation curves and SIDM’s
feedback on baryons. I also acknowledge the IllustrisTNG collaborations for
making their simulation data publicly available. 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. The opinions expressed in this
publication are those of the authors and do not necessarily reflect the views
of the John Templeton Foundation.

##  Appendix A The DM-only limit of the Cored-DZ model

Figure 11:  The inner density evolution (top) and the DM density profiles
(bottom) based on the parametric model for the $\beta 4$ profile (dotted) and
the Cored-DZ profile (dashed).

As a validation, we show in this appendix that the Cored-DZ model reduces to
the $\beta 4$ model in the DM-only limit for a large range in the density
profile. For this purpose, we follow the equations in Ref. Freundlich  _et
al._ (2020) to compute the $a$ and $c$ parameters of the DZ profile through
two other parameters $s_{1}$ and $c_{2}$ as

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{1.5s_{1}-2\left(3.5-s_{1}\right)\left(r_{1}/R_{\mathrm{cut}}\right)^{1/2}c_{2}^{1/2}}{1.5-\left(3.5-s_{1}\right)\left(r_{1}/R_{\mathrm{cut}}\right)^{1/2}c_{2}^{1/2}}$
(24) $\displaystyle c$ $\displaystyle=$
$\displaystyle\left(\frac{s_{1}-2}{\left(3.5-s_{1}\right)\left(r_{1}/R_{\mathrm{cut}}\right)^{1/2}-1.5c_{2}^{-1/2}}\right)^{2}$

where $s_{1}$ represents the logarithmic density slope at $0.01R_{\rm cut}$
and $c_{2}=R_{\rm cut}/r_{-2}$ is determined by the radius $r_{-2}$, at which
the logarithmic density slope equals $-2$.

For DM-only halos, evaluate

$\displaystyle s_{1,\rm DMO}$ $\displaystyle=$
$\displaystyle-\left(r_{1}^{2}\rho_{s}^{-1}r_{s}^{-1}\left(-2\rho_{s}r_{1}^{-1}(1+r_{1}/r_{s})^{-3}\right.\right.$
$\displaystyle\left.\left.-\rho_{s}r_{s}r_{1}^{-2}(1+r_{1}/r_{s})^{-2}\right)(1+r_{1}/r_{s})^{2}\right),$
$\displaystyle c_{\rm DMO}$ $\displaystyle=$ $\displaystyle R_{\rm
cut}/r_{s},$

using the NFW parameters $\rho_{s}$, $r_{s}$ and halo mass $M_{h}$ computed
within $R_{\rm cut}$ yields a DZ profile that largely mimics the NFW profile
in the range $0.05r_{s}$ to $5r_{s}$. See the red curves in Fig. 11 for a
comparison of the two.

##  Appendix B Testing the hybrid approach

Figure 12:  Distributions of maximum circular velocity ($V_{\rm max,DM}$) and
radius at which this velocity occurs ($R_{\rm max,DM}$) for sampled halos,
considering the dark matter component under the influence of baryons. The
upper panels display results from the hybrid approach, while the lower panels
provide results from the integral approach. The predictions for three SIDM
models are provided. From left to right, the results correspond to cross
sections per mass of $\sigma/m=1~{}\rm cm^{2}/g$, $\sigma/m=10~{}\rm
cm^{2}/g$, and a velocity dependent cross section with $\sigma_{0}/m=100~{}\rm
cm^{2}/g$ and $w=100~{}\rm km/s$, respectively. All halos were uniformly
sampled using a flat halo mass function across the mass range
$M_{h}\in[10^{9},10^{14}]\rm\ M_{\odot}$.

As discussed in the main text, the hybrid approach can significantly reduce
the computational complexity. Here, we present three examples to compare the
performances of the hybrid and the integral approaches: two with constant
cross sections per mass, $\sigma/m=1~{}\rm cm^{2}/g$, $\sigma/m=10~{}\rm
cm^{2}/g$, and the other with a velocity dependent cross section specified
with $\sigma_{0}/m=100~{}\rm cm^{2}/g$ and $w=100~{}\rm km/s$.

Figure 12 and Fig. 13 illustrate the distributions of maximum circular
velocity and radius at which it occurs, for the dark matter-only contribution
and the dark matter plus baryon contributions, respectively. These figures
demonstrate that the hybrid approach effectively captures the overall
characteristics of SIDM, particularly the spread and how it diverges from the
CDM scenario. While there are notable differences in the detailed shapes,
these are largely mitigated when the baryon contributions are considered in
the rotation curves. Depending on the problem’s complexity and the required
level of precision, readers may select the appropriate method to employ.

##  Appendix C Adiabatic contraction in different models

Figure 13:  Same as Fig. 12, but plot the total circular velocity ($V_{\rm
max}$) and radius at which this velocity occurs ($R_{\rm max}$), considering
both the dark matter and baryon contributions. The upper panels display
results from the hybrid approach, while the lower panels provide results from
the integral approach. The predictions for three SIDM models are provided.
From left to right, the results correspond to cross sections per mass of
$\sigma/m=1~{}\rm cm^{2}/g$, $\sigma/m=10~{}\rm cm^{2}/g$, and a velocity
dependent cross section with $\sigma_{0}/m=100~{}\rm cm^{2}/g$ and
$w=100~{}\rm km/s$, respectively. All halos were uniformly sampled using a
flat halo mass function across the mass range $M_{h}\in[10^{9},10^{14}]\rm\
M_{\odot}$.

In this appendix, we perform a comparative analysis of three models that
address baryon-induced adiabatic contraction in different ways. Alongside the
Cored-DZ and Contracted $\beta 4$ models, we examine an additional model that
incorporates supplementary contraction effects on top of the Contracted $\beta
4$ model, drawing inspiration from the approach introduced in Ref. Gnedin  _et
al._ (2004).

By calibrating the Contracted $\beta 4$ model to N-body simulations, the
effect of contraction is incorporated to some levels. However, the capability
of the modeling is limited by the parametric form of the density profile. In
Fig. 15, we present the Contracted $\beta 4$ model results in a format
consistent with that of Fig. 5, reproduced in Fig. 14 for ease of comparison.
In Fig. 15, a comparison between the model-predicted curves and their
simulated counterparts reveals a diminishing agreement towards inner halo
regions of halos equal or heavier than $10^{12}~{}\rm M_{\odot}$, indicating
inadequacies in modeling the contraction effect.

Given this result, one could hypothesize that the variation from the
Contracted $\beta 4$ profile is fully induced by the part of the contraction
effect that is independent of SIDM. The SIDM relevant parts are effectively
accounted for in the Contracted $\beta 4$ model, including the core collapse
time, core size, and the transition from core to outer NFW profile.

To assess the degree of such a separation, we adopt the method outlined in
Ref. Gnedin  _et al._ (2004), incorporating an additional baryonic contraction
effect into the Contracted $\beta 4$ profile. This method employs the
adiabatic invariant $M(\bar{r})r$, where $\bar{r}=0.85(r/R_{\rm
200})^{0.8}R_{\rm 200}$ represents an orbitally averaged radius. This radius
is calibrated against hydrodynamical simulations in Ref. Gnedin  _et al._
(2004). The choice between using $R_{\rm 200}$ and $R_{\rm vir}$ for
calculating $\bar{r}$ has a negligible impact for our purposes, and we
consistently use $R_{\rm 200}$ in this analysis. Starting with a tabulated
initial total mass profile $M_{i}(r)$, we calculate the contracted radial
position $r_{f,j}$ using the formula

$\frac{r_{j}}{r_{f,j}}=1-f_{b}+\frac{M_{b}(\bar{r}_{f,j})}{M_{i}(\bar{r}_{j})},$

where $r_{j}$ represents each radial position, $f_{b}=M_{b}/M_{\rm 200}$ is
the ratio of baryon mass to total halo mass, and $M_{b}(r)$ is the baryon mass
profile modeled using the Hernquist profile. The contracted dark matter mass
profile at each $r_{f,j}$ is then calculated as
$M_{h}(r_{f,j})=(1-f_{b})M_{i}(r_{j})$.

The results are presented in Fig. 16. Upon comparison with Fig. 15, a notable
enhancement in agreement is observed, particularly evident in the inner
regions for halo masses exceeding $10^{12}~{}\rm M_{\odot}$. However, for
halos with masses $10^{11}~{}\rm M_{\odot}$, the performance of the Contracted
$\beta 4$ model is already satisfactory, and the additional contraction
appears to overcompensate. Notably, the over-estimation around the core
transitions in the DM13+baryon2 and DM13 extreme cases persists after
accounting for the inner contraction in Fig. 16. Therefore, the presumed
separation does not work well, and the Cored-DZ model outperforms the other
two models.

Figure 14:  The same figure as Fig. 5, duplicated here for an easier
comparison of the other cases. The simulated (colored bands) and model
predicted (colored curves) halo density profiles are presented at three
representative gravothermal phases: $t/t_{c}\approx 0,0.2$, and $1$. The DM12
and DM13 scenarios use a contracted CDM profile as the initial condition,
whereas the DM11 scenarios commence with an instant insertion method. In the
left panel ($t/t_{c}\approx 0$), the DM11 cases are depicted at $t=0.25$ Gyr
to allow some initial evolution away from the original NFW profile. At
$t/t_{c}\approx 1$, the core collapse time, as calculated using Eq. (III.2),
is found to be 10% (30%) shorter than the simulated DM13+baryon2 (DM13
extreme). To align the profiles for equivalent gravothermal phases, we adjust
the timing of the simulated curves accordingly in these specific cases.

Figure 15:  Similar to Fig. 14, but with the colored curves now representing
predictions from the $\beta 4$ model.

Figure 16:  Similar to Fig.14, but here the colored curves depict predictions
from the modified $\beta 4$ model, which now includes an additional
contraction effect following the method introduced in Ref.Gnedin  _et al._
(2004). This enhancement in the model leads to improved alignment between the
model predictions and simulation results, except for the DM11 scenarios at
$t/t_{c}\approx 0.2$, where the additional contraction effects overlap with
adjustments made in Eq. (16). Despite these improvements, the overall quality
of the model-to-simulation comparison only occasionally matches that of the
Cored-DZ model.

## References

  * Yang  _et al._ (2024a) Daneng Yang, Ethan O. Nadler, Hai-Bo Yu, and Yi-Ming Zhong, “A parametric model for self-interacting dark matter halos,” JCAP 02, 032 (2024a), arXiv:2305.16176 [astro-ph.CO] .
  * Flores and Primack (1994) Ricardo A. Flores and Joel R. Primack, “Observational and Theoretical Constraints on Singular Dark Matter Halos,” ApJ Letters 427, L1 (1994), arXiv:astro-ph/9402004 [astro-ph] .
  * Spergel and Steinhardt (2000) David N. Spergel and Paul J. Steinhardt, “Observational Evidence for Self-Interacting Cold Dark Matter,” Phys. Rev. Lett. 84, 3760–3763 (2000), arXiv:astro-ph/9909386 [astro-ph] .
  * de Blok  _et al._ (2001) W. J. G. de Blok, Stacy S. McGaugh, Albert Bosma, and Vera C. Rubin, “Mass density profiles of LSB galaxies,” Astrophys. J. Lett. 552, L23–L26 (2001), arXiv:astro-ph/0103102 .
  * Kuzio de Naray  _et al._ (2010) Rachel Kuzio de Naray, Gregory D. Martinez, James S. Bullock, and Manoj Kaplinghat, “The Case Against Warm or Self-Interacting Dark Matter as Explanations for Cores in Low Surface Brightness Galaxies,” ApJ Letters 710, L161–L166 (2010), arXiv:0912.3518 [astro-ph.CO] .
  * Oh  _et al._ (2011) Se-Heon Oh, W. J. G. de Blok, Elias Brinks, Fabian Walter, and Jr. Kennicutt, Robert C., “Dark and Luminous Matter in THINGS Dwarf Galaxies,” Astrophys. J. 141, 193 (2011), arXiv:1011.0899 [astro-ph.CO] .
  * Oman  _et al._ (2015) Kyle A. Oman, Julio F. Navarro, Azadeh Fattahi, Carlos S. Frenk, Till Sawala, Simon D. M. White, Richard Bower, Robert A. Crain, Michelle Furlong, Matthieu Schaller, Joop Schaye, and Tom Theuns, “The unexpected diversity of dwarf galaxy rotation curves,” MNRAS 452, 3650–3665 (2015), arXiv:1504.01437 [astro-ph.GA] .
  * Navarro  _et al._ (1996) Julio F. Navarro, Vincent R. Eke, and Carlos S. Frenk, “The cores of dwarf galaxy halos,” Mon. Not. Roy. Astron. Soc. 283, L72–L78 (1996), arXiv:astro-ph/9610187 .
  * Gnedin and Zhao (2002) Oleg Y. Gnedin and HongSheng Zhao, “Maximum feedback and dark matter profiles of dwarf galaxies,” Mon. Not. Roy. Astron. Soc. 333, 299 (2002), arXiv:astro-ph/0108108 .
  * Read and Gilmore (2005) Justin I. Read and G. Gilmore, “Mass loss from dwarf spheroidal galaxies: The Origins of shallow dark matter cores and exponential surface brightness profiles,” Mon. Not. Roy. Astron. Soc. 356, 107–124 (2005), arXiv:astro-ph/0409565 .
  * Governato  _et al._ (2010) F. Governato, C. Brook, L. Mayer, A. Brooks, G. Rhee, J. Wadsley, P. Jonsson, B. Willman, G. Stinson, T. Quinn, and P. Madau, “Bulgeless dwarf galaxies and dark matter cores from supernova-driven outflows,” Nature (London) 463, 203–206 (2010), arXiv:0911.2237 [astro-ph.CO] .
  * Di Cintio  _et al._ (2014) Arianna Di Cintio, Chris B. Brook, Aaron A. Dutton, Andrea V. Macciò, Greg S. Stinson, and Alexander Knebe, “A mass-dependent density profile for dark matter haloes including the influence of galaxy formation,” Mon. Not. Roy. Astron. Soc. 441, 2986–2995 (2014), arXiv:1404.5959 [astro-ph.CO] .
  * Chan  _et al._ (2015) T. K. Chan, D. Kereš, J. Oñorbe, P. F. Hopkins, A. L. Muratov, C. A. Faucher-Giguère, and E. Quataert, “The impact of baryonic physics on the structure of dark matter haloes: the view from the FIRE cosmological simulations,” Mon. Not. Roy. Astron. Soc. 454, 2981–3001 (2015), arXiv:1507.02282 [astro-ph.GA] .
  * Tollet  _et al._ (2016) Edouard Tollet, Andrea V. Macciò, Aaron A. Dutton, Greg S. Stinson, Liang Wang, Camilla Penzo, Thales A. Gutcke, Tobias Buck, Xi Kang, Chris Brook, Arianna Di Cintio, Ben W. Keller, and James Wadsley, “NIHAO - IV: core creation and destruction in dark matter density profiles across cosmic time,” MNRAS 456, 3542–3552 (2016), arXiv:1507.03590 [astro-ph.GA] .
  * Lazar  _et al._ (2020) Alexandres Lazar, James S. Bullock, Michael Boylan-Kolchin, T. K. Chan, Philip F. Hopkins, Andrew S. Graus, Andrew Wetzel, Kareem El-Badry, Coral Wheeler, Maria C. Straight, Dušan Kereš, Claude-André Faucher-Giguère, Alex Fitts, and Shea Garrison-Kimmel, “A dark matter profile to model diverse feedback-induced core sizes of $\Lambda$CDM haloes,” MNRAS 497, 2393–2417 (2020), arXiv:2004.10817 [astro-ph.GA] .
  * Macciò  _et al._ (2020) Andrea V. Macciò, Samuele Crespi, Marvin Blank, and Xi Kang, “NIHAO – XXIII. Dark matter density shaped by black hole feedback,” Mon. Not. Roy. Astron. Soc. 495, L46–L50 (2020), arXiv:2004.03817 [astro-ph.GA] .
  * Santos-Santos  _et al._ (2020a) Isabel M. E. Santos-Santos, Julio F. Navarro, Andrew Robertson, Alejandro Benítez-Llambay, Kyle A. Oman, Mark R. Lovell, Carlos S. Frenk, Aaron D. Ludlow, Azadeh Fattahi, and Adam Ritz, “Baryonic clues to the puzzling diversity of dwarf galaxy rotation curves,” MNRAS 495, 58–77 (2020a), arXiv:1911.09116 [astro-ph.GA] .
  * Tulin and Yu (2018) Sean Tulin and Hai-Bo Yu, “Dark matter self-interactions and small scale structure,” Physics Reports 730, 1–57 (2018), arXiv:1705.02358 [hep-ph] .
  * Adhikari  _et al._ (2022) Susmita Adhikari, Arka Banerjee, Kimberly K. Boddy, Francis-Yan Cyr-Racine, Harry Desmond, Cora Dvorkin, Bhuvnesh Jain, Felix Kahlhoefer, Manoj Kaplinghat, Anna Nierenberg, Annika H. G. Peter, Andrew Robertson, Jeremy Sakstein, and Jesús Zavala, “Astrophysical Tests of Dark Matter Self-Interactions,” arXiv e-prints , arXiv:2207.10638 (2022), arXiv:2207.10638 [astro-ph.CO] .
  * Spergel and Steinhardt (2000) David N. Spergel and Paul J. Steinhardt, “Observational evidence for selfinteracting cold dark matter,” Phys. Rev. Lett. 84, 3760–3763 (2000), arXiv:astro-ph/9909386 .
  * Kochanek and White (2000) C. S. Kochanek and Martin J. White, “A Quantitative study of interacting dark matter in halos,” Astrophys. J. 543, 514 (2000), arXiv:astro-ph/0003483 .
  * Kamada  _et al._ (2017) Ayuki Kamada, Manoj Kaplinghat, Andrew B. Pace, and Hai-Bo Yu, “How the Self-Interacting Dark Matter Model Explains the Diverse Galactic Rotation Curves,” Phys. Rev. Lett. 119, 111102 (2017), arXiv:1611.02716 [astro-ph.GA] .
  * Ren  _et al._ (2019) Tao Ren, Anna Kwa, Manoj Kaplinghat, and Hai-Bo Yu, “Reconciling the Diversity and Uniformity of Galactic Rotation Curves with Self-Interacting Dark Matter,” Physical Review X 9, 031020 (2019), arXiv:1808.05695 [astro-ph.GA] .
  * Santos-Santos  _et al._ (2020b) Isabel M. E. Santos-Santos, Julio F. Navarro, Andrew Robertson, Alejandro Benítez-Llambay, Kyle A. Oman, Mark R. Lovell, Carlos S. Frenk, Aaron D. Ludlow, Azadeh Fattahi, and Adam Ritz, “Baryonic clues to the puzzling diversity of dwarf galaxy rotation curves,” MNRAS 495, 58–77 (2020b), arXiv:1911.09116 [astro-ph.GA] .
  * Correa  _et al._ (2022) Camila A. Correa, Matthieu Schaller, Sylvia Ploeckinger, Noemi Anau Montel, Christoph Weniger, and Shinichiro Ando, “TangoSIDM: Tantalizing models of Self-Interacting Dark Matter,” Mon. Not. Roy. Astron. Soc. 517, 3045–3063 (2022), arXiv:2206.11298 [astro-ph.GA] .
  * Yang  _et al._ (2023a) Daneng Yang, Ethan O. Nadler, and Hai-Bo Yu, “Strong Dark Matter Self-interactions Diversify Halo Populations within and surrounding the Milky Way,” Astrophys. J. 949, 67 (2023a), arXiv:2211.13768 [astro-ph.GA] .
  * Balberg  _et al._ (2002) Shmuel Balberg, Stuart L. Shapiro, and Shogo Inagaki, “Self-Interacting Dark Matter Halos and the Gravothermal Catastrophe,” Astrophys. J. 568, 475–487 (2002), arXiv:astro-ph/0110561 [astro-ph] .
  * Balberg  _et al._ (2002) Shmuel Balberg, Stuart L. Shapiro, and Shogo Inagaki, “Selfinteracting dark matter halos and the gravothermal catastrophe,” Astrophys. J. 568, 475–487 (2002), arXiv:astro-ph/0110561 .
  * Koda and Shapiro (2011) Jun Koda and Paul R. Shapiro, “Gravothermal collapse of isolated self-interacting dark matter haloes: N-body simulation versus the fluid model,” MNRAS 415, 1125–1137 (2011), arXiv:1101.3097 [astro-ph.CO] .
  * Nishikawa  _et al._ (2020) Hiroya Nishikawa, Kimberly K. Boddy, and Manoj Kaplinghat, “Accelerated core collapse in tidally stripped self-interacting dark matter halos,” Phys. Rev. D 101, 063009 (2020), arXiv:1901.00499 [astro-ph.GA] .
  * Sameie  _et al._ (2020) Omid Sameie, Hai-Bo Yu, Laura V. Sales, Mark Vogelsberger, and Jesús Zavala, “Self-Interacting Dark Matter Subhalos in the Milky Way’s Tides,” Phys. Rev. Lett. 124, 141102 (2020), arXiv:1904.07872 [astro-ph.GA] .
  * Zavala  _et al._ (2019) Jesús Zavala, Mark R. Lovell, Mark Vogelsberger, and Jan D. Burger, “Diverse dark matter density at sub-kiloparsec scales in Milky Way satellites: Implications for the nature of dark matter,” Phys. Rev. D 100, 063007 (2019), arXiv:1904.09998 [astro-ph.GA] .
  * Kaplinghat  _et al._ (2019) Manoj Kaplinghat, Mauro Valli, and Hai-Bo Yu, “Too Big To Fail in Light of Gaia,” Mon. Not. Roy. Astron. Soc. 490, 231–242 (2019), arXiv:1904.04939 [astro-ph.GA] .
  * Kahlhoefer  _et al._ (2019) Felix Kahlhoefer, Manoj Kaplinghat, Tracy R. Slatyer, and Chih-Liang Wu, “Diversity in density profiles of self-interacting dark matter satellite halos,” JCAP 12, 010 (2019), arXiv:1904.10539 [astro-ph.GA] .
  * Correa (2021) Camila A. Correa, “Constraining velocity-dependent self-interacting dark matter with the Milky Way’s dwarf spheroidal galaxies,” Mon. Not. Roy. Astron. Soc. 503, 920–937 (2021), arXiv:2007.02958 [astro-ph.GA] .
  * Turner  _et al._ (2021) Hannah C. Turner, Mark R. Lovell, Jesús Zavala, and Mark Vogelsberger, “The onset of gravothermal core collapse in velocity-dependent self-interacting dark matter subhaloes,” Mon. Not. Roy. Astron. Soc. 505, 5327–5339 (2021), arXiv:2010.02924 [astro-ph.GA] .
  * Slone  _et al._ (2023) Oren Slone, Fangzhou Jiang, Mariangela Lisanti, and Manoj Kaplinghat, “Orbital evolution of satellite galaxies in self-interacting dark matter models,” Phys. Rev. D 107, 043014 (2023), arXiv:2108.03243 [astro-ph.CO] .
  * Silverman  _et al._ (2022) Maya Silverman, James S. Bullock, Manoj Kaplinghat, Victor H. Robles, and Mauro Valli, “Motivations for a large self-interacting dark matter cross-section from Milky Way satellites,” Mon. Not. Roy. Astron. Soc. 518, 2418–2435 (2022), arXiv:2203.10104 [astro-ph.GA] .
  * Minor  _et al._ (2021) Quinn E. Minor, Sophia Gad-Nasr, Manoj Kaplinghat, and Simona Vegetti, “An unexpected high concentration for the dark substructure in the gravitational lens SDSSJ0946+1006,” Mon. Not. Roy. Astron. Soc. 507, 1662–1683 (2021), arXiv:2011.10627 [astro-ph.GA] .
  * Yang and Yu (2021) Daneng Yang and Hai-Bo Yu, “Self-interacting dark matter and small-scale gravitational lenses in galaxy clusters,” Phys. Rev. D 104, 103031 (2021), arXiv:2102.02375 [astro-ph.GA] .
  * Gilman  _et al._ (2021) Daniel Gilman, Jo Bovy, Tommaso Treu, Anna Nierenberg, Simon Birrer, Andrew Benson, and Omid Sameie, “Strong lensing signatures of self-interacting dark matter in low-mass haloes,” Mon. Not. Roy. Astron. Soc. 507, 2432–2447 (2021), arXiv:2105.05259 [astro-ph.CO] .
  * Gilman  _et al._ (2023) Daniel Gilman, Yi-Ming Zhong, and Jo Bovy, “Constraining resonant dark matter self-interactions with strong gravitational lenses,” Phys. Rev. D 107, 103008 (2023), arXiv:2207.13111 [astro-ph.CO] .
  * Nadler  _et al._ (2023) Ethan O. Nadler, Daneng Yang, and Hai-Bo Yu, “A Self-interacting Dark Matter Solution to the Extreme Diversity of Low-mass Halo Properties,” Astrophys. J. Lett. 958, L39 (2023), arXiv:2306.01830 [astro-ph.GA] .
  * Pollack  _et al._ (2015) Jason Pollack, David N. Spergel, and Paul J. Steinhardt, “Supermassive Black Holes from Ultra-Strongly Self-Interacting Dark Matter,” Astrophys. J. 804, 131 (2015), arXiv:1501.00017 [astro-ph.CO] .
  * Choquette  _et al._ (2019) Jeremie Choquette, James M. Cline, and Jonathan M. Cornell, “Early formation of supermassive black holes via dark matter self-interactions,” JCAP 07, 036 (2019), arXiv:1812.05088 [astro-ph.CO] .
  * Feng  _et al._ (2021) Wei-Xiang Feng, Hai-Bo Yu, and Yi-Ming Zhong, “Seeding Supermassive Black Holes with Self-interacting Dark Matter: A Unified Scenario with Baryons,” Astrophys. J. Lett. 914, L26 (2021), arXiv:2010.15132 [astro-ph.CO] .
  * Feng  _et al._ (2022) Wei-Xiang Feng, Hai-Bo Yu, and Yi-Ming Zhong, “Dynamical instability of collapsed dark matter halos,” JCAP 2022, 036 (2022), arXiv:2108.11967 [astro-ph.CO] .
  * Meshveliani  _et al._ (2023) Tamar Meshveliani, Jesús Zavala, and Mark R. Lovell, “Gravothermal collapse of self-interacting dark matter halos as the origin of intermediate mass black holes in Milky Way satellites,” Phys. Rev. D 107, 083010 (2023), arXiv:2210.01817 [astro-ph.GA] .
  * Fischer  _et al._ (2024) Moritz S. Fischer, Lenard Kasselmann, Marcus Brüggen, Klaus Dolag, Felix Kahlhoefer, Antonio Ragagnin, Andrew Robertson, and Kai Schmidt-Hoberg, “Cosmological and idealized simulations of dark matter haloes with velocity-dependent, rare and frequent self-interactions,” Mon. Not. Roy. Astron. Soc. 529, 2327–2348 (2024), arXiv:2310.07750 [astro-ph.CO] .
  * Dhanasingham  _et al._ (2023) Birendra Dhanasingham, Francis-Yan Cyr-Racine, Charlie Mace, Annika H. G. Peter, and Andrew Benson, “Anisotropic strong lensing as a probe of dark matter self-interactions,” Mon. Not. Roy. Astron. Soc. 526, 5455–5470 (2023), arXiv:2306.10109 [astro-ph.CO] .
  * Zhang  _et al._ (2024) Xingyu Zhang, Hai-Bo Yu, Daneng Yang, and Haipeng An, “Self-interacting dark matter interpretation of Crater II,” (2024), arXiv:2401.04985 [astro-ph.GA] .
  * Vogelsberger  _et al._ (2014) Mark Vogelsberger, Jesus Zavala, Christine Simpson, and Adrian Jenkins, “Dwarf galaxies in CDM and SIDM with baryons: observational probes of the nature of dark matter,” Mon. Not. Roy. Astron. Soc. 444, 3684–3698 (2014), arXiv:1405.5216 [astro-ph.CO] .
  * Robertson  _et al._ (2019) Andrew Robertson, David Harvey, Richard Massey, Vincent Eke, Ian G. McCarthy, Mathilde Jauzac, Baojiu Li, and Joop Schaye, “Observable tests of self-interacting dark matter in galaxy clusters: cosmological simulations with SIDM and baryons,” Mon. Not. Roy. Astron. Soc. 488, 3646–3662 (2019), arXiv:1810.05649 [astro-ph.CO] .
  * Robertson  _et al._ (2021) Andrew Robertson, Richard Massey, Vincent Eke, Joop Schaye, and Tom Theuns, “The surprising accuracy of isothermal Jeans modelling of self-interacting dark matter density profiles,” Mon. Not. Roy. Astron. Soc. 501, 4610–4634 (2021), arXiv:2009.07844 [astro-ph.CO] .
  * Rose  _et al._ (2023) Jonah C. Rose, Paul Torrey, Mark Vogelsberger, and Stephanie O’Neil, “Unravelling the interplay between SIDM and baryons in MW haloes: defining where baryons dictate heat transfer,” Mon. Not. Roy. Astron. Soc. 519, 5623–5636 (2023), arXiv:2206.14830 [astro-ph.GA] .
  * Mastromarino  _et al._ (2023) Claudio Mastromarino, Giulia Despali, Lauro Moscardini, Andrew Robertson, Massimo Meneghetti, and Matteo Maturi, “Properties and observables of massive galaxies in self-interacting dark matter cosmologies,” Mon. Not. Roy. Astron. Soc. 524, 1515–1528 (2023), arXiv:2212.01403 [astro-ph.CO] .
  * Rahimi  _et al._ (2023) Elham Rahimi, Evan Vienneau, Nassim Bozorgnia, and Andrew Robertson, “The local dark matter distribution in self-interacting dark matter halos,” JCAP 02, 040 (2023), arXiv:2210.06498 [astro-ph.CO] .
  * Zhong  _et al._ (2023) Yi-Ming Zhong, Daneng Yang, and Hai-Bo Yu, “The impact of baryonic potentials on the gravothermal evolution of self-interacting dark matter haloes,” Mon. Not. Roy. Astron. Soc. 526, 758–770 (2023), arXiv:2306.08028 [astro-ph.CO] .
  * Creasey  _et al._ (2017) Peter Creasey, Omid Sameie, Laura V. Sales, Hai-Bo Yu, Mark Vogelsberger, and Jesús Zavala, “Spreading out and staying sharp - creating diverse rotation curves via baryonic and self-interaction effects,” MNRAS 468, 2283–2295 (2017), arXiv:1612.03903 [astro-ph.GA] .
  * Fry  _et al._ (2015) A. Bastidas Fry, F. Governato, A. Pontzen, T. Quinn, M. Tremmel, L. Anderson, H. Menon, A. M. Brooks, and J. Wadsley, “All about baryons: revisiting SIDM predictions at small halo masses,” Mon. Not. Roy. Astron. Soc. 452, 1468–1479 (2015), arXiv:1501.00497 [astro-ph.CO] .
  * Kong  _et al._ (2024) Demao Kong, Daneng Yang, and Hai-Bo Yu, “Cold Dark Matter and Self-interacting Dark Matter Interpretations of the Strong Gravitational Lensing Object JWST-ER1,” Astrophys. J. Lett. 965, L19 (2024), arXiv:2402.15840 [astro-ph.GA] .
  * Fischer  _et al._ (2023) Moritz S. Fischer, Nils-Henrik Durke, Katharina Hollingshausen, Claudius Hammer, Marcus Brüggen, and Klaus Dolag, “The role of baryons in self-interacting dark matter mergers,” Mon. Not. Roy. Astron. Soc. 523, 5915–5933 (2023), arXiv:2302.07882 [astro-ph.CO] .
  * Robertson  _et al._ (2018) Andrew Robertson _et al._ , “The diverse density profiles of galaxy clusters with self-interacting dark matter plus baryons,” Mon. Not. Roy. Astron. Soc. 476, L20–L24 (2018), arXiv:1711.09096 [astro-ph.CO] .
  * Robles  _et al._ (2017) Victor H. Robles, James S. Bullock, Oliver D. Elbert, Alex Fitts, Alejandro González-Samaniego, Michael Boylan-Kolchin, Philip F. Hopkins, Claude-André Faucher-Giguère, Dušan Kereš, and Christopher C. Hayward, “SIDM on FIRE: hydrodynamical self-interacting dark matter simulations of low-mass dwarf galaxies,” MNRAS 472, 2945–2954 (2017), arXiv:1706.07514 [astro-ph.GA] .
  * Elbert  _et al._ (2018) Oliver D. Elbert, James S. Bullock, Manoj Kaplinghat, Shea Garrison-Kimmel, Andrew S. Graus, and Miguel Rocha, “A Testable Conspiracy: Simulating Baryonic Effects on Self-interacting Dark Matter Halos,” Astrophys. J. 853, 109 (2018), arXiv:1609.08626 [astro-ph.GA] .
  * Despali  _et al._ (2019) Giulia Despali, Martin Sparre, Simona Vegetti, Mark Vogelsberger, Jesús Zavala, and Federico Marinacci, “The interplay of Self-Interacting Dark Matter and baryons in shaping the halo evolution,” Mon. Not. Roy. Astron. Soc. 484, 4563 (2019), arXiv:1811.02569 [astro-ph.GA] .
  * Essig  _et al._ (2019) Rouven Essig, Samuel D. Mcdermott, Hai-Bo Yu, and Yi-Ming Zhong, “Constraining Dissipative Dark Matter Self-Interactions,” Phys. Rev. Lett. 123, 121102 (2019), arXiv:1809.01144 [hep-ph] .
  * Yang  _et al._ (2023b) Shengqi Yang, Xiaolong Du, Zhichao Carton Zeng, Andrew Benson, Fangzhou Jiang, Ethan O. Nadler, and Annika H. G. Peter, “Gravothermal Solutions of SIDM Halos: Mapping from Constant to Velocity-dependent Cross Section,” Astrophys. J. 946, 47 (2023b), arXiv:2205.02957 [astro-ph.CO] .
  * Jiang  _et al._ (2023) Fangzhou Jiang _et al._ , “A semi-analytic study of self-interacting dark-matter haloes with baryons,” Mon. Not. Roy. Astron. Soc. 521, 4630–4644 (2023), arXiv:2206.12425 [astro-ph.CO] .
  * Yang  _et al._ (2023c) Shengqi Yang, Fangzhou Jiang, Andrew Benson, Yi-Ming Zhong, Charlie Mace, Xiaolong Du, Zhichao Carton Zeng, Annika H. G. Peter, and Moritz S. Fischer, “A quantitative comparison between velocity dependent SIDM cross sections constrained by the gravothermal and isothermal models,” (2023c), arXiv:2305.05067 [astro-ph.CO] .
  * Elbert  _et al._ (2018) Oliver D. Elbert, James S. Bullock, Manoj Kaplinghat, Shea Garrison-Kimmel, Andrew S. Graus, and Miguel Rocha, “A Testable Conspiracy: Simulating Baryonic Effects on Self-Interacting Dark Matter Halos,” Astrophys. J. 853, 109 (2018), arXiv:1609.08626 [astro-ph.GA] .
  * Sameie  _et al._ (2018) Omid Sameie, Peter Creasey, Hai-Bo Yu, Laura V. Sales, Mark Vogelsberger, and Jesus Zavala, “The impact of baryonic discs on the shapes and profiles of self-interacting dark matter haloes,” Mon. Not. Roy. Astron. Soc. 479, 359–367 (2018), arXiv:1801.09682 [astro-ph.GA] .
  * Outmezguine  _et al._ (2023) Nadav Joseph Outmezguine, Kimberly K. Boddy, Sophia Gad-Nasr, Manoj Kaplinghat, and Laura Sagunski, “Universal gravothermal evolution of isolated self-interacting dark matter halos for velocity-dependent cross-sections,” Mon. Not. Roy. Astron. Soc. 523, 4786–4800 (2023), arXiv:2204.06568 [astro-ph.GA] .
  * Yang and Yu (2022) Daneng Yang and Hai-Bo Yu, “Gravothermal evolution of dark matter halos with differential elastic scattering,” JCAP 2022, 077 (2022), arXiv:2205.03392 [astro-ph.CO] .
  * Navarro  _et al._ (1997) Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White, “A Universal Density Profile from Hierarchical Clustering,” Astrophys. J. 490, 493–508 (1997), arXiv:astro-ph/9611107 [astro-ph] .
  * Zhao (1996) HongSheng Zhao, “Analytical models for galactic nuclei,” Mon. Not. Roy. Astron. Soc. 278, 488–496 (1996), arXiv:astro-ph/9509122 .
  * Dekel  _et al._ (2017) Avishai Dekel, Guy Ishai, Aaron A. Dutton, and Andrea V. Maccio, “Dark-matter halo profiles of a general cusp/core with analytic velocity and potential,” MNRAS 468, 1005–1022 (2017), arXiv:1610.00916 [astro-ph.GA] .
  * Freundlich  _et al._ (2020) Jonathan Freundlich, Fangzhou Jiang, Avishai Dekel, Nicolas Cornuault, Omry Ginzburg, Rémy Koskas, Sharon Lapiner, Aaron Dutton, and Andrea V. Macciò, “The Dekel-Zhao profile: a mass-dependent dark-matter density profile with flexible inner slope and analytic potential, velocity dispersion, and lensing properties,” MNRAS 499, 2912–2933 (2020), arXiv:2004.08395 [astro-ph.GA] .
  * Nelson  _et al._ (2019a) Dylan Nelson _et al._ , “The IllustrisTNG simulations: public data release,” Comput. Astrophys. Cosmol. 6, 2 (2019a), arXiv:1812.05609 [astro-ph.GA] .
  * Nelson  _et al._ (2019b) Dylan Nelson, Annalisa Pillepich, Volker Springel, Ruediger Pakmor, Rainer Weinberger, Shy Genel, Paul Torrey, Mark Vogelsberger, Federico Marinacci, and Lars Hernquist, “First Results from the TNG50 Simulation: Galactic outflows driven by supernovae and black hole feedback,” Mon. Not. Roy. Astron. Soc. 490, 3234–3261 (2019b), arXiv:1902.05554 [astro-ph.GA] .
  * Pillepich  _et al._ (2019) Annalisa Pillepich _et al._ , “First results from the TNG50 simulation: the evolution of stellar and gaseous discs across cosmic time,” Mon. Not. Roy. Astron. Soc. 490, 3196–3233 (2019), arXiv:1902.05553 [astro-ph.GA] .
  * Hernquist (1990) Lars Hernquist, “An Analytical Model for Spherical Galaxies and Bulges,” Astrophys. J. 356, 359 (1990).
  * Lynden-Bell and Eggleton (1980) D. Lynden-Bell and P. P. Eggleton, “On the consequences of the gravothermal catastrophe,” Mon. Not. Roy. Astron. Soc. 191, 483–498 (1980).
  * Yang  _et al._ (2024b) Daneng Yang, Ethan O. Nadler, and Hai-Bo Yu, “Testing the parametric self-interacting dark matter model using matched halos in cosmological simulations,” (2024b), arXiv:24XX.XXXX [astro-ph.GA] .
  * Gnedin  _et al._ (2004) Oleg Y. Gnedin, Andrey V. Kravtsov, Anatoly A. Klypin, and Daisuke Nagai, “Response of dark matter halos to condensation of baryons: Cosmological simulations and improved adiabatic contraction model,” Astrophys. J. 616, 16–26 (2004), arXiv:astro-ph/0406247 .
  * Ryden and Gunn (1987) Barbara S. Ryden and James E. Gunn, “Galaxy Formation by Gravitational Collapse,” Astrophys. J. 318, 15 (1987).
  * Ando  _et al._ (2024) Shin’ichiro Ando, Shunichi Horigome, Ethan O. Nadler, Daneng Yang, and Hai-Bo Yu, “SASHIMI-SIDM: Semi-analytical subhalo modelling for self-interacting dark matter at sub-galactic scales,” (2024), arXiv:2403.16633 [astro-ph.CO] .
  * Dutton and Macciò (2014) Aaron A. Dutton and Andrea V. Macciò, “Cold dark matter haloes in the Planck era: evolution of structural parameters for Einasto and NFW profiles,” Mon. Not. Roy. Astron. Soc. 441, 3359–3374 (2014), arXiv:1402.7073 [astro-ph.CO] .
  * Moster  _et al._ (2013) Benjamin P. Moster, Thorsten Naab, and Simon D. M. White, “Galactic star formation and accretion histories from matching galaxies to dark matter haloes,” MNRAS 428, 3121–3138 (2013), arXiv:1205.5807 [astro-ph.CO] .
  * Carleton  _et al._ (2019) Timothy Carleton, Raphaël Errani, Michael Cooper, Manoj Kaplinghat, Jorge Peñarrubia, and Yicheng Guo, “The formation of ultra-diffuse galaxies in cored dark matter haloes through tidal stripping and heating,” MNRAS 485, 382–395 (2019), arXiv:1805.06896 [astro-ph.GA] .
  * Behroozi  _et al._ (2013) Peter S. Behroozi, Risa H. Wechsler, and Charlie Conroy, “The Average Star Formation Histories of Galaxies in Dark Matter Halos from z = 0-8,” Astrophys. J. 770, 57 (2013), arXiv:1207.6105 [astro-ph.CO] .
  * van der Wel  _et al._ (2014) A. van der Wel _et al._ , “3D-HST+CANDELS: The Evolution of the Galaxy Size-Mass Distribution since $z=3$,” Astrophys. J. 788, 28 (2014), arXiv:1404.2844 [astro-ph.GA] .
  * Ludlow  _et al._ (2013) Aaron D. Ludlow, Julio F. Navarro, Michael Boylan-Kolchin, Philip E. Bett, Raul E. Angulo, Ming Li, Simon D. M. White, Carlos Frenk, and Volker Springel, “The Mass Profile and Accretion History of Cold Dark Matter Halos,” Mon. Not. Roy. Astron. Soc. 432, 1103 (2013), arXiv:1302.0288 [astro-ph.CO] .
  * Correa  _et al._ (2015) Camila A. Correa, J. Stuart B. Wyithe, Joop Schaye, and Alan R. Duffy, “The accretion history of dark matter haloes – III. A physical model for the concentration–mass relation,” Mon. Not. Roy. Astron. Soc. 452, 1217–1232 (2015), arXiv:1502.00391 [astro-ph.CO] .