Constraining Dissipative Dark Matter Self-Interactions

We study the gravothermal evolution of dark matter halos in the presence of dissipative dark matter self-interactions. Dissipative interactions are present in many particle-physics realizations of the dark-sector paradigm and can significantly accelerate the gravothermal collapse of halos compared to purely elastic dark matter self-interactions. This is the case even when the dissipative interaction timescale is longer than the free-fall time of the halo. Using a semi-analytical fluid model calibrated with isolated and cosmological $N$-body simulations, we calculate the evolution of the halo properties -- including its density profile and velocity dispersion profile -- as well as the core-collapse time as a function of the particle model parameters that describe the interactions. A key property is that the inner density profile at late times becomes cuspy again. Using 21 dwarf galaxies that exhibit a core-like dark matter density profile, we derive constraints on the strength of the dissipative interactions and the energy loss per collision.

I. Introduction. The elusive nature of dark matter in terrestrial experiments combined with hints for nontrivial dynamics from astrophysical systems has led to the dark sector paradigm: the dark matter may be connected to a plethora of hidden particles with their own interactions, see [1][2][3] for overviews. These dark-sector interactions may modify the formation and evolution of dark matter halos and alter their inner structure. Astrophysical observations can in turn provide important tests on the microscopic physics in the dark sector.
In this Letter, we explore observational consequences of a generic dark-sector model, where dark matter particles have both elastic and dissipative self-interactions. Self-interacting dark matter (SIDM) has been motived to solve long-standing issues of the prevailing cold dark matter model on galactic scales, see [4] for a review. Most SIDM studies focus on the elastic scattering limit. However, in many particle physics realizations of SIDM [5][6][7][8][9][10][11][12][13][14], dark matter particles also have dissipative collisions. We show that observations of constant dark matter density cores in many dwarf galaxies can be used to test dissipative dark matter self-interactions.
We focus on the "mild cooling regime", in which the cooling timescale is longer than the free-fall time of the halo. In this case, the halo mostly stays in hydrostatic equilibrium and contracts as a whole without fragmentation, as opposed to situations with strong cooling [36][37][38][39].
After introducing a physical model to capture the bulk cooling, we perform numerical simulations to trace the evolution of the halo and calibrate the results against both isolated and cosmological N -body simulations. Finally, we derive strong limits on the strength of dissipative interactions in the dark sector. In the Supplemental Material, we provide additional details and results to further support our main text.
II. Methodology. To understand halo evolution in the presence of dissipative interactions, we employ a semi-analytical fluid model, which has been used to study globular clusters [40][41][42] and halos consisting of SIDM without dissipation [43][44][45][46][47]. Since this method is computationally inexpensive, we are able to scan a wide range of parameter space. Moreover, it can resolve the very inner regions of the simulated halo.
For an isolated halo, we assume spherical symmetry and use the following set of transport equations to describe the gravothermal evolution in the radial direction where M (r, t) is the fluid mass enclosed within radius r at a time t, ρ(r, t) is the local density, ν(r, t) is the one-dimensional velocity dispersion, L(r, t) is the luminosity, C(r, t) is the volumetric bulk cooling rate, G is the gravitational constant, and (∂ t ) M denotes the Lagrangian time derivative. The temperature is related to The evolution of the density profiles for an SIDM halo, LSB F583-1 (see text for halo parameters), assuming purely elastic dark matter self-interactions ("no cooling", blue) and self-interactions with an additional dissipative interaction ("with cooling", red) as indicated. Numbers show the Knudsen number for the innermost shell at a given time. We project the evolution of the density of the innermost shell on to the ρ-t plane and mark the separation of stages 1 → 2 and 2 → 3 with and , respectively. We set β = 0.60. Right: The evolution of the 1D velocity dispersion profile with respect to the enclosed mass. The setup and labels are the same as the left panel.
ν as mν 2 = k B T , where k B is the Boltzmann constant. We assume the dark matter particle is monatomic and set the adiabatic index γ = 5/3. The elastic and dissipative interactions are encoded in the conduction ∂L/∂r and the cooling term C, respectively. In this work, we assume both the elastic and inelastic cross sections are independent of the dark matter velocity. Dark matter elastic self-scattering allows radial heat conduction. This can be characterized by comparing the mean free path λ = 1/nσ, where n is the local number density and σ is the cross section, to the scale height H = ν 2 /4πGρ. The ratio of λ to H is the Knudsen number, Kn ≡ λ/H, which indicates the importance of heat conduction induced by elastic scattering. We refer to regions with Kn > 1 (Kn < 1) as long-mean-free-path (short-mean-free-path) regions. Note that Kn ≈ t r /t dy , where t r ≈ λ/ν is the local relaxation time for the elastic scattering and t dy = H/ν is the dynamical time of the halo. The luminosity L is a function of the temperature gradient where the conductivity κ = (κ −1 lmfp + κ −1 smfp ) −1 reduces to the conductivity of the long-mean-free-path (κ lmfp ) and short-mean-free-path (κ smfp ) regions in the appropriate limits, i.e., κ lmfp = (3β/2)nH 2 k B /t r 0.27βnν 3 σk B /(Gm), and κ smfp = (75π/256)nλ 2 k B /t r 2.1νk B /σ [44][45][46][47]. In κ lmfp , the numerical factor β = 0.60 (0.45) is set by calibrating the fluid model to isolated (cosmological) N -body simulations; see Supplemental Material for details.
Since we assume the energy released during the dissipative collision is not reabsorbed by dark matter particles in the halo, the cooling rate C appears as a bulk term in Eq. (1), which can be written as a function of the model parameters, where ν loss ≡ E loss /m is the "velocity loss" that parameterizes the energy loss per collision; t r and σ are the relaxation time and the cross section of the dissipative interaction, respectively; and we restrict scattering to particles whose kinetic energy exceeds E loss . This model of cooling captures the essential features of dissipative interactions and can be applied to various cases. For instance, in an atomic dark matter model [5,7,11], dark atoms could release dark radiation through collisional excitations followed by spontaneous emission of a massless dark photon, which escapes the halo without re-absorption. We solve Eq. (1) with the boundary conditions at t = 0 of M = L = 0 for the inner boundary and L = 0 for the outer boundary. We assume the initial halo mass distribution follows an NFW profile [48], ρ(r) = ρ s r 3 s /r(r + r s ) 2 , where ρ s and r s are the scale density and radius, respectively. In our simulations, we reformulate Eq. (1) in terms of a set of dimensionless variables based on r s and ρ s [45,47]. Our numerical procedure is similar to [45,47]. We assume the halo is always in hydrostatic equilibrium and divide the halo into radial Lagrangian zones that have a fixed mass. At each time Dimensionless collapse timetc ≡ tc(4πGρs) 1/2 as a function ofσ when cooling is absent. We show the isolated N -body simulation result from [46] and a recast of the fluid model of [47] with β = 0.60. The gray line shows tc with β = 0.45, calibrated from cosmological N -body simulations [26]. In setting the constraints, we focus onσ ≤ 10 −1 (solid). Right: Ratio of the collapse time, ξ ≡ t c /tc, as a function ofν loss for different values of σ /σ. Note that forσ = 10 −4 − 10 −1 , we find that ξ has a very weakσ dependence.
step, we allow heat conduction and cooling across these zones, followed by hydrostatic relaxation. III. Gravothermal Evolution. To illustrate the effect of the dissipative interactions, we consider a dwarf halo with mass 8 × 10 10 M and characteristic halo parameters r s = 6.5 kpc and ρ s = 1.28 × 10 7 M /kpc 3 . This is motivated by the SIDM fit to the galactic rotation curve of LSB F583-1 [32], which exhibits a cored density profile. We take σ/m = 3 cm 2 /g as in [32] and consider σ = 0 as well as σ = σ and ν loss = 13 km/s. Fig. 1 shows the density vs radius (left) and the velocity dispersion vs enclosed mass (right) over time with (red) and without (blue) bulk cooling. Each curve in Fig. 1 is labeled with a Kn value for the innermost simulated shell. From the density evolution, we see the process can be divided into three stages: (1) Core expansion. Heat conduction is inwards (L < 0) and Kn 1. The halo evolves quickly to a quasi-isothermal state. (2) Selfsimilar collapse. Heat is conducted outwards (L > 0) and Kn slowly decreases. The self-similar collapse results in a cuspy density profile and with log-slope of approximately −2, a characteristic feature if the cooling is absent or mild. (3) Post-self-similar collapse. Here Kn < 1 at the center and the inner density suddenly begins to increase dramatically. In Fig. 1, the symbol denotes the 1 → 2 transition, when the innermost shell is at its least dense and its luminosity vanishes; the symbol denotes the 2 → 3 transition, when Kn = 1. From the velocitydispersion evolution, we observe that the cold inner halo (r < r s ) quickly thermalizes with the maximum velocitydispersion of the initial NFW profile, which is about ν ∼ 0.3 (4πGρ s r 2 s ) 1/2 (orν ≡ ν/(4πGρ s r 2 s ) 1/2 ∼ 0.3 in normalized quantities). The inner halo maintains this velocity dispersion during much of stage 2, before increasing sharply at very late times.
For concreteness, we define a collapse time as the time at which Kn = 0.1 for the innermost shell, and we denote the collapse time with (without) inelastic cooling as t c (t c ). Since the evolution of the third stage is very fast, t c and t c are largely determined by the time of the first two stages. The most important effect of the dissipative interaction is to significantly speed up the collapse time, t c < t c . For LSB F583-1 with the model parameters chosen in Fig. 1, the collapse time with cooling is shortened by about a factor of 20, resulting in t c ≈ 8.5 Gyr. This amount of cooling is disfavored, because the final density profile is too steep to be consistent with the observed profile of LSB F583-1 [32].
We perform a suite of simulations, varying the model parameters within the following range of values in dimensionless units:σ ≡ (σ/m)ρ s r s = 10 −4 -10 3 , σ /σ = 10 −3 -1, andν loss ≡ ν loss /(4πGρ s r 2 s ) 1/2 = 0-5 with evenly log-spaced steps. In Fig. 2 (left), we show results for the halo evolution with pure elastic self-scattering and no cooling. Forσ 1, there is a simple scaling relation betweent c ≡ (4πGρ s ) 1/2 t c andσ, namelŷ t c ≈ (150/β)σ, which can be expressed as In this regime, a largeσ speeds up the thermal evolution of the halo and shortens the collapse timescale. However, asσ 1, the inverse proportionality is lost because the mean free path is too short and heat conduction is actually suppressed [43,49]. Below, when setting constraints on dissipative dark matter, we restrict toσ ≤ 0.1, along with σ ≤ σ, so that the initial halo is in the optically-thin regime. For the parameters shown in Fig. 1,σ = 0.1 corresponds to σ/m = 0.1/r s ρ s = 5.8 cm 2 /g, so the choice of σ/m = 3 cm 2 /g satisfies the condition.
Our results are in good agreement with [47], wherê σ = 0.088 and β = 0.75 were chosen. To compare with cosmological N -body simulations of dwarf halos in [26], we take the Pippin halo parameters, r s = 2.7 kpc and ρ s = 1.73 × 10 7 M /kpc 3 , and apply Eq. (4). The estimated core-collapse time is t c ≈ 80 Gyr for σ/m = 10 cm 2 /g and t c ≈ 16 Gyr for σ/m = 50 cm 2 /g for β = 0.60, consistent with the absence of core collapse and the presence of a mild collapse, respectively, observed in the simulations. It is remarkable that although our method does not take into account mergers in halo formation, it does reliably capture the overall features of its thermal evolution. A careful calibration directly to the cosmological N -body simulations yields even better agreement for β = 0.45. Fig. 2 (right) shows the reduction of the collapse time, ξ ≡ t c /t c , from dissipative interactions. Althoughσ sets the overall timescale of t c , we find that ξ is sensitive to σ /σ andν loss , but notσ itself. Overall, the maximal reduction is achieved whenν loss ≈ 0.3 for a wide range of σ /σ. The origin of this scale is clear from Fig. 1 and the surrounding discussion, where we showed thatν of the inner halo is near 0.3 for most of the halo's evolution. Forν loss 0.3, the energy loss is small per collision, while forν loss 0.3, inelastic scattering can only occur among particles on the high-velocity tail or very late in the halo evolution (stage 3). The collapse time in the presence of cooling is then t c = ξ(σ /σ,ν loss )t c , where t c is given by Eq. (4) and ξ can be read from Fig. 2 (right). Forν loss < 0.2, we find an approximate formula ξ ≈ exp −ν loss (σ /σ) 1/2 /0.027(4πGρ s r 2 s ) 1/2 . We see that the collapse time can be reduced as much as a factor of 10 3 , indicating that dissipative scattering can be important for the evolution of the SIDM halo.
IV. Astrophysical Implications. The dark matter distribution in a halo can be inferred from stellar kinematics. For many dwarf and LSB galaxies, dark matter dominates the dynamics, and the stars and gas particles trace the gravitational potential well of the halo. Ref. [32] analyzed the rotation curve data of 30 spiral galaxies and found that they can be fitted with an SIDM model with elastic self-scattering cross section σ/m = 3 cm 2 /g. In fact, 21 of the galaxies have a constant density core and none of them show evidence of core collapse. Thus, we use this sample to constrain the dissipation parameters σ and ν loss by demanding that the collapse timescale, t c , is longer than 10 Gyr, the overall age of galaxies. Fig. 3 shows regions (shaded) where core collapse occurs in less than 10 Gyr for individual galaxies, taking halo parameters r s and ρ s from [32] as input. In these regions, the inner density profiles of the associated halos at 10 Gyr are much steeper than inferred from the stellar kinematics [32]; see Supplemental Material for more details on determining the exclusion limits. In solid (dashed) purple, we show the boundary from all galaxies imposing the constraint t c < 10 Gyr with calibration : Disfavored parameter space (σ /σ, ν loss ) from the absence of core collapse in dwarf and low surface brightness galaxies assuming σ/m = 3 cm 2 /g. We take the fitted dark matter halo parameters of the galaxies with a low-baryon content from [32] (listed in Tab. II and labelled for the outer ones) and require t c ≥ 10 Gyr in solid (dashed) purple for fluid simulations with calibration parameter β = 0.60 (0.45). We highlight LSB F583-1 with t c ≥ 10 Gyr and the four benchmark points discussed in Sec. IV.
It is useful to find scaling relations between the collapse time and other halo parameters, such as mass 200 . Thus, for fixed M 200 , more highly concentrated halos collapse much faster than halos with low concentration, while the dependence of t c on the halo mass is mild. In addition, t c is sensitive to ν loss =ν loss (4πGρ s r 2 We have studied the gravothermal evolution of dark matter halos in the presence of dissipative dark matter self-interactions. After introducing a simple but well-motivated model to capture the cooling effect, we performed numerical simulations based on a semi-analytical approach and obtained numerical relations between the core-collapse time and the model parameters, which can be easily adapted for specific particle physics realizations of dissipative dark matter. Utilizing the density cores inferred from the rotation curves of 21 dwarf galaxies, we put strong constraints on the inelastic scattering cross section and the energy loss per collision. It is of interest to generalize our analysis to include velocity-dependent cross sections, which we leave for future work. Our formalism can be extended also to other scenarios, e.g., those proposed in [51][52][53], where dark matter particles are heated from energy release due to dark-sector interactions. hgal, Scott Tremaine, and Sean Tulin for useful discussions.

Supplemental Material
We provide additional information and results, including the derivation of the cooling rate (A), the procedure for numerical simulations (B), the condition for the mildcooling regime (C), calibration with isolated and cosmological N -body simulations (D), the halo parameters of 21 dwarf and LSB galaxies taken from [32], which were used to derive the bounds shown in Fig. 3 (E), snapshots of density profile for the benchmark points of LSB F583-1 (F), and additional details on our limit-setting procedure (G).

A. The Cooling Rate
The volumetric cooling rate C is defined as the energy loss per unit volume per unit time, which can be calculated as where . . . denotes the thermal average with respect to the relative velocity of the two incoming particles Assuming the dark matter velocity follows the Boltzmann distribution the thermal average over the relative velocity for quantity X is given by In computing the cooling rate, we impose a lower cut-off on v rel to take into account the fact that the energy loss during an inelastic collision should be smaller than the total kinetic energy of the two colliding particles in the center-of-mass frame, where µ = m/2 is the reduced mass. Thus, the cooling term can be written as If σ is velocity-independent, we have The cooling effect is small when ν ν loss , and it becomes significant when ν > ν loss .

B. Numerical Simulations
In performing our numerical simulations, we have used dimensionless variables by taking the ratio of a physical quantity (x) to its corresponding fiducial value, i.e.,x ≡ x/x 0 , and then rewrote Eq. (1) into dimensionless form. The fiducial quantities are built from the halo parameters ρ s and r s , as shown in Tab. I.  In the fluid model, a self-gravitating halo is assumed to be in hydrostatic equilibrium. The conservation of mass, momentum, and energy resembles the zeroth, first, and the second moments of the Boltzmann equation. Under the assumption that the velocity dispersion is spherically symmetric, they form a closed set of equations that truncates the Bogoliubov hierarchy. We segmented the halo into N = 150 evenly log-spaced concentric shells in radius {r 1 ,r 2 , · · · ,r N }. Following the treatment in Ref. [47], we take values of the extensive quantities (M i , L i ) at the radius of i-th sphere and those of the intensive ones (ρ i ,ν i ,Ĉ i ) as the average between i-th and (i−1)-th spheres. The algorithm uses Lagrangian zones: the radiî r i are allowed to change, but the mass in each shell is fixed.
Each step of the evolution is separated into two stages: thermal energy is exchanged by conduction and/or cooling, after which hydrostatic relaxation brings the system back to equilibrium. We assume that the entropy, s ∝ ln(ν 3 /ρ), is conserved during hydrostatic relaxation. The workflow is as follows: 1. Compute the luminosity and the cooling profileŝ L i andĈ i based on profile inputr i ,ρ i , andν i , and particle physics inputσ,σ , andν loss .
2. Allow a small passage of time ∆t and compute the specific energy change ∆û i (û i ≡ 3ν 2 i /2) by conduction and cooling, assuming fixed density. Eq. (1) gives and we updateû i by ∆û i . The time step ∆t is sufficiently small (|∆û i /û i | < 10 −3 ) such that the linear approximations used in step 3 below are valid.
3. Upon updatingû i , the i-th shell is no longer virialized. To return to hydrostatic equilibrium, we perturbr i ,ρ i , andν i while keeping the massM i and entropyŝ i = ln(ν 3 i /ρ i ) of the shell fixed. We treat mass conservation, hydrostatic equilibrium relations, and energy conservation, shown in Eq. (1), at the linear order and solve them for all shells simultaneously. For numerical accuracy, we iteratively perform the perturbation 10 times until hydrostatic equilibrium is established everywhere.
4. Re-establishing hydrostatic equilibrium gives new values forr i ,ρ i , andν i . We return to step 1 and update the luminosityL i and cooling profileĈ i .

5.
Track the Kn for the innermost shell. The evolution is terminated when Kn drops below 0.1 (stage 3).
The above procedure is coded in C++ with the eigen 3.2.8 library for linear algebra [54].

C. Strong Cooling vs. Mild Cooling
Cooling is strong if the kinetic energy gain of the infalling particles can be efficiently removed from the halo on a time scale smaller than the free-fall time, t ff = 3π/(32Gρ), which is close to t dy . This is similar to the "isothermal collapse" or "free-fall collapse" of the evolution of protostars (see e.g. [42,55]). A halo under isothermal collapse can fragment into multiple dark clumps and may also lead to the formation of a dark disk or a dark bulge. The strong cooling condition is set by t cool t ff , where the cooling timescale is t cool = (3/2)ρν 2 /C. In the strong cooling regime, the inner part of a dwarf-size halo will collapse and fragment within a time of t ff ∼ 0.1 Gyr for typical densities of ρ s = 10 6 − 10 7 M /kpc 3 . Since t ff 10 Gyr, the corresponding parameters for strong cooling can be excluded by dwarf galaxy observations. By definition, in the strong cooling regime the timescale of the thermal energy change (t cool ) is smaller than that of hydrostatic relaxation (t dy ). This invalidates the assumptions in the fluid simulations that the system stays in hydrostatic equilibrium. Hence, we should only trust our simulations for parameters that provide t cool t dy . As we show now, mild cooling indeed describes much of our parameter space.
The ratio between t cool and t ff is given by We take the initial density profile to be given by the NFW profile andσ ≤ 0.1. Forσ = 0.1, the ratio t cool /t ff at the inner most region (r = 0.01r s ) is ≥ 3 for all possiblê ν loss and ≥ 10 for eitherν loss ≤ 0.07 or ≥ 0.3. At larger radii (r ≥ 0.01r s ), Kn monotonically increases while ν monotonically decreases. Both factors contribute to a fast growth of t cool /t ff . Reducingσ increases t cool /t ff . Consequently, this extends the range ofν loss that permits t cool /t ff 1 and for which our analysis is valid. We explicitly checked that the parameter space for the 21 galaxies that saturate t c = 10 Gyr with β = 0.60 and 0.45 all have t cool /t ff > 30 1, i.e., the halo collapse begins when the cooling is still mild.
As the evolution proceeds and enters stage 3, Kn for the innermost region drops below 1. This means that at this point of the evolution, t cool /t ff ∝ Kn for the inner part of the halo also drops below 1, and the inner core may begin to fragment as it isothermally collapses. Nevertheless, given that the pre-factor (σ/σ )[(ν 2 loss /ν 2 )(1+ν 2 loss /ν 2 ) exp(−ν 2 loss /ν 2 )] −1 is greater than 1, the isothermal collapse is unlikely to happen before entering stage 3. Thus, it does not significantly affect our fluid-model-based collapse time estimation, nor does it affect the conclusion of a cuspy inner density profile. On the other hand, different techniques are necessary to understand the ultimate fate of the matter in this region and whether or not gravitational collapse to a massive black hole is possible.

D. Calibration and Cross-check with N -body Simulations
Unlike the conductivity in the small mean-free-path (smfp) region that is well-described by kinetic theory [56], the conductivity in the long mean-free-path (lmfp) region, κ lmfp , cannot be fully determined from first principles. The existence of the free coefficient β in the expression of κ lmfp reflects this fact (see discussion below Eq. (2)). Its value can be fixed by comparing the evolution of the inner density profile from the fluid model to that of the N -body simulations. This is possible, since κ lmfp ∝ β determines the collapse time of stage 2, which takes the bulk of the entire evolution time.
Various calibrations have been presented in the literature. Ref. [45,46] used N -body simulations of an isolated halo with an initial profile given by a self-similar core-like profile (a.k.a., the Balberg-Shapiro-Inagaki profile) and find β 0.75. However, Ref. [46] also showed that when applying β = 0.75 to a halo with an initial profile given by the NFW profile, the collapse time inferred from a fluid model is about ∼ 20% shorter than that from the N -body simulations. We re-examine this calibration by using our own analysis of the fluid model, which we compare to the N -body simulation data of Ref. [46]. The parameters for the simulated halo are ρ s = 1.49 × 10 6 M /kpc 3 and r s = 11.1 kpc, with an elastic cross section σ/m = 25.44 cm 2 /g [57]. We find that 0.59 ≤ β ≤ 0.61 agrees with will all three stages of the evolution, see Fig. 4. Note that the difference between β = (0.59, 0.60, 0.61) only becomes evident near the end of stage 2. We take the calibrated β (0.59 ≤ β ≤ 0.61) and cross check the fluid simulation results with that of the cosmological N -body simulations in Ref. [26]. Ref. [26] showed the snapshots of density and 3D velocity dispersion profiles for today for a dwarf galaxy halo that they call "Pippin". They are shown as solid lines in the left and right panels of Fig. 5 (note that we translate the 3D-velocity dispersion to a 1D-velocity dispersion by ν = ν 3D / √ 3). The parameters for Pippin [26] are M vir = 9 × 10 9 M , r s = 2.7 kpc, and ρ s = 1.73 × 10 7 M /kpc 3 , where the subscript "vir" indicates that the halo boundaries are set by their virial radius. On the fluid simulation side, we assume the initial profile to be an NFW profile for σ/m = 0 today. We also assume the evolution time of the halo is 10 Gyr. We then take the Pippin halo setup and run the fluid simulations for σ/m = 0.1, 0.5, 1, 5, 10, and 50 cm 2 /g. The simulations are truncated at 10 Gyr and snapshots for the density and velocity dispersion profiles with β = 0.60 are shown as dashed lines in the left and right panel of Fig. 5, respectively. Varying β to 0.59 from 0.61 does not bring noticeable changes to the snapshots. Generically, the fluid model overestimates the density and underestimates the velocity dispersion in comparison with the values provided by the cosmological N -body simulations. Nevertheless, the discrepancies in the density and velocity-dispersion snapshot are within a factor of 2 and 1.2 respectively.
For σ/m = 50 cm 2 /g, Pippin sees a "moderate collapse". This is also observed in the fluid model, where we see the halo evolves close to the end of stage 2 at t = 10 Gyr with σ/m = 50 cm 2 /g. Nevertheless, the Pippin halo exhibit a slower collapse than the fluid simulated halo with β = 0.60. To better see this, we notice that κ lmfp ∝ σ/m, i.e., a SIDM halo with σ/m = 5 cm 2 /g σ/m [cm 2 /g] evolves 5 times faster than that of σ/m = 1 cm 2 /g during stage 2, the bulk of the evolution. This proportionality induces a degeneracy: the snapshot for a Pippin halo with σ/m = 5 cm 2 /g at 10 Gyr is identical to that with σ/m = 1 cm 2 /g at 50 Gyr. Adopting this degeneracy, we translate the density snapshots for Pippin halo with σ/m = 0.1, 0.5, 1, 5, 10, and 50 cm 2 /g at t = 10 Gyr into snapshots with σ/m = 1 cm 2 /g at t = 1, 5, 10, 50, 100, and 500 Gyr respectively. The averaged inner densities for each snapshot are shown as colored circles in Fig. 6. In comparison with the fluid simulations with β = 0.60 (dashed line), the Pippin simulations show a slower collapse between 100 Gyr and 500 Gyr. This slower collapse is better captured by the fluid simulations with β = 0.45 (dotted line), since that a smaller β implies a slower collapse. In Fig. 5, we thus also show snapshots with β = 0.45 for various cross section strengths at 10 Gyr as dotted lines.
The change of β from 0.60 to 0.45 can be translated to a lengthening of evolution timescales by a factor 4/3. Thus, the limits in Fig. 3 that saturate t c = 10 Gyr with β = 0.45 can be interpreted as the limits saturate a shorter collapse timescale t c = 7.4 Gyr with β = 0.60.
Our fluid simulations do not take into account mergers and environmental effects of the continuous infall of the background matter and the presence of baryons. A major merger can rebuild cuspy profiles and reset the evolution clock, delaying the thermal evolution. The selfsimilar accretion shock heating can prevent gravothermal collapse of an SIDM halo if the mass accretion rate is high [58]. However, cosmological N -body simulations show that both major merger and mass accretion rates decrease sharply towards low redshifts, and they become negligible for halos at present [59,60]. This may explain the good agreement between the simulations and the semi-analytical estimates for the Pippin halo. Thus, we expect our results are robust for near-field galaxies at redshift 0. Since our constraints on the dissipation parameters are based on the dwarf galaxies with low baryon content, baryons do not play an important role in setting the limits. It is plausible that the presence of baryons would speed up core collapse [61,62].

E. Data of Dwarf and LSB Galaxies
In Tab. II, we list all dwarf/LSB galaxies with low baryon-content fitted in [32]. We translated the dex on the concentration c 200 into the value of c 200 according to the concentration function used in [32] (taken from [63]). The subscript "200" indicates that the halo boundaries are set by the radius where the averaged density is 200 times of the critical density of the Universe. Based on the concentration c 200 and total mass M 200 , we determined r s and ρ s and further construct fiducial quantities listed in Tab. I for each halo. In Sec. IV, we use LSB F5831-1 as an example and illustrate the estimates of the collapse time t c , t c for five particle physics benchmarks (summarized in Tab. III for the reader's convenience). In Fig. 7, we show the snapshots of the density evolutions for each benchmark to  Fig. 7 represent density snapshots at different times. The black dashed line shows the density profile inferred from the rotation curve fit. The 3D views of the density evolution of benchmark (a) and (b) are shown in Fig. 1 left.

G. Procedure for Setting Limits on Dark Matter Parameter Space
In this subsection, we discuss our procedure for determining the excluded region of Fig. 3 in further detail.
Our criterion for excluding a point in the dissipative dark matter parameter space is equivalent to ruling out those parameters for which the dwarf and LSB halos reach stage 3 of their gravothermal evolution in less than 10 Gyr. This criterion is conservative, since it requires that a cusp has been reestablished through the entire inner core region. Before our criterion is achieved, several other deviations from the cored profile may become apparent. For example, the central density of the halo will begin to exceed its minimum value well before the runaway collapse of stage 3 begins. Likewise, the log-slope of the density becomes steeper than the cored prediction at some point during stage 2.
However, these alternative criteria are somewhat more difficult to observe rigorously. As discussed in Sec. D, the core density from a fluid simulation is up to a factor of 2 greater than that from the cosmological N -body simulation. We should regard this difference as a systematic uncertainty of the fluid model. To claim a robust exclusion based on the central density, a very large difference in the central density compared to the cored profile is necessary.
To illustrate our exclusion criterion, consider again Fig. 7.
As seen from panel (a), when the cooling is absent (benchmark (a)), the evolution already develops a large core around 2 Gyr for σ/m = 3 cm 2 /g. Collapse eventually happens, but takes about 173 Gyr. The density profile is almost static for most of the halo's existence. The prediction of the central density of the fluid model is about 1.3 times larger than the inferred density profile from the rotation curve fit of LSB F583-1, but this difference should be tolerated given the uncertainties in the rotation curve measurement and the systematic uncertainties in the fluid model.
Next we look at panels (b) and (c). The corresponding benchmark particle physics parameters (b) and (c) of Tab. III are inside the excluded region for LSB F583-1 but near the boundary, as shown in Fig. 3. As is apparent in these panels, the inner core density is already 10 times larger than from the fit at an age of ∼ 7 Gyr. After the halo is ∼ 8 Gyr old, the entire inner halo is cuspy and the evolution reaches t c . It is also worth noting that the core density (core size) at the maximal expansion is higher (smaller) than that of benchmark (a). This is because at stage 1, bulk cooling counters the inward heat flow, and hence limits the development of the core. As a result, the particles in the inner halo experience fewer collisions during stage 1. This yields a smaller core with a higher core density.
As we go deep into the excluded parameter space, the collapse time t c becomes even smaller, and the core size (density) at the maximal expansion becomes even smaller (higher) as shown in panel (d). Finally, when we reach benchmark (e) (panel (e)), we are near the strong cooling region and the collapse time is 223 Myr, which is near the free-fall time for the inner core. As shown in panel (e), the self-interactions never produce a large core beyond 0.1r s . This yields an obvious contradiction with the observed density profile.