Gravothermal Solutions of SIDM Halos: Mapping from Constant to Velocity-dependent Cross Section

Self-interacting dark matter (SIDM) is a promising class of candidate to explain the cored dark matter (DM) halos implied by observations of local dwarf galaxies (de Blok et al. 2008; Oh et al. 2015), while maintaining the accurate reproduction of large-scale cosmic properties provided by the cold dark matter (CDM) paradigm (Peebles 1982; Planck Collaboration et al. 2011, 2014; Anderson et al. 2012). In the SIDM model, first proposed by Spergel & Steinhardt (2000), in addition to the gravitational interactions, DM particles also scatter with each other. The effect of scattering is negligible on large scales due to the low DM particle number density and the resulting long relaxation time. However, the collisional mean free path of DM particles can become smaller than the orbital Jeans scale near the halo center, causing the formation of an isothermal core and even the gravitational catastrophe (Lynden-Bell et al. 1968) at small scales.

The gravothermal fluid formalism originally developed for studying the gravitational catastrophe of stellar clusters (Lynden-Bell et al. 1968; Eggleton 1971) has been extended to isolated SIDM halos and provides rich physical insights into the nature of SIDM halo evolution (e.g., Lynden-Bell & Eggleton 1980; Balberg et al. 2002; Balberg & Shapiro 2002; Shapiro 2018; Essig et al. 2019; Nishikawa et al. 2020). Although it makes multiple simplifying assumptions, numerical solutions to the gravothermal fluid formalism show surprisingly good agreement with idealized N-body simulations and also match reasonably well with cosmological SIDM zoom-in simulations (Koda & Shapiro 2011; Elbert et al. 2015; Essig et al. 2019; Outmezguine et al. 2022). A nice property of the gravothermal fluid formalism is that the fluid equations themselves can be rewritten into a form that is independent of the halo density and radius scales. Therefore, although computationally expensive, the fluid formalism only needs to be numerically solved once and can be rescaled to SIDM halos with arbitrary masses, sizes, and concentrations.

All previous comparisons between the gravothermal solutions and numerical simulations have assumed a constant DM particle scattering cross section. However, if SIDM is the explanation for the apparent deviations from CDM behavior on small scales, then recent observations indicate that the cross section, σ, for DM particle self-interaction is likely varying among systems of different scales. Specifically, σ is constrained to be of the order of 0.1 cm² g^–1, 1 cm² g^–1, and 10 cm² g^–1 for galaxy cluster, galactic, and dwarf-galaxy halos, respectively, where the velocity dispersion of DM particles decreases from ∼1000 to ∼1 km s⁻¹ (e.g., Kaplinghat et al. 2016; Elbert et al. 2018; Sagunski et al. 2021; Andrade et al. 2022). A velocity-dependent cross section σ ∝ v⁻⁴ at high velocity is also motivated by particle physics, assuming a Yukawa or Coulomb potential for the self-interaction (Feng et al. 2009; Ackerman et al. 2009; Ibe & Yu 2010; Kummer et al. 2018). Motivated by this, more recent SIDM numerical simulations have started to focus on velocity-dependent cross-section models (e.g., Nadler et al. 2020). Testing whether the gravothermal fluid solutions are still in good agreement with numerical simulations assuming a velocity-dependent cross-section model becomes crucial in developing better understanding of the SIDM N-body simulations, as well as justifying SIDM cross-section constraints derived by gravothermal fluid formalism-based methods.

In this work, we show using an analytic gravothermal fluid formalism that assuming different cross-section values/models effectively rescales SIDM halo evolution rates. In particular, there exists a one-on-one time mapping between any two SIDM halos simulated by the gravothermal fluid approach such that their density, velocity, and all the other physical properties are identical up to some rescaling factors. Based on this finding, we propose a fast mapping method that transfers the gravothermal numerical solutions assuming a constant SIDM cross section to those with arbitrary velocity-dependent cross-section models. To test the performance of this mapping method we perform multiple idealized N-body simulations assuming velocity-dependent cross-section scenarios with Arepo (Springel 2010). We find excellent agreement between the mapped gravothermal solutions and N-body simulation results. This mapping approach is powerful because it eliminates the degrees of freedom introduced by assumptions in the SIDM cross-section model, such that one only needs to solve the gravothermal fluid equations once and one can then rescale the solution to halos of arbitrary sizes, masses, and cross-section models. Due to the high computational efficiency, this method is useful for exploring the continuous parameter space and constraining SIDM cross sections. It is also applicable to semianalytic galaxy formation models to efficiently generate SIDM halo populations. To demonstrate the utility of this approach, we combine the mapping method with rotation curve measurements of the low-surface-brightness (LSB) galaxy UGC 128 (van der Hulst et al. 1993) and illustrate the parameters of the Navarro–Frenk–White (NFW) density profile (Navarro et al. 1996) and velocity-dependent cross-section models favored by observations.

As this work was nearing completion, Outmezguine et al. (2022; hereafter O22) introduced a mapping method with a similar spirit, but which differs in detail from our work. We find that O22 only compared the mapped gravothermal solutions with idealized N-body simulations that assume constant cross sections (Koda & Shapiro 2011). Example applications of their mapping method are also not discussed in detail. We have therefore expanded this work to include comparison of our approach with that of O22.

Calibrating the gravothermal fluid formalism to cosmological zoom-in simulations is an even more interesting topic to explore. However, cosmological simulations generally contain more complicated physics and can break multiple assumptions made in the current fluid model. As a first step, in this work we will only compare gravothermal solutions with idealized SIDM N-body simulations for isolated halos with constant or velocity-dependent cross sections. We leave the more complicated comparisons of gravothermal solution versus cosmological simulation to future works.

The plan of this paper is as follows. In Section 2 we briefly introduce our idealized SIDM halo simulation sets. Section 3 is an introduction to the gravothermal fluid formalism. We briefly review the long-mean-free-path (lmfp) gravothermal solution mapping method among different constant-cross-section scenarios introduced by Balberg et al. (2002), Koda & Shapiro (2011), and Essig et al. (2019) in Section 4. We then explain the methodology of mapping the lmfp gravothermal solutions from constant cross section σ to velocity-dependent cross section σ(v) in Section 5 and comparing the mapped gravothermal solutions with N-body simulation results as well as O22 in Section 6. We briefly discuss the generalization of this lmfp mapping method to the short-mean-free-path (smfp) and intermediate regime in Section 7. In Section 8 we combine the mapping method with observational data and show parameter constraints of the velocity-dependent cross-section models. We conclude in Section 9.

In this work we perform idealized N-body simulations of an isolated Milky Way-sized halo with NFW initial density profile:

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{{\rm{s}}}}{\tfrac{r}{{r}_{{\rm{s}}}}{\left(1+\tfrac{r}{{r}_{{\rm{s}}}}\right)}^{2}},\end{eqnarray} \tag{ 1 }$$

where ρ_s = 4.2 × 10⁶ M_⊙ kpc^–3 and r_s = 24.54 kpc. The halo virial radius and mass, defined with overdensity Δ_vir = 99.2 with respect to the critical density, are r_vir = 280.6 kpc and M_vir = 10^12.1 M_⊙. To guarantee a finite halo mass, we apply an exponential cutoff to the halo density profile at r_vir. Convergence test results show that such an exponential cutoff at r ≳ 10r_s has negligible influence on the gravothermal evolution of the SIDM halo. The properties of this halo are chosen to be consistent with those of the host halo from the CDM zoom-in simulation presented in Nadler et al. (2020), which has also been resimulated in velocity-dependent SIDM models. ⁶

However, our selection of these halo properties is an arbitrary choice for testing the performance of the gravothermal solution mapping method introduced in Section 5. Since the gravothermal fluid formalism is scale-free (i.e., independent of ρ_s and r_s besides some rescaling factors), comparisons between the gravothermal solutions and isolated SIDM halo simulations of arbitrary sizes and masses are equivalent.

To avoid extremely long collapsing timescales for halo cores, and therefore expensive simulations, we select a large constant-cross-section factor σ_c = 30 cm² g^–1 and consider a simple velocity-dependent cross-section model introduced in Gilman et al. (2021),

$$\begin{eqnarray}&&\sigma (v)=\displaystyle \frac{{\sigma }_{c}}{{\left(1+{v}^{2}/{\omega }^{2}\right)}^{2}},\end{eqnarray} \tag{ 2 }$$

where ω is a characteristic velocity beyond which the cross section drops as v⁻⁴. This asymptotic behavior is motivated by Coulomb or Yukawa interactions (Feng et al. 2009; Ackerman et al. 2009; Ibe & Yu 2010; Kummer et al. 2018). The factor 1 in the denominator is used to avoid a divergent cross section for low-velocity particles in the numerical simulations.

We consider four different cross-section models where ω = ∞ (i.e., σ = σ_c = 30 cm² g^–1 is a constant) and ω = 1000, 500, and 400 km s⁻¹. Notice that since the 1D rms velocity of DM particles at the halo center and at the moment of core formation (when the halo central density reaches minimum) is about 100 km s⁻¹, corresponding to a 3D rms relative velocity of about 250 km s⁻¹, those four scenarios correspond to constant cross section, and weak, intermediate, and strong cross-section velocity dependence, respectively.

We use SpherIC (Garrison-Kimmel et al. 2013) to generate five initial conditions with identical resolutions and halo properties, but different random seeds. Through this we can quantify the effects of discreteness noise in the initial conditions on the evolution of the halo in the N-body simulations. We set the particle mass to 10⁶ M_⊙ so the halo in each simulation is resolved with 1.7 × 10⁶ particles. For each cross-section scenario we simulate five halos with the Arepo SIDM module, and we direct the reader to Vogelsberger et al. (2012, 2014, 2019) for a detailed introduction to the numerical method. This SIDM simulation pipeline has been adopted by many previous works (e.g., Zeng et al. 2022) and has been proved to generate valid simulation results for cross sections up to 60 cm² g^–1.

Consider a spherically symmetric and isolated halo. Assuming gravothermal equilibrium is achieved at every radius and every moment, the halo evolution can be described by the following four equations (see also Balberg et al. 2002; Essig et al. 2019):

$$\begin{eqnarray}&&\displaystyle \frac{\partial M}{\partial r}=4\pi {r}^{2}\rho ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {v}_{\mathrm{rms}}^{2})}{\partial r}=-\displaystyle \frac{{GM}\rho }{{r}^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{L}{4\pi {r}^{2}}=-\kappa \displaystyle \frac{\partial T}{\partial r},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\rho {v}_{\mathrm{rms}}^{2}}{\gamma -1}{\left(\displaystyle \frac{\partial }{\partial t}\right)}_{M}\mathrm{ln}\displaystyle \frac{{v}_{\mathrm{rms}}^{2}}{{\rho }^{\gamma -1}}=-\displaystyle \frac{1}{4\pi {r}^{2}}\displaystyle \frac{\partial L}{\partial r}.\end{eqnarray} \tag{ 6 }$$

Throughout this paper we use M(r) to denote halo enclosed mass at radius less than r. Then, ρ(r), v_rms(r), L(r), κ(r), and $T(r)={{mv}}_{\mathrm{rms}}^{2}(r)/{k}_{{\rm{B}}}$ correspond to halo density, 1D rms velocity averaged over the Maxwell–Boltzmann (MB) distribution, luminosity, thermal conductivity, and temperature at radius r respectively. Here, m is the DM particle mass and k_B is the Boltzmann constant. Since we do not consider the movement of the halo center throughout this work, the rms velocity and velocity dispersion among DM particles are identical at all radii. Assuming the DM particles act as a monatomic ideal gas, we set the adiabatic index γ = 5/3. The Jeans equation (Equation (4)) forces gravothermal equilibrium where the outward pressure $\partial (\rho {v}_{\mathrm{rms}}^{2})/\partial r$ is balanced by the gravitational force. Equation (5) defines the heat flux and Equation (6) expresses the first law of thermodynamics. Together Equations (5) and (6) govern the quasi-static time evolution of an SIDM halo.

For ideal-gas environments where the heat conduction is solely realized by the thermal motion and collisions among gas particles, the kinetic theory of gases gives heat conductivity as

$$\begin{eqnarray}&&{\kappa }_{\mathrm{smfp}}\propto \displaystyle \frac{n{\lambda }^{2}{k}_{{\rm{B}}}}{{t}_{r}},\end{eqnarray} \tag{ 7 }$$

and it is conventional to define κ_smfp = 2.1v_rms k_B/σ, derived from integration over the MB velocity distribution. Here, n = ρ/m is the local DM particle number density, λ = 1/(n σ) is the mean free path, and t_r = λ/v₁₂ is the relaxation time. Here we use σ to denote DM particle scattering cross section, which is either a constant or a velocity-dependent quantity. We use v₁₂ to distinguish particle relative velocity from the 1D rms velocity v_rms used for the temperature definition.

Equation (7) does not apply to gravothermal systems of low density and low self-interaction cross section, where most of the particles tend to orbit several times before colliding with one another. Using the Jeans length $H=\sqrt{{v}_{\mathrm{rms}}^{2}/(4\pi G\rho )}$ to approximate particles' characteristic orbital radii, Equation (7) is only accurate when λ ≪ H (also referred as the smfp regime). In this work we only care about halo evolution processes that satisfy λ ≫ H, also referred as the lmfp regime. To derive the heat conductivity in the lmfp regime one just needs to replace the mean free path λ with the scale height H (Spitzer 1987):

$$\begin{eqnarray}&&{\kappa }_{\mathrm{lmfp}}\propto \displaystyle \frac{{{nH}}^{2}{k}_{{\rm{B}}}}{{t}_{r}}.\end{eqnarray} \tag{ 8 }$$

It is conventional to assume ${\kappa }_{{\rm{lmfp}}}=(3\beta /2){{nH}}^{2}{k}_{{\rm{B}}}/{t}_{r}\,=0.27\beta {{nv}}_{{\rm{rms}}}^{3}\sigma {k}_{{\rm{B}}}/(Gm)$, where β is an adjustable parameter to match the fluid numerical solutions with N-body simulations. The value of β cannot be derived from first principles and it is therefore treated as a free parameter to be calibrated against N-body simulations. It is known that calibration to N-body simulations results in β of order unity, but with some variation in the range 0.5–1.5 depending on the specific simulation and cross-section parameters (Essig et al. 2019). We will show in Section 6 that in the lmfp regime, it is valid to absorb constant factors resulting from integration over the MB velocity distributions into β when σ is constant. However, we will demonstrate that a more careful treatment is required for cases of relative velocity-dependent cross section σ(v₁₂).

A nice feature of this gravothermal fluid formalism is that, given the halo initial density profile as NFW with characteristic radius r_s and density ρ_s, the set of Equations (3)–(6) can be rewritten into the scale-free form such that the unitless gravothermal solutions can be applied to an arbitrary halo by rescaling. Specifically, all physical quantities in the equation set x in the fluid formalism can be separated into a unit factor x₀ and a scale-free factor $\hat{x}$ (i.e., $x={x}_{0}\hat{x}$) such that Equations (3)–(6) can be written into the unitless form:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \hat{M}}{\partial \hat{r}} & = & \hat{\rho }{\hat{r}}^{2},\\ \displaystyle \frac{\partial (\hat{\rho }{\hat{v}}_{\mathrm{rms}}^{2})}{\partial \hat{r}} & = & -\displaystyle \frac{\hat{M}\hat{\rho }}{{\hat{r}}^{2}},\\ \displaystyle \frac{\hat{L}}{{\hat{r}}^{2}} & = & -\hat{\kappa }\displaystyle \frac{\partial {\hat{v}}_{\mathrm{rms}}^{2}}{\partial \hat{r}},\\ \hat{\rho }{\hat{v}}_{\mathrm{rms}}^{2}{\left(\displaystyle \frac{\partial }{\partial \hat{t}}\right)}_{\hat{M}}\mathrm{ln}\displaystyle \frac{{\hat{v}}_{\mathrm{rms}}^{3}}{\hat{\rho }} & = & -\displaystyle \frac{1}{{\hat{r}}^{2}}\displaystyle \frac{\partial \hat{L}}{\partial \hat{r}}.\end{array}\end{eqnarray} \tag{ 9 }$$

The unit variables are (Essig et al. 2019)

$$\begin{eqnarray}\begin{array}{rcl}{r}_{0} & = & {r}_{s},\,{\rho }_{0}={\rho }_{s},\,{M}_{0}=4\pi {r}_{s}^{3}{\rho }^{s},\\ {v}_{{\rm{rms}},0} & = & \sqrt{4\pi G{\rho }_{s}}{r}_{s},\,{t}_{0}=1/\sqrt{4\pi G{\rho }_{s}},\\ {L}_{0} & = & {\left(4\pi \right)}^{5/2}{G}^{3/2}{\rho }_{s}^{5/2}{r}_{s}^{5},\,{(\sigma /m)}_{0}=1/({\rho }_{s}{r}_{s}),\end{array}\end{eqnarray} \tag{ 10 }$$

and it can be derived that ${\hat{\kappa }}_{\mathrm{lmfp}}=0.27\times 4\pi \beta \hat{\rho }{\hat{v}}_{\mathrm{rms}}^{3}(\hat{\sigma /m})$ and ${\hat{\kappa }}_{\mathrm{smfp}}=2.1{\hat{v}}_{\mathrm{rms}}/(\hat{\sigma /m})$. Hereafter we will abbreviate $(\hat{\sigma /m})$ as $\hat{\sigma }$ for simplicity. We assume $\hat{\kappa }=1/({\hat{\kappa }}_{\mathrm{lmfp}}^{-1}+{\hat{\kappa }}_{\mathrm{smfp}}^{-1})$ to join the lmfp and smfp regimes smoothly.

The steps to solve Equation (9) numerically are (Essig et al.2019)

• 1.  
  Determine input parameters β and $\hat{\sigma }$.
• 2.  
  Grid the halo with NFW density profile $\hat{\rho }=(1/\hat{r})/{\left(1+\hat{r}\right)}^{2}$ into 150 log radial bins spread uniformly in the range $-2\leqslant \mathrm{log}\hat{r}\leqslant 3$. ^{7 , 8}  We assume boundary conditions $\hat{M}=\hat{L}=0$ at the inner radius and $\hat{L}=0$ at the outer radius to solve for the radial profile of ${\hat{v}}_{\mathrm{rms}}$ through the Jeans equation. Following Pollack et al. (2015) and Essig et al. (2019), we take values of ${\hat{M}}_{i}$ and ${\hat{L}}_{i}$ of the ith shell at the log radial grid ${\hat{r}}_{i}$, while ${\hat{\rho }}_{i}$, ${\hat{v}}_{\mathrm{rms},i}$ are taken as the averaged values between the (i − 1)th and the ith shells.
• 3.  
  Let the fluid code take a small time step $d\hat{t}$ such that the specific energy $\hat{u}=3{\hat{v}}_{\mathrm{rms}}^{2}/2$ changes by a factor of no more than 10⁻³ among all the shells:

  $$\begin{eqnarray}&&\displaystyle \frac{d{\hat{u}}_{i}}{d\hat{t}}=-\displaystyle \frac{{\rm{\partial }}{\hat{L}}_{i}}{{\rm{\partial }}{\hat{M}}_{i}}.\end{eqnarray} \tag{ 11 }$$
• 4.  
  After taking one time step forward, perturb $\hat{r}$, $\hat{\rho }$, and ${\hat{v}}_{\mathrm{rms}}$ of each shell 10 times such that the gravothermal equilibrium (Equation (4)) is established again throughout the halo. During the perturbation, the mass, specific energy, and entropy $s=\mathrm{ln}({\hat{v}}_{\mathrm{rms}}^{3}/\hat{\rho })$ of each shell are conserved.
• 5.  
  Repeat steps 2, 3, and 4 iteratively until core collapse.

We checked that the basic assumptions made in solving Equation (9), i.e., (i) the halo is spherically symmetric, and (ii) the halo mass and scale do not vary with time, are satisfied in the idealized N-body simulations for isolated halos. However, those assumptions may be not valid for cosmological simulations due to halo interactions, mergers, and mass accretion.

Notice that most of the current gravothermal fluid studies are interested in SIDM models with moderate cross section and β values near unity. Such a setup results in the halo remaining in the lmfp regime for most of its life until core collapse occurs. Moreover, σ decreases drastically with v_rms at v_rms ≳ ω under the scenario of a velocity-dependent cross section, leaving the heat conduction more likely to be dominated by the lmfp term. We therefore focus on the gravothermal solutions in the lmfp regime in this work. Specifically, we choose input parameters $\hat{\sigma }=0.01,\beta =0.6$ in step 1 and solve Equation (9) numerically to derive the gravothermal solutions, which will be used by the mapping method introduced in Sections 4–6. We stop the evolution of the fluid code when $\beta \hat{\sigma }\hat{t}$ reaches 173. At the last time step, the halo is still in the lmfp regime and has entered the core-collapsing stage, where its central density increases rapidly. The time evolution of the ratio between ${\hat{\kappa }}_{\mathrm{lmfp}}$ and ${\hat{\kappa }}_{\mathrm{smfp}}$ at the halo center is shown in Figure 1 by the blue solid curve.

Zoom In Zoom Out Reset image size

We notice that the halos simulated in Section 2 adopt a much greater cross section of $\hat{\sigma }=0.65$. To check the validity of the lmfp assumption we numerically solve Equation (9) for another set of input parameters, $\{\beta =0.82,\hat{\sigma }=0.65\}$, and present the time evolution of the halo central density as the red dashed curve of Figure 1. This halo still evolves in the lmfp regime during most of its lifetime before core collapse, although the increase in cross section boosts the halo to be closer to the smfp regime since ${\hat{\kappa }}_{\mathrm{lmfp}}/{\hat{\kappa }}_{\mathrm{smfp}}\propto \beta {\hat{\sigma }}^{2}$. Ideally we would like to simulate a halo with smaller cross section to better match the assumption made in this work. However, this will delay the onset of core collapse and lead to more expensive and unstable idealized N-body simulations. We will show later that even for this large constant cross section $\hat{\sigma }=0.65$ the lmfp assumption still holds, and the simulated halo evolution matches the lmfp gravothermal solution well. The smfp tension will be further alleviated when we consider the velocity-dependent cross-section scenarios, where the velocity dependence will largely decrease the effective cross section and bring the halo further away from the smfp regime.

The gravothermal evolution of SIDM halos described by Equation (9) is shown in Figure 2. Specifically, at the NFW initial state (shown by the red solid curves) where $\hat{\rho }=(1/\hat{r})/{\left(1+\hat{r}\right)}^{2}$, the halo temperature is highest at $\hat{r}\approx 0.8$. The heat therefore flows to both the halo center and outskirts. The halo center keeps absorbing heat until local thermal equilibrium is reached. At this point a cored region is formed where the density and the rms velocity are uniform, as shown by the green dashed curves. DM particles in the halo cored region continuously lose energy due to outward heat transfer. In this process their orbital scale decreases and the kinetic temperature increases. As shown by the blue dashed–dotted curves, the halo core density and temperature will keep increasing and will move further away from thermal equilibrium with the halo outskirts. Such a regime is also referred to as gravothermal catastrophe or core collapse.

Zoom In Zoom Out Reset image size

While the gravothermal model is computationally much more efficient than running an N-body simulation it is still too slow to apply to large ensembles of halos as might be required in a semianalytic model of subhalo evolution, or when generating realizations of subhalo populations for analysis of gravitational lensing systems (e.g., Benson 2012; Gilman et al. 2021). For example, the calculations shown in this section required around 24 CPU hours to compute. While this is trivially manageable for a single halo, it is impractical to apply to ensembles of millions of halos. As such it is useful to be able to map a single solution of the gravothermal model to more general cases. We will explore this possibility in the remainder of this paper.

Additionally, in the Appendix we provide a triple power-law empirical model calibrated to this set of fully lmfp solutions with high precision. The advantage of this analytical model is that it captures the halo core size, core collapsing time, and the radius where the cored density profile joins smoothly to the halo NFW outskirts with simple formulas, and can be easily implemented into semianalytic models of galaxy formation.

In this section we briefly review the mapping method discussed in Balberg et al. (2002), Koda & Shapiro (2011), and Essig et al. (2019) that transfers the gravothermal solutions under a constant cross section ${\hat{\sigma }}_{1}$ to another constant cross section ${\hat{\sigma }}_{2}$, assuming the halo is mostly evolving in the lmfp regime.

The third and fourth equations in Equation (9) show that the halo evolution rate is degenerated with radius-independent factors in the heat conductivity $\hat{\kappa }$. In the lmfp regime ${\hat{\kappa }}_{\mathrm{lmfp}}\propto \beta \hat{\sigma }$, therefore $\beta \hat{\sigma }$ and the inverse of the halo evolution time $\hat{t}$ are degenerate such that the time evolution of the gravothermal solution can be parameterized by $\beta \hat{\sigma }\hat{t}$. In other words, there exists a one-to-one mapping ${\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}{\hat{t}}_{{\rm{a}}}\leftrightarrow {\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}{\hat{t}}_{{\rm{b}}}$ between two sets of gravothermal solutions with parameters $\{{\beta }_{{\rm{a}}},{\hat{\sigma }}_{{\rm{a}}}\}$ and $\{{\beta }_{{\rm{b}}},{\hat{\sigma }}_{{\rm{b}}}\}$ such that

$$\begin{eqnarray}&&{\hat{x}}_{{\rm{a}}}(\hat{r},{\hat{t}}_{{\rm{a}}})={\hat{x}}_{{\rm{b}}}\left(\hat{r},{\hat{t}}_{{\rm{a}}}\times \displaystyle \frac{{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}}{{\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}}\right),\end{eqnarray} \tag{ 12 }$$

where $\hat{x}(\hat{r},\hat{t})$ stands for an arbitrary scale-free halo property at radius $\hat{r}$ and time $\hat{t}$, including $\hat{\rho }$, ${\hat{v}}_{\mathrm{rms}}$, $\hat{M}$, $\hat{L}$, etc. We show the scale-free halo central density $\hat{{\rho }_{c}}$ versus $\beta \hat{\sigma }\hat{t}$ for two sets of gravothermal solutions with very different input parameters in Figure 3. The time evolution of these two halos is almost identical after the first few time steps, where ${\hat{\kappa }}_{\mathrm{lmfp}}$ is comparable to ${\hat{\kappa }}_{\mathrm{smfp}}$ for the $\hat{\sigma }=0.65$ case, and the lmfp assumption is not very accurate (see Figure 1).

Zoom In Zoom Out Reset image size

Although the gravothermal fluid formalism considers a very simplified picture, it provides numerical solutions in surprisingly good agreement with idealized N-body simulations of isolated halos up to a slightly varying β factor. As an example, in Figure 4 we present comparisons of the evolution of the halo central density over time between the gravothermal solutions (red solid curve) and the averaged Arepo N-body simulations (black solid curve) for a Milky Way-sized SIDM halo with ρ_s = 4.2 × 10⁶ M_⊙ kpc^–3, r_s = 24.54 kpc, and σ/m = 30 cm² g^–1 (see Section 2 for simulation details). To test convergence, we reduce the particle number by a factor of 4 in each N-body realization and present the averaged simulation results as the cyan dashed curve. We find the rescaled lmfp gravothermal solution with β = 0.82 captures the overall halo evolution and reproduces the core collapsing time in the Arepo simulation results. We notice that the fluid numerical solution has a halo central density that drops faster than the N-body results prior to isothermal core formation (0 ≤ t/Gyr ≤ 2). The central density also grows slightly faster than N-body at t ≳ 18 Gyr. Those differences are, again, caused by the fact that we have selected a very large cross section for the N-body simulations. The basic assumption of this mapping method, that the halo is always in the lmfp regime, is not of high accuracy. Since our focus is to extend the lmfp gravothermal fluid formalism to the velocity-dependent cross-section scenario and capture the overall halo evolution in the N-body simulations, in this work we will accept this slight disagreement. Moreover, the smfp problem is alleviated when we simulate velocity-dependent cross-section models.

Zoom In Zoom Out Reset image size

In this section we introduce a mapping method to transfer the gravothermal solutions under a constant cross section $\hat{\sigma }$ to those under a velocity-dependent cross-section model σ(v). In numerical simulations the velocity used to determine the scattering probability between two particles is generally the particle relative velocity v₁₂, but in this section we will consider a simpler case where the cross section is determined by the particle rms velocity. We will discuss this mapping method under the more practical σ(v₁₂) scenario in Section 6.

In this paper we will use the cross-section model introduced by Equation (2) as an example, but this mapping method is applicable to an arbitrary σ(v) model.

In Section 4 we have shown that there exists a one-to-one mapping ${\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}{\hat{t}}_{{\rm{a}}}\leftrightarrow {\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}{\hat{t}}_{{\rm{b}}}$ between two sets of gravothermal solutions with different β and constant cross sections. Now consider two sets of gravothermal solutions, one with constant cross section $\{{\beta }_{{\rm{a}}},{\hat{\sigma }}_{{\rm{a}}}\}$ and the other with a velocity-dependent cross section $\{{\beta }_{{\rm{b}}},{\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{\mathrm{rms},\ {\rm{b}}})\}$. The linear time mapping does not hold anymore since ${\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{\mathrm{rms},\ {\rm{b}}})$ is a function of radius and time rather than a constant. However, if in every small time interval we may assume a characteristic halo-wide cross section (or, equivalently, a characteristic velocity) such that we can drop the radial dependence of the cross section $\hat{\sigma }({\hat{v}}_{\mathrm{rms},\ {\rm{b}}}(\hat{r},{\hat{t}}_{{\rm{b}}}))\approx \hat{\sigma }({\hat{v}}_{\mathrm{rms},\ {\rm{b}}}^{c}({\hat{t}}_{{\rm{b}}}))$, then there exists the following one-to-one mapping for every small time step:

$$\begin{eqnarray}\begin{array}{ccc}{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}d{\hat{t}}_{a} & = & {\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{{\rm{r}}{\rm{m}}{\rm{s}},\,{\rm{b}}}^{c}({\hat{t}}_{{\rm{b}}}))d{\hat{t}}_{{\rm{b}}}\\ & = & {\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{{\rm{r}}{\rm{m}}{\rm{s}},\,{\rm{a}}}^{c}({\hat{t}}_{{\rm{a}}}))d{\hat{t}}_{{\rm{b}}}.\end{array}\end{eqnarray} \tag{ 13 }$$

Here the second equation holds because alternating a time-dependent yet radius-independent SIDM cross section has no influence on the gravothermal fluid solution besides nonlinearly stretching the time axis. In other words, there still exists the following one-to-one mapping between these two sets of gravothermal solutions:

$$\begin{eqnarray}\begin{array}{rcl}{\hat{x}}_{{\rm{a}}}(\hat{r},{\hat{t}}_{{\rm{a}}}) & = & {\hat{x}}_{{\rm{b}}}(\hat{r},{\hat{t}}_{{\rm{b}}})\\ & = & {\hat{x}}_{{\rm{b}}}\left(\hat{r},{\displaystyle \int }_{0}^{{\hat{t}}_{{\rm{a}}}}\displaystyle \frac{{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}}{\left.\left.{\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{{\rm{r}}{\rm{m}}{\rm{s}},\,{\rm{a}}}^{c}(\hat{t}{{\rm{{\prime} }}}_{{\rm{a}}}\right)\right)}d\hat{t}{{\rm{{\prime} }}}_{{\rm{a}}}\right).\end{array}\end{eqnarray} \tag{ 14 }$$

In practice one has the constant-cross-section gravothermal solution that specifies the radial profiles of halo density and rms velocity at an array of time ${\hat{t}}_{{\rm{a}}}=\{{\hat{t}}_{{\rm{a}}}^{0},{\hat{t}}_{{\rm{a}}}^{1},\ldots ,{\hat{t}}_{{\rm{a}}}^{N}\}$ (here ${\hat{t}}_{{\rm{a}}}^{i}$ is the time of the ith step), and aims to map it to some velocity-dependent cross-section scenario. To achieve this one will need to first compute the time step $\delta {\hat{t}}_{{\rm{a}}}=\{{\hat{t}}_{{\rm{a}}}^{0},{\hat{t}}_{{\rm{a}}}^{1}-{\hat{t}}_{{\rm{a}}}^{0},\ldots ,{\hat{t}}_{{\rm{a}}}^{N}-{\hat{t}}_{{\rm{a}}}^{N-1}\}$, and then use Equation (13) to compute the mapped time step $\delta {\hat{t}}_{{\rm{b}}}=[{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}/{\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}({\hat{v}}_{\mathrm{rms},\ {\rm{a}}}^{c}({\hat{t}}_{{\rm{a}}}))]\delta {\hat{t}}_{{\rm{a}}}$. Summing $\delta {\hat{t}}_{{\rm{b}}}$ up, a new time axis ${\hat{t}}_{{\rm{b}}}$ is derived such that ${\hat{x}}_{{\rm{b}}}(\hat{r},{\hat{t}}_{{\rm{b}}})={\hat{x}}_{{\rm{a}}}(\hat{r},{\hat{t}}_{{\rm{a}}})$.

The remaining problem is to determine how to choose the characteristic rms velocity ${\hat{v}}_{\mathrm{rms}}^{c}$ at each time. Since the SIDM core evolves the most with time, while the halo outskirts show very little time evolution (see Figure 2), we find that simply using the halo central rms velocity shows very good mapping performance:

$$\begin{eqnarray}&&{\hat{v}}_{\mathrm{rms}}^{c}(\hat{t})={\hat{v}}_{\mathrm{rms}}(\hat{r}={\hat{r}}_{\min },\hat{t}).\end{eqnarray} \tag{ 15 }$$

Here ${\hat{r}}_{\min }={10}^{-2}$ is the innermost radial bin of the fluid formalism.

To test this mapping method, we assume ${\hat{\sigma }}_{{\rm{c}}}=0.01$, β = 0.6, and numerically solve Equation (9) for two cross-section models. For the first test we assume a constant cross section $\hat{\sigma }={\hat{\sigma }}_{{\rm{c}}}$, while for the second test we assume $\hat{\sigma }\,={\hat{\sigma }}_{c}/{\left(1+{\hat{v}}_{{\rm{rms}}}^{2}/{\hat{\omega }}^{2}\right)}^{2}$. We consider four velocity-dependent cross-section models with $\hat{\omega }=0.01$, 0.1, 0.3, and 1.0, corresponding to very strong, strong, intermediate, and weak velocity dependence. Figure 5 presents the unitless time of the constant-cross-section gravothermal solution rescaled by this mapping method. Notice that the 1D velocity dispersion of dark matter particles ${\hat{v}}_{\mathrm{rms}}$ varies from 0.1 to 0.4 near the halo center (see Figure 2 bottom panel), which is small compared to 1.0. The halo SIDM cross section for the $\hat{\omega }=1.0$ case is therefore almost a constant $\hat{\sigma }={\hat{\sigma }}_{c}$, corresponding to a nearly trivial time mapping relation shown by the magenta dotted curve. Decreasing $\hat{\omega }$ suppresses the total cross-section value $\hat{\sigma }$ and also boosts the velocity dependence of the particle scattering probabilities. As a result, the rescaled $\hat{t}$ enlarges for the lower $\hat{\omega }$ cases, and the time mapping nonlinearity becomes more significant.

Zoom In Zoom Out Reset image size

Comparisons between the mapped constant-cross-section gravothermal solutions (green dashed lines) and the actual velocity-dependent gravothermal solutions (black solid curves) are shown in Figure 6. We find this simple remapping method can successfully reproduce the gravothermal solutions in velocity-dependent cross-section models in all the cases, although it overestimates the core collapsing time by 14%, 13%, 7.5%, and 1.3%, respectively. This overestimation is caused by Equation (15), where we set the halo central rms velocity as a characteristic velocity that determines the overall DM particle cross section at every moment. A true characteristic v_rms is generally lower than the central velocity, and therefore our method slightly underestimates the effective cross section at each time and will further overestimate the core collapsing time.

Zoom In Zoom Out Reset image size

We also check the performance of the mapping method introduced in O22. Instead of using a time-varying characteristic rms velocity ${\hat{v}}_{\mathrm{rms}}^{{\rm{c}}}$, O22 makes a more simplified assumption to use the rms velocity at the halo center and at the instance of maximum core size as the characteristic value. Effectively the characteristic rms velocity is assumed to be

$$\begin{eqnarray}&&{\hat{v}}_{\mathrm{rms}}^{{\rm{c}}}={\hat{v}}_{\mathrm{rms}}(\hat{r}={\hat{r}}_{\min },\hat{t}={\hat{t}}_{{\hat{\rho }}_{c,\min }}).\end{eqnarray} \tag{ 16 }$$

We show the time evolution of the halo central density mapped by O22 as magenta dotted curves in Figure 6. We show that the O22 method only underestimates the core collapsing time by 6%, 5%, 2.3%, and 0.3% for the $\hat{\omega }=0.01$, 0.1, 0.3, and 1.0 scenarios, but it cannot correctly capture the halo evolution prior to the instance of maximum core size.

The accuracy of our mapping method can be further boosted if we compute the characteristic ${\hat{v}}_{\mathrm{rms}}^{{\rm{c}}}$, or, equivalently, the characteristic cross section, at each time step more carefully. Specifically, at every time step we can define the characteristic cross section as the averaged cross section weighted by the mass of multiple shells:

$$\begin{eqnarray}&&{\hat{\sigma }}^{{\rm{c}}}=\displaystyle \frac{{{\rm{\Sigma }}}_{0}^{{i}_{\mathrm{out}}}4\pi {\hat{r}}_{i}^{2}\delta {\hat{r}}_{i}{\hat{\rho }}_{i}{\hat{\sigma }}_{b,c}({\hat{v}}_{\mathrm{rms},{\rm{i}}})}{{{\rm{\Sigma }}}_{0}^{{i}_{\mathrm{out}}}4\pi {\hat{r}}_{i}^{2}\delta {\hat{r}}_{i}{\hat{\rho }}_{i}},\end{eqnarray} \tag{ 17 }$$

where i is the shell index, and i_out corresponds to the shell that joins smoothly to the halo NFW outskirts, beyond which the DM particle self-interaction is negligible. The radius of shell i_out is given by Equation (A4) and Equation (A8) in the Appendix. After this step we can map the time step of the constant-cross-section gravothermal solution as $\delta {\hat{t}}_{{\rm{b}}}=[{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{a}}}/{\beta }_{{\rm{b}}}{\hat{\sigma }}_{{\rm{b}}}^{c}({\hat{t}}_{{\rm{a}}})]\delta {\hat{t}}_{{\rm{a}}}$. We find that this more careful way of computing the characteristic cross section effectively decreases the fractional error of the mapping method in predicting the core collapsing moment to 6.4%, 5.6%, 2.9%, and 0.5% for $\hat{\omega }=0.01$, 0.1, 0.3, and 1.0 cases, providing performance equally as good as O22 while maintaining a good match to the actual gravothermal numerical solutions before the instance of maximum core size. However, we do not use this more complicated treatment as the default method because it can be computationally expensive when applied to a large number of halos. Furthermore, in practice the mapping error in the core collapsing time of ∼10% at most is small compared to the uncertainties in the β factor.

To calibrate the gravothermal solution remapping method introduced in Section 5 to N-body simulations, we consider cross-section models that are determined by the relative velocity v₁₂ between two scattering particles. Such velocity-dependent models will influence the heat conductivity averaged over the MB velocity distribution 〈κ_lmfp〉, resulting in additional complexity.

It is worth recalling how 〈κ_lmfp〉 is computed for the case of constant cross section. By convention,

$$\begin{eqnarray}&&{\kappa }_{\mathrm{lmfp}}({v}_{12})=\displaystyle \frac{3}{2}\beta {{nH}}^{2}\displaystyle \frac{{k}_{{\rm{B}}}}{{t}_{r}}=\displaystyle \frac{3\beta {k}_{{\rm{B}}}}{8\pi {\rm{G}}}\displaystyle \frac{n\sigma }{m}{v}_{\mathrm{rms}}^{2}{v}_{12},\end{eqnarray} \tag{ 18 }$$

and the relative velocity can be defined as ${v}_{12}\,=\sqrt{{v}_{1}^{2}+{v}_{2}^{2}-2{v}_{1}{v}_{2}\cos {\theta }_{12}}$, where θ₁₂ is the angle between the velocity vectors v ₁ of particle 1 and v ₂ of particle 2. To compute the averaged heat conductivity over the MB velocity distribution, previous works compute the averaged v₁₂. Specifically, let us define the MB distribution and the averaged v₁₂ as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{\rm{M}}{\rm{B}}}(v){d}^{3}v & = & {f}_{{\rm{M}}{\rm{B}}}(v){v}^{2}dv\sin (\theta )d\theta d\phi ={v}^{2}{\left(\displaystyle \frac{m}{2\pi {k}_{{\rm{B}}}T}\right)}^{3/2}\\ & & \times \,\exp \left(-\displaystyle \frac{{{mv}}^{2}}{2{k}_{{\rm{B}}}T}\right)\sin (\theta )dvd\theta d\phi ,\\ \langle {v}_{12}\rangle & = & {\displaystyle \int }_{0}^{\pi }\sin {\theta }_{1}d{\theta }_{1}{\displaystyle \int }_{0}^{2\pi }d{\phi }_{1}{\displaystyle \int }_{0}^{2\pi }d{\phi }_{12}\\ & & \times \,{\displaystyle \int }_{0}^{\infty }{v}_{1}^{2}{f}_{{\rm{M}}{\rm{B}}}({v}_{1})d{v}_{1}\\ & & \times \,{\displaystyle \int }_{0}^{\infty }{v}_{2}^{2}{f}_{{\rm{M}}{\rm{B}}}({v}_{2})d{v}_{2}{\displaystyle \int }_{0}^{\pi }{v}_{12}\sin ({\theta }_{12})d{\theta }_{12}\\ & = & \displaystyle \frac{4}{\sqrt{\pi }}\sqrt{\displaystyle \frac{{k}_{{\rm{B}}}T}{m}}=\displaystyle \frac{4}{\sqrt{\pi }}{v}_{{\rm{r}}{\rm{m}}{\rm{s}}}=2.26{v}_{{\rm{r}}{\rm{m}}{\rm{s}}}.\end{array}\end{eqnarray} \tag{ 19 }$$

Using this result for 〈v₁₂〉 in Equation (18) gives the definition assuming a constant DM particle cross section:

$$\begin{eqnarray}&&\langle {\kappa }_{{\rm{lmfp}}}\rangle =\displaystyle \frac{3\beta {k}_{{\rm{B}}}}{8\pi G}\displaystyle \frac{n\sigma }{m}{v}_{{\rm{rms}}}^{2}\langle {v}_{12}\rangle =0.27\beta {{nv}}_{{\rm{rms}}}^{3}\sigma {k}_{{\rm{B}}}/(Gm).\end{eqnarray} \tag{ 20 }$$

One may argue that a better definition of the averaged heat conductivity might be 〈κ_lmfp〉 ∝ 〈v₁ v₂ v₁₂〉, or $\langle {v}_{12}^{3}\rangle$, or $\langle {v}^{n}\rangle /{v}_{\mathrm{rms}}^{n-3}$, where n is an arbitrary number. The specific choice does not matter for constant-cross-section scenarios because those different integrations will only cause a difference of a constant factor in Equation (20), which can be absorbed into the free parameter β. This is no longer true if a velocity-dependent cross section is assumed. Specifically, the definition $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}^{n}\sigma ({v}_{12})\rangle /({v}_{\mathrm{rms}}^{n-3}{\sigma }_{c})$ results in $\langle {v}^{n}\rangle /{v}_{\mathrm{rms}}^{n-3}$ for ω → ∞ . However, with finite ω the vⁿ term serves as a weight in the integration. Larger n biases the integration over the MB distribution toward higher velocities, corresponding to a smaller effective cross section and longer halo core collapse timescale. Therefore it is necessary to define 〈κ_lmfp〉 in a physically motivated way that matches N-body results.

Notice that the heat flux is determined by the product of the conductivity and the temperature gradient: H∂T/∂r ∼ ΔT ∼ ΔE is the characteristic temperature variation and energy transferred within a particle's orbit. Here temperature is also a statistical property averaged over the MB distribution $T\propto \langle {v}_{12}^{2}\rangle$, so the energy transfer rate is proportional to ${v}_{12}^{3}\sigma ({v}_{12})$ (Colquhoun et al. 2021). We therefore define the averaged heat conductivity as

$$\begin{eqnarray}&&\langle {\kappa }_{{\rm{lmfp}}}\rangle =0.27\beta {{nv}}_{{\rm{rms}}}^{3}{k}_{{\rm{B}}}/(Gm)\displaystyle \frac{\langle {v}_{1}{v}_{2}{v}_{12}^{3}\sigma ({v}_{12})\rangle }{\langle {v}_{1}{v}_{2}{v}_{12}^{3}\rangle }.\end{eqnarray} \tag{ 21 }$$

Here the v₁ v₂ term in the numerator comes from the scale height H². One v₁₂ term is contributed by 1/t_r, and the other ${v}_{12}^{2}$ term comes from the temperature gradient. The numerator is defined such that Equation (21) reduces to Equation (20) at ω → ∞ . This definition is slightly different from the ansatz assumed in O22, where $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{12}^{3}\sigma ({v}_{12})\rangle$. We notice that the triple integral $\langle {v}_{1}{v}_{2}{v}_{12}^{3}\sigma ({v}_{12})\rangle$ needs to be computed for every time interval, which significantly slows down the mapping calculations. However, it can be shown that the integration results are only sensitive to $\omega /\sqrt{T}\propto \omega /{v}_{\mathrm{rms}}$ times β σ_c , we therefore create a lookup table for the integration result spanning the parameter range 0.01 ≤ ω/v_rms ≤ 100 and use the linearly interpolated integration results to achieve fast mapping.

To map a set of constant-cross-section gravothermal solutions with parameters $\{{\beta }_{{\rm{a}}},{\hat{\sigma }}_{c,a}\}$ to those assuming a velocity-dependent cross-section model with parameters $\{{\beta }_{{\rm{b}}},{\hat{\sigma }}_{c,b},\hat{\omega }\}$, we perform the following steps:

• 1.  
  Compute the time step $\delta {\hat{t}}_{{\rm{a}}}=\{{\hat{t}}_{{\rm{a}}}^{0},{\hat{t}}_{{\rm{a}}}^{1}-{\hat{t}}_{{\rm{a}}}^{0},\ldots ,{\hat{t}}_{{\rm{a}}}^{N}-{\hat{t}}_{{\rm{a}}}^{N-1}\}$.
• 2.  
  Compute the mapped time step $\delta {\hat{t}}_{{\rm{b}}}$ =$[{\beta }_{{\rm{a}}}{\hat{\sigma }}_{{\rm{c}},{\rm{a}}}/{\beta }_{{\rm{b}}}]$ $[\langle {v}_{1}{v}_{2}{v}_{12}^{3}{\rangle }_{{\hat{t}}_{{\rm{a}}}}/\langle {v}_{1}{v}_{2}{v}_{12}^{3}{\hat{\sigma }}_{b}({v}_{12}){\rangle }_{{\hat{t}}_{{\rm{a}}}}]\delta {\hat{t}}_{{\rm{a}}}$. Here $\langle x{\rangle }_{{\hat{t}}_{{\rm{a}}}}$ implies computing the averaged x over an MB distribution where the temperature is given by the constant-cross-section gravothermal solutions at the innermost radius and time ${\hat{t}}_{{\rm{a}}}$. We choose the halo central temperature as a characteristic value at every time step with the same argument as discussed in Section 5.
• 3.  
  Sum over all $\delta {\hat{t}}_{{\rm{b}}}$ to obtain a new time axis such that ${\hat{x}}_{{\rm{b}}}(\hat{r},{\hat{t}}_{{\rm{b}}})={\hat{x}}_{{\rm{a}}}(\hat{r},{\hat{t}}_{{\rm{a}}})$.

Assuming that in the idealized N-body simulation β does not vary with ω, we compare the mapped gravothermal solutions and Arepo N-body simulations in velocity-dependent cross-section models introduced by Equation (2) with different ω. We find excellent matches, as presented in Figure 7. We notice that this comparison is only performed in a limited ω ≥ 400 km s⁻¹ range for this Milky Way-sized halo. Ideally, we would extend this comparison to cross-section models with stronger velocity dependence, i.e., with even lower characteristic velocity ω. However, a halo collapses more slowly as ω decreases, which makes the N-body simulation less stable and computationally more expensive. Furthermore, for the Milky Way-sized halo we simulate in this work ω ≤ 400 km s⁻¹ is a less interesting scenario, where the halo collapsing time is much longer than the age of the universe.

Zoom In Zoom Out Reset image size

We also show the gravothermal solution mapped by O22 by the magenta dotted curves in Figure 7. Specifically, O22 suggests to define the averaged heat conductivity as

$$\begin{eqnarray}&&\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto {K}_{n}={\left.\displaystyle \frac{\langle {v}_{12}^{n}\sigma ({v}_{12})\rangle }{\langle {v}_{12}^{n}\rangle }\right|}_{{r}_{\min },{t}_{{\rho }_{c,\min }}},n=3.\end{eqnarray} \tag{ 22 }$$

It is emphasized in O22 that the assumption 〈κ_lmfp〉 ∝ K₃ should be further tested by N-body simulations. We find this assumption overestimates the effective cross section as well as the heat conductivity, leading to a faster core collapse than the N-body results. However, we find assuming 〈κ_lmfp〉 ∝ K₅ and a linear time rescaling gives mapping performance as precise as Equation (21). This good match is also proved independently by Yang & Yu (2022) through N-body simulations.

The lmfp gravothermal solution mapping method introduced in Sections 4–6 can be easily generalized to the smfp regime as well as to the intermediate regime. However, in this work we do not test the validity of the smfp and intermediate-regime mapping methods introduced in this section due to difficulties in both computation and physics. Neither will we apply these mappings in the following sections. From the computational aspect, we find the gravothermal fluid formalism introduced in Section 3 runs into numerical issues when solving systems with very large cross section. Furthermore, the Arepo idealized SIDM simulation pipeline introduced in Section 2 is not tested for isolated halos with cross section greater than 60 cm² g^–1. From the aspect of physics, a very large SIDM cross section may change the formation of structure from that of CDM. The initial density profile of the SIDM halo with large cross section is also uncertain and unlikely to be NFW (Ahn & Shapiro2003).

Let us first revisit Equation (9) in the fully smfp regime. Now, since we switch focus from ${\hat{\kappa }}_{\mathrm{lmfp}}\propto \beta \hat{\sigma }$ to ${\hat{\kappa }}_{\mathrm{smfp}}\propto 1/\hat{\sigma }$, it becomes clear that the time mapping relation changes from ${\beta }_{a}{\hat{\sigma }}_{a}{\hat{t}}_{a}\leftrightarrow {\beta }_{b}{\hat{\sigma }}_{b}{\hat{t}}_{b}$ to ${\hat{t}}_{a}/{\hat{\sigma }}_{a}\leftrightarrow {\hat{t}}_{b}/{\hat{\sigma }}_{b}$ for the models with constant SIDM cross section. For the relative velocity-dependent cross-section scenarios one only needs to modify the time interval mapping as

$$\begin{eqnarray}&&\delta {\hat{t}}_{b}=\displaystyle \frac{\langle {v}_{12}^{3}{\rangle }_{{\hat{t}}_{a}}/{\hat{\sigma }}_{c,a}}{\langle {v}_{12}^{3}/{\hat{\sigma }}_{b}({v}_{12}){\rangle }_{{\hat{t}}_{a}}}\delta {\hat{t}}_{a},\end{eqnarray} \tag{ 23 }$$

to map gravothermal solutions computed under constant cross sections ${\hat{\sigma }}_{c,a}$ to those assuming velocity-dependent cross-section models ${\hat{\sigma }}_{b}({v}_{12})$. Here one v₁₂ term in the MB integration comes from the mean free path λ and another ${v}_{12}^{2}$ term comes from the energy transfer as has been discussed in Section 6. The two key differences between the lmfp and smfp mapping methods are: (1) The smfp heat conductivity and therefore the gravothermal solutions are no longer dependent on the free parameter β. (2) In the lmfp regime a halo collapses faster with a larger cross section, but the situation is reversed in the smfp regime.

To consider the intermediate regime where the mean free path and the scale height are comparable, the time mapping relation among constant-cross-section models has the form ${\hat{\kappa }}_{a}{\hat{t}}_{a}\leftrightarrow {\hat{\kappa }}_{b}{\hat{t}}_{b}$, where $\hat{\kappa }=1/({\hat{\kappa }}_{\mathrm{lmfp}}^{-1}+{\hat{\kappa }}_{\mathrm{smfp}}^{-1})$. For the relative velocity-dependent cross-section scenarios the time interval mapping relation is

$$\begin{eqnarray}&&\delta {\hat{t}}_{b}=\displaystyle \frac{\langle {\hat{\kappa }}_{a}{v}_{12}^{2}{\rangle }_{{\hat{t}}_{a}}}{\langle {\hat{\kappa }}_{b}{v}_{12}^{2}{\rangle }_{{\hat{t}}_{a}}}\delta {\hat{t}}_{a},\end{eqnarray} \tag{ 24 }$$

where the v₁₂ terms in the MB integration again come from the energy transfer. The challenge of computing Equation (24) compared to the asymptotic lmfp and smfp scenarios is that constant parameters β, ${\hat{\sigma }}_{c}$, as well as the scale-free density $\hat{\rho }$ can no longer be pulled out from the integration due to this more complicated expression for heat conductivity. It is therefore necessary to prepare a 4D lookup table that simultaneously covers the variations of $\{\beta ,{\hat{\sigma }}_{c},\hat{\omega }/{\hat{v}}_{\mathrm{rms}},\hat{\rho }\}$ to achieve fast mapping calculations.

Both the gravothermal fluid formalism and idealized N-body simulations are computationally expensive in the sense that they cannot be directly used to explore the continuous SIDM particle parameter space. The strength of the lmfp mapping method introduced in this work is that one needs to numerically solve for the gravothermal fluid equations only once, and can then extend this solution to halos with arbitrary cross-section models, as long as the halo is mostly evolving in the lmfp regime. It therefore serves as a powerful tool for deriving constraints on SIDM cross-section model parameters.

Halos hosting dwarfs and LSB galaxies are promising laboratories for constraining SIDM models because these small galaxies serve as powerful tracers of the DM halo density profile, while their shallow baryonic potentials have little influence on the DM dynamics (although see Pontzen & Governato 2012; Read et al. 2016, for example). As a demonstration of the utility of this work, we constrain the cross-section model parameter space favored by the measurement of the rotation curve of the LSB dwarf galaxy UGC 128 (van der Hulst et al. 1993; Kamada et al. 2017). We select UGC 128 because in this system the circular velocity is mostly determined by the dark matter component, while the stellar disk has negligible effect (Kamada et al. 2017). The fact that our model does not capture baryonic physics is therefore not expected to significantly influence the SIDM parameters inferred from UGC 128's rotation curve. We notice that the rotation curves of dozens of dwarf/LSB galaxies have been measured carefully (e.g., Kuzio de Naray et al. 2008; Oh et al. 2015; Essig et al. 2019). Some of the measurements have been used to derive stringent upper bounds on cross-section models (e.g., Read et al. 2018; Jiang et al. 2021). While a careful and thorough exploration of the SIDM cross-section model parameter space using all the available dwarf/LSB observational data will be an important step to take, in this work we choose to show just an example application of the gravothermal solution mapping method. We therefore illustrate regions of parameter space favored by the data on UGC 128 alone and defer a more complete analysis of a larger sample of dwarf/LSB galaxies to a subsequent work. In this work we trust the measurement of UGC 128's rotation curve and ignore potential underestimation of the error caused by beam smearing, coherent turbulence, the lack of determination of pressure support, and the treatment of inclination errors (Oman et al. 2019; Sellwood et al. 2021). The derived constraints can therefore be overtightened.

We assume β = 0.5 and a halo evolution time of 10 Gyr. We then use emcee (Foreman-Mackey et al. 2013) to conduct the Markov Chain Monte Carlo (MCMC) sampling over the 4D parameter space $-2.5\leqslant {\rm{log}}({\sigma }_{c}/{{\rm{cm}}}^{2}\,{{\rm{g}}}^{-1})\leqslant 4.0$, $0.0\leqslant {\rm{log}}(\omega /{\rm{km}}\,{{\rm{s}}}^{-1})\leqslant 4.0$, $5.0\leqslant {\rm{log}}({\rho }_{s}/{M}_{\odot }\,{{\rm{kpc}}}^{-3})\leqslant 8.0$, and $-1.0\leqslant {\rm{log}}({r}_{s}/{\rm{kpc}})\leqslant 2.0$, assuming uniform priors on these logarithmic variables. Here σ_c and ω are parameters introduced by the velocity-dependent cross-section model $\sigma ({v}_{12})={\sigma }_{c}/{\left(1+{v}_{12}^{2}/{\omega }^{2}\right)}^{2}$, and {ρ_s , r_s } are the NFW profile scaling parameters. In each MCMC step we use the gravothermal mapping method introduced in Section 6 to compute the halo density profile at t = 10 Gyr, and transfer the dark matter density profile to the rotation curve through ${v}_{\mathrm{circ}}(r)=\sqrt{{GM}(\lt r)/r}$. Here v_circ(r) is the circular velocity at radius r. We adopt a Gaussian likelihood function to direct 20 MCMC random walkers with the measurement of UGC 128's rotation curve. To get a quantitative sense of how strongly the baryonic gravitational potential in this system may influence the parameter fitting results, we perform two independent MCMC fits for the total rotation curve and the rotation curve of the dark matter halo, respectively. Contributions to the rotation curve from stars and gas are fit by Kamada et al. (2017). A triangle plot for the parameter constraints is presented in Figure 8(left panel).

Zoom In Zoom Out Reset image size

The MCMC fitting results show that the NFW parameters ρ_s and r_s are more sensitive to the differences between the total and DM-only rotation curves, while the $\mathrm{log}{\sigma }_{c}-\mathrm{log}{r}_{s}$ parameter constraints show little response to this data variation. This is fortunate as we are generally more interested in learning about the cross-section parameters rather than fitting NFW parameters of different halos with high accuracy. However, the generality of this feature should be further tested using a larger observational sample.

The MCMC results also show that the NFW scaling parameters ρ_s and r_s are degenerate. This is expected since the increase in amplitude of the density profile caused by increasing ρ_s can be compensated by decreasing r_s , which effectively shifts the halo density profile to smaller radii. Despite this parameter degeneracy we are still able to constrain {ρ_s , r_s } within a range much narrower than the flat priors. The best-fit NFW parameters for this halo fitted from the total rotation curve are ${\rm{log}}({\rho }_{s}/{M}_{\odot }\,{{\rm{kpc}}}^{-3})={6.49}_{-0.16}^{+0.17}$ and ${\rm{log}}({r}_{s}/{\rm{kpc}})=1.34\pm 0.09$ (1σ errors). For the rotation curve with the fitted stellar disk and gas components subtracted, the best-fit NFW scales are ${\rm{log}}({\rho }_{s}/{M}_{\odot }\,{{\rm{kpc}}}^{-3})=6.40\pm 0.17$ and ${\rm{log}}({r}_{s}/{\rm{kpc}})=1.35\pm 0.10$. There is a "L"-shaped degeneracy between the cross-section parameters $\mathrm{log}{\sigma }_{c}$ and $\mathrm{log}\omega$. One can also see such a degeneracy from the asymptotic behavior of the cross-section model. Specifically, the cross-section model reduces to the constant factor σ_c for high ω, and the parameter constraints become independent of ω. On the other hand, lowering ω effectively decreases the cross section, resulting in a weaker upper bound on σ_c. The turnover of this "L"-shaped band occurs on a scale comparable to the characteristic rms velocity of the halo. Therefore, the degeneracy between σ_c and ω can be at least partially broken by combining similar parameter constraints for systems of different sizes, masses, and characteristic velocity dispersions, and we defer a more careful analysis to a future work. We also emphasize that in this fit we have assumed β = 0.5 and a halo evolution time of 10 Gyr. Due to the $\beta \hat{\sigma }\hat{t}$ degeneracy introduced in Sections 4–6, assuming alternative values for β and halo evolution time would rescale the σ_c constraints by a factor of 1/(β t). The constraints on parameters {ω, ρ_s , r_s } are insensitive to the assumptions on β and halo age.

The single "L"-shaped band in the $\mathrm{log}{\sigma }_{c}-\mathrm{log}\omega$ parameter space favored by the observed rotation curve was a surprise to us. Specifically, we have assumed an NFW initial density profile, for which the halo central density first decreases, reaching a maximum core size, and then increases with time. There will be in general two β σ t instances where the cored halo density profiles are very similar to each other (unless the halo is measured to be precisely at the moment of maximum core size, or the halo has entered the late core collapsing stage where its central density is higher than in the initial state): one from the core formation process and another one from the core collapsing process. The existence of two gravothermal solutions with similar halo density profiles but very different $\beta \hat{\sigma }\hat{t}$ is one of the causes of the huge dispersion in SIDM cross-section constraints among different works (Kamada et al. 2020). For example, SIDM cross sections of the Milky Way's dwarf satellite galaxies constrained by isothermal models such as those of Kamada et al. (2017) and Valli & Yu (2018) only pick out the low-cross-section solution because the isothermal models by design only work for the SIDM halo core formation process. Cross sections constrained by isothermal-based methods are therefore relatively low, ranging from 0.1 to 40 cm² g^–1. On the other hand, previous studies such as Correa (2021) that use the gravothermal fluid formalism to constrain SIDM cross sections may only focus on the halo core collapsing process, and give higher constraints of 30–200 cm² g^–1 for similar systems. We find the "L"-shaped band favored by UGC 128's rotation curve points to one solution during the core formation process. As an example, we pick one point in the best-fit parameter space σ_c = 1 cm² g^–1, ω = 10⁴ km s⁻¹, and compare the mapped halo density profile as well as the rotation curve at t = 10 Gyr with observations in Figure 9. We show the best-fit gravothermal solutions to the total and DM-only rotation curves as black and blue solid curves, respectively. We find fixing the best-fit ρ_s , r_s , and ω values, but enlarging σ_c to about 100 cm² g^–1 gives very similar mapped halo density profile at t = 10 Gyr, although under such a large cross section the halo has started core collapse. The gravothermal solutions with identical halo central density to the best-fit solution, but a much larger cross section, are shown by dashed lines in Figure 9. The red/green dashed line corresponds to the case of total/DM-only rotation curve. We show in the right panel of Figure 9 that, for either the total or the DM-only case, rotation curve measurements of UGC 128 at radii less than about 5 kpc do not favor either one of those two gravothermal solutions, owning to their similar halo central densities. To clarify the origin of this tight cross-section constraint from UGC 128, in the left panel of Figure 9 we compare these two sets of gravothermal solutions to NFW profiles with the best-fit {ρ_s , r_s } parameters to the total/DM-only rotation curve, shown by the black/blue dotted curve. During the halo core formation process where σ_c = 1 cm² g^–1, the halo density profile is only suppressed in the cored region (r ≲ 5 kpc) compared to the NFW initial condition. However, for the σ_c ≈ 100 cm² g^–1 case the halo density is altered by DM self-interactions over a larger radial range (r ≲ 100 kpc). As a result the halo density profile outside the cored region becomes steeper than the initial NFW profile. This subtle difference is distinguishable by the rotation curve measurements outside the isothermal core. However, we emphasize that all the above analysis relies on the assumption that the data on UGC 128's rotation curve and its measurement uncertainties are statistically valid. The error bars on the rotation curve can be underestimated, especially in the central region of the galaxy, due to center offsets, noncircular motions, inclination, the effects of pressure support, and many other observational uncertainties (see de Blok 2010 for a review). While a careful reprocessing of the data on UGC 128's rotation curve is beyond the scope of this work, we attempt to account for potentially underestimated uncertainties by increasing the error bars on UGC 128's rotation curve of the innermost three radius bins, where the radius is less than 5 kpc, by a factor of 2. The MCMC fitting results are shown in Figure 8(right panel). Assuming larger measurement uncertainties near the galaxy center, the "L"-shaped bands in the $\mathrm{log}{\sigma }_{c}-\mathrm{log}\omega$ space become significantly wider such that cross sections as large as σ_c ∼ 100 cm² g^–1 can be preferred by the rotation curve measurement at 95% confidence level.

Zoom In Zoom Out Reset image size

We note that dissipative cooling processes (Essig et al. 2019), the effects of a baryonic disk (Jiang et al. 2022), and tidal effects (Zeng et al. 2022) can all speed up the onset of halo core collapse. On the other hand, Meshveliani et al. (2022) suggests that halo mass accretion and merger can delay the core collapse process, and proposes a correction term to the halo core collapsing time in order to account for this effect. However, the goal of this work is simply to introduce a convenient nonlinear method for rescaling time for the gravothermal fluid formalism without considering any of the above effects. Nevertheless, such a preliminary constraint already shows the advantage of our mapping method in the sense that it can rapidly explore the continuous parameter space and clearly present the parameter degeneracies.

In this work we introduce a fast mapping method that transfers the gravothermal fluid formalism solutions from the constant-cross-section scenario to arbitrary velocity-dependent cross-section models. Such a mapping is intrinsically built into the gravothermal fluid formalism. Specifically, when a halo is evolved fully in either the lmfp or the smfp regime, the cross section only influences the heat conductivity of the isolated halo at each radius and moment, which further determines the halo evolution rate. Varying the assumptions of the cross-section model therefore only remaps the time axis of the halo evolution, and has no effect on the radial profile amplitudes/shapes of halo density, velocity, enclosed mass, etc. Such a time remapping is linear between two different constant cross sections, but becomes nonlinear between constant and velocity-dependent cross-section models.

The thermal conductivity of an SIDM halo varies at different radii and times. To achieve fast mapping we use the thermal conductivity at the halo center at each instance as a characteristic value. We test the performance of this mapping method through comparing the mapped gravothermal solutions, assuming a constant cross section, with actual gravothermal numerical solutions assuming different velocity-dependent cross-section models. We find that the fast mapping method introduced in this work nicely captures the full time evolution of an SIDM halo, although it slightly overestimates the core collapsing time due to the simplified assumption about the instantaneous heat conductivity of the halo. We show that the accuracy of this mapping method can be further boosted through a more careful treatment of the calculation of instantaneous heat conductivity, but argue that since the fractional error of the default mapping is small compared to the uncertainties of halo evolution time and the free parameter β, the more complicated mapping method is unnecessary.

Calibrating gravothermal solutions with N-body simulations under the assumption of a velocity-dependent cross section is more challenging since in numerical simulations the scattering probabilities of particles are generally determined by their relative velocities, while the constant-cross-section gravothermal fluid formalism only cares about the particle rms velocity/velocity dispersion in each radial shell. Assuming the SIDM particle velocity in each halo shell follows the Maxwell–Boltzmann distribution, we define the averaged thermal conductivity as $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{1}{v}_{2}{v}_{12}^{3}\sigma ({v}_{12})\rangle$, where v₁, v₂ are the velocities of two particles, and v₁₂ is their relative velocity. The velocity terms in the averaging bracket of the MB distribution are contributed by particle orbiting scales, collisional relaxation times, and energy transfer. We show that this physically motivated definition of the averaged heat conductivity provides an excellent match between the mapped gravothermal solutions and Arepo idealized N-body simulation results, assuming multiple different velocity-dependent cross-section models.

The advantage of this mapping method lies in its high computational efficiency. Specifically, it takes about 24 CPU hours of computation with our gravothermal fluid code to reach the core-collapse regime for an SIDM halo in the lmfp regime, and introducing a velocity-dependent cross-section model brings at least two extra degrees of freedom. It is therefore impractical to continuously explore the cross-section model parameter space directly with the gravothermal fluid code. Luckily, by adopting the mapping method introduced in this work one only needs to solve the fluid code once, and one can fast map this set of numerical solutions to halos with arbitrary assumptions on cross section. Besides deriving parameter constraints, this mapping method is also suitable to be implemented into semianalytic galaxy formation models in order to achieve fast generation of SIDM halo populations.

As an example application, we select the LSB system UGC 128 where the baryonic disk potential has negligible influence on the galaxy rotation curve, and use the mapping method to constrain four free parameters in the NFW and velocity-dependent cross-section models. We have assumed β = 0.5 and a halo evolution time of 10 Gyr, but we show analytically that the β t assumption will only rescale constraints on the cross-section scaling parameter σ_c. The other three parameters {ω, ρ_s, r_s} are not sensitive to the assumed values for β and halo age. To get a quantitative sense of how much the baryonic gravitational potential in this system may influence the parameter fitting results, we perform two independent MCMC fits for the total rotation curve and the dark matter halo rotation curve, respectively, finding that the SIDM parameters of interest are unaffected. We also show that the accurately measured rotation curve of UGC 128 can distinguish two very similar gravothermal solutions during the core formation and core collapsing phases, and favors the low-cross-section solution. We notice that the constraints derived in this work do not account for other possible nonlinear effects that can influence the halo core collapsing time, including cooling, an additional gravitational potential contributed by a baryonic disk or black hole, tidal effects, halo mass accretion, and merger. The constraints can be too tight since we have trusted the measurements of UGC 128's rotation curve, and have ignored possible underestimation of the error. Nevertheless the advantage of this fast gravothermal solution mapping method is recognized in the sense that it can quickly explore the continuous 4D parameter space with the first-principles gravothermal solutions, and effectively capture degeneracies among the SIDM halo parameters of interest.

We also compare the performance of this mapping method with a similar strategy proposed in O22. We find the simpler method introduced in O22, where the characteristic heat conductivity is assumed to be a time- and radius-independent constant, outperforms the method introduced in this work in reproducing the core collapse time when mapping between gravothermal solutions with v_rms-dependent cross-section models. However, this work performs better in mapping the full halo evolution process, while O22 fails to capture the evolution of the halo central density before the instant of maximum core size if the cross-section model has strong velocity dependence. Correctly mapping the SIDM halo core formation process is as important as capturing the core collapse time because: (1) most current SIDM cosmological simulations evolve halos with mild self-interaction cross section for only around 10 Gyr, such that the simulation stops before an isothermal core has been completely formed (e.g., Elbert et al. 2015; Nadler et al. 2020), and (2) observations may prefer an SIDM model with small cross section. If this is the case, most halos hosting observable galaxy clusters/galaxies/dwarfs are still approaching their maximum core size. O22 propose to estimate the averaged thermal conductivity in a particle relative velocity-dependent scenario as $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{12}^{n}\sigma ({v}_{12})\rangle$, where n is to be calibrated to N-body simulations. We find that a physically motivated definition $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{1}{v}_{2}{v}_{12}^{3}\sigma ({v}_{12})\rangle$ outperforms the ansatz assumed in O22, where n = 3. However, we show that combining the linear time rescaling method introduced by O22 and the definition $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{12}^{5}\sigma ({v}_{12})\rangle$, also suggested by Yang & Yu (2022), gives an equally good match to the N-body simulation results as $\langle {\kappa }_{\mathrm{lmfp}}\rangle \propto \langle {v}_{1}{v}_{2}{v}_{12}^{3}\sigma ({v}_{12})\rangle$.

We thank Yiming Zhong for sharing his gravothermal fluid code and for many helpful discussions. We thank Mark Vogelsberger for sharing the Arepo N-body simulation package including the SIDM module. Computing resources used in this work were made available by a generous grant from the Ahmanson Foundation. This work was supported in part by the NASA Astrophysics Theory Program under grant No. 80NSSC18K1014.

In this appendix we introduce a convenient empirical model for the density profile of SIDM halos, calibrated to the lmfp gravothermal solution used in this work.

Balberg et al. (2002) has proved that the gravothermal density profile evolves with time self-similarly, and the density in the halo outskirts falls with radius as $\hat{\rho }\propto {\hat{r}}^{-2.19}$. This derivation is elegant, but not directly applicable to most of the numerical SIDM simulations. Specifically, in most idealized N-body as well as some cosmological simulations the halo starts from an NFW initial state where $\hat{\rho }\propto {\hat{r}}^{-3}$ at large radii. As a result, when a core has been formed the halo density still drops as $\hat{\rho }\propto {\hat{r}}^{-2.19}$ slightly outside the core radius. At even larger radii the density profile joins smoothly with the NFW initial condition, and drops as $\hat{\rho }\propto {\hat{r}}^{-3}$. At this stage the halo density profile can be captured by a triple power law.

We first define the instant of halo core collapse as the time at which the halo central density diverges:

$$\begin{eqnarray}&&\mathrm{log}(\beta \hat{\sigma }\hat{t})=E=2.238,\end{eqnarray} \tag{ A1 }$$

and the instant of maximum core size as the time when the halo central density reaches its minimum:

$$\begin{eqnarray}&&\mathrm{log}(\beta \hat{\sigma }\hat{t})=F=1.341.\end{eqnarray} \tag{ A2 }$$

We find at the early stage of the halo evolution, i.e., $\mathrm{log}(\beta \hat{\sigma }\hat{t})\lt F$, the halo density profile can still be described by the NFW profile multiplied by a tanh central cut:

$$\begin{eqnarray}&&\hat{\rho }=\displaystyle \frac{\tanh (\hat{r}/{r}_{\mathrm{core}})}{\hat{r}{\left(1+\hat{r}\right)}^{2}},\end{eqnarray} \tag{ A3 }$$

where we use free parameter r_core to capture the core size. We find the following empirical model described the time evolution of the halo core size at $\mathrm{log}(\beta \hat{\sigma }\hat{t})\leqslant 1.341$ with high accuracy:

$$\begin{eqnarray}&&\mathrm{log}{r}_{\mathrm{core}}={A}_{{r}_{\mathrm{core}}}{\left(\mathrm{log}(\beta \hat{\sigma }\hat{t})\right)}^{2}+{B}_{{r}_{\mathrm{core}}}\mathrm{log}(\beta \hat{\sigma }\hat{t})+{C}_{{r}_{\mathrm{core}}}.\end{eqnarray} \tag{ A4 }$$

Here ${A}_{{r}_{\mathrm{core}}}=-0.1078$, ${B}_{{r}_{\mathrm{core}}}=0.3737$, and ${C}_{{r}_{\mathrm{core}}}=-0.7720$ are calibrated to the gravothermal solutions.

As mentioned before, in later phases of halo evolution, $\mathrm{log}(\beta \hat{\sigma }\hat{t})\gt 1.341$, when a core has been formed, the density profile can be described by triple power law:

$$\begin{eqnarray}&&\hat{\rho }(\hat{r})=\displaystyle \frac{{\rho }_{{\rm{core}}}}{1+{\left(\hat{r}/{r}_{{\rm{core}}}\right)}^{s}{\left(1+\hat{r}/{r}_{{\rm{out}}}\right)}^{3-s}}.\end{eqnarray} \tag{ A5 }$$

Here s = 2.19 is provided by the gravothermal solution self-similar analysis of Balberg et al. (2002). Again we use free parameter r_core to trace the characteristic halo core size at each time step. Similarly, the parameter ${\rho }_{\mathrm{core}}$ captures the halo central density, and r_out traces the radius where the density profile influenced by self-interaction joins smoothly to the NFW outskirts. We develop the following empirical model for ${\rho }_{\mathrm{core}}$, r_core, and r_out, calibrated to the lmfp gravothermal solution:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{\rho }_{\mathrm{core}} & = & {A}_{{\rho }_{\mathrm{core}}}{\left(\mathrm{log}\beta \hat{\sigma }\hat{t}-(E+3)\right)}^{2}+{C}_{{\rho }_{\mathrm{core}}}\\ & & +\,\displaystyle \frac{{D}_{{\rho }_{\mathrm{core}}}}{{\left(E+0.0001-\mathrm{log}\beta \hat{\sigma }\hat{t}\right)}^{0.02}},\end{array}\end{eqnarray} \tag{ A6 }$$

where ${A}_{{\rho }_{\mathrm{core}}}=0.05771,{C}_{{\rho }_{\mathrm{core}}}=-21.64,{D}_{{\rho }_{\mathrm{core}}}=21.11$,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{r}_{\mathrm{core}} & = & {A}_{{r}_{\mathrm{core}}}{\left(\mathrm{log}\beta \hat{\sigma }\hat{t}-(E+2)\right)}^{2}+{C}_{{r}_{\mathrm{core}}}\\ & & +\displaystyle \frac{{D}_{{r}_{\mathrm{core}}}}{{\left(E+0.0001-\mathrm{log}\beta \hat{\sigma }\hat{t}\right)}^{0.005}},\end{array}\end{eqnarray} \tag{ A7 }$$

where ${A}_{{r}_{\mathrm{core}}}=-0.04049,{C}_{{r}_{\mathrm{core}}}=43.07,{D}_{{r}_{\mathrm{core}}}=-43.07$, and

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{r}_{\mathrm{out}} & = & {A}_{{r}_{\mathrm{out}}}{\left(\mathrm{log}\beta \hat{\sigma }\hat{t}-(E+2)\right)}^{2}+{C}_{{r}_{\mathrm{out}}}\\ & & +\displaystyle \frac{{D}_{{r}_{\mathrm{out}}}}{{\left(E+0.04-\mathrm{log}\beta \hat{\sigma }\hat{t}\right)}^{0.005}},\end{array}\end{eqnarray} \tag{ A8 }$$

where ${A}_{{r}_{\mathrm{out}}}=0.02403,{C}_{{r}_{\mathrm{out}}}=-4.724,{D}_{{r}_{\mathrm{out}}}=5.011$.

We compare the halo density profile given by the gravothermal solution and this empirical model at five characteristic moments in Figure 10. This simple empirical model accurately captures the overall time evolution and radial trend of the gravothermal solutions throughout the halo evolution process. We also show the halo central density using the exact result from the gravothermal fluid model and our empirical model (along with their fractional difference) in Figure 11. The fractional error of this empirical model is better than 10% at most times. Even at the very late core collapsing stage its accuracy is still better than 50%. This degree of accuracy is sufficient for many applications, considering the fact that the random fluctuations in N-body simulation results can cause the halo central density to fluctuate at a similar level. The fractional error shows noncontinuous behavior at the maximum core moment $\beta \hat{\sigma }\hat{t}=1.341$ because we switch the empirical description of the density profile from a tanh cutoff to the triple power law here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This analytic empirical model has wide applications. For example, it can be implemented into semianalytic galaxy formation models for fast SIDM halo simulation. As another example, we show in Section 5 that the characteristic radius r_out provided by this model can be used to boost the accuracy of the mapping method introduced in this work.