Accelerated core collapse in tidally stripped self-interacting dark matter halos

TL;DR

This paper investigates how self-interacting dark matter halos evolve under tidal forces, showing that core collapse can be accelerated, potentially leading to black hole formation in dwarf galaxies.

Contribution

It introduces a semi-analytic model calibrated with simulations to study core collapse in tidally stripped dark matter halos, highlighting the impact of tidal effects on collapse timescales.

Findings

Core collapse times are highly sensitive to outer density profiles.

Tidal stripping accelerates core collapse in satellite galaxies.

Potential formation of intermediate mass black holes from accelerated collapse.

Abstract

We use a semianalytic approach that is calibrated to N-body simulations to study the evolution of self-interacting dark matter cores in galaxies. We demarcate the regime where the temporal evolution of the core density follows a well-defined track set by the initial halo parameters and the cross section. Along this track, the central density reaches a minimum value set by the initial halo density. Further evolution leads to an outward heat transfer, inducing gravothermal core collapse such that the core shrinks as its density increases. We show that the time scale for the core collapse is highly sensitive to the outer radial density profile. Satellite galaxies with significant mass loss due to tidal stripping should have larger central densities and significantly faster core collapse compared to isolated halos. Such a scenario could explain the dense and compact cores of dwarf galaxies in the Local Group like Tucana (isolated from the Milky Way), the classical Milky Way satellite Draco, and some of the ultrafaint satellites. If the ultimate fate of core collapse is black hole formation, then the accelerated time scale provides a new mechanism for creating intermediate mass black holes.

Tables

Table 1. Table 1: The core mass and 3D velocity dispersion for various profiles with σ m = 5 ​ cm 2 / g subscript 𝜎 𝑚 5 superscript cm 2 g \sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g} at the time of maximum SMFP-core formation t ~ f subscript ~ 𝑡 𝑓 \tilde{t}_{f} , roughly when the condition κ SMFP ≤ κ LMFP subscript 𝜅 SMFP subscript 𝜅 LMFP \kappa_{\mathrm{SMFP}}\leq\kappa_{\mathrm{LMFP}} becomes true over the entire core region. Note that this definition of the core mass M SMFP subscript 𝑀 SMFP M_{\mathrm{SMFP}} differs from the core mass M core subscript 𝑀 core M_{\textrm{core}} used in Fig.  3 . We use M SMFP subscript 𝑀 SMFP M_{\mathrm{SMFP}} when extrapolating the empirical relation found in Ref.  Balberg  et al. ( 2002 )

Profile	${\tilde{t}}_f$	$t_f\left[\mathrm{Gyr}\right]$	$M_{\mathrm{SMFP}}\left({\tilde{t}}_f\right)\left[{10}^7{M}_{\odot }\right]$	$v_{3\text{D},c}\left({\tilde{t}}_f\right)\left[\mathrm{km}/s\right]$
NFW	374.56	95.66	5.7	83
TNFW $r_t=3r_s$	183.39	46.84	4.6	81
TNFW $r_t=r_s$ after $3\mathrm{Gyr}$	60.10	15.34	3.0	59
TNFW $r_t=r_s$	52.99	13.53	2.6	61

Equations

$$\frac{L}{4\pi r^2}=-\kappa \nabla T,$$

$$\begin{array}{c}\frac{\partial \tilde{L}}{\partial \tilde{r}}=-{\tilde{r}}^2\tilde{\rho }{\tilde{v}}^2{\left(\frac{\partial }{\partial \tilde{t}}\right)}_{\tilde{M}}\ln \left(\frac{{\tilde{v}}^3}{\tilde{\rho }}\right)\\ \tilde{L}=-\frac{3}{2}{\tilde{r}}^2\tilde{v}{\left(\frac{a}{b}{{\tilde{\sigma }}_m}^2+\frac{1}{C}\frac{1}{\tilde{\rho }{\tilde{v}}^2}\right)}^{-1}\frac{\partial {\tilde{v}}^2}{\partial \tilde{r}},\end{array}$$

$$\frac{\partial \tilde{M}}{\partial \tilde{r}}={\tilde{r}}^2\tilde{\rho }$$

$$\frac{\partial (\tilde{\rho }{\tilde{v}}^2)}{\partial \tilde{r}}=-\frac{\tilde{M}\tilde{\rho }}{{\tilde{r}}^2}$$

$$\frac{\partial \tilde{L}}{\partial \tilde{r}}=-{\tilde{r}}^2\tilde{\rho }{\tilde{v}}^2{\left(\frac{\partial }{\partial \tilde{t}}\right)}_{\tilde{M}}\ln \left(\frac{{\tilde{v}}^3}{\tilde{\rho }}\right)$$

$$\tilde{L}=-\frac{3}{2}{\tilde{r}}^2\tilde{v}{\left(\frac{a}{b}{{\tilde{\sigma }}_m}^2+\frac{1}{C}\frac{1}{\tilde{\rho }{\tilde{v}}^2}\right)}^{-1}\frac{\partial {\tilde{v}}^2}{\partial \tilde{r}},$$

$$\delta \tilde{u}\simeq -\left(\frac{\partial \tilde{L}}{\partial \tilde{M}}\right)\delta \tilde{t},$$

$$\frac{{\tilde{p}}_{i+1}^{\prime }-{\tilde{p}}_i^{\prime }}{({\tilde{r}}_{i+1}^{\prime }-{\tilde{r}}_{i-1}^{\prime })/2}=-\frac{{\tilde{M}}_i}{{\tilde{r}}_i^{\prime 2}}\frac{{\tilde{\rho }}_{i+1}^{\prime }+{\tilde{\rho }}_i^{\prime }}{2}$$

$$\tilde{L}=-\frac{3C}{2}{\tilde{r}}^2\tilde{\rho }{\tilde{v}}^3{(1+x^{\alpha })}^{-1/\alpha }\frac{\partial {\tilde{v}}^2}{\partial \tilde{r}},$$

Full text

I Introduction

Self-interacting dark matter (SIDM) Spergel and Steinhardt (2000); Firmani et al. (2000); Kaplinghat et al. (2016) is a compelling framework for explaining the small-scale structure formation puzzles Bullock and Boylan-Kolchin (2017). SIDM simulations show that elastic collisions of dark matter particles transfer heat towards the colder central regions of dark matter halos, lowering central densities and creating constant density cores Vogelsberger et al. (2012); Rocha et al. (2013); Peter et al. (2013); Zavala et al. (2013); Vogelsberger et al. (2014); Fry et al. (2015); Fitts et al. (2018), which has long been recognized as a means of alleviating the longstanding core-cusp problem de Blok (2010); Salucci et al. (2012).

The core-cusp issue is tricky, because galaxies exhibit a large diversity of rotation curves Kuzio de Naray et al. (2010); Oman et al. (2015), while following an approximate radial acceleration relation McGaugh et al. (2016). Recent work has shown that SIDM models with a cross section of few $\mathrm{cm}^{2}/\mathrm{g}$ at low velocities can solve these interconnected problems for spiral galaxies Kamada et al. (2017); Ren et al. (2018).

In this paper, we investigate the evolution of a dark matter halo in the presence of self-interactions using a semianalytic method, originally developed to study gravothermal collapse in globular clusters Lynden-Bell and Eggleton (1980); Quinlan (1996) and later applied to isolated SIDM halos Balberg et al. (2002); Koda and Shapiro (2011). This method allows us to track the full halo evolution at scales $\leq 100~{}\mathrm{pc}$, which are expensive to achieve with N-body simulations. We first characterize the full process of SIDM halo evolution, whereby large cores can be created today for cross sections of $\sim\!1~{}\mathrm{cm}^{2}/\mathrm{g}$ and the gravothermal collapse phase sets in for larger cross sections, as seen in simulations. We show that this temporal evolution is accelerated if the outer region of the halo is stripped, resulting in a higher density core today.

Our key result is that a large cross section leads to central densities that are higher than for field halos if the halo has experienced tidal stripping early in its formation history. Our results are applicable to local group galaxies that have interacted with the Milky Way (MW) in the past (such as Tucana) and satellites of the MW (including ultrafaint galaxies). The timescale for halo collapse depends on characteristic density of the initial halo and the details of the truncation process; thus, gravothermal collapse contributes to the diversity of MW satellites, for similar cross sections that can also explain the diversity of spiral galaxies Kamada et al. (2017); Ren et al. (2018). We also discuss the possible consequences of core collapse of the SIDM halo for intermediate-mass black hole (BH) formation.

II Gravothermal fluid model

We consider a spherical halo with a density profile $\rho(r,t)$ and an enclosed mass of $M(r,t)$ at radius $r$ and time $t$. We assume that this halo structure is set very early on ($\sim\!\mathrm{Gyr}$), as is relevant for the small mass halos that we focus on in this work. The halo is assumed to be isotropic and spherically symmetric, and dark matter is modeled as a thermally conducting fluid in quasistatic virial equilibrium. The dark matter particles with mass $m$ self-interact with a cross section per mass $\sigma_{m}$. They have a one-dimensional (1D) velocity dispersion $v(r,t)$ with a corresponding three-dimensional (3D) velocity dispersion $v_{3\textrm{D}}=\sqrt{3}v$. The relaxation time is defined as $t_{r}\equiv\lambda/(av)$, where $\lambda=1/(\rho\sigma_{m})$ is the mean-free path and $a=\sqrt{16/\pi}$ is a coefficient relevant for hard-sphere scattering of particles with a Maxwell-Boltzmann velocity distribution Balberg et al. (2002).

Through self-interactions, heat can flow from one region of the halo to another with a luminosity $L(r,t)$ through a spherical shell located at radius $r$. Heat transfer is given by Fourier’s law of thermal conduction,

[TABLE]

where $\kappa$ is the thermal conductivity. Strictly speaking, this equation is valid in the short mean-free path (SMFP) regime in which $\lambda\ll H$, where $H\equiv\sqrt{v^{2}/(4\pi G\rho)}$ is the gravitational scale height or Jean’s length. In this case, $\kappa_{\textrm{SMFP}}=(3/2)b\rho\lambda^{2}/(amt_{r})$ Balberg et al. (2002), where $b=25\sqrt{\pi}/32\approx 1.38$ is the effective impact parameter, calculated in Chapman-Enskog theory Chapman et al. (1990). In the long mean-free path (LMFP) regime in which $\lambda\gg H$, the thermal conduction formula empirically well-describes the gravothermal collapse of globular clusters with $\kappa_{\textrm{LMFP}}=(3/2)C\rho H^{2}/(mt_{r})$ Lynden-Bell and Eggleton (1980), where $C$ is a calibration parameter, described below. We interpolate between these regimes via $\kappa^{-1}=\kappa^{-1}_{\textrm{SMFP}}+\kappa^{-1}_{\textrm{LMFP}}$ Balberg et al. (2002).

Using the above ansatz, we can write the heat conduction equation in dimensionless variables as:

[TABLE]

where $\tilde{r}\equiv r/r_{s}$, $\tilde{v}\equiv v/v_{0}$, $\tilde{M}\equiv M/M_{0}$, $\tilde{t}\equiv t/t_{0}$, $\tilde{\sigma}_{m}\equiv\sigma_{m}/(4\pi r_{s}^{2}M_{0}^{-1})$, $\tilde{\rho}\equiv\rho/\rho_{s}$, and $\tilde{L}\equiv L/[(GM_{0}^{2}){(r_{s}t_{0})}^{-1}]$. We use the mass scale $M_{0}=4\pi r_{s}^{3}\rho_{s}$, velocity scale $v_{0}=\sqrt{GM_{0}/r_{s}}$, and timescale $t_{0}^{-1}=a\sigma_{m}v_{0}\rho_{s}$.

We numerically solve Eq. (2) along with the equations for mass conservation and hydrostatic equilibrium, using the techniques described in Ref. Pollack et al. (2015). See Appendix A for details of our numerical implementation.

We calibrate $C$ using the Pippin SIDM halos in the LMFP regime Elbert et al. (2015). We solve the gravothermal fluid equations using the Pippin CDM halo directly [as opposed to a Navarro-Frenk-White (NFW) fit to it] as our initial condition. As shown in Fig. 1, we find that $C=0.75$ produces halo density profiles that are in reasonable agreement with those from Pippin for $\sigma_{m}\leq 10~{}\mathrm{cm}^{2}/\mathrm{g}$, while the 3D velocity dispersions are systematically lower. The semianalytic approach enforces hydrostatic equilibrium, which does not hold in the simulation, particularly for small cross sections. Nonetheless, the velocity dispersions match qualitative features towards the center of the halo. In particular, for the largest cross sections shown, the halo is isothermal and the velocity dispersion has a negative radial gradient.

We note that the calibration depends on the simulation (e.g., isolated halo Koda and Shapiro (2011) versus cosmological Elbert et al. (2015)). As $C$ controls the importance of the LMFP regime, increasing (decreasing) its value causes a delay (advance) in the time of core collapse. For the range of $\sigma_{m}$ under consideration, however, mildly changing the calibration does not strongly affect the agreement in Fig. 1; see Appendix B. Additionally, deep in the LMFP regime, a change in $C$ can be compensated by a change in $\sigma_{m}$. Thus, while the precise values of $\sigma_{m}$ in our subsequent analyses are subject to systematic uncertainties of the calibration, our main results comparing halos with and without truncation (as we now discuss) still hold.

III Tidally stripped halos

We proceed to investigate the consequences of self-interactions on tidally stripped halos. Tidal effects impact both MW satellites (e.g., Tucana III Li et al. (2018); Simon et al. (2017) and Triangulum II Kirby et al. (2015)) and field dwarfs (e.g., Tucana Teyssier et al. (2012)). This process strips away stars and the outer regions of the subhalo and is thus potentially applicable to the dwarf galaxies on sufficiently radial orbits that experience significant mass loss.

For a profile $\rho_{\mathrm{NFW}}(r)=\rho_{s}/[(r/r_{s}){(1+(r/r_{s}))}^{2}]$, we model the effect of tidal stripping by changing the profile for $r>r_{t}$ to $\rho_{\mathrm{NFW}}(r_{t})\times(r_{t}/r)^{p_{t}}$, where $p_{t}=5$ Penarrubia et al. (2010). This simple model is meant to illustrate the effects of mass loss from the outer region, not to represent a realistic description of SIDM subhalo orbital history. We consider truncations at $r_{t}=r_{s}$ and $3\,r_{s}$, corresponding to a halo mass loss of $\sim\!90\%$ and $\sim\!70\%$, respectively, for a typical concentration of $20$. A halo truncated at the same $r_{t}$ with a shallower outer slope of $p_{t}=4$ undergoes core collapse slightly later, but the effect is qualitatively similar to $p_{t}=5$. We choose $\rho_{s}=0.019~{}\mathrm{M}_{\odot}/\mathrm{pc}^{3}$ and $r_{s}=2.59~{}\mathrm{kpc}$, obtained from fitting the Pippin CDM profile to an NFW profile using the Colossus cosmology code Diemer (2018).

IV Temporal evolution of central density

We focus on the central density $\rho_{c}$, defined as the average density of the innermost region of the halo specified in our analysis (at $\tilde{r}<0.01$, as described in Appendix A). We show the evolution of $\rho_{c}$ for initially NFW and initially TNFW profiles in the left panel of Fig. 2. For the cross sections of interest, the halos remain in the LMFP limit well into the runaway phase of collapse (the nearly vertical portion of the curves). In the LMFP regime, Eq. (2) is independent of the value of $\tilde{\sigma}_{m}$ and thus the gravothermal evolution shown holds for all $\sigma_{m}\leq 10~{}\mathrm{cm}^{2}/\mathrm{g}$.

For all three profiles, the central density drops as the core of the halo forms and rises again as the core begins to collapse. The minimum core density occurs approximately when the luminosity $L$ in the central region of the halo transitions from being negative (positive temperature gradient) to being positive (negative temperature gradient). The minimum density for halos with no tidal stripping is about $3\,\rho_{s}$, independent of the cross section. For the cases with tidal stripping, the collapse time becomes shorter and the minimum core density is higher. For the NFW profile, only for cross sections $\sigma_{m}\gtrsim 4~{}\mathrm{cm}^{2}/\mathrm{g}$ do we find the central density rise as the core begins to collapse within the lifetime of the Universe. However, for the TNFW profile with truncation at $r_{t}=r_{s}$ and $r_{t}=3\,r_{s}$, cross sections of $\sigma_{m}\gtrsim 0.4~{}\mathrm{cm}^{2}/\mathrm{g}$ and $2~{}\mathrm{cm}^{2}/\mathrm{g}$, respectively, have started to collapse by today.

To roughly gauge the impact of infall time, we allow the halo to evolve as before for a period of time before abruptly truncating it. We neglect the impact of multiple pericenter passages in our simplified analysis. We show the central density evolution of a halo truncated at $r_{t}=r_{s}$ after $3~{}\mathrm{Gyr}$ (or $z\simeq 2$) for an SIDM cross section of $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$ in the right panel of Fig. 2. Truncation times of $3-6~{}\mathrm{Gyr}$ are appropriate for the closest MW dwarfs Rocha et al. (2012). Such an extreme tidal stripping event leads to almost 2 orders of magnitude increase in the central density, showing the importance of this effect for nearby satellites. Note that an initial truncation at $r_{t}=3\,r_{s}$ barely alters the central density evolution away from the NFW case within the lifetime of the Universe, as seen from the left panel of Fig. 2 for $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$.

After truncation, hydrostatic equilibrium significantly lowers the pressure of the halo beyond the point of truncation, where most of the mass is lost, causing the velocity dispersion (and thus temperature) to decrease substantially. Heat flows towards the colder truncated part of the halo from the region near the scale radius, where the temperature is highest. As a result, heat is diffused more quickly within the truncated halo, leading to a faster formation of the isothermal core and thus an accelerated evolution for core collapse. We provide a more detailed description of the halo evolution in Appendix C.

V Observational consequences for the local group

Simulated CDM halos have circular velocities that are systematically higher than those observed for field dwarf galaxies Garrison-Kimmel et al. (2014) and MW satellites Boylan-Kolchin et al. (2012), while dark matter self-interactions reduce circular velocities to the observed range of $10-20~{}\mathrm{km}/\mathrm{s}$ Elbert et al. (2015). A notable exception is the field dwarf Tucana Fraternali et al. (2009); Kirby et al. (2014), which has been isolated for a long time ($\sim\!10~{}\mathrm{Gyr}$) but is observed to have a high circular velocity $v_{\mathrm{circ}}\simeq 33.7~{}\mathrm{km}/\mathrm{s}$ at $r_{1/2}\simeq 0.28~{}\mathrm{kpc}$ Kirby et al. (2014); Garrison-Kimmel et al. (2014), even for CDM Fitts et al. (2018); Garrison-Kimmel et al. (2014).

Accounting for truncation from tidal stripping, we have shown that the central density of a halo is larger than expected from CDM (if it has begun the process of core collapse), with a correspondingly larger circular velocity. Indeed, the TNFW halo with a delayed truncation shown in the right panel of Fig. 2 produces a circular velocity of $v_{\mathrm{circ}}=34~{}\mathrm{km}/\mathrm{s}$ at the same radius, in line with observations of Tucana. Our suggestion that Tucana may have experienced large tidal effects early on is consistent with a proposal to explain the isolation of Tucana Sales et al. (2007), wherein Tucana came in with a companion and was ejected due to three-body interactions, while the companion became bound to the MW (and perhaps fully disrupted). This early ejection scenario also seems consistent with the kinematics of Tucana Fraternali et al. (2009) and the early cutoff in star formation Saviane and Held (1996); Monelli et al. (2010).

Similar to Tucana, the MW satellite Draco (dwarf spheroidal) has a large central density Read et al. (2018), which may be difficult to explain with cross sections $\sigma_{m}>1~{}\mathrm{cm}^{2}/\mathrm{g}$, assuming evolution like a field halo Valli and Yu (2018). Draco most likely had a close pericenter passage according to Gaia data Fritz et al. (2018) and hence our arguments would suggest it could also have enhanced its central density compared to the field SIDM halo evolution. Note that the presence of the disk will increase tidal stripping and enhance the effects we have discussed D’Onghia et al. (2010); Garrison-Kimmel et al. (2017); Kelley et al. (2018).

The ultrafaint satellites of the MW are likely to have even smaller pericenters Simon (2018), owing to their large velocities and closer distances. These galaxies would certainly have been impacted by the shortened timescales and higher core densities predicted by our semianalytic method. A subhalo with $V_{\max}=10~{}\mathrm{km}/\mathrm{s}$ that is tidally truncated at $r=r_{s}$ at around $z=2$ would evolve to have a smaller $r_{\rm max}$ and a denser core by up to an order of magnitude (compared to a similar field SIDM halo) around $30-50~{}\mathrm{pc}$ radius (typical half-light radii for ultrafaints). The denser core may have a bearing on the survival of these subhalos Penarrubia et al. (2010), and the adiabatic contraction of the stars (as the SIDM density increases) would impact the compactness of the stellar distribution.

VI Intermediate-mass black holes

We note that the mass in the central region of the halo remains unaffected by a tidal stripping, which should only affect the halo at distances near or beyond the scale radius. Defining the core mass $M_{\textrm{core}}$ to be the mass within a core radius $r_{\textrm{core}}$ at which $\rho(r_{\textrm{core}})=0.1\,\rho_{c}$, we show the relation between $M_{\textrm{core}}$ and $v_{3\textrm{D},c}$ in Fig. 3 after the core has begun to collapse. This relation should hold for both truncated and field halos.

When the halo is in the LMFP regime, $\log(M_{\textrm{core}})$ decreases as $\log({v_{3\textrm{D},c}}^{2})$ increases, with a slope of $-4.17$, in agreement with the self-similar collapse scenario Balberg et al. (2002). Eventually, the SMFP heat transfer in the core begins to take over and the core becomes optically thick, inhibiting the diffusion of heat and mass from the core and thus resulting in a lower rate of mass loss.

An interesting consequence of the runaway gravothermal collapse in the SMFP regime is the possible formation of a BH Balberg and Shapiro (2002); Pollack et al. (2015); Choquette et al. (2018), which may exist today due to accelerated collapse from tidal stripping. If we define the core mass $M_{\mathrm{SMFP}}$ to be the mass within the radius at which $\kappa_{\textrm{SMFP}}=\kappa_{\textrm{LMFP}}$, then we can apply previous work that found $\mathrm{d}\log(M_{\mathrm{SMFP}})/\mathrm{d}\log({v_{3\textrm{D},c}}^{2})\simeq-0.85$ asymptotically Balberg and Shapiro (2002); Balberg et al. (2002). If a BH should form at the onset of relativistic instability at $v_{3\textrm{D},c}\approx c/3$ Balberg and Shapiro (2002); Balberg et al. (2002), then $M_{\mathrm{SMFP}}=10^{7}-10^{8}\,\mathrm{M}_{\odot}$ at $v_{3\textrm{D},c}=60~{}\mathrm{km}/\mathrm{s}$ (see Table 1 in Appendix C) would create BHs with mass $M_{\mathrm{BH}}=30-300\,\mathrm{M}_{\odot}$. This estimate involves a large extrapolation of our results and should be considered as motivation for further work rather than a concrete prediction.

Depending on the mass accretion history, these black holes may grow further. Finally, for rare halos that collapse early (with high concomitant densities), it is possible that this rapid gravothermal collapse may result in seeds for supermassive BHs, without assuming $\sigma_{m}\gg 1$ for a subdominant dark matter fraction Pollack et al. (2015); Choquette et al. (2018).

The exact properties of the SIDM cross section (e.g., velocity dependence, dissipative interactions, etc.) may become relevant as the extrapolation spans over orders of magnitude in $v_{c}$, and such interactions could remove particles from the collapsing core Chu and Garcia-Cely (2018); Vogelsberger et al. (2018); Essig et al. (2018); Choquette et al. (2018). Consideration of the angular momentum of SIDM particles may also require a closer examination of the final stage of collapse.

VII Conclusions

We applied the gravothermal fluid model, calibrated to cosmological SIDM simulations Elbert et al. (2015), to tidally stripped halos. We showed that for $\sigma_{m}<10~{}\mathrm{cm}^{2}/\mathrm{g}$, the halo core density follows a track defined by the outer halo parameters ($\rho_{s}$ and $r_{s}$), attaining a minimum value of $\sim 3\,\rho_{s}$. We found that mass loss from the outer regions can shorten the timescale of this evolution and increase the minimum SIDM density.

The accelerated evolution due to tidal stripping opens up the possibility of intermediate-mass BH formation in the rapid gravothermal collapse phase in dwarf halos for $\sigma_{m}\gtrsim 5~{}\mathrm{cm}^{2}/\mathrm{g}$. The increase in the core density due to gravothermal collapse could explain the high central density of the isolated Local Group dwarf Tucana, assuming it came close to the MW about $10~{}\mathrm{Gyr}$ ago. The same arguments suggest that the central dark matter density of satellites that have come close to the MW disk, like Draco and the ultrafaint dwarfs, would be enhanced for large cross sections. In SIDM models with a cross section per mass of a few $\mathrm{cm}^{2}/\mathrm{g}$, which are favored by fits to the rotation curves of field galaxies Ren et al. (2018), our results suggest that the structural properties of satellite galaxies should be correlated with their orbital histories.

Acknowledgements.

We thank Anna Kwa and Annika Peter, who contributed in the early stages of this work. We also thank Jason Pollack for useful conversations. HN was supported by NSF Grant No. 1519353, and MK was supported by NSF Grant No. PHY1620638.

Appendix A Numerical implementation

We describe our procedure, which follows that in Ref. Pollack et al. (2015), for solving the gravothermal evolution equations

[TABLE]

which describe mass conservation, hydrostatic equilibrium, the first law of thermodynamics, and heat conduction, respectively, in dimensionless variables. For an initial density profile of the NFW form, the characteristic scales of the halo are the NFW scale density $\rho_{s}$ and radius $r_{s}$. With these quantities, we write the dimensionless radius $\tilde{r}\equiv r/r_{s}$, 1D velocity dispersion $\tilde{v}\equiv v/v_{0}$, cross section per mass $\tilde{\sigma}_{m}\equiv\sigma_{m}/(4\pi r_{s}^{2}M_{0}^{-1})$, density $\tilde{\rho}\equiv\rho/\rho_{s}$, and luminosity $\tilde{L}\equiv L/[(GM_{0}^{2}){(r_{s}t_{0})}^{-1}]$. We use the mass scale $M_{0}=4\pi r_{s}^{3}\rho_{s}$, velocity scale $v_{0}=\sqrt{GM_{0}/r_{s}}$, and timescale $t_{0}^{-1}=a\sigma_{m}v_{0}\rho_{s}$. There are thus two main inputs: the initial density profile $\rho_{\mathrm{init}}$ and the SIDM cross section per mass $\sigma_{m}$.

The first step is the discretization of the spherical halo. We divide the halo into $N$ concentric shells, and the outer radii of the shells are logarithmically spaced between $\tilde{r}_{\min}=\tilde{r}_{1}$ and $\tilde{r}_{\max}=\tilde{r}_{N}$. In this work, we choose $\tilde{r}_{\min}=0.01$, $\tilde{r}_{\max}=100$, and $N=400$.

Extensive variables associated with the $i$-th shell are defined at $\tilde{r}_{i}$: $\tilde{M}_{i}$ and $\tilde{V}_{i}$ are the dimensionless mass and volume contained within the radius $\tilde{r}_{i}$, and $\tilde{L}_{i}$ is the dimensionless luminosity at $\tilde{r}_{i}$.

Intensive variables of the $i$-th shell are the dimensionless density $\tilde{\rho}_{i}$, pressure $\tilde{p}_{i}$, specific energy $\tilde{u}_{i}$, and velocity dispersion $\tilde{v}_{i}$; for concreteness, these quantities are defined at the midpoint of the $i$-th shell, $(\tilde{r}_{i}+\tilde{r}_{i-1})/2$. They are also related to one another through the equipartition theorem and ideal gas law: $\tilde{p}_{i}=\tilde{\rho}_{i}\tilde{v}_{i}^{2}=\frac{2}{3}\tilde{\rho}_{i}\tilde{u}_{i}$.

From the input initial density, it is straightforward to obtain $\tilde{M}_{i}$ and $\tilde{p}_{i}$ from Eqs. (3a) and (3b), respectively. Given a value of $\tilde{\sigma}_{m}$, we can determine $\tilde{L}_{i}$ from Eq. (3c). With all quantities initialized, we evolve the halo by taking a small heat conduction time evolution step, followed by adjustments to maintain hydrostatic equilibrium.

To perform the heat conduction time step, we rewrite Eq. (3c) as

[TABLE]

where the change in density is taken to be negligible (such that the mass within the shells remains constant) for a sufficiently small time step $\delta\tilde{t}$. This condition is ensured throughout the calculation by imposing $\delta\tilde{t}$ to be smaller than the minimum local relaxation time of the halo by $\epsilon_{t}(\ll 1)$. In particular, we require $\delta\tilde{t}=\epsilon_{t}\min_{0\leq i\leq N}(t_{r,i}/t_{0})$, where $t_{r,i}$ is the local relaxation time of the $i$-th shell, and we set $\epsilon_{t}=10^{-4}$.

The heat conduction causes an increase in the specific energy, which pushes the halo out of hydrostatic equilibrium. With the assumption of negligible change in $\tilde{\rho}_{i}$ in this step, the increase in $\tilde{u}_{i}$ is reflected in the following step through $\tilde{p}_{i}=\frac{2}{3}\tilde{\rho}_{i}\tilde{u}_{i}$. To return the halo to hydrostatic equilibrium, we adjust the shell location $\tilde{r}_{i}\to\tilde{r}^{\prime}_{i}\equiv\tilde{r}_{i}+d\tilde{r}_{i}$, pressure $\tilde{p}_{i}\to\tilde{p}^{\prime}_{i}\equiv\tilde{p}_{i}+d\tilde{p}_{i}$, and density $\tilde{\rho}_{i}\to\tilde{\rho}^{\prime}_{i}\equiv\tilde{\rho}_{i}+d\tilde{\rho}_{i}$.

The conservation of mass $M_{i}$ enforces $d\tilde{\rho}_{i}=-\tilde{\rho}_{i}(d\tilde{V}_{i}/\tilde{V}_{i})$, where $d\tilde{V}_{i}$ is the change in volume associated with the shift in radii. The conservation of entropy requires $\tilde{p}_{i}\tilde{V}^{5/3}_{i}$ to be constant after the adjustment. Through this process, the discretized form of the hydrostatic condition becomes

[TABLE]

for $i=1,\dots,N-1$. We define a fixed inner boundary shell $\tilde{r}_{0}=0$ with $\tilde{L}_{0}=0$ (i.e., no source of heat at the center of the halo), and we fix the boundary of the outermost shell such that $d\tilde{r}_{N}=0$.

When linearized, Eq. (5) takes up the form of $a_{i}\,d\tilde{r}_{i-1}+b_{i}\,d\tilde{r}_{i}+c_{i}\,d\tilde{r}_{i+1}=d_{i}$, where $a_{i}$, $b_{i}$, $c_{i}$, and $d_{i}$ are constants. The result is a set of $N-1$ tridiagonal equations for $N-1$ variables $d\tilde{r}_{1},\dots,d\tilde{r}_{N-1}$. After solving this system of equations, we update all other variables. The hydrostatic adjustment is repeated in between every heat conduction step until $\max_{0<i<N}|d\tilde{r}_{i}/\tilde{r}_{i}|<\epsilon_{r}$ is satisfied, where we set $\epsilon_{r}=10^{-3}$.

Appendix B Calibration

In the gravothermal fluid model, the two parameters $C$ and $b$ need to be calibrated against N-body simulations, as they adjust the efficiency of the LMFP and SMFP heat transfer process, respectively, in Eq. (3d). Reference Koda and Shapiro (2011) matched the collapse time of pure LMFP evolution (i.e., ${\sigma_{m}}^{2}/b\to 0$) to first calibrate $C$ for the case of an initially self-similar solution of an isolated halo. The value of $b$ was then decreased from $b=1.38$ to $0.25$ to match the collapse times of self-similar halos with larger values of $\sigma_{m}$. Instead of matching the collapse times, we determine the values of $C$ and $b$ by directly comparing the calculated density profiles after $13~{}\mathrm{Gyr}$ of evolution against the corresponding profiles of the Pippin halo, generated from a cosmological simulation Elbert et al. (2015).

Figure 4 shows the results of SIDM evolutions with different $(C,b)$ values: $S_{0}=(0.75,1.38)$, $S_{C}=(0.45,1.38)$, and $S_{b}=(0.75,0.003)$. Used in Ref. Essig et al. (2018), $S_{C}$ provides a better match of the density profile for $\sigma_{m}=50~{}\mathrm{cm}^{2}/\mathrm{g}$ but worsens agreement for lower $\sigma_{m}$. However, the Pippin halo with $\sigma_{m}=50~{}\mathrm{cm}^{2}/\mathrm{g}$ exhibits signs of core collapse, so we do not use it to calibrate $C$, which controls LMFP heat transfer. Instead, we calibrate $C$ against the Pippin halo for cross sections $\sigma_{m}<50~{}\mathrm{cm}^{2}/\mathrm{g}$, keeping $b=1.38$ fixed. The profiles with $S_{0}$ in Fig. 4 show good agreement for $\sigma_{m}\leq 10~{}\mathrm{cm}^{2}/\mathrm{g}$, but a clear discrepancy for $\sigma_{m}=50~{}\mathrm{cm}^{2}/\mathrm{g}$. We note that the results for $S_{0}$ and $S_{C}$ are essentially the same for the key points in this paper.

In the spirit of Ref. Koda and Shapiro (2011), we allow for adjustments in $b$ and find that $b=0.003$ lowers the density profile with $\sigma_{m}=50~{}\mathrm{cm}^{2}/\mathrm{g}$. This solution, however, never reduces the density as low as the $10~{}\mathrm{cm}^{2}/\mathrm{g}$ case requires, and it pushes the early evolution into the SMFP regime. In addition, the dispersion profile for $10~{}\mathrm{cm}^{2}/\mathrm{g}$ (as well as $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$) does not have the negative radial gradient at $t=13~{}\mathrm{Gyr}$ that the simulation exhibits.

In light of the discussion above, we see no compelling reason to adopt the $S_{b}$ solution. Instead, the key issue is likely that the heat conduction $\kappa^{-1}={\kappa_{\textrm{SMFP}}}^{-1}+{\kappa_{\textrm{LMFP}}}^{-1}$ does not correctly capture the physics of the transition between the SMFP and LMFP regimes. We could generalize this interpolation to $\kappa^{-\alpha}={\kappa_{\textrm{LMFP}}}^{-\alpha}+{\kappa_{\textrm{SMFP}}}^{-\alpha}$, such that

[TABLE]

where $x\equiv(aC/b){\tilde{\sigma}_{m}}^{2}\tilde{\rho}{\tilde{v}}^{2}$. This generalization would allow us to increase $C$ and tune $\alpha$ to obtain a faster evolution, while decreasing the impact of SMFP regime. We leave such an investigation for future work.

Appendix C Halo evolution

In Fig. 5, we show the time evolution of the halo central density with various values of $\sigma_{m}$ for an initially NFW profile. We use the dimensionless time $\tilde{t}=t/t_{0}$, so larger values of $\sigma_{m}$ correspond to longer dimensionless times. For the cross section $\sigma_{m}=10~{}\mathrm{cm}^{2}/\mathrm{g}$, we evolve the central density beyond $13~{}\mathrm{Gyr}$ to show the gravothermal catastrophe as the central density sharply rises at $\tilde{t}\gtrsim 350$.

There are a few important features to note. First, the curves for $\sigma\leq 5~{}\mathrm{cm}^{2}/\mathrm{g}$ align with the $\sigma_{m}=10~{}\mathrm{cm}^{2}/\mathrm{g}$ curve when extended, indicating that the evolution for $\sigma_{m}\leq 10~{}\mathrm{cm}^{2}/\mathrm{g}$ occurs in the LMFP regime at times of interest and is insensitive to the value of $\sigma_{m}$. Second, the deviations of the $50$ and $100~{}\mathrm{cm}^{2}/\mathrm{g}$ curves away from the lower $\sigma_{m}$ curves at small values of $\tilde{t}$ indicate that SMFP physics is not negligible at these early times. At late times, the large $\sigma_{m}$ curves both bifurcate from the self-similar solution. Finally, for all values of $\sigma_{m}$ we tested, there is a minimal central density of $\tilde{\rho}_{c}=2.8$ that the SIDM halo achieves at $\tilde{t}\simeq 45$.

To demonstrate the formation of the core and its subsequent collapse, we consider the evolution of an initially NFW halo with $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$. Figure 6 shows snapshots of the halo density profile, 3D velocity dispersion, and luminosity at 5 different times. We choose the same calibration $S_{0}=(C,b)=(0.75,1.38)$ and halo parameters ($\rho_{s}=0.019~{}\mathrm{M}_{\odot}/\mathrm{pc}^{3}$ and $r_{s}=2.59~{}\mathrm{kpc}$) used in the main text. The state of the halo in each panel of Fig. 6 is as follows:

$\tilde{t}=0$: Initial NFW profile. The initial profile satisfies $\kappa_{\mathrm{SMFP}}\geq\kappa_{\mathrm{LMFP}}$ throughout the halo for $\sigma_{m}\lesssim 10~{}\mathrm{cm}^{2}/\mathrm{g}$ (cf. Fig. 5). In fact, heat transfer remains dominated by the LMFP physics until the core of the halo becomes sufficiently dense at $\tilde{t}\simeq 374$. 2. 2.

$\tilde{t}=1$: Core expansion. Self-interactions thermalize the inner region of the halo, where the mean-free path $\lambda=1/(\rho\sigma_{m})$ is short. As a result, the temperature ($\propto{\tilde{v}}^{2}$) near the center begins to increase, and the particles are pushed towards $r\sim r_{s}$. This leads to the formation of a gradually expanding core. 3. 3.

$\tilde{t}=53$: Formation of the isothermal core. The luminosity becomes entirely positive (shortly after a minimum central density of $\tilde{\rho}_{c}=2.8\,\rho_{s}$ is achieved at $\tilde{t}=45$) after a period of core expansion, resulting from the thermalization within the innermost region of the halo. The core becomes isothermal, as $\tilde{v}$ is now nearly constant within this region. This time marks the turnaround of the central density evolution, seen in the left panel of Figs. 2 and 5. 4. 4.

$\tilde{t}=351$: LMFP core collapse. The constant-density core becomes denser as it shrinks in size, while the density profile maintains its shape since the evolution is self-similar. The temperature uniformly increases within the core such that it remains isothermal. This phase corresponds to the solid lines in Fig. 3. 5. 5.

$\tilde{t}=374.56$: SMFP core collapse. The condition $\kappa_{\mathrm{SMFP}}\geq\kappa_{\mathrm{LMFP}}$ is satisfied entirely within the core, and the halo central density deviates from the self-similar solution. At the time shown, the SMFP core has attained its maximum mass. The density and temperature begins to increase more rapidly within the SMFP core as the core becomes optically thick to self-interactions. This final phase of evolution is shown as dashed lines in Fig. 3; the core retains a substantial amount of its mass until it reaches the relativistic instability to form a BH.

Figure 7 shows analogous snapshots for the evolution of an initially TNFW profile at $r_{t}=r_{s}$ with $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$. We describe the state of the halo in each panel in the following:

$\tilde{t}=0$: Initial TNFW profile. The outer region of the halo ($r\geq r_{s}$) is assumed to be tidally stripped upon formation of an NFW halo, and the density is suppressed by a factor of $(r/r_{t})^{5}$ beyond the truncation radius. The truncation reduces the pressure of the halo beyond the truncation radius, due to enforcing hydrostatic equilibrium. Although the density is decreased beyond the truncation radius as well, the drop in pressure is more drastic, causing the velocity dispersion (or equivalently, temperature) profile to decrease as well. Thus, the outer part of the halo gets significantly colder, and the temperature peak around $r\sim r_{s}$ becomes smaller and narrower compared to the NFW profile without truncation. 2. 2.

$\tilde{t}=1$: Core expansion. Compared to the initially NFW halo in Fig. 6 at the same time, the core is closer to being fully thermalized due to the less pronounced peak in velocity dispersion. 3. 3.

$\tilde{t}=5.2$: Formation of the isothermal core. The luminosity becomes entirely positive, and the core expansion halts at a much earlier time than seen in Fig. 6. As a result, the isothermal core is more concentrated: its size is smaller, and its density is higher ($\rho_{c}\simeq 7.1\,\rho_{s}$). 4. 4.

$\tilde{t}=20$: LMFP core collapse. The core contracts slowly in the same way as seen for the initially NFW halo in Fig. 6 through the LMFP heat transport. 5. 5.

$\tilde{t}=52.99$: SMFP core collapse. The halo reaches the same state as the initially NFW halo at $\tilde{t}=374.56$. Collapse occurs $13.5~{}\mathrm{Gyr}$ after halo formation, implying that such a truncation causes a gravothermal catastrophe with $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$ on timescales relevant for what we may observe today.

Lastly, we show in Fig. 8 analogous snapshots for the evolution of an initially NFW halo that is truncated at $r_{t}=r_{s}$ after $3~{}\mathrm{Gyr}$ with $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$. We describe the state of the halo in each panel in the following:

$\tilde{t}=11.71$: Before truncation. The halo follows the evolution of the initially NFW profile (the first two panels in Fig. 6 are the same for this case). The negative luminosity within $r\leq r_{s}$ indicates that the core is gradually expanding. 2. 2.

$\tilde{t}=11.75$: After truncation. Compared to the initially TNFW halo in Fig. 7, the (instantaneous) late truncation at $\tilde{t}=11.75$ ($t=3~{}\mathrm{Gyr}$) results in a slightly different halo configuration. When the late truncation lowers the temperature of the halo beyond the truncation radius, the already well-thermalized core immediately becomes hotter than the $r\sim r_{s}$ region where the temperature used to peak. This flips the sign of the luminosity within the core from negative to positive and triggers the core contraction phase seen in the initially NFW halo after $\tilde{t}=53$ in Fig. 6. 3. 3.

$\tilde{t}=20$: $2~{}\mathrm{Gyr}$ after truncation. The halo properties have smoothed out after the sudden truncation, and the core is gradually contracting under the LMFP evolution. Although the central density is lower than the initially TNFW profile at $\tilde{t}=20$, the central velocity dispersion is higher due to the initial thermalization phase prior to truncation. 4. 4.

$\tilde{t}=51$: $13~{}\mathrm{Gyr}$ after the initial formation of the NFW halo. While the initially NFW halo reaches the minimum core density around the same time (panel 3 of Fig. 6), the core of the truncated halo has a density that is $\sim\!30$ times higher. 5. 5.

$\tilde{t}=60.10$: SMFP core collapse. At $15.3~{}\mathrm{Gyr}$, the innermost region of the halo is undergoing the same collapse process that its NFW counterpart reaches at $95.7~{}\mathrm{Gyr}$ after its formation.

Finally, in Table 1, we show the values of the mass of the SMFP core $M_{\textrm{SMFP}}$ and the 3D velocity dispersion at the center of the halo $v_{3\mathrm{D},c}$ when $M_{\textrm{SMFP}}$ reaches the maximum value. In the main text, we used these values in the empirical relation Balberg et al. (2002) to estimate the mass of BH that forms in the gravothermal collapse in the SMFP regime for the case of $\sigma_{m}=5~{}\mathrm{cm}^{2}/\mathrm{g}$. We note that the last panels of Figs. 6, 7, and 8 show the halo profiles at these times, the central regions of which have all deviated from the self-similar evolution.