Core-collapse, evaporation and tidal effects: the life story of a self-interacting dark matter subhalo

Self-interacting dark matter (SIDM) cosmologies admit an enormous diversity of dark matter (DM) halo density profiles, from low-density cores to high-density core-collapsed cusps. The possibility of the growth of high central density in low-mass halos, accelerated if halos are subhalos of larger systems, has intriguing consequences for small-halo searches with substructure lensing. However, following the evolution of $\lesssim 10^8 M_\odot$ subhalos in lens-mass systems ($\sim 10^{13}M_\odot$) is computationally expensive with traditional N-body simulations. In this work, we develop a new hybrid semi-analytical + N-body method to study the evolution of SIDM subhalos with high fidelity, from core formation to core-collapse, in staged simulations. Our method works best for small subhalos ($\lesssim 1/1000$ host mass), for which the error caused by dynamical friction is minimal. We are able to capture the evaporation of subhalo particles by interactions with host halo particles, an effect that has not yet been fully explored in the context of subhalo core-collapse. We find three main processes drive subhalo evolution: subhalo internal heat outflow, host-subhalo evaporation, and tidal effects. The subhalo central density grows only when the heat outflow outweighs the energy gain from evaporation and tidal heating. Thus, evaporation delays or even disrupts subhalo core-collapse. We map out the parameter space for subhalos to core-collapse, finding that it is nearly impossible to drive core-collapse in subhalos in SIDM models with constant cross sections. Any discovery of ultra-compact dark substructures with future substructure lensing observations favors additional degrees of freedom, such as velocity-dependence, in the cross section.


INTRODUCTION
Even after decades of being widely accepted as the most plausible answer to astrophysical anomalies such as the mass deficit in galaxy rotation curves (Rubin et al. 1978) and the spatial offset of the total mass center from the luminous matter center during cluster collisions (e.g., the Bullet Cluster; Clowe et al. 2006), dark matter's physical nature remains essentially unknown to us. The standard cold dark matter (CDM) paradigm has been well-tested on large, cosmological scales (Planck Collaboration et al. 2020;Alam et al. 2021), but is challenged on galactic scales (Bullock & Boylan-Kolchin 2017;Buckley & Peter 2018). The tensions between CDM predic-★ E-mail: zeng.408@osu.edu tions and galactic observations can be sorted into three categories. First, there is the 2000's-era mismatch between the small number of observed satellite galaxies and the large number of galaxy-scale subhalos in simulation -the missing satellites problem (Klypin et al. 1999;Moore et al. 1999). Recently, the discovery of many new satellites, a careful accounting to survey selection functions, and insights from simulation work instead now hint that, CDM may underpredict the abundance of luminous satellites (Willman et al. 2005;Zucker et al. 2006;Tollerud et al. 2008;Koposov et al. 2009;Walsh et al. 2009;Belokurov et al. 2009;Drlica-Wagner et al. 2015;Laevens et al. 2015;Torrealba et al. 2016;Garrison-Kimmel et al. 2017;Newton et al. 2018;Kim et al. 2018;Jethwa et al. 2018;Torrealba et al. 2018;Kelley et al. 2019;Homma et al. 2019;Nadler et al. 2020a;Kim & Peter 2021). Second, the dynamical mass in the inner region of dwarf galaxies (in both the dwarf spheroidal satellites of the Milky Way as well as rotation-supported field dwarfs) is low compared to CDM simulations without baryons, which is classically understood as an issue of core vs. cusp in the shape (slope) of density profiles (Moore 1994;Adams et al. 2014;Oh et al. 2015). Recently, this problem has expanded, with discoveries that observed dwarf galaxies display a large diversity in rotation curves as well as stellar dynamics (diversity problem; Kuzio de Naray & McGaugh 2014;Oman et al. 2015;Errani et al. 2018;Read et al. 2019;Relatores et al. 2019;Zavala et al. 2019;Santos-Santos et al. 2020;Li et al. 2020;Hayashi et al. 2020). The two former problems converge in the third problem, in the dynamical properties of the most luminous satellites of the Milky Way, showing that CDM halos whose central densities match the bright satellites are much lower in mass than expected from abundance matching (too-big-to-fail problem; Boylan-Kolchin et al. 2011, 2012Tollerud et al. 2014;Jiang & van den Bosch 2015;Kaplinghat et al. 2019).
Introducing baryons into simulations has been shown to have the potential of addressing the small scale problems within the CDM framework. Baryon physics determines which halos contain luminous baryons at all (White & Rees 1978;Gnedin 2000;Benson et al. 2002;Okamoto et al. 2008). Various baryonic feedback processes, such as stellar feedback and supernova explosion, help with the redistribution of dark matter and core formation in the halo center (Governato et al. 2012;Brooks & Zolotov 2014;Chan et al. 2015;Read et al. 2016;Wetzel et al. 2016;Fitts et al. 2017). The shallower potential, as a consequence, also leads to more satellite galaxies being tidally disrupted, decreasing the abundance of satellites predicted by CDM (D'Onghia et al. 2010;Brooks et al. 2013;Zhu et al. 2016;Garrison-Kimmel et al. 2017;Despali & Vegetti 2017). On the other hand, if the mass fraction of stellar objects is high (∼ 10%), this extra baryonic potential could possibly compensate for the feedback and make a cuspy central density, potentially fitting in the cuspy side of the diversity problem (Di Cintio et al. 2014;Tollet et al. 2016;Hopkins et al. 2018;Lazar et al. 2020).
Alternatively, modifications to CDM have been proposed to alleviate the tension on small scales while preserving the large scale success. A large category of such modifications is hidden-sector models, in which non-gravitational beyond-standard-model interactions among dark matter particles are generic features (Spergel & Steinhardt 2000;Battaglieri et al. 2017). In these self-interacting dark matter models (SIDM) the scattering, often (but not always) elastic, leads to momentum exchange between dark matter particles, thermalization of particles in halos, and thus the formation of cored density profiles in the center of the dark halo. This halo phenomenology has been demonstrated with numerical N-body simulations with only dark matter (Davé et al. 2001;Colín et al. 2002;Vogelsberger et al. 2012;Rocha et al. 2013;Zavala et al. 2013;Elbert et al. 2015;Vogelsberger et al. 2019), and analytically via isothermal Jeans equations Valli & Yu 2018;Robertson et al. 2021, Jiang et al in prep.). The shallower gravitational potential wells of subhalos result in them being more vulnerable to tidal stripping, potentially alleviating the tension in satellite abundance (Colín et al. 2002;Vogelsberger et al. 2012;Robles et al. 2019;Dooley et al. 2016), or otherwise causing this tension to grow with the satellite completeness correction ("toomany-satellites" problem; Kim & Peter 2021;Bhattacharyya et al. 2021;Nadler et al. 2020b).
Although low-density cores were long considered the trademark signature of SIDM, a dramatically different phenomenology may take over the central density of the SIDM halo, long after the core-formation phase. SIDM scatterings lead to heat outflow from the inner, hotter region to the outer, cooler layers of the halo. Such heat loss in the inner region results in the infall of dark matter particles to more bound orbits, where they become even hotter than before. Due to this negative heat capacity of the self-gravitating system, the trend of dark matter in-fall accelerates itself as the negative heat gradient steepens, leading to a denser and more cuspy central density, known as the core-collapse process or gravothermal catastrophe (Lynden-Bell & Wood 1968;Kochanek & White 2000;Colín et al. 2002;Balberg et al. 2002;Koda & Shapiro 2011;Pollack et al. 2015;Zavala et al. 2019;Essig et al. 2019;Nishikawa et al. 2020;Sameie et al. 2020;Turner et al. 2021;Correa 2021). Indeed, while low-density cores may be formed in other non-CDM cosmologies, core-collapse is a unique signature of SIDM, distinguishing it from other alternative DM models. The two phases of core-formation and core-collapse suggest that dark matter halos can have a large diversity of central densities and profiles in an SIDM Universe, thus potentially capable of solving the diversity problem of galaxy rotation curves (Kamada et al. 2017;Zavala et al. 2019;Sameie et al. 2020;Kaplinghat et al. 2020). The high density also makes subhalos robust to destruction, avoiding a "too many satellites" problem (Kim et al. 2018;Kelley et al. 2019;Kim & Peter 2021). However, the time scale for SIDM core-collapse to happen is much longer than the age of the Universe unless the cross section is large (Balberg et al. 2002;Koda & Shapiro 2011;Elbert et al. 2015;Essig et al. 2019, see also Sec. 4.1 in this work). Excitingly, it has been recently reported that tidal stripping by the host halo can noticeably accelerate core-collapse of the subhalo, since the removal of dark matter from the outer layers helps with the formation of a negative temperature gradient and makes heat outflow more efficient (Nishikawa et al. 2020;Correa 2021, also Sec. 4.2 in this work). This tidal acceleration of core-collapse has been used to explain the diversity of the Milky Way's dwarf spheroidal galaxies Sameie et al. 2020).
Apart from observations of the Milky Way's satellite galaxies, another way to detect core-collapse is through substructure lensing.
Substructure lensing is the perturbation to a strongly lensed system, due to the existence of dark matter substructures/subhalos in the foreground lens (with in-field halos along line-of-sight as a systematic; Mao & Schneider 1998;Dalal & Kochanek 2002;Chiba 2002). Properties of subhalos can be inferred from the perturbed lensing image, either by flux anomalies of lensed quasars (Dalal & Kochanek 2002;Xu et al. 2009;Gilman et al. 2018Gilman et al. , 2020a, or distorted images of lensed galaxies (Koopmans 2005;Vegetti & Koopmans 2009;Vegetti et al. 2010Vegetti et al. , 2014Minor et al. 2020), thereby setting constraints on dark matter physics. Observed substructure lensing systems so far have shown subhalo mass function in consistency with CDM (Dalal & Kochanek 2002;Vegetti et al. 2014;Hsueh et al. 2020;Gilman et al. 2020a). Core-collapsed subhalos, if they do exist, are promising targets because their highly concentrated mass distributions leave greater distortions on strongly lensed images than cored or CDM substructures (Gilman et al. 2021;Yang & Yu 2021).
A problem for theoretical predictions, though, is that substructure lensing can be used to probe subhalos down to 10 6 -a factor of 10 7 smaller than typical host halos as main lenses-which is too computationally expensive to resolve in cosmological simulations, especially since ensembles of simulated systems are required. We are thus motivated to build a hybrid semi-analytical + N-body method capable of simulating the evolution of such small SIDM subhalos in this work, where we implement the host halo in the form of an analytic density profile, while tracking the subhalo with N-body particles. Compared to previous works with a similar approach Sameie et al. 2020), our method includes not only gravitational interaction between the host and subhalo, but also the dark matter self-interaction. The latter, which we denote as "evaporation" for the rest of this paper (also known as "SIDM ram pressure" in Kummer et al. 2018, Nadler et al. 2020b, leads to energy gain and mass loss in the subhalo, and we show to be significant for the onset of the subhalo core-collapse (Sec. 3.1 and Sec. 4.3).
Unlike cosmological simulations (or zoom-ins of them) with SIDM which mostly concentrate on the host halo, or statistics of subhalos such as the mass function and radial distribution (Robles et al. 2019;Banerjee et al. 2020), in this study we mainly focus on the evolution and possible core-collapse of individual subhalos, using idealized ensembles of one-host-one-subhalo systems. Our work is structured as follows: in Sec. 2 we introduce our semi-analytic method for subhalo evaporation. In Sec. 3 we present our code validation procedure, which includes comparing to live host simulations for SIDM scenarios of constant and velocity-dependent cross sections, and for subhalos that remain cored or are core-collapsing. Sec. 4 is the heart of the paper, in which we present results from a series of production runs with our method. The major theme of this section is tracking the detailed evolution of SIDM subhalo central density, as a sign of whether and when the subhalo core-collapses, and describing how the relevant physical processes imprint themselves on the evolution of the central density. We also map out the parameter space needed for subhalo core-collapse, which generally requires an initial subhalo concentration much higher than cosmological values, when the SIDM model is velocity-independent. In Sec. 5 we summarize and discuss our work.

METHOD
There are two major effects from the host halo on the subhalo: the gravitational/tidal effects (including tidal heating and stripping), and the host-subhalo evaporation, rooted from the self-interaction between the host and subhalo dark matter. The latter, which matters because of the large gap between the velocities of host and subhalo dark matter particles in a paired two-body scattering, has not been included in previous semi-analytical approaches Sameie et al. 2020), nor systematically studied in the context of subhalo core-collapse. The evaporation effect blows away particles throughout the whole subhalo (note: we use the term 'evaporation' to denote the host-subhalo dark matter interaction throughout this paper, because typically most of the subhalo particles will become unbound -but not necessarily all of them; see Sec. 4.3 and Appendix D for more detailed discussion), leading to its additional mass loss and halting the process of core-collapse.
In this work, we use a quasi-analytic approach to model this evaporation effect. For each subhalo N-body particle, we evaluate the evaporation probability based on knowledge of the particle's position and velocity with respect to the host center. This approach, with any user-specified analytic host profile, can also be easily extended to embed additional analytic potentials, such as those due to baryons. The probability that a certain subhalo particle scatters off a virtual host particle ('virtual' because the simulation does not really include any host particles, we only sample them when evaluating the host-sub scatterings) is evaluated as where is the time between the previous and the current timestep, / is the self-interaction transfer cross section per unit mass (since SIDM particles are identical in the simulation, both forward and backward scatterings are suppressed in momentum transfer, thus we consider the definition of = ∫ Ω (1 − cos 2 ) Ω; see Tulin et al. 2013a;Cline et al. 2014;Boddy et al. 2016;Tulin & Yu 2018), ℎ is the local SIDM mass density of the host at the position of the target subhalo particle, and |v ℎ − v | is the mean relative velocity between this subhalo particle and virtual host particles in its vicinity. Since the scattering probability between any two DM particles depends linearly on their relative velocity, the expected evaporation probability ℎ also scales linearly with the mean relative velocity |v ℎ − v |. Here v is the velocity of the well-tracked subhalo particle. Since we do not have real host particles in the simulation, we evaluate this mean relative velocity with a statistical approach. We have explicitly checked the isotropy of host particle velocities with host-only runs, and found that ≈ ≈ ≡ ℎ and that ≈¯≈¯≈ 0. Thus we assume v ℎ to be a Gaussian distribution centered at 0 with standard deviation value ℎ , directly measure the radial distribution ℎ ( ) in these host-only runs, and load it in a tabulated form in the simulation to parameterize the velocity of the host halo particles. For convenience of calculation, we build the coordinate system such that -axis is aligned with v , finding where and ℎ are the norm of v and v ℎ , and Erf refers to the Gauss error function.
We incorporate our semi-analytical treatment of the evaporation process as a patch to the Arepo code (Springel 2010), which already has a well-tested SIDM module built in (Vogelsberger et al. 2012;Vogelsberger et al. 2014). At each timestep in our simulation, for each gravitationally active subhalo particle, it is probabilistically determined whether it scatters with another subhalo particle (probability = ), is 'evaporated' via scattering with a virtual host particle ( = ℎ ), or does not scatter at all ( = 1 − − ℎ ). To determine this, we compare the corresponding probabilities with a random number ∈ (0, 1). We evaluate of a subhalo particle labelled with the default SIDM module of Arepo, by summing up the two-body scattering probability over its neighbors. ℎ is evaluated according to Eqs. (1) and (2). In terms of the choice of the scattering partner of particle , if < , the first neighboring subhalo particle that satisfies ≥ is assigned to be paired with particle (Vogelsberger et al. 2012). If < < + ℎ , we sample the 3D velocity of a virtual host particle with a Gaussian distribution centered at 0 and has the standard deviation ℎ measured from our host-only SIDM simulation, to be the partner of subhalo particle . During the two-body scattering process, both particles have their velocities redistributed spherically symmetrically relative to the center of mass of this two-body system. After the scattering, the sampled virtual host particle velocity is discarded, therefore the analytic host is always static in our simulation and there is no feedback from the subhalo on the host, which is a reasonable approach when we choose a subhalo that is small compared to the host. In this work, we have only considered elastic, isotropic, and single-state scattering processes.
Because our goal is to determine when (sub)halos start corecollapsing, we must determine a time-stepping algorithm to resolve the core-collapse, and define a criterion for stopping the simulation before the core goes too far into the fluid regime. The central density of a halo going into the core-collapse phase can be high, and so is the scattering probability per unit time. Thus we dynamically confine the simulation timestep of each particle to make sure that its total scattering probability is always smaller than 1. The default SIDM module of Arepo has such a setup built in, evaluating the SIDM timestep (sidm) and gravity timestep (grav), and selecting the smaller of them to be . But it also limits how small (sidm) can be, relative to (grav), preventing the timestep from going arbitrarily small. In this work we effectively disable this lower limit on (sidm) and terminate the simulation if the halo central density grows to 100 times the initial central density. This termination is necessary for these particle-type simulations for two reasons. First, although we can currently afford the reduction of the timesteps by 3 to 4 orders of magnitude compared to a similar-scale CDM simulation, the computational cost is formidable when the central density gets much higher. Second, when the central density is too high, a larger number of particles are located near or under the force resolution limit, characterized by ∼ 3 times the softening length, as we will show in the code validation part, and also discussed in Turner et al. 2021. We argue that additional analytical methods or hydrodynamic simulation techniques (such as in Kummer et al. 2019) should be used to better understand the physics in this ultradense region, which is beyond the scope of this work.
To track the process of SIDM (sub)halo core-collapse, we measure the central density of the (sub)halo by marking the 50 simulation particles that have the greatest local densities, and defining the average as cen50 . This set of the top 50 particles ranked by local SIDM cross sections T / = 3 cm 2 /g, 6 cm 2 /g, 10 cm 2 /g, 30 cm 2 /g are shown in green, yellow, red and pink colors. One specific choice of velocitydependent SIDM model, the atomic dark matter model (ADM) from Cline et al. 2014, is shown in orange. The fine-tuned parameters of the ADM model are described in Sec. 3.2. Bottom: The corresponding host density profiles after evolution for 5 Gyr. Colors are the same as above. The CDM counterpart is shown by black dots. density is updated at each simulation timestep, not each snapshot output time. For our current core-collapse simulations, we define a tentative collapse criterion of where we terminate the simulation. Compared to fully live host runs, our semi-analytic method reduces the computational time by over two orders of magnitude for cored subhalos with 1/1000 the mass of the host, and by a factor of a few tens for core-collapsing subhalos.

CODE VALIDATION
In this section, we validate our evaporation model against fully live host simulations, testing both constant cross sections and a velocity-dependent cross section. A subhalo is injected into either an analytical or live host at its orbit apocenter. The evolution of its bound mass and its detailed density/mass profile are tracked and compared. Throughout this work, we use 200 to define the boundary of a halo, within which the mean density is 200 times the mean matter density of the Universe. We assume a flat ΛCDM cosmology with Hubble parameter = 70 km/s/Mpc and Ω = 0.3. We as well fix the host halo mass host = 10 13.5 ( 200 ∼ 1000 kpc), which is typical of a strong gravitational lens system (Birrer et al. 2017).
To implement the analytic host consistently, for each SIDM cross section, we measure the density profile of the live host halo after it reaches equilibrium with DM self-interaction and has a cored host density profile (see Sec. 3.1 for details). Some of the SIDM cross sections we use in this paper and corresponding host density profiles are presented in Fig. 1. We then load these host density profiles in tabulated form to Arepo, and log-linearly interpolate (i.e. log[ ℎ ] scales linearly with log[ ]) between these measured data points, as the input ℎ ( ) in Eq. 1. This external input of ℎ ( ) also accounts for the host gravitational field. To test SIDM models with a velocity-dependent cross section (which we plan to explore in detail in future work), we select an atomic dark matter model (hereafter 'ADM') from Cline et al. 2014 for the purpose of validation. We present its T / − rel relation in Fig. 1a and the host density profile as the orange scattered dots in Fig. 1b.
We note that in this work, the analytic host density profile remains static in time, so effectively we are simulating the evolution of late-accreted subhalos, after the thermal equilibrium within the host center has been reached. A time-dependent analytic host density will be used for a more realistic modelling of SIDM subhalo population in the future.
As we show in the next two sections, our semi-analytic evaporation method is in good agreement with the live host simulations, in both the low-density-cored and core-collapsing subhalos. However, we find two subtle systematics that, while not affecting our overall results, are important to highlight and describe here before we show our validation results. First, in the live host simulations, the subhalo accretes dark matter from the host as it travels along its orbit. However, we are able to distinguish the captured host dark matter from the original subhalo dark matter when measuring the mass profile of the subhalo. This systematic of host-particle-capturing accounts for only ∼ 1% for subhalos with mass 1/1000 of the host.
Second, dynamical friction in the live host runs results in the shrinking of subhalo orbits, thus affecting both evaporation and tides. When the whole host is modelled as an analytical field, dynamical friction cannot be easily represented without introducing more complex semi-analytical models (Peñarrubia & Benson 2005;Boylan-Kolchin et al. 2008;Petts et al. 2015) that would add to the complexity of this work. Dynamical friction slows down the bulk motion of the subhalo as it orbits around the host, leading to a stronger orbit decay and earlier arrival at pericenter, where both evaporation and tidal effects are strongest. However, its importance is significantly reduced when the subhalo mass is small relative to the host. We present a comparison between sub-host mass ratios of 1:1000 and 1:100 in the next section and show that the discrepancy is much smaller for the former case, and that the accuracy of the 1 : 100 case is noticeably improved with inclusion of the dynamical friction model from Petts et al. 2015 (see Appendix A for details).
Thus we find that our approach of including an analytic host works best for the evolution of relatively small subhalos. This is the regime that is computationally challenging for simulations with live hosts, and modelling these small halos is a main reason behind our analytic approach. In short, we show that the systematic errors in our approach become vanishingly small for our target use cases.

SIDM with constant cross section
In this section, we validate our method of subhalo evolution for SIDM with a constant cross section T / = 6 cm 2 /g. The simulated system is composed of a host halo with mass of 200 = 10 13.5 , and a subhalo with mass ratio 1:1000 (or 1:100 for a comparison case) with respect to the host. This choice of host-subhalo mass ratio is aimed to reduce dynamical friction and accretion of host dark matter as much as possible, while keeping the live host simulations computationally inexpensive.
Both the host and the subhalos are generated with initial NFW profiles using the SpherIC code . The concentration 200m of the host is 6.5, while that of the subhalo is set to be 14.8, 40 and 80 for three control groups. As a reference, the cosmological value of 200m for the subhalo at redshift = 0 is around 24 (Duffy et al. 2008;Komatsu et al. 2009). Thus our three initial concentrations of subhalos include one less concentrated than the cosmological mean, one more concentrated, and the extremely high concentrations of 80 for the subhalo is artificially set, meant to test our code in the context where core-collapse of the SIDM subhalo may occur for relatively low values of the self-interaction cross section. We evolve the host halo in isolation with dark matter self-interactions for 5 Gyr before we insert the subhalo, so that the host comes to equilibrium, and its density profile is evolved from the initial cuspy NFW to a cored state. The density profile and velocity dispersion of the host are then measured and used as host model for the analytic evaporation modelling. The SIDM subhalos are not preevolved prior to infall in the same way, because the timescales for core-formation and core-collapse are different for subhalos with different concentrations (Essig et al. 2019). Such pre-evolution could introduce bias between control groups. In future work, we will treat this point more holistically in a semi-cosmological context, introducing different infall times and thus various pre-evolution times of subhalos.
The subhalo is initially placed at its apocenter at apo = 0.7 200 from the center of the host, following the subhalo population model of Peñarrubia et al. 2010. The magnitude of its velocity is determined such that the distance between the pericenter of its orbit and the host center, peri , is 1/10 of its initial displacement apo . In later sections of our production runs, we will fix this apocenter distance apo and vary the pericenter-apocenter ratio, peri / apo , to characterize different subhalo orbits.
The mass of each simulation particle is set to be = 10 6 for our validation runs. For the force resolution, we follow the criteria of Equation (20) and (21)  .
where is the initial number of particles in a halo, is the concentration and = 200 / is the scale radius. Because our main goal is to make sure that the subhalo is well resolved, these , and are all drawn from the initial condition of the subhalo instead of the host. For example, the simulations with a subhalo of 10 10.5 with initial = 14.8 have = 0.294 kpc while for the subhalos with = 80 we have = 0.035 kpc. To test the robustness of our evaporation mechanism in terms of the particle mass resolution, we have also prepared corresponding runs with 10 times more particles, denoted as '10 ' hereafter, where the particle mass is set to = 10 5 instead of the fiducial choice of 10 6 . Note that we only run the high-resolution tests for the analytical potential cases, and not the live host simulation runs.
As a (sub)halo core-collapses, its central region contracts to become an ultra-dense structure, of which the length scale might be comparable to, or even smaller than, the gravitational softening length. A typical/average particle distance within the dense core of a core-collapsing (sub)halo, for example for the subhalo with [ sub = 10 10.5 , = 80] and simulation particle mass = 10 5 , is about 7 pc, while the softening length of gravity is = 14 pc as set by Eq. (4). Generally any structure at the subhalo center is not numerically robust given that the average particle separation is below the gravitational softening length. We prepared a companion run with a smaller softening length by an order of magnitude, = 1.4 pc, to test the convergence of our force resolution. We find that the evolution of the central densities of these two runs are nearly indistinguishable, with only < 3% difference in their core-collapse times, suggesting that the default force resolution is reasonable up to the termination point. This is probably because the (sub)halo only enters such a dense regime after the onset of core-collapse, which is a quick, runaway process and discrepancy of force resolution does not have the time to accumulate yet. Therefore, we suggest that our default force resolution set by Eq. (4) should be generally good enough for SIDM subhalo evolution (see more discussion in Sec. 5).
In this section, we fix the cross section of the dark matter selfinteraction to be T / = 6 cm 2 /g, which is much higher than the constraint of T / 1 cm 2 /g for group-scale halos such as our hosts (Yoshida et al. 2000;Natarajan et al. 2002;Rocha et al. 2013;Peter et al. 2013;Tulin & Yu 2018;Elbert et al. 2018), but plausible to produce core-collapsed subhalos (Nishikawa et al. 2020;Sameie et al. 2020;Correa 2021). We will explore a wider range of cross sections in Sec. 4.
The mass-loss histories of the control groups of subhalos, the discrepancy between the mass loss of subhalos with live hosts and with our evaporation model, and the corresponding host-subhalo separations are plotted in Fig. 2a. We use Amiga Halo Finder (AHF Knollmann & Knebe 2009) to track the subhalo evolution in each snapshot. We compare our evaporation modelling against live host simulations (dashed lines vs scattered dots) with cases of sub-host mass ratio 1 : 1000 and initial concentration = 14.8, 40, 80 (red, orange, blue), and mass ratio 1 : 100 with = 40 (pink). The = 80 group simulations (blue) all end well before the time limit of 10 Gyr because we terminate them when the halo central density reaches the termination criterion described in Eq. (3), just as the subhalo approaches the fluid phase of core-collapse. We find that for the mass ratio 1 : 1000 cases, the discrepancy is generally smaller than 10%. The spike-like features in the discrepancy panel at each pericenter are clear signatures of missing dynamical friction in our model, that the insufficient orbit decay results in an overall time delay in the subhalo mass loss history in our modelling.
We validate this hypothesis by increasing the subhalo mass by a factor of ten, so that the mass ratio is 1 : 100 (pink, dots vs dashed line). Compared to 1 : 1000 subhalos, this 1 : 100 subhalo suffers from a much larger, yet also diverging, discrepancy in the massloss history, due to the stronger orbital decay caused by enhanced dynamical friction. We then add a tentative correction for dynamical friction dyn on top of our model, based on Petts et al. 2015 (see Appendix A for details), as shown in the pink dash-dotted line.  Figure 2. The mass-loss history and orbital distance between subhalo and host for different control groups. a) Top: Remaining mass fraction of the subhalos vs. time. Subhalos with initial concentrations of = 14.8, 40, 80 are in red, orange and blue colors separately, and an extra group for = 40 but with a 10 times more massive subhalo is shown in pink, to highlight the importance of dynamical friction dyn . The blue data set ends much earlier, because the = 80 subhalo core-collapses and the simulation is terminated. The live host simulations are shown in scattered dots, while our evaporation modelling in dashed lines. Middle: The ratio between the remaining mass fraction of subhalos within our evaporation model and the ones from live host runs. Bottom: The host-subhalo separation as a function of time. b) Comparison between cases with evaporation model turned on (dashed lines) and off (dotted lines), and also the convergence test over the particle mass resolution (dashed lines for default resolution vs. dot-dashed lines for 10 times more particles.) Figure 3. The density profiles of subhalos with T / = 6 cm 2 /g, sub-host mass ratio 1:1000 and peri : apo = 1 : 10, for six selected time stamps. Subhalos with initial = 14.8, 40, 80 are in red, orange and blue colors. Different control groups with a live host; with an analytical host plus our evaporation model; with an analytical host but without evaporation; and with our evaporation model and also with particle mass resolution increased by a factor of 10, are shown in scattered dots, dashed lines, dotted lines and dash-dotted lines respectively. The '+' data sets for = 14.8 and = 80 subhalos in the (< ) vs panels (may need to zoom in to see) represent the mass profile of subhalos after removing particles that were bound to the host at the beginning but are captured by the subhalo in later times. The six times are: a) the beginning of the simulation; b) near the core-collapse time for the = 80, evaporation-off case (blue-dotted line); c) after the first pericenter, halted core-collapse process for other = 80 cases with evaporation (note that central decreases by nearly an order of magnitude compared to panel b); d) the core-collapse time for remaining = 80 cases; e) near the second pericenter; f) the end of the simulation. Readers are referred to Appendix B for a full evolution movie.
With the inclusion of this dynamical friction model, our evolution model for the 1 : 100 subhalo is noticeably improved in its accuracy, in terms of both the orbit and mass loss history. However, we do not implement this model of dynamical friction as a regular part of our semi-analytical model, because dynamical friction remains small for our main targets, the small (mass ratio 1:1000 and smaller) subhalos, where our evolution model has achieved satisfactory accuracy. Furthermore, our implementation of this dynamical friction model itself is not perfect (see Appendix A for more discussion), and potentially introduces additional uncertainties, such as the orbit over-decay we already see in the 1 : 100 subhalo (see the orbit panel of Fig. 2a). Thus we only use the dynamical friction model as a demonstration here, and note that it is not included in the rest part of this work. Another source of discrepancy between our subhalo evolution model and a live host simulation is the capture of host dark matter by the subhalo. Its initial offset is shown in the middle, Δ sub , panel of Fig. 2a, at time 0, denoting the host particles that immediately become bound to the subhalo when it is placed in the host. We can see that for subhalos with sub-host mass ratio 1 : 1000, captured host dark matter only counts for around 1% of its total mass, while it is ∼ 10% for the 1 : 100 subhalo. Thus we expect that for higher mass ratios, greater than 1:1000 (e.g. we consider the mass ratio of 1:100 to be 'smaller' than 1:1000), which is the physical case driving the creation of this semi-analytic evaporation model, both dynamical friction and the accretion of host dark matter would be a smaller factor, and the discrepancy between our code and a full live host SIDM simulation would be further reduced. Note that the subhalo continues to capture dark matter from the host while it orbits, but this still only counts for a small portion of the subhalo mass, as we will see later in the subhalo's density profiles in Fig. 3.
In Fig. 2b we show convergence tests for the particle mass resolution of our evaporation code, comparing the default resolution cases with ones with 10 times more particles (dashed vs dashdotted). The level of convergence is remarkable for all the groups of simulations. The mass loss histories of the higher resolution runs are nearly indistinguishable from the ones with lower resolution, for = 40 and 80. We observe differences of around 10% for the subhalo with lowest concentration at late times, when the subhalo has already lost more than 95% of its initial mass. The convergence regarding particle mass resolution can be also seen later in the detailed evolution of the subhalo's density profile in Fig. 3.
In Fig. 2b, we also show how the mass-loss history of the subhalo changes if our evaporation model is turned off (dotted) compared to our fiducial case in which both tidal and evaporative mass loss are included (dashed lines). The differences in mass-loss history start near the first pericenter, where both tidal mass loss and evaporation are at their highest rates, regardless of the subhalo's initial concentration. The mass-loss histories diverge with increasing time. For our particular choice of orbit and cross section, the simulations that include evaporation indicate that the subhalos are much less massive than in the tidal-only cases at the end of the simulation (10 Gyr). The largest difference at the end of the simulation is for the lowest concentration halo-a difference of nearly an order of magnitude. This highlights the necessity of including the physics of host-subhalo evaporation properly.
We explore the evolution of the density profile in Fig. 3. We present the detailed evolution of the subhalo density profiles and mass profiles for six selected snapshots: a) the beginning of the simulations; b) near the core-collapse time of the 'evaporation off' semi-analytical case with = 80 (blue dotted line); c) after the first pericenter passage, the disruption of core-collapse of other = 80 cases by the host-sub evaporation (their central densities decrease by about an order of magnitude); d) core-collapse of other = 80 cases (blue dots, dashed, dash-dotted); e) the second pericenter passage; f) the end of the simulation. Readers may refer to Appendix B for a full evolution movie of density/mass profiles.
As can be seen in Fig. 3, we find that the density/mass profiles of subhalos in our semi-analytical modelling of the evaporation are in overall good agreement with the live host simulations (dashed lines vs dots). Our model is able to capture the fast growth of central density of the core-collapsing subhalos, as well as the evolution of cored subhalos. By contrast with the tidal-field-only cases (dashed vs dotted lines), we can see the necessity of including evaporation. In this figure we show again the robustness of our method in terms of particle mass resolution, where the control groups with higher and lower resolutions (dash-dotted vs dashed lines) show no systematic disagreement. We notice some discrepancy between our model and the live host runs in Fig. 3f, at late times in the simulation, especially for the lowest concentration subhalos. One reason for this is that, as shown in the orbit panel of Fig. 2a, there is small but non-negligible extra orbit decay induced by dynamical friction in the live host runs, which results in stronger evaporation at this end time.
Overall we note that our method shows good agreement with live host simulations in terms of both the total mass and detailed density/mass profile, for SIDM with constant cross sections.

SIDM with velocity-dependent cross section
Although we use constant cross sections for most parts of this work, we show that our analytic evaporation model also extends to velocity-dependent cross sections. This extension will be important for our future work because, as we show later in Sec. 4, SIDM with a constant cross section is unlikely to produce core-collapsed subhalos, even in the absence of known cross section constraints. And, it is challenging to fit observations from clusters to satellite galaxies with a single constant cross section (Tulin & Yu 2018). Moreover, most theoretical models from particle physics predict velocity-dependent cross sections (Feng et al. 2009;Boddy et al. 2014;Tulin & Yu 2018). Therefore, we validate our code with a specific model of velocity-dependent SIDM in this sub-section.
We choose an atomic dark matter (ADM) velocity-dependent cross section model from Cline et al. 2014 (see also Goldberg & Hall 1986;Kaplan et al. 2010Kaplan et al. , 2011Cyr-Racine & Sigurdson 2013) where dark matter consists of 'dark atoms' composed of a dark proton and a dark electron. We manually fine-tune the input parameters, the dark atom fine structure constant = 0.05, dark atom mass d = 50 GeV and the dark proton-to-electron mass ratio d = 400, in order to generate a velocity-dependent cross section profile that is of order 10 cm 2 /g for velocities close to the subhalo velocity dispersion of O (10) km/s, and of order 1 cm 2 /g near the orbital velocity of O (100) km/s (see the orange line of Fig.  1a). The subhalo velocity dispersion is the key velocity scale for subhalo core-collapse, while the subhalo orbital velocity is key to host-subhalo evaporation strength (Nadler et al. 2020b;Jiang et al. 2021). This is meant to test our method in a scenario where the selfinteraction among subhalo dark matter is efficient, thus plausible for driving core-collapse, yet also with a relatively weak but nonnegligible evaporation from the host.
The setup of ADM subhalos is nearly the same as in Sec. 3.1, only that we set peri / apo = 1/20 for the orbit because the evaporation effect is stronger at a closer pericenter, and that the ADM host halo has a smaller core than hosts with constant cross sections (see Fig. 1b). In Fig. 4, we present the density evolution of ADM subhalos for selected time stamps: a) the snapshot right before the core-collapse of = 80 subhalos; b) the second pericenter encounter, where both evaporation and tidal effects are strongest; c) the last snapshot of the simulation. As expected from the low cross section around the typical bulk velocity of the subhalo (of order 100 km/s), the evaporation is relatively low in the ADM cases compared to those in Sec. 3.1, thus including it or not (dashed vs dotted lines) only makes a small difference in both the corecollapsing and cored subhalos for our choice of initial conditions. For this particular set of initial conditions with high concentration, the weak evaporation does not disrupt the core-collapse for the = 80 subhalo, compared to the constant cross section case T = 6 cm 2 /g above ( Fig. 3b to Fig. 3d). However, for the cored subhalo of = 14.8, including evaporation or not leads to a ∼ 50% difference in its final mass in the ADM case.
This set of validation tests with velocity-dependent cross section further demonstrates the remarkable accuracy of our method, with an even smaller discrepancy from the live host simulation compared to the constant cross sections.

RESULTS
In this section, we present results from our production simulation runs, with which we study how different physical processes drive the evolution of a (sub)halo central density, including the intrinsic corecollapse of an SIDM halo due to its internal heat outflow, the tidal effects of a host potential, and host-subhalo evaporation. We then scan the relevant parameter space to show how these processes depend on SIDM subhalo properties, and eventually determine whether a subhalo core-collapses or not. We map out the critical boundary for subhalo core-collapse in the parameter space. All simulations in this section are conducted with the semi-analytical treatment described in the previous two sections, hence we note that analysis on subhalos below are most accurate for small subhalos (< 1:1000 host mass). We summarize our simulation suite in Table 1. In Sec. 4.1, we simulate SIDM halos in isolation and evolve them until corecollapse, in order to provide a numerical scaling relation of halos' intrinsic core-collapse time c as a function of input parameters. Our time scaling is compared to a similar one from Essig et al. 2019. This is an additional way to validate our simulations, and to explore a regime that has previously been missed in similar numerical fits for scaling relations for core-collapse times.
We scrutinize the acceleration of subhalo core-collapse by tidal effects Nishikawa et al. 2020;Sameie et al. 2020;Correa 2021) in Sec. 4.2, including tidal stripping and tidal heating. Tidal stripping refers to the process by which outlying dark matter is removed from the subhalo by tidal forces from the host halo. Tidal heating, which is, to first order, the work done by the tidal force on the subhalo dark matter relative to the subhalo center of mass, heats up the subhalo as it moves along its orbit (Gnedin et al. 1999a,b;Gnedin & Ostriker 1999;Taylor & Babul 2001;Pullen et al. 2014; Van den Bosch et al. 2018;Yang et al. 2020). The strength of tidal heating at different layers of the subhalo increases with distance from the subhalo center. We show that the acceleration of subhalo core-collapse by tidal effects is significant for subhalos with relatively low initial concentrations, but reverses (i.e., tides decelerate core-collapse) at ultra-high concentrations ( 75 for our choice of sub , T / and orbit). We argue that the net impact of the tidal field arises from the competition between c and the orbital time.
Sec. 4.3 serves two purposes. First, we contrast the subhalos with full evaporation against ones with the same tidal field but no host-subhalo evaporation. This allows us to highlight the importance of properly modelling evaporation for SIDM subhalo evolution. Second, by putting these subhalos on circular orbits in this subsection, we make the interaction among different physical processes clearer by maintaining nearly time-invariant strength of both evaporation and tides, before we dive into more realistic and complex cases later. We demonstrate three relevant processes as heating/cooling (energy gain/loss) terms: internal heat outflow from the subhalo core as a cooling term, which is driven by the self-interaction of the subhalo dark matter; and evaporation and tidal heating as heating terms. We correlate the evolution of these heating/cooling terms with that of the central density cen50 , the indicator of whether or when the subhalo core-collapses, and come to our primary conclusion that the subhalo central density grows only when there is a net cooling (energy loss) within the subhalo core.
We generalize our analysis in Sec. 4.4, generating subhalos with four varying parameters: SIDM cross section T / , the subhalo mass sub , the subhalo initial concentration , and its orbit characterized by peri : apo . We confirm our conclusion from the circular orbit runs, that the evolution of the subhalo central density is driven by the net heating/cooling within its core. We show that whether an SIDM subhalo can eventually core-collapse is extremely sensitive to its initial parameters. Therefore, in Sec. 4.5 we scan the multi-dimensional parameter space with a few hundred simulations, and map a boundary of critical parameters for SIDM subhalos with constant cross sections to survive the host-sub evaporation and eventually core-collapse.

Core-collapse time scaling for an isolated halo
In this section, we provide fits for a scaling relation of the corecollapse time c of an isolated SIDM halo using our particle-based simulations, and compare them with similar ones from Essig et al.
where s and s are the scale radius and normalization density of NFW profile, and is a constant factor. Alternatively, the collapse time can be described by virial quantities, In this work, we prepare a set of similar isolated SIDM halos by sampling T / , 200c and 200c to check with the power-law exponents of Eq. (6). Note that our overall definition of halo mass 200m throughout this work is counted in the region within which the average density is 200 times cosmological mean mass density. For this section only, we use a halo mass definition of 200 , the '200 times the critical density' definition of Essig et al. 2019, for a fair comparison.
As specified in Sec. 2, we measure the central density of the halo cen50 at each simulation timestep to track the process of SIDM halo core-collapse. We define the halo core-collapse state when cen50 grows by two orders of magnitude, and terminate the simulation. This criterion is approximately comparable to the ∼ 1 core-collapse state in Essig et al. 2019 (see their Fig. 1), and being the 1D number density, velocity dispersion and mass density. In Essig et al. 2019, the halo core-collapse state is defined at ∼ 0.1. However, since the time difference between the two states ∼ 1 and ∼ 0.1 is negligible relative to c , as the core-collapse process exponentially accelerates itself once initiated, our criterion of halo core-collapse is comparable to that of Essig et al. 2019.
We sample isolated halo simulations in three groups, within each fixing two parameters of [ / , 200c , 200c ] and varying the other one. The scaling relations of c against these three parameters are presented in We can confirm that the scaling of c follows almost perfectly a power law as observed in Essig et al. 2019, with the slopes of all our parameters slightly smaller in magnitude than their fitting results. The differences in our slopes are physical in origin. , because we need higher concentrations and more frequent scattering to drive core-collapse within the simulation time limit of 10 Gyr. As can be seen in Fig. 2 of Essig et al. 2019, the slope of c indeed gradually decreases withˆ, until reaching the minimum of at about = 5. This is because the heat conduction is actually suppressed when the dark matter scattering enters the highly-frequent regime and thus becomes more localized (Balberg et al. 2002;Agrawal et al. 2017;Essig et al. 2019). Other than this we conclude that our scaling relation of c is in agreement with that of Essig et al. 2019.
We combine the c of all the halos in the three subplots of Core-collapse occurs in a shorter time for more massive (thus larger velocity dispersion), more concentrated (thus larger inner density) halos with higher cross sections, because the heat conduction is boosted by the self-interaction before it is suppressed when the mean-free-path of self-interaction becomes too small. The reversal in the high-frequency scattering regime actually hints that the ultimate fate of a core-collapsing (sub)halo may not be collapse to a singularity, but instead the formation of a dense core in thermal equilibrium that has lost thermal contact with the outer parts of the (sub)halo, unless the core is so dense that the relativistic instability takes over (Shapiro & Teukolsky 1986;Pollack et al. 2015;Feng et al. 2021b,a). We also refer readers to more recent works on SIDM and the formation of black holes, such as Feng et al. 2021b, Choquette et al. 2019, Shapiro 2018, Cruz et al. 2021and Di Cintio et al. 2017.
As we argue later in the following sub-sections of Sec. 4, this trend in core-collapse time with concentration, cross section, and halo mass observed in isolated halos is preserved to some extent for subhalos as well.

Tidal effects on subhalo core-collapse
As has been observed in both semi-analytical treatments of truncating the subhalo (Nishikawa et al. 2020;Correa 2021) and N-body simulations of SIDM subhalos with an analytical host potential Sameie et al. 2020), tidal effects that detach dark matter from the outskirts of an SIDM subhalo enhance the formation of the negative temperature gradient, and thus accelerate the core-collapse process (see also Sec. 1 for more context). In this section we run simulations as in Sameie et al. 2020, with a similar setup of host and subhalo as in Sec. 3.1 ( host = 10 13.5 , sub = 10 10.5 , T / = 6 cm 2 /g, r peri : r apo = 1 : 10). Our goal is to understand the role that the tidal field plays in subhalo core-collapse, which is more complex than shown in prior work, and which depends on the timescales of orbiting and core-collapse.
We use the (sub)halo central density cen50 as the indicator  of whether and when the (sub)halo core-collapses, and show the comparison of subhalos and the same halos in isolation with black and pink colors in Fig. 6. cen50 ( ) is updated at each simulation timestep, and saved for output when a ±20% change is detected, i.e.
(a) (b) (c) Figure 7. The density profiles and radial velocity dispersion profiles of SIDM halos in isolation (pink lines) or in a host potential (black lines), with initial concentration = 60 (dashed lines) or 75 (solid lines), and fixed cross section T / = 6 cm 2 /g and sub = 10 10.5 (with analytic host mass 10 13.5 ). A same CDM halo in the tidal field is shown in the yellow dotted line as a comparison. The three time stamps are: a) core-collapse time of the = 75 isolated halo, which is earlier than the core-collapse of the = 75 subhalo in a tidal field; b) after the first pericenter, both = 60 halos are relaxed again; c) core-collapse of the = 60 subhalo in a tidal field, which is earlier than its counterpart for the same halo in isolation. Readers are referred to Appendix C for a full evolution movie. the data points are equally spaced along the -axis but not -axis for the top panel in order to capture the feature where the central density exponentially grows in a short time. In this figure, we consider four sets of simulations with initial concentrations = 45, 60, 75, 90. Our results confirm the tidal acceleration on subhalo core-collapse for relatively low , in our cases of = 45 or = 60. However, the net tidal effect on subhalo core-collapse is reversed at high concentration ( = 75), and the effect of tides nearly vanishes at our highest concentration ( = 90). This trend was not observed in previous studies. To understand the tidal deceleration on corecollapse that we first saw in the = 75 case, we select the two groups = 60 and = 75 and scrutinize the detailed evolution of their densities and radial velocity dispersion as shown in Fig. 7. We also include CDM simulations for comparison.
The leftmost figure of Fig. 7 presents the first pericenter encounter, where the tidal effects are the strongest. For isolated halos (pink), their velocity dispersion is high at the center and dropping outward, with a flat core at the center for the lower concentration case and a steep slope for the higher concentration, a sign that corecollapse has started for the latter. However, for the subhalos (black lines for SIDM, orange for CDM), we observe a rise in the velocity dispersion from the center to a peak at large radii due to tidal heating, for both = 60 and = 75. The sharp drop of radial velocity dispersion beyond that peak, at large , is due to the heated particles being stripped away by the tidal force. The small peak-like feature at around = 10 kpc is moving outward with time (see the animation in Appendix C for details), showing the outflow of inner dark matter particles that have been heated up by tidal heating, which is also shared by the CDM subhalo (yellow line). This well-synchronised movement of the peak among CDM and SIDM models indicates that this outflow is indeed in the form of particles rather than heat transfer via dark matter self-interaction.
Tidal heating affects the evolution of core-collapse in subhalos. The core-collapse of SIDM halos relies on the formation of the negative temperature gradient, which makes core-collapse a selfaccelerating process. But when tidal heating is strong, the negative temperature gradient of the subhalo is disrupted, because dark matter at larger radius of the subhalo are heated more (see Sec. 4.3 for a more quantitative description). This blocks the heat from being transferred outward. Thus the core-collapse is delayed for the = 75 halo when it is in the tidal field as we see in its density profile (upper panel of Fig. 7a). However, tidal effects have a net acceleration on the core-collapse of the = 60 subhalo in the same potential field, because after the pericenter passage the whole subhalo regains virial equilibrium and forms a steeper velocity gradient due to the loss of dark matter heated up and stripped by tidal forces, as seen in Fig. 7b, and eventually core-collapses earlier than its isolated counterpart as seen in Fig. 7c. We can see that while tidal stripping accelerates the core-collapse of a subhalo, tidal heating plays a mixed role, that it delays the core-collapse when the subhalo approaches the pericenter, but becomes an accelerant after the pericenter. But, tidal heating and stripping, both rooted from the same tidal force, are not completely separable. Thus, tidal effects as a whole have a mixed impact on the core-collapse of SIDM subhalos, and simple truncation of the subhalo density profile does not capture all of it, as tidal heating is left out.
As shown in Fig. 6, for the plunging orbit of peri : apo = 1/10, the mixed impact of the tidal field on the central density is a result of the competition between time scales of the subhalo intrinsic (when in isolation) core-collapse c and its orbital motion period orb .

Tidal-field-only vs. full evaporation effects on subhalo evolution: a case study with circular orbits
In this section, we simulate a set of subhalos ( = 60, T / = 6 cm 2 /g, M sub = 10 10.5 M ) on circular orbits. By simulating circular orbits, we remove the periodic behavior of the host-subhalo evaporation and tidal effects, such that the comparison between different physical processes is clearer. A more complete and realistic study of subhalos with varying parameters will be shown in the next two sections. Here we also compare the cases with only tidal field (evaporation off) against those with full evaporation.
In the top panel of Fig. 8, we show the central density cen50 of different subhalos as a function of time, where we use dashed lines to denote simulations with evaporation turned off, and solid lines for cases with full evaporation. The three colors mark subhalos on three different circular orbits. It can be seen that the host-subhalo evaporation always delays the subhalo core-collapse, with an increasing impact as the radius of the circular orbit cir decreases. For our orbit with the largest radius of cir = 0.7 200m , evaporation has only a minimal effect; while at the orbit with the smallest radius of cir = 0.1 200m , evaporation is so strong that the whole subhalo is gradually evaporated and never core-collapses.
We further explore the interplay between the two orbit-related effects and SIDM-induced heat flow in the lower panels of Fig. 8.
Here, we present supplementary information about how relevant physical processes affect subhalo evolution, from the top to bottom: the orbital distance between the host and subhalo, the host-subhalo evaporation heating rate (Eq. D2), the average internal cooling rate within the surface of maximum luminosity of the subhalo (luminosity as in Eq. 10), and the tidal heating rate (Eq. 11). All of these heating/cooling rates are in units of energy per unit mass per unit time, as indicated in the luminosity panel.
The evaporation heating rate is the product of the event rate of a subhalo particle scattering off a (virtual) host halo particle and the average energy transferred to the subhalo during one such scattering event. Note that the energy gain of the subhalo per event may not equal the energy transferred to the evaporated subhalo particle. In fact, it is very possible that the after-evaporation velocity of a subhalo particle exceeds the escape velocity, leading to its unbinding and mass loss of the subhalo. Thus the heat transferred during the host-subhalo particle scattering is partially carried away by the escaping subhalo particle. A more careful analysis of the net effect of evaporation on the subhalo should involve at least the mass loss, the adiabatic expansion induced by the mass loss, heating that does not lead to unbinding, and re-scatterings between escaping particles and the rest of the subhalo, which is complicated and beyond the scope of this work. We propose two limits for the evaporation heating rate: the scenario where all evaporated subhalo particles remain bounded to the subhalo, as the maximal heating limit; and the scenario where all evaporated subhalo particles are instantaneously expelled from the subhalo, as the minimal heating limit (complete thermalization limit vs. instantaneous expulsion limit; see Appendix D for a detailed discussion). Because of the large difference between the typical velocity of the subhalo particles and that of the host particles, most of the host-subhalo evaporation scatterings should lead to expulsion of the subhalo particles, thus here we choose the minimal heating case to evaluate the evaporation heating rate. The event rate (probability per unit time) is given by ℎ / of Eq. (1), and the heat transferred to the subhalo per such  Evolution of subhalos on different circular orbits, from top to bottom: the central density of the subhalo cen50 , the host-subhalo separation, the evaporative heating rate in the subhalo, the heat outflow at the subhalo velocity core boundary as a cooling term, and the average tidal heating rate within the subhalo core. Circular orbits with three radii are chosen: brown lines for cir = 0.7 200m , magenta lines for cir = 0.3 200m , and blue lines for cir = 0.1 200m ; two semi-analytical models for SIDM subhalo evolution are shown, dashed lines for tidal-field only and solid lines for cases with full evaporation. The corresponding core-collapse time for a same halo, but in isolation without a host potential, is shown as the brown circle in the top panel. Note that the bottom three panels of heating/cooling terms have the same units. event per unit mass equals the change in the binding energy within the subhalo core due to the mass loss, < − − 2 /2 >, where and are the gravitational potential and velocity of subhalo particles relative to the subhalo center, averaged within the subhalo core radius.
We follow the formalism of Essig et al. 2019 to trace the internal heat flow of the subhalo, where the luminosity ( ) as a function of the temperature gradient is the key factor. The luminosity is given by where lmfp ≈ 0.27 3 T / and smfp ≈ 2.1 / T are the heat conductivity for the long-mean-free-path and short-meanfree-path regimes, is the one dimensional velocity dispersion, = / is the number density of dark matter, B is the Boltzmann constant, and = 2 / B is the temperature, assuming dark matter follows monatomic thermal dynamics. ( ) is the total heat outflow at the surface with radius per unit time, thus we take ( )/ (< ) to be the average cooling rate on all dark matter within the surface . We trace the time evolution of ( )/ (< ) evaluated at the radius of maximal , which generally coincides with the size of the velocity core of the subhalo (not the density core, see Fig. 6 of Nishikawa et al. 2020 for an insightful set of plots illustrating the difference), as the strength of the cooling of the subhalo core.
For the tidal heating rate, we follow the semi-analytical evaluation of Pullen et al. 2014 (their Eq. 15), which accounts for the breakdown of the impulse approximation of the gravitational encounter when the orbit and encounter time scales are comparable, and also includes corrections for higher order heating effects: where shock = orb / orb and orb are the gravitational shock and orbit timescales, with orb and orb being the orbital distance and velocity of the subhalo, and is the radius of the subhalo at which the differential tidal heating rate is calculated. We consider two tensors: , is the tidal tensor (Eq. 10 of Pullen et al. 2014); and , is the time integral of , : with = orb being the orbital distance between the subhalo and the host center. We can see that the strength of tidal heating grows as 2 , meaning that although it heats up the inner part of the subhalo and delays the central density growth, it heats up the outer layers even more. The stronger heating at large subhalo radius could either directly hinder the formation of the negative temperature gradient (larger cooling by ) or help with it by boosting stripping efficiency, at different phases of an orbital period, as we have observed in Sec. 4.2. But this complex, mixed role of the tidal field on subhalo corecollapse via effects on subhalo outliers has been absorbed in the luminosity term we show above. Thus, in Fig. 8 and similar figures below, we evaluate the average tidal heating rate within the radius of max , the same radius as the cooling term / , to highlight the direct heating term of the tidal field within the subhalo central core: With all the three heating/cooling terms explained and evaluated at the same radius of the subhalo, the velocity core radius, we are able to compare them and their effects on subhalo evolution, as shown in the three panels of Fig. 8. The evaporation rates are flat on circular orbits (except for the blue subhalo, whose average binding energy within the core decreases), because they are only affected by SIDM cross section, the orbital parameters, the host properties and the average binding energy within the subhalo core to the first order. The heat outflow rates at the core radius are flat most of the time as well, until the late time exponential growth, which corresponds to the steep temperature gradient when the subhalo is core-collapsing. The / panel suffers from overall greater fluctuations than the evaporation panel, mainly because is the first-order derivative of dark matter 1-D velocity dispersion, which is already a noisy term itself, limited by the simulation resolution and also the resolution of radial binning when we measure the subhalo velocity profile. The tidal heating rates are also relatively flat with time because of the circular orbit, with a slight trend of decreasing near the core-collapse time, as max where we measure the tidal heating rate shrinks. The wiggles in the tidal heating panel are similarly due to the finite resolution in measuring , since the tidal heating rate is measured at the same radius of maximal . Tidal heating is also less significant than the other two terms for the parameters of subhalos that we are considering, except for the case with the smallest orbital radius cir = 0.1 200 (blue color), where although the data we present yield to noise, we can see that the tidal heating is comparable to the other two heating/cooling terms at order of magnitude level.
These three subhalos on circular orbits serve as representative examples to analyze the interactions of different physical processes and the subhalo central density. For the subhalo with the largest orbital radius vir = 0.7 200 (brown color), both the evaporation and tidal effects are relatively weak on it-outweighed by the subhalo internal cooling term / by one and two orders of magnitude respectively, and thus the evolution of its central density is dominated by its internal heat flow. As a result, the core-collapse time of the subhalo in brown color, either with full evaporation or only with the tidal field, is close to the of the same halo in isolation (the small circle in the cen50 panel of Fig. 8). For the subhalo with the smallest orbital radius vir = 0.1 200 (blue color), the total heating outweighs cooling by about one order of magnitude from the beginning. This leads to the subhalo in blue color having a monotonically decreasing central density, and it is completely disrupted at around 4 Gyr. The subhalo with the intermediate orbital radius vir = 0.3 200 (magenta color) is an interesting, marginal case. The cooling term / and the evaporation heating are close in magnitude, with tidal heating being insignificant. Therefore, we see that the central density of the subhalo in magenta color grows slowly for a long time after its core-formation phase, and the subhalo core-collapses long after the of its counterpart halo in isolation. This is the result of the close competition between the total heating and cooling processes.
Furthermore, we find additional evidence of the mixed effects of tidal heating and tidal stripping. In the 1D model of Nishikawa et al. 2020, in the absence of evaporation, core-collapse is increasingly accelerated for decreasing tidal radius. In the case of circular orbits, we expect the smallest orbits to have the largest amount of tidal stripping, and the smallest tidal radii. In the dashed lines in Fig. 8, which show the central density evolution without evaporation, core-collapse proceeds nearly identically for our two innermost orbits, despite the fact that tidal stripping is most severe for the innermost (blue) orbit. The tidal heating rate, however, is also significantly stronger for this orbit than our middle (magenta) orbit. This has potentially important implications for the close-to-host-center subhalos that are especially relevant to substructure lensing.
In this section, our primary conclusion is that there is a strong link between the evolution of the subhalo central density and the heating/cooling terms, that the central density grows (and finally undergoes core-collapse) when the cooling term outweighs the total heating terms inside the velocity core, and grows faster when the net cooling is stronger. We discuss more general cases with eccentric orbits in the next sub-section, where periodical behavior of relevant heating/cooling terms are added back.

Full simulation suite including evaporation
In the code validation runs (Sec. 3.1), we deliberately selected a parameter set [ sub = 10 10.5 , T / = 6 cm 2 /g, peri : apo = 1/10, = 80] such that this particular subhalo is right on the edge of core-collapse. In other words, the unphysically high = 80 is close to the minimal (critical) concentration needed for subhalo corecollapse with other parameters specifying the subhalo and its orbit held fixed, as we shall see later. In this section we vary these four parameters around the critical concentration to explore the transition cases. We find that whether a subhalo core-collapses or not is highly sensitive to these parameters. Similar to the previous section, to fully understand the physics behind subhalo core-collapse, we track the evolution of the subhalo central density cen50 , together with a breakdown of the strength of relevant physical processes: hostsubhalo evaporation and tidal heating as heating terms, and heat outflow at the velocity core as a cooling term. In Fig. 9 we show the cases where and T / are varied, in Fig. 10 and sub , and in Fig. 11 and peri : apo .
In Fig. 9, we present four subhalo simulations: SIDM cross sections T / = 6 cm 2 /g (orange) or 10 cm 2 /g (green), and initial concentrations = 75 (thin lines) and = 80 (thick lines), with the same fixed subhalo mass 10 10.5 (i.e. subhalo-host mass ratio 1:1000) and peri : apo = 1 : 10 (apocenter fixed at 0.7 200 ≈ 700 kpc, unless otherwise specified). We choose = 75 or 80 because these subhalos are right on the edge of undergoing core-collapse within 10 Gyrs. For example, in the T = 10 cm 2 /g with = 80 case, we clearly see the stages of the subhalo coreformation first, where cen50 drops by one order of magnitude, then core-collapses, during which the central density increases from the lowest-density core by three orders of magnitude in a short time. However, for the subhalo with the same SIDM cross section T = 10 cm 2 /g but slightly smaller initial concentration = 75, the core-collapse process happens more slowly and less violently, and is halted before the first pericenter. This hints that whether a subhalo will finally core-collapse or not on cosmological timescales may have a sharp transition in parameter space. A similar disruption of the core-collapse is observed for the cases of the smaller cross section T / = 6 cm 2 /g (orange) near the first pericenter, but compared to the higher cross section cases, the subhalo survives the strongest evaporation and its core-collapse process is resumed at late times for = 80. For the subhalo with T = 6 cm 2 /g and = 75, the core-collapse process is further disrupted again at the second pericenter, where its central density has already grown by a factor of 10, and similarly resumed later, reaching our core-collapse criterion at the very end of 10 Gyr. Thus we record = 75 as the critical concentration for core-collapse for T / = 6 cm 2 /g and = 80 for T / = 10 cm 2 /g, when the mass ratio is 1:1000 and the orbit is characterized by peri : apo = 1 : 10. This seemingly counter-intuitive behavior that the subhalo with a lower cross section could have an easier core-collapse, contrary to what we have seen  The self-interaction cross section is T / = 6 cm 2 /g (orange) or 10 cm 2 /g (green), with chosen to be 75 (thinner line) or 80 (thicker line), near the critical concentration for subhalo core-collapse. The corresponding intrinsic core-collapse time scales c calculated by Eq. (9), the expected core-collapse time when these halos are in the field rather than subhalos, are shown in corresponding small circles in the top panel. Note that we have highlighted the range of ( / ) max of the orange thin line as the light orange shadow in the three lower panels, such that the comparison of the heating/cooling energy scales becomes visually easier. for isolated halos, highlights the power of host-subhalo evaporation in disrupting the core-collapse process of the subhalo.
The expected core-collapse time for corresponding halos in isolation from the host are shown in small circles of Fig. 9, calculated from Eq. (9). The subhalo with T = 10 cm 2 /g and = 80 core-collapses at about the same time we expect for the same halo in isolation, because it happens well before the first pericenter encounter, and the host-induced effects have not become significant yet. For other cases though, their core-collapse time scales are much longer than the isolated halos, showing the impact of the host halo on SIDM subhalo evolution.
All three heating/cooling terms in the lower three panels of Fig.  9 are, as expected, well-synchronized with the central density and orbital distance plots. The evaporative heating rate, proportional to T / , anti-correlates with the orbital distance, since both the host density and the bulk velocity of the subhalo increase as the orbital distance shrinks (see Eq. (1)). The internal cooling term on the subhalo core ( / )| max tracks the central density, which soars when the subhalo tries to core-collapse, and is disrupted together with the central density when the evaporation is strong. Tidal heating also follows the corresponding orbit period, but with a clear extra drop when the subhalo is core-collapsing, because of the shrinking of max at which the tidal heating rate is evaluated. The net cooling/heating of these three terms again determines the growth/decay of the central density cen50 , consistent with our primary conclusion from the previous section with circular orbits, with the evaporative heating and internal cooling being more dominant than tidal heating.
Apart from the heating/cooling perspective, although not as accurate and complete, it is also straightforward to understand the core-collapse of an SIDM subhalo in terms of the competition of time scales. When the subhalo intrinsic core-collapse time c , i.e. the isolated case fit by Eq. (9), is smaller than half the orbit period, the subhalo may eventually core-collapse, but when its c is longer than orb /2, its core-collapse process may either be permanently disrupted (the green thin line of T / = 10 cm 2 /g in Fig. 9), or temporarily halted near the pericenter and resumed at late times near the apocenter (the orange thick and thin lines of T / = 6 cm 2 /g), depending on the strength of evaporation. The latter scenario predicts a large diversity of subhalo central density.
We show similar plots in Fig. 10, varying the subhalo mass sub instead of T / . The group of = 10 8.5 with the same = 75 or 80 as before is not close to core-collapse, compared to the group of = 10 10.5 (blue lines vs orange), due to the long c in isolation (see blue circles vs orange) for low-mass halos. The long intrinsic core-collapse time scale c of a less massive halo is due to the shallower potential well. This yields a smaller velocity dispersion and lower DM self-interaction rate within it (see Eq. (9) and discussions in Sec. 4.1). Another way to understand it is through the / panel of Fig. 10, in which we clearly see that the heat outflow of the 10 8.5 subhalo is weaker than that of the 10 10.5 subhalo by about two orders of magnitude. This is because in the long-mean-free-path regime for cored SIDM halos, the luminosity scales approximately with 5 (see Sec. 4.3). The evaporation heating rate is also lower for smaller subhalos of 10 8.5 , because the average binding energy is lower, but only by about one order of magnitude. As a result, unlike the case with the 10 10.5 subhalo, the cooling by / for the 10 8.5 subhalo is visibly smaller than the total heating for a big chunk of the orbit. Therefore, the critical for less massive subhalos to have core-collapse at fixed cross section will be larger, as we will show in the next section.
Similar results with varying orbits are presented in Fig. 11, where, compared to the fiducial group, we fix the apocenter of the subhalo orbit and alter peri : apo between the default 1 : 10 and the variant 1 : 5. We can see that on this less plunging orbit, both the evaporation and tidal effects are weakened, with evaporation still being the dominant heating term. The net disruption on the subhalo central density near the pericenter is weaker for the less plunging orbit (see the sharp drop of cen50 in the top panel, purple vs orange). As a result, the critical needed for core-collapse on the 1 : 5 orbit becomes 60, while its counterpart is 80 for the 1 : 10 orbit, for a fixed cross section of T / = 6 cm 2 /g.
In this section we see that, as we concluded in the previous section, the core-collapse behavior of an SIDM subhalo on realistic orbits is also a result of the competition between the three relevant physical processes: internal heat transport, host-subhalo evaporation and tidal effects. Because of this, the evolution of a subhalo is very sensitive to its initial parameters, and shows a large diversity of subhalo central density, especially when near the critical parameter sets of ending up in core-collapse. In the next sub-section we map out a boundary in the parameter space for subhalo core-collapse. to core-collapse (defined by central density growing by two orders of magnitude, see Sec. 2) within 10 Gyrs. Different orbits peri : apo = 1 : 2, 1 : 5, 1 : 10 are shown in green, blue and red lines, while isolated halos are in black. The cases with evaporation switched on or off are shown in solid and dashed lines. b) Same as a), but for sub = 10 8.5 .

Mapping parameter space for subhalo core-collapse
From the previous sections, we find that the critical initial concentration for driving subhalo core-collapse is sensitive to other parameters: sub , T / and peri : apo (or, more generally, to the orbit). In this section, we map out the core-collapse criterion for isolated halos, subhalos in the absence of SIDM evaporation, and subhalos with the full evaporation model in the parameter space of [ , sub , T / , peri : apo ].
Properties of the host, the subhalo's apocenter and the subhalo's evolution prior to infall matter as well, but we defer a more careful analysis of those effects to future work. The host halo mass host is relevant because it determines the velocity dispersion of host particles, which affects the evaporation rate, and mediates tidal heating and stripping. However, in this study, we fix host = 10 13.5 as a representative system for substructure lensing (Birrer et al. 2017), which is promising for constraining the 18 Z. Zeng et al. likelihood of finding core-collapsed subhalos in observations. The apocenter distance apo controls the orbital time and thus is another relevant parameter, is fixed at apo = 0.7 200 ∼ 700 kpc, following the subhalo population model of Peñarrubia et al. 2010. The evolution of the subhalo with dark matter self-interaction prior to its inall, which we denote as the 'pre-evolution', determines the actual initial state (e.g. NFW, cored, or early stage of core-collapse) of the subhalo at infall. In this work we effectively set the pre-evolution time pre to be 0 Gyr, since we include subhalos with extremely large cross sections or initial concentrations, of which the core-collapse time scales in isolation are small. We include a simple demonstration of the effect of different pre in Appendix E. host , apo and pre should as well be realistically sampled in future works of population modelling of SIDM subhalos, but in this work we fix them for the purpose of reducing computational burden.
The mass resolution in these simulations is chosen to be p = 10 5 for sub = 10 10.5 and p = 10 3 for sub = 10 8.5 , ten times better resolution than the one used in the validation simulations, with robustness of the resolution tested in Sec. 3.1. This higher resolution is meant to reduce the uncertainty introduced by different realizations of the subhalo. With about 10 5.5 simulation particles, the uncertainty in (sub)halo central density (< 1 kpc for a 10 10.5 halo, most relevant for core-collapse) of different realizations is smaller than the difference in central densities of (sub)halos with 10% difference in initial concentration. Thus we claim a 10% error bar on the critical plots in Fig. 12 (not shown explicitly for plot readability), and believe that this estimate is conservative. Of these simulations, the non-core-collapsed subhalos cost a few tens of CPU hours while the core-collapsed ones cost up to a thousand or more, due to the extreme refinement of the simulation timestep (four to five orders of magnitude relative to cored subhalos) on account of the high central density.
In Fig. 12 we show the critical needed for subhalos to have core-collapsed in 10 Gyr, against the self-interaction cross section T / , for subhalo mass of 10 10.5 in Fig. 12a and 10 8.5 in Fig. 12b. Subhalos on different orbits peri : apo = 1 : 2, 1 : 5, 1 : 10 are shown in green, blue and red colors, while same halos in isolation are in black lines. Cases with or without our evaporation model are plotted in solid or dashed lines.
As we can see in both figures of Fig. 12, for isolated halos (black), critical monotonically decreases as sub and T / increase, as expected from Sec. 4.1. For subhalos with only tidal effects (dashed lines), we find that their critical concentrations are always smaller or equal to those of the isolated ones, with a decreasing trend of critical as orbit gets more plunging. This is due to the tidal acceleration of core-collapse. The 'tidal deceleration' we find in Sec. 4.2 is not reflected in this comparison between the isolated halos and tidal-only subhalos, because it happens in specific cases with c orb /2, while for the boundary cases with c ∼ 10 Gyr, all subhalos have experienced at least a full orbital period. For the same reason as isolated halos, the critical of subhalos in a tidal field monotonically decreases with SIDM cross section and subhalo mass.
The cases get more complicated with full evaporation. On top of the analysis of tidal-only subhalos above, the strength of evaporation increases with SIDM cross section and peri : apo . Thus larger critical is required for more plunging orbits, with an only exception at the high cross section end / = 60 cm 2 /g for the less massive subhalos sub = 10 8.5 , which is because the initial concentration is high enough to counter the evaporation before pericenters, but tidal effects have come into play as an acceleration factor to core-collapse. But critical does not monotonically scale with T / as in the cases of tidal-only or isolated (sub)halos, because both the subhalo intrinsic core-collapse and evaporation are boosted with larger T / . Thus we only numerically map out critical boundaries in the parameter space of subhalo core-collapse with evaporation, which, as we can see in Fig. 12, are more complex than the other two models.
The most important take-away message from Fig. 12, however, is that core collapse is not feasible for subhalos in constant cross section SIDM cosmologies. Core-collapse only occurs within the Hubble time for extremely high initial concentrations (compared to the cosmological CDM subhalo concentration of 200 ∼ 24 for 200 = 10 10.5 and 200 ∼ 40 for 200 = 10 8.5 , see Duffy et al. 2008, Diemer & Joyce 2019, Wang et al. 2020. These concentrations are unphysically high for standard inflationary power spectra, although on the small halos' end they may be resulted from non-standard cosmologies such as a distorted primordial power spectrum (Gosenca et al. 2017;Delos et al. 2018b,a) or an early matter-dominated era (Erickcek 2015;Delos et al. 2019;Barenboim et al. 2021). However, without evaporation, core-collapse can occur with typical halo concentrations if the cross section at velocity scales typical of particles within the subhalo is high. This suggests that a strongly velocity-dependent cross section -such that the cross section relevant to evaporation is low even as the cross section relevant to heat transport within the subhalo is high (Turner et al. 2021;Correa 2021;Jiang et al. 2021) -may allow for the development of core-collapsed subhalos in a strong-lens-mass host. Alternatively, turning on another degree of freedom in the cross section, such as allowing energy dissipation in scattering (Essig et al. 2019;Huo et al. 2020), may allow our constraints here to be evaded. We will explore these extra degrees of freedom, especially in the context of popular velocity-dependent models, in future work.

SUMMARY AND DISCUSSION
Small dark matter substructures are promising targets for constraining microphysical models and cosmological phenomenology of dark matter, and can be directly studied with substructure lensing techniques. Self-interacting dark matter, a class of dark matter model that allows for scattering between dark-matter particles, is even more promising to be constrained by substructure lensing, because core-collapsed (sub)halos may form as a result of the thermalization and heat outflow within the (sub)halo (Gilman et al. 2021;Yang & Yu 2021). These core-collapsed subhalos make excellent lenses due to their highly concentrated central densities and convergence, bringing larger distortions to the strongly lensed images of the host halos. However, for SIDM, because of the computational expense, traditional cosmological simulations and staged simulations with live host halos are barely able to resolve small subhalos with mass 10 8 , which is within the sensitivity scope of coming substructure lensing observations (Vegetti et al. 2014;Hezaveh et al. 2016;Nierenberg et al. 2014Nierenberg et al. , 2017Gilman et al. 2018Gilman et al. , 2020aHsueh et al. 2020). This is even more costly for SIDM subhalos that are core-collapsing.
In this work, we present a hybrid semi-analytical + N-body method for staged simulations of SIDM subhalo evolution with the host halo modeled as an analytical static density distribution, where for the first time the host-subhalo evaporation is properly incorporated to ensure full consistency. In principle, we are able to trace the detailed evolution of arbitrarily small SIDM subhalos with this new method. Our method is implemented as a patch to the SIDM module of the hydro-gravity code Arepo (Springel 2010;Vogelsberger et al. 2012Vogelsberger et al. , 2014, and the evolution of subhalos in such analytical hosts is accelerated by orders of magnitude compared to a full live host simulation. We validate the method, according to both the subhalo mass loss history and density/mass profile evolution, with live host runs. The mass-loss history of the subhalos in our method is consistent with the live host runs at the level of ±10% for sub-host mass ratio 1:1000. We show that the lingering discrepancy can be mostly explained by a small amount of dynamical friction that is missing in our analytic host model. Thus our method works best for relatively small subhalos ( 1 : 1000 host halo mass; for larger subhalos, simulation with a live host is the most accurate, and inexpensive as well) and we expect that the accuracy of our method increases for even smaller subhalos ( 10 10 ), which are the main targets of substructure lensing in the near future. The detailed evolution of the subhalo density/mass profiles also show that our method is able to track both cored and core-collapsing subhalos at a remarkable level. We tested constant and velocity-dependent cross section models for code validation. Both are in good agreement with the fully live host simulations. We will explore the physics of core-collapse of velocity-dependent SIDM subhalos in substructure lenses in future work.
To explore the conditions for SIDM subhalos to core-collapse and study the relevant physical processes that drive such central density evolution, we ran a set of simulations scanning the subhalo parameter space of SIDM constant cross section, subhalo mass, peri-center distance and subhalo initial concentration, [ / , sub , peri : apo , ]. We fixed the host mass to a typical strong lens system (Birrer et al. 2017), and subhalo orbit apocenter at a typical infall radius (Peñarrubia et al. 2010). We use the subhalo central density cen50 as the indicator of whether and when the subhalo core-collapses, and find that cen50 grows if and only if there is a net cooling (energy loss) out of the three cooling/heating terms: the heat outflow from the subhalo core as a cooling term, and host-subhalo evaporation and tidal heating as heating terms. The SIDM subhalo central density decreases first during the coreformation era, then increases as a result of the heat outflow from the core (Colín et al. 2002). However, this increase is disrupted when the subhalo comes near its pericenter, where both the evaporation and tidal heating are the strongest. Thus we have observed three possible fates of an SIDM subhalo in our simulations: 1. that it either core-collapses before the first peri-center, 2. is largely evaporated at the first peri-center and can never core-collapse, 3. or is disrupted but not fatally enough such that its central density may resume the increase again after the peri-center passage. In the last case, we may expect multiple cycles of growth and loss in the central density evolution, suggesting a large diversity of subhalo density profiles.
We find that whether a subhalo core-collapses in a Hubble time, defined by its central density growing by two orders of magnitude (comparable to the criterion of Essig et al. 2019), is highly sensitive to the subhalo initial parameters, especially the concentration . Even a small change in may lead to completely different evolution histories. Therefore, we sample a few hundred subhalos from the four-dimensional parameter space of [ / , sub , peri : apo , ] to map the boundaries of critical conditions for subhalos to have corecollapse. We find that when host-subhalo evaporation is included, the critical/minimum is always higher by a factor of a few compared to the cosmological ( 60 for most / and orbits we consider vs. ∼ 24 for sub = 10 10.5 ), no matter what values we choose for the other parameters. However, for field halos or subhalos without the evaporation, the gap between the critical and the cosmological ones is greatly reduced (∼ 30 vs ∼ 24, for sufficiently high cross sections / 30 cm 2 /g), in agreement with Correa 2021. This suggests that it is nearly impossible to have completely core-collapsed subhalos in the simplest case of SIDM with constant cross sections. More degrees of freedom, such as velocity dependence of T / (Vogelsberger et al. 2012;Zavala et al. 2013;Nadler et al. 2020b;Banerjee et al. 2020;Turner et al. 2021;Correa 2021) or energy dissipation during the two-body scattering (Schutz & Slatyer 2015;Essig et al. 2019;Vogelsberger et al. 2019;Huo et al. 2020), must be introduced to achieve subhalo core-collapse. Subhalos with constant cross sections and only tidal effects (but not evaporation) can be seen as being the limit of an extremely velocity-dependent cross section, in the limit where the relative speed between the subhalo and host is significantly different from the subhalo internal velocity dispersion. Because the relevant heating/cooling processes are sensitive to different velocity scales, the core-collapse rate of (sub)halos can be used to constrain velocity-dependence of SIDM models as well. We plan to include such advanced models of SIDM in future work, which are simple extensions to our current method (as shown in the validation for velocity-dependent SIDM in Sec.

3.2).
Our work suggests additional pathways for future progress to connect substructure lensing observations with predictions: • Particle simulation near the hydrodynamic limit. By sampling dark matter as macro simulation particles, it is straightforward to implement the two-body scattering processes in particle-based simulations such as Arepo, but they are limited by the finite mass/force resolution when the central density increases exponentially, when the cores transition to the short-mean-free-path limit in core-collapse. Our core-collapse criterion (Eq. 3) is comparable to the short meanfree-path criterion of Knudsen number ∼ 1 (i.e. mean free path comparable to gravity scale height) in Essig et al. 2019. We terminate the simulation at this point because the timestep required to follow the scattering accurately becomes too small. Thus we suggest that transitioning to a hydrodynamical simulation model or analytical treatment should be implemented when an SIDM halo is approaching the core-collapse stage.
• Softening length As we have mentioned in Sec. 3.1, although the length scale of an ultra-dense core of a core-collapsing (sub)halo could become smaller than the gravitational softening length, only a small numerical error is introduced near the core-collapse time, since the core-collapse is a fast, runaway process. However, we also note that this small error could possibly have a significant effect on when and whether a subhalo core-collapses for certain rare marginal cases, such as the one with [ sub = 10 10.5 , T / = 6 cm 2 /g, peri : apo = 1/10, = 75] (thin orange line in Fig. 9), where the cooling and heating processes have a close competition right at the second pericenter, and a small offset along the timeaxis could lead to dramatic difference. This could happen when the core-collapse time of the subhalo is well-synchronised with its peri-center arrival, thus should be only on rare occasions. The force resolution for CDM subhalos should be generally good enough for SIDM subhalos.
• Baryonic effects. The inclusion of baryons in simulations can dramatically change predictions of halo properties relative to those from dark-matter-only simulations. When baryonic physics is included in SIDM simulations, the halo central density may be either more cuspy (Kamada et al. 2017;Sameie et al. 2018;Robertson et al. 2019;Sameie et al. 2021), or more cored (Robles et al. 2017), or display a variety of cusps and cores (Creasey et al. 2017;Robertson et al. 2018;Robles et al. 2019;Despali et al. 2019;Burger et al. 2021), potentially depending on the halo mass, environment (Rob-20 Z. Zeng et al. les et al. 2019), and SIDM cross section ). This will affect the effective concentration of both host halos and subhalos. Higher (lower) concentrations of host halos would strengthen (weaken) the tidal effects (further strengthened with a central galaxy potential) and evaporation, and thus disrupt more (less) subhalos before they can core-collapse. In contrast, higher (lower) effective concentrations of subhalos lead to more (less) core-collapse in them (Feng et al. 2021b). The net effect of baryons can thus be complicated and diverse, and systematic studies of hydrodynamic simulations with SIDM are needed to unpack it. For relatively small subhalos within the sensitivity scope of upcoming substructure lensing observations, the baryon fraction will be small enough to be safely ignored, but the baryonic effects on the host remain significant. Since substructure lensing is most sensitive to subhalo lenses concentrated towards the main lens, modelling baryonic effects correctly and consistently becomes even more important.
• A diversity of (sub)halo central densities. Core-collapse of SIDM subhalos has been proposed as a plausible solution to multiple small-scale problems of CDM, such as the cusp-core problem and the too-big-to-fail problem. In this paper we show that fully core-collapsed subhalos are nearly impossible to form for SIDM with constant cross sections, given the cosmological (sub)halo concentrations, but this does not rule out SIDM as a solution to the small-scale problems. With periodical ups and downs of the central density as we have observed (e.g. the green and orange thin lines in Fig. 9), the subhalos could be even more diversified and interesting, even though subhalos that once entered early core-collapse stage but are evaporation-disrupted later are not counted as completely core-collapsed in our critical parameter space. This might be able to explain the mystery in Turner et al. 2021, that two out of nine subhalos in a Milky Way analogue host halo have cored, decreasing central densities, while the other seven show clear sign of entering core-collapsing phase (see also Correa 2021). To fit the diversity of observed subhalos into SIDM, either with constant cross sections or velocity-dependent ones, a proper population modelling of the subhalos with realistic initial parameters is needed. We plan to sample multiple subhalos from realistic merger trees as an extension to our single subhalo evolution model in this work (Yang et al. 2020), where the possible interactions between the subhalos are also naturally included.
This work presents our first step in studying the evolution of SIDM subhalos, including the core-collapsing scenario, in the context of strong lensing systems. In the near future, we can extend our work in multiple directions, such as: a more holistic study of SIDM subhalo population with realistic assembly histories -potentially generating mocks of substructure lensing observations as well; implementing baryonic effects, either with hydrodynamic simulations or with analytical approximations; and further simplifying our method to a purely semi-analytical approach. These are all necessary steps to generate large ensembles of mock observations to constrain the nature of dark matter with future substructure lensing observations. , T / = 6 cm 2 /g, same as the violet data sets in Fig. 2a. agreement with the live host simulation than without the inclusion of this dynamical friction term, but still find an under-decay of the orbit in the semi-analytic framework. Instead, if we set hm = 0, the discrepancy between orbits becomes smaller, with the modelling orbits having an over-decay (see Fig. A1 for a comparison). In neither case can we get dynamical friction perfectly correct and the live host orbit lies between these two cases. Thus the 'true' mass loss history from our method, assuming we track the orbit decay perfectly, should lie between the dashed and dotted lines in Fig. A1 as well. Therefore, we only use the modelled dynamical friction as a demonstration, showing that the discrepancy in subhalo mass loss history is indeed reduced, if the orbits are in better agreement.
We conclude that the missing dynamical friction is a primary systematic in our semi-analytical method of subhalo evolution. Thus, our model works best for relatively small subhalos, with a mass 1/1000 of the host, for which the dynamical friction is small, as we have shown in the main text. And in fact, such small halos are the reason for our semi-analytic model development in the first place.

APPENDIX B: ANIMATION: DENSITY PROFILE EVOLUTION IN CODE VALIDATION
This animation shows the full evolution of density/mass profiles for the subhalo simulation from Sec. 3.1, corresponding to Fig.  3. Adobe Acrobat Reader is required to make the animation run.

APPENDIX C: ANIMATION: TIDAL HEATING
This animation shows the disturbed velocity dispersion profiles of subhalos in a tidal field, when these subhalos approach pericenter and tidal heating comes into effect. Readers are referred to Sec. 4.2, Fig. 6 and Fig. 7 for more details. Adobe Acrobat Reader is required to make the animation run.

APPENDIX D: EVAPORATION HEATING: COMPLETE THERMALIZATION LIMIT VS INSTANTANEOUS EXPULSION LIMIT
The net effect of host-subhalo evaporation is complicated, since after scattering with a host particle, a subhalo particle could have a velocity that is either larger or smaller than the escape velocity of the subhalo (see also Kummer et al. 2018 for a quantitative discussion). When a subhalo particle's after-evaporation velocity is smaller than the escape velocity, the heat gained from this scattering event is completely absorbed by the subhalo. When its velocity is greater than the escape velocity, the major effect on the subhalo would be the adiabatic expansion caused by this instantaneous mass loss, with possible re-scatterings with other subhalo particles on its path-ofescape as a secondary effect, leading to partial thermalization.
Thus we propose two natural limits when considering the heating of evaporation (within the subhalo core): the complete thermalization where all evaporated subhalo particles remain bound to the subhalo, thus all heat transferred to the evaporated particles are absorbed by the subhalo; and the complete instantaneous explusion limit, where all evaporated particles immediately leave the subhalo, with no re-scatterings. The complete thermalization scenario is the upper limit of evaporation heating, while the instantaneous expulsion is the lower limit. Complete thermalization Instantaneous expulsion Figure D1. The evaporation heating rate, in the complete thermalization limit (solid lines) and instantaneous expulsion limit (dashed lines). The subhalo in orange color has a parameter set of [ / = 6 cm 2 /g, sub = 10 10.5 , = 75, peri : apo = 1 : 10, apo = 0.7 200 ], same as the thin orange line in Fig. 9; and the subhalo in magenta color has [ / = 6 cm 2 /g, sub = 10 10.5 , = 60, peri : apo = 1 : 1, apo = 0.3 200 ], same as the magenta solid line in Fig. 8.
In the complete thermalization limit, the evaporation heating per unit mass per unit time is given by the product of the event rate (probability per unit time) and the average change in subhalo particles' kinetic energy during one evaporation event: where ℎ / is given by Eq.
(1), and and are the initial and final velocities before and after the host-subhalo particle scattering, relative to the subhalo center. We evaluate the average < 2 − 2 > over the 50 particles that have the maximum local densities, to represent the center particles of the subhalo.
In the instantaneous expulsion limit, the evaporation heating per unit mass per unit time is given by the product of the event rate and the increase in subhalo center's binding energy per unit mass: where is the gravitational potential of a subhalo particle and is its velocity relative to the subhalo center. The explusion leads to a subhalo particle's binding energy increases from ( + 2 /2) to 0. Note that for the expulsion limit we evaluate the averaged binding energy within the subhalo core, unlike in the thermalization limit where we fix the number of particles to take average over. This is because in the thermalization limit, the kinetic energy change is nearly universal subhalo-wide, because the velocity difference of the host-subhalo particle pair is much larger than the velocity dispersion of the subhalo. But in the expulsion limit, a subhalo particle's binding energy is sensitive to its radius from the subhalo center. Thus we must choose the core radius, Lmax , as the boundary within which we evaluate the averaged binding energy.
As we can see from the comparison of these two limits in Fig.  D1, the maximal heating of complete thermalization limit is indeed always larger than the minimal heating in the instantaneous limit. The difference ranges between a few times to slightly larger than one dex. For eccentric orbits, the difference between these two limits peaks at the pericenter. This is because for the maximal heating case, the change in kinetic energy per evaporation event, < 2 − 2 > /2 in Eq. (D1), peaks at the pericenter where the subhalo's orbiting velocity is at maximum; but for the minimal heating case, the change in binding energy, < − − 2 2 > in Eq. (D2), is not much affected by the orbit, or even slightly decreases at the pericenter, since the subhalo's gravitational potential gets shallower with more particles evaporated.
We expect the exact evaporation heating rate to lie between the complete thermalization limit and instantaneous expulsion limit as shown above, but it should be closer to the minimal heating of instantaneous expulsion scenario. Because of the large gap between the typical velocity of the subhalo particles and that of the host halo particles, most of the host-subhalo evaporation scatterings should lead to expulsion of the subhalo particles. Re-scatterings with other fellow subhalo particles on the escape-path of these evaporated ones should be relatively rare, as the outskirts of the subhalo are always in the long-mean-free-path regime of SIDM.

APPENDIX E: PRE-EVOLUTION OF SUBHALOS
In this work, all subhalos are initialized with NFW profiles according to the input concentrations and subhalo masses at the infall time. However, in reality, subhalos would have evolved with dark matter self-interaction for some time prior to infall, which we denote as the 'pre-evolution'. This pre-evolution time is another relevant free parameter that sets the 'realistic' initial condition of the subhalo at infall, but due to the computational cost we are not able to do a thorough analysis over another degree of freedom in this work. Below we show a simple demonstration of the effect that different pre-evolution time have on the subhalo.
In Fig. E1 we show the central density as a function of time for the two subhalos with [ sub = 10 10.5 , / = 6 cm 2 /g, c = 75, r peri : r apo = 1 : 10] and [ sub = 10 10.5 , / = 10 cm 2 /g, c = 75, r peri : r apo = 1 : 10] (the orange thin line and green thin line in Fig. 9), with pre-evolution time of 0 Gyr, 0.5 Gyr and 1 Gyr separately. We can see that for the / = 10 cm 2 /g subhalo, pre-evolving the subhalo for 0.5 or 1 Gyr is very different from not pre-evolving it -the latter is completely disrupted while core-collapse of the cases with pre-evolution is not quite delayed. For the / = 6 cm 2 /g subhalo, pre-evolving it for 0.5 or 1 Gyr still sees delay in the core-collapse time, but not as much as the case without pre-evolution (∼1 orbital period vs. ∼2 orbital periods).
With realistic pre-evolution time, the critical concentration for subhalo core-collapse could be reduced, especially for SIDM with large cross sections, since the intrinsic core-collapse time for halos in isolation is short. our next paper (ongoing research) will be focused on simulating an ensemble/population of subhalos, with realistic assembly histories into the host. Thus the pre-evolution problem will be naturally dealt with in our next work. This paper has been typeset from a T E X/L A T E X file prepared by the author.