Constraining Velocity-dependent Self-Interacting Dark Matter with the Milky Way's dwarf spheroidal galaxies

The observed anti-correlation between the central dark matter (DM) densities of the bright Milky Way (MW) dwarf spheroidal galaxies (dSphs) and their orbital pericenter distances poses a potential signature of self-interacting dark matter (SIDM). In this work we investigate this possibility by analysing the range of SIDM scattering cross section per unit mass, $\sigma/m_{\chi}$, able to explain such anti-correlation. We simulate the orbital evolution of dSphs subhaloes around the MW assuming an analytical form for the gravitational potential, adopting the proper motions from the Gaia mission and including a consistent characterization of gravitational tidal stripping. The evolution of the subhaloes density profile is modelled using the gravothermal fluid formalism, where DM particle collisions induce thermal conduction that depends on $\sigma/m_{\chi}$. We find that models of dSphs, such as Carina and Fornax, reproduce the observed central DM densities with fixed $\sigma/m_{\chi}$ ranging between $30$ and $50$ cm$^{2}$g$^{-1}$, whereas other dSphs prefer larger values ranging between $70$ and $100$ cm$^{2}$g$^{-1}$. These cross sections correlate with the average collision velocity of DM particles within each subhalo's core, so that systems modelled with large cross sections have lower collision velocities. We fit the cross section-velocity correlation with a SIDM particle model, where a DM particle of mass $m_{\chi}=0.648\pm 0.154$ GeV interacts under the exchange of a light mediator of mass $m_{\phi}=0.636\pm 0.055$ MeV, with the self-interactions being described by a Yukawa potential. The outcome is a cross section-velocity relation that explains the diverse DM profiles of MW dSph satellites and is consistent with observational constraints on larger scales.


INTRODUCTION
The standard cosmological paradigm of Λ collisionless cold dark matter (ΛCDM) accurately predicts the large-scale structure of the Universe (Springel et al. 2006), however its success is less certain over the small scales, the regime relevant to the substructure within galactic haloes. The deficit of observed low-mass satellite galaxies within the Local Group relative to predictions from analytical theory (Press & Schechter 1974) and CDM N-body simulations (Klypin et al. 1999), along with the discrepancy between the low dark matter densities of some galaxies and the cuspy and dense subhaloes predicted by simulations (e.g. Flores & Primack 1994;Moore 1994;Moore et al. 1999), motivated to question the cold and collisionless nature of dark matter (DM). The alternative was that DM might have non-negligible selfinteractions (Spergel & Steinhardt 2000).
Self-interacting dark matter (hereafter SIDM) assumes E-mail: camila.correa@uva.nl that DM particles experience collisions with each other, these collisions transfer heat towards the colder central regions of DM haloes, lowering central densities and creating constant density cores (e.g. Davé et al. 2001;Colín et al. 2002;Vogelsberger et al. 2012;Rocha et al. 2013;Dooley et al. 2016;Vogelsberger et al. 2019;Robles et al. 2019). DM particle collisions also lead to less concentrated subhaloes that are more prone to tidal disruption, as well as to the evaporation of subhaloes via ram pressure stripping exerted by the larger host. In this manner various SIDM models have produced halo mass functions that are significantly suppressed on small scales (e.g. Buckley et al. 2014;Cyr-Racine et al. 2016;Vogelsberger et al. 2016, who also assumed a suppression in the power spectrum, and e.g. Nadler et al. 2020 that assumed standard SIDM). An important disagreement between the prediction of pure CDM simulations and observations is the so-called toobig-to-fail problem (hereafter TBTF), which states that the most massive subhaloes in CDM simulations are too dense in the centre to host the observed satellites of the Milky Way (Boylan-Kolchin et al. 2011. The TBTF problem, although solved by invoking baryonic physics and environmental effects (e.g. Sawala et al. 2016;Dutton et al. 2016;Wetzel et al. 2016;Fattahi et al. 2016), has been revisited in the SIDM paradigm due to its prevalence in other galaxies besides the Milky Way (e.g. Garrison-Kimmel et al. 2014;Papastergis et al. 2015).
The SIDM models that alleviate the TBTF problem require a DM interaction cross section per unit mass, σ/m χ , to be larger than σ/m χ >1 cm 2 g −1 (Zavala et al. 2013). Despite its success in solving the TBTF problem and other small-scales discrepancies (e.g. Rocha et al. 2013;Kamada et al. 2017;Ren et al. 2019), the excitement caused by SIDM diminished due to the strong constraints set by X-ray and lensing observations of galaxy clusters, that set an upperlimit on the scattering cross section of the order of 1 cm 2 g −1 (Miralda-Escudé 2002;Peter et al. 2013). This upper-limit was later supported by observations of bounds from major mergers (Randall et al. 2008;Harvey et al. 2015;Wittman et al. 2018) and bright central galaxy wobbles (Kim et al. 2017;Harvey et al. 2019).
An exciting alternative for SIDM, however, is that the cross section may depend on the relative velocity of DM particles, in such a way that DM behaves as a collisional fluid on small scales while it is essentially collisionless over large scales. Such velocity-dependence is supported by many theoretical models that argue that DM exists in a 'hidden sector', where forces between DM particles are mediated by analogues to electroweak or strong forces (e.g. Pospelov et al. 2008;Arkani-Hamed et al. 2009;Buckley & Fox 2010;Feng et al. 2010;Boddy et al. 2014;Tulin & Yu 2018). Velocitydependent SIDM models have been explored on galaxy cluster scales (e.g., Robertson et al. 2017Robertson et al. , 2019Banerjee et al. 2020) and Milky Way (MW)-mass systems (e.g. Vogelsberger et al. 2012;Zavala et al. 2013;Nadler et al. 2020).
The works of Read et al. (2018) and Valli & Yu (2018) analysed the density profiles of the MW dwarf spheroidal galaxies (dSphs), aiming to place constraints on σ/m χ on dwarf galaxy scales. Read et al. (2018) focused on Draco and claimed that its high central density gives an upper bound on the SIDM cross section of σ/m χ <0.57 cm 2 g −1 . Valli & Yu (2018) used a similar methodology as Read et al. (2018), although they also found that Draco's high central density could be described by a σ/m χ ∼0.4 cm 2 g −1 , they concluded that the remaining dSphs probed different cross sections, ranging between 0.1 and 40 cm 2 g −1 .
Recently, Kaplinghat et al. (2019) have revisited the TBTF problem of MW dSphs. They have reported an interesting anti-correlation between the central DM densities of the bright dSphs and their orbital pericenter distances, so that the dSphs that have come closer to the MW centre are more dense in DM than those that have not come so close. Read et al. (2019) proposes that the anti-correlation is the result of baryonic effects. The gas expelled by stellar feedback 'heats' the surrounding DM, lowering the haloes' central density (Navarro et al. 1996). If the effect repeats over several cycles of star formation, it accumulates, leading eventually to transform the DM cusp into a core (e.g. Read & Gilmore 2005;Pontzen & Governato 2012). However, for dwarf galaxies, such as Draco and Ursa Minor, that are DM-dominated and stopped forming stars long ago (∼ 10 Gyr), the DM cusp is formed again, thus explaining the high DM central densities. Read et al. (2019) modelled the stellar kinematics to infer the DM distribution of the MW dSphs, and showed that all dwarfs except for Fornax, are well fitted by a cuspy profile. Even if the 'lack' of baryonic effects is responsible for the high DM densities in the dSphs, it does not explain the origin of the anti-correlation with pericenter distance.
SIDM, on the contrary, can potentially explain the observed anti-correlation. DM collisions lead to an outward heat transfer that induces gravothermal core collapse, i.e. the central density increases with time (Balberg et al. 2002;Elbert et al. 2015). This gravothermal collapse would be accelerated by mass loss via tidal stripping (Nishikawa et al. 2020) and correlate with how close the satellite galaxies come to the centre of the MW. According to this scenario, the anti-correlation could not only be a potential signature of SIDM Nishikawa et al. 2020;Sameie et al. 2020), but it would also invalidate the upper limits on σ/m χ from Read et al. (2018) and Valli & Yu (2018), since their analysis does not include the effects of gravothermal collapse.
The goals of this study are to investigate the possibility that the anti-correlation is a signature of SIDM, and analyse the ranges of σ/m χ that could potentially explain it. To do so, we simulate the orbital evolution of dSphs subhaloes around the MW by adopting the proper motions from the Gaia mission (Helmi et al. 2018;Fritz et al. 2018), assuming an analytical form for the MW gravitational potential and including a consistent characterization of gravitational tidal stripping. The evolution of the density profile of SIDM subhaloes is simulated using the gravothermal fluid formalism, which was originally developed to study the evolution of globular clusters (Lynden-Bell & Eggleton 1980), but it has also been applied to isolated SIDM haloes (Balberg et al. 2002;Koda & Shapiro 2011;Shapiro 2018). This method allows us to track the DM halo evolution within scales smaller than 100 pc, which are largely expensive to resolve with Nbody simulations, as well as to easily cover a wide range of parameter space.
This work is organised as follows. Section 2 outlines our model setup. We present our results in Section 3 and compare to observational data. We discuss constraints on the cross section-velocity relation in Section 3.3. Comparison with previous works, as well as the challenges of the model and impact of initial conditions are discussed in Section 4. Finally, Section 5 summarises our key results.

SIDM HALO MODEL
The SIDM halo model derived in this work connects the gravothermal fluid approximation, with orbit integration and tidal stripping modelling. Subsection 2.1 describes the gravothermal fluid formalism, the orbital evolution is introduced in Subsection 2.2 and Subsection 2.3 describes the gravitational tidal stripping modelling. We summarise the complete model and describe the initial conditions in Subsection 2.4.

Gravothermal collapse
We consider a spherical halo in isolation and in quasi-static virial equilibrium, with a density profile ρ(r, t) and an enclosed mass of m(<r, t) at radius r and time t. We assume that the ensemble of gravitating particles is well approximated by a fluid-like description, where the effective temperature is identified with the square of the one-dimensional velocity dispersion, v(r, t), and thermal heat conduction is employed to reflect the manner in which the close-encounter large-angle (hard-sphere) scatterings combine to transfer energy in the system. The quasi-static approximation means that, while the fluid evolves thermally, it always satisfies hydrostatic equilibrium at each moment.
The fundamental equations of the model are mass conservation, hydrostatic equilibrium, energy flux equation, and the first law of thermodynamics, where L(r) the luminosity through a sphere at r. The flux equation (eq. 3) can be written into a single expression that considers both the cases where the the mean free path between collisions is significantly shorter (known as SMFP regime) or larger (LMFP regime) than the system size (Balberg et al. 2002), as follows where H ≡ v 2 /(4πGρ) is the gravitational scale height of the system, λ is the collisional scale for the mean free path given by λ = 1/(ρσ m ), with σ m = σ/m χ the cross section per unit mass, and t r ≡ λ/(av) is the relaxation time, with a = 16/π for hard-sphere scattering of particles with a Maxwell-Boltzmann velocity distribution (Balberg et al. 2002).
In the SMFP regime the continuum assumption applies and the collection of particles can accurately be treated as a continuous fluid. From this regime, the effective impact parameter b = 25 √ π/32 ≈ 1.38 in eq. (5) can be derived from first principles (e.g. Chapman et al. 1990). In the LMFP regime the model needs to be calibrated using N-body simulations. Previous studies have tested the parameter C that determines the radial heat conduction for isolated (Balberg et al. 2002;Koda & Shapiro 2011;Essig et al. 2019) and cosmological N-body simulations (Elbert et al. 2015;Nishikawa et al. 2020;Essig et al. 2019) with purely elastic DM selfinteractions. Essig et al. (2019) and Nishikawa et al. (2020) assumed spherical symmetry but not isolation. Essig et al. (2019) used the cosmological simulation from Elbert et al. (2015) and showed that C = 0.45−0.6 closely reproduces the density and velocity dispersion from the simulation. Nishikawa et al. (2020) also used Elbert et al. (2015) simulation and con-cluded that for σ/m χ > 10 cm 2 g −1 , C = 0.75 and b = 0.003 is needed to reproduce the subhalo's DM density. Both studies also reported differences with respect to the simulation in the subhalo density of up to a factor of 2 for σ/m χ = 50 cm 2 g −1 .
We take this discrepancy into account when we constrain the cross section in Section 3.3, and adopt C = 0.75 as reported by Balberg et al. (2002), who assumed spherical symmetry and isolation for the modelling of SIDM haloes. However, we also discuss how the different values of C = 0.45 − 0.6 and b = 0.003 − 1.38 impact on our key results in Subsection 4.2.

Orbital evolution of MW spheroidal galaxies
Throughout this work we model the internal dynamics and orbital evolution of the nine most luminous MW dSph galaxies. These include Ursa Minor (hereafter UM), Draco, Sculptor, Sextans, Fornax, Carina, LeoII and LeoI, and Canes Venatici I (hereafter CVnI). We focus on these systems because they have the highest quality kinematic data and the largest samples of spectroscopically confirmed member stars to resolve the dynamics at small radii.
The second data released by the Gaia mission (Brown et al. 2018;Helmi et al. 2018) has largely increased the precision and amount of astrometric data of Galactic stars, making possible the determination of spatial motions of many dSphs orbiting the MW halo (Helmi et al. 2018;Fritz et al. 2018). Using the proper motions determined by Fritz et al. (2018), and the publicly available code galpy 1 (Bovy 2015), we integrate the orbits of the dSphs adopting the static MW-Potential14 model, which has been shown to be consistent with various observations (see Bovy 2015 for details). For the MW dark matter halo mass we assume a virial mass of M 200 = 10 12 M , defined as the total within R 200 , radius within which the mean density is equal to 200 times the critical density of the Universe, ρ crit . In Appendix D we show that assuming a lighter MW halo mass of 0.8 × 10 12 M or a heavier model of 1.6 × 10 12 M does not largely affect our key results. Fig. 1 shows the time evolution of the galactocentric distance of the MW dSph galaxies. The black dashed line indicates the time evolution of the MW's virial radius, R 200 , calculated using the halo accretion history model from Correa et al. (2015a,c) 2 . It can be seen from Fig. 1 that the dSphs of UM, Draco, Sculptor and Carina became satellites of the MW nearly 8-9 Gyrs ago and since then they have completed many orbits, whereas LeoII, Sextans, CVnI and Fornax have completed two orbits around the MW. LeoI crossed the MW's virial radius for the first time roughly 2 Gyrs ago.

Gravitational tidal stripping
An important aspect in the evolution of small subhaloes relative to haloes in the field, is that when subhaloes are accreted by a larger halo, they begin to lose mass due to the We calculate the rate of mass loss, dm/dt, of the subhaloes hosting the dSphs as they orbit around the MW by adopting the following tidal stripping rate where m(>r t ) is the subhalo mass outside the instantaneous tidal radius r t , τ orb = 2π/ω with ω the instantaneous angular velocity of the subhalo, and α = 1 (see van den Bosch et al. 2018). The tidal radius is calculated as which corresponds to the scenario where a subhalo of mass m is on a circular orbit of radius R, with angular speed Ω(= V circ (R)/R), around a halo of mass M (e.g., King 1962;Tollet et al. 2017). Eq. (7) gives the approximate amount of mass stripped from the subhalo over a short time-step, but it does not indicate how the density profile is modified by it. To model the truncation in the density profile, we employ the fitting functions of , that follow the structural evolution of a tidally truncated subhalo and solely depend on the fraction of mass stripped.
Green & van den Bosch (2019) used the DASH library (Ogiya et al. 2019) of high-resolution, idealized dark matter only collisionless N-body simulations that follow the evolution of an individual subhalo as it orbits within the fixed, analytical potential of its host halo. Both the fixed host halo and the initial subhalo are spherically symmetric, each with a Navarro-Frenk-White (hereafter NFW, Navarro et al. 1997) density profile, where ρ s and r s are the scale density and radius where the logarithmic density slope is equal to −2. The ratio between the halo's virial radius R 200 and the scale radius, r s , defines the concentration parameter c 200 = R 200 /r s of the profile.
Green & van den Bosch (2019) provide the best-fit parameters for the transfer function, H(r, t, f b , c 200 ), defined as the ratio of the evolved subhalo density profile relative to the initial profile, H(r, t, f b , c 200 ) = ρ(r, t)/ρ(r, t = 0), with f b the bound fraction (mass that remains bound to the subhalo while it is tidally stripped) and c 200 subhalo concentration parameter at t = 0 (see Green & van den Bosch 2019 for more details).
In our model, however, the subhalo DM density profile depends exclusively on the density profile at the previous time-step, not on the initial profile. We therefore assume that the density profile at time-step t n , ρ(r, t n ), is calculated from the density profile at a previous time-step t n−1 , ρ(r, t n−1 ), via the transfer function as where H depends on the bound fraction, defined as the fraction of mass that remains bound after it lost dM mass between t n and t n−1 , and c 200 (t n−1 ) the concentration parameter of the density profile at t n−1 . Although we calculate the amount of mass loss during each t n , we do not apply eq. (10) at the end of every time step. This is because the time-step size can become very small, making dM a negligible quantity. Instead, we apply eq. (10) and truncate the density profile every 250 Myr. We have found that during this period of time, the cumulative mass loss of subhaloes reaches on average 1 − 2% of their total mass. In Appendix E we show that truncating the density profile every 350 Myr, instead of 250 Myr, slightly decelerates gravothermal collapse, conversely a more frequent truncation of the density accelerates gravothermal collapse. Section 5 discusses the impact of the truncation time parameter on our key results. Note that eq. (10) is a strong variation of the Green & van den Bosch (2019) model. We compare the outcome by evolving a 10 8.84 M subhalo that has an initial NFW profile with concentration 15.7, it follows the orbit of UM (shown in top panels in Fig. 1) and loses roughly 40% of its initial mass. We evolve the subhalo using the Green & van den Bosch (2019) model, as well as the modified model shown in eq. (10).
The top panels of Fig. 2 show the density (top-left) and enclosed mass (top-right) as a function of radius of a NFW subhalo. The dashed lines in the panels correspond to the initial density and mass profile, whereas the solid lines correspond to the final profiles after 10 Gyr of evolution. From the top-left panel it can be seen that the density profile is largely truncated when we apply eq. (10) in comparison to the Green & van den Bosch (2019) model, however the large difference only occurs in the outskirts of the halo, beyond the virial radius. The top-right panel shows that the final masses agree, indicating that the truncation imposed by eq. (10) closely follows the rate of mass loss of the  model.
An important caveat to consider when we apply eq. (10) is that the transfer function H was fitted according to the structural evolution of an NFW CDM subhalo in the outer regions. Differently, the density profile of SIDM subhaloes largely deviates from the NFW shape in the inner regions, but not in the outer regions, since there the rate of DM-DM particle interactions is low. Given that the transfer function mostly affects the outer regions of the density profile (as shown in the top panels of Fig. 2), we test the evolution of a SIDM subhalo with an initial cored-shape density profile. The bottom panels of Fig. 2 show the same as the top but for the SIDM subhalo. After 10 Gyrs of evolution, the subhalo loses roughly 43% of its initial mass, in close agreement with the evolution of the CDM subhalo. Like in the top panels, the density profile is largely truncated when we apply eq. (10) in comparison to the Green & van den Bosch (2019) model, but mostly beyond the virial radius.
Given the good agreement between the rate of mass loss and truncated profiles of SIDM and CDM subhaloes in the outer regions, we assume that applying eq. (10) to the truncation of SIDM subhaloes is a good approximation. Note, however, that the rate of mass loss given by eq. (7) should be taken as a lower limit. SIDM subhaloes lose more mass due to ram-pressure stripping and the presence of baryons, effects that we do not model in this work. Ram-pressure stripping is caused by DM self-interactions with the host halo particles, that drive material out of subhaloes, to the extend of being able to completely evaporate subhaloes (e.g. see Vogelsberger et al. 2019). Baryons will not affect the gravothermal SIDM modelling presented in Subsection 2.1, but will enhance the effect of tidal stripping due to the presence of a Galactic disk.

Integration of the equations & initial conditions
The gravothermal model comprises eqs.
(1-5) that govern the evolution of the subhalo's density profile given the cross section per unit mass, σ/m χ , and the initial subhalo density profile, ρ init . We set ρ init to follow an NFW profile, and divide the spherical subhalo into 150 logarithmically spaced concentric shells, ranging between r min = 10 −2 kpc and r min = 10 2 kpc.
We solve the gravothermal model by re-writing the equations into non-dimensional form, to do so we introduce a characteristic mass, density and radius, and numerically integrate the equations over a time-step ∆t. For each iteration, ∆t is restricted to be where "∼" denotes the variables non-dimensional form, and ∆r the radial spacing of the density profile. We further describe the non-dimensional terms and solution of the gravothermal equations in Appendix A.
We combine the gravothermal model with the orbit integration model, so that at each time-step we calculate the subhalo's distance to the galaxy centre, d GC , as well as the instantaneous angular velocity, ω. We use these quantities Top panels: Density (left) and enclosed mass (right) as a function of radius for a subhalo that has an initial virial mass of 10 8.84 M , an initial NFW profile with concentration 15.7, it follows the orbit of UM over 10 Gyr as it loses 40% of its initial mass ( f b = 0.6). The dashed line corresponds to the initial density and mass profile, whereas the solid lines correspond to the final profiles after 10 Gyr. Bottom panels: same as top panels but for a subhalo with an initial virial mass of 10 9.07 M and an initial cored profile. From the left panels it can be seen that the density profiles are largely truncated when we apply eq. (10) in comparison to the Green & van den Bosch (2019) model, however the large difference only occurs in the outskirts of the halo, beyond the virial radius. The right panels show that the enclosed mass of the two models agree when initializing with NFW-shape or cored-shape profile.
to determine the amount of mass lost between two consecutive time-steps using eqs. (7) and (8). At each time step we also evolve the MW's virial mass following eq. (6), as well as the MW's halo density profile, which we assume to be an NFW profile that follows the concentration-mass relation of the form from Correa et al. (2015c). Every 250 Myrs, we apply eq. (10) in order to truncate the density profile according to the amount of mass lost during that period.
Subhaloes hosting the MW dSphs are initialised with the orbital parameters taken from Fritz et al. (2018), namely the distance to the MW centre, d GC , radial velocities, v R and tangential velocities, v T , which are used to calculate the orbital evolution of each subhalo. Fritz et al. (2018), as well as Helmi et al. (2018), reported uncertainties in the dwarf distances and radial velocities of the order of 7−8%. We have found that such errors do not change the results presented in the following section. The errors of the tangential velocities, however, are larger, of the order of 20% for all dwarfs except CVnI, LeoI and LeoII, which range between 60 and 110%. These errors can change the orbits of the dwarfs to a great extent, thus altering the rate of mass loss and affecting the period of gravothermal collapse. In Section 4.5 we discuss how the errors in v T impact on our results.
Besides the orbital parameters, each subhalo is initialised with two free parameters: the cross section per unit mass, σ/m χ , and the initial subhalo virial mass, M 200,init . The latter is used to estimate the DM halo concentration parameter at the initial cosmic time t = 3.5 Gyr (redshift z = 1.87), using the concentration-mass relation of Correa et al. (2015c), which in turn is used to initialise ρ init . Table 1 lists the orbital parameters taken from Fritz et al. (2018), the initial virial mass and set of parameters that describe the initial NFW profile for each subhalo. Note that the initial concentrations for the dSphs are quite low (c 200 ∼ 6 − 7), in opposite to typical z = 0 values for 10 9 M systems of c 200 ∼ 15−20. This is because of two reasons. First c 200,init is set by assuming that 10 Gyrs ago, before the dSph galaxies became MW's satellites, those galaxies were hosted by field haloes that followed the median concentration-mass relation for z = 1.87. Secondly, the concentration-mass relation from Correa et al. (2015c) is not a best-fit extrapolation from cosmological simulations, it is a semi-analytic model that combines an analytic model for the halo mass accretion history, based on extended Press Schechter (EPS) theory (Press & Schechter 1974), with an empirical relation between concentration and formation time (Correa et al. 2015b). Because the semi-analytic model is based on EPS theory, it can be applied to wide ranges in mass, redshift and cosmology. Throughout this work we assume Planck13 cosmology (Planck Collaboration et al. 2014) with Ω m , Ω Λ , h, σ 8 , n s equal to 0.307, 0.693, 0.6777, 0.8288, 0.9611, respectively. In Section 4.4 we discussed how changing the initial values of the concentration parameter impacts on our results.
We run the SIDM halo model for the nine systems hosting the most massive MW dSphs. The evolution begins when the Universe is 3.5 Gyr old, at a point when none of the systems have yet crossed the MW's virial radius, and it finishes at present time, covering 10.2 Gyrs of evolution. In Section 4.3 we show that using an NFW profile for the initial density profile of subhaloes and MW, rather than a cored profile, does not modify our results.

RESULTS
In this Section we present the results obtained with the SIDM halo analytic model. We begin by illustrating how the joint framework of gravothermal evolution and gravitational tidal stripping shapes the evolution of the density profile of each system. We next describe the range of values of the free parameters, σ/m χ and M 200,init , that reproduce the central DM densities reported in Kaplinghat et al. (2019). Finally we show that there is a promising range of velocity-dependent cross section models that explain the anti-correlation of central density and pericenter distance for the nine most massive MW dSphs.

SIDM halo evolution
In this Section we show the evolution of the subhalo that hosts the galaxy Carina. For this system the model was initialised with a cross section of σ/m χ = 40 cm 2 g −1 . We as-sume that 10.2 Gyrs ago (redshift z = 1.87), Carina had a virial mass of M 200 = 2 × 10 9 M , and its density followed the NFW profile with a scale density and radius of 4.2 × 10 6 M kpc −3 and 2.09 kpc, respectively. Fig. 3 shows the density (left panels) and velocity (right panels) profiles at different times. The velocity corresponds to the particles' average collision velocity, v = (4/ √ π)v (for a Maxwellian distribution), with v the 1-D velocity dispersion. Each line in the figure is coloured according to the lookback time as shown by the colour bars on the right. In the top panels, the dark blue lines correspond to the density and velocity of the system when it begins to evolve, 10.2 Gyrs ago, whereas the orange lines correspond to the evolution of the system between 7.4 and 8.8 Gyrs ago. It can be seen that after a few time-steps the cusp in the central region disappears, the central density decreases and a core of roughly constant density begins to form. The top right panel shows that while the density profile rapidly forms a central core, the particles' velocity increases. The density in the inner regions decreases due to DM-DM collisions, that expel some DM particles from the central regions into further out orbits, at the same time the velocity increases because collisions increase the mean velocity of particles.
DM-DM particle collisions add energy to the core, causing particles to move into larger orbits at lower velocities. Through collisions the subhalo's core becomes a system with negative heat capacity, where adding energy cools down the system, while the extended subhalo becomes a large thermal reservoir that absorbs the core's energy. The subhalo stabilises as it forms a high temperature core with negative heat capacity. This is the period between the end of the core expansion phase and the beginning of the gravothermal collapse-phase, that began roughly 9 Gyrs ago for this system and lasts for roughly 3 Gyrs. The bottom panels of Fig. 3 show the evolution of the subhalo in the gravothermalcollapse phase. During this phase, the high temperature, negative-heat capacity core in contact with the cold extended halo gives up heat, getting hotter rather than colder. The hot core then contracts, the central density increases, leading to the gravothermal collapse phase. In the case of Carina, the subhalo reaches a stable central density of 10 7 M kpc −3 . During the last 7 Gyrs, as it goes into the gravothermal collapse phase, its density increases an order of magnitude, reaching a value of 2 × 10 8 M kpc −3 .
An important aspect in the evolution of Carina is the result of the joint gravothermal evolution and gravitational tidal stripping modelling. Differently from previous studies, the contraction phase of the core is not followed by an increase in the particle's velocity as it would be expected. Instead the particles' velocity decreases during the last 4 Gyr of evolution, this is because during this period hydrostatic equilibrium significantly lowers the pressure of the subhalo when it loses mass, causing the velocity to decrease. Since heat flows towards the colder extended halo, heat is diffused at a faster rate when mass is tidally stripped, leading to a faster formation of the isothermal core and thus an accelerated evolution for core collapse.
All systems undergo a similar evolution to the one described in this section, with the only difference that a few systems reach a higher or lower central density during the gravothermal collapse phase, and others lose more or less  Table 1. Form left to right: list of orbital parameters and initial conditions. The first column indicates the name of the dSph galaxy that corresponds to the observational estimates for the galactocentric distance, d GC , radial and tangential velocities, v R and v T , taken from Fritz et al. (2018). The fifth and sixth columns from the left correspond to the initial virial mass and concentration, M 200,init and c 200,init , each subhalo is initialised at cosmic time 3.5 Gyr (z = 1.87) before infalling onto the MW system. The seventh and eighth columns indicate the respective scale density and radius, ρ s and r s , of the initial NFW density profile, ρ init . mass as they orbit around the MW. The following section describes the dependence of central density evolution on the scattering cross section.

Central density evolution
The evolution of the central DM density of the subhalo, along with its mass loss rate, largely depends on the scattering cross section. At fixed initial mass, a large cross section leads to a larger rate of DM-DM collisions that produce a shallower and lower density core. Similarly, the larger rate of DM-DM collisions leads to less concentrated subhaloes, making them more prone to tidal disruption and mass loss.
This dependency on the cross section can be seen in   The left panel of Fig. 4 shows that the central density quickly drops when the core of the subhalo forms, and it rises again as the core begins to collapse. For both cases, with or without tidal stripping, the central density reaches a minimum stable value, roughly independent of the cross section. For the model that includes mass loss from tidal stripping, the collapse time becomes shorter than the age of the Universe (as also shown by e.g. Nishikawa et al. 2020), and the central density reaches higher values for a higher cross section.
The right panel shows that for the case of no tidal stripping, the subhalo's virial mass slightly increases during its evolution, this is because M 200 is calculated by integrating the subhalo's density profile until the mean density is 200 × ρ crit , with ρ crit (z) decreasing with decreasing lookback time (decreasing redshift). For the model with tidal stripping, the panel shows that during the last 2.5 Gyr, the rate of mass loss is higher for higher cross sections. Fig. 4 shows that cross sections ranging between σ/m χ = 32 cm 2 g −1 and σ/m χ = 40 cm 2 g −1 are able to explain the observed central DM density of Carina (within the uncertainty), as well as the final virial mass of the system. However, this range of parameters seems to only apply to Carina. Fig. 5 shows the evolution of DM density at 150 pc and virial mass for the remaining subhaloes hosting the MW dSphs: UM, Draco, Sculptor, Sextans, Fornax, LeoII, LeoI and CVnI. From the figure it can be seen that there is a large variety of cross sections that reproduce the observed central densities. Draco, for instance, prefers lower values of σ/m χ , ranging between σ/m χ = 23 cm 2 g −1 and σ/m χ = 25 cm 2 g −1 , whereas LeoII requires σ/m χ ranging between σ/m χ = 120 cm 2 g −1 and σ/m χ = 90 cm 2 g −1 .
An interesting case shown in Fig. 5 is that of LeoI, since this subhalo crossed the MW's virial radius roughly 2 Gyrs ago, it has not lost a large amount of mass from tidal interactions with the MW. Therefore the models with and without mass loss agree. and large values (∼ 10 cm 2 g −1 ) on the scale of dwarf galaxies. They showed that the velocity-dependent SIDM model resolved the TBTF problem, and it provided a particular good match to the spread in dwarf satellite central densities seen around the MW. Elbert et al. (2015) showed that even larger cross sections, e.g. σ/m χ = 50 cm 2 g −1 , also alleviate the TBTF problem and produce constant density cores of size 300-1000 pc, comparable to the half-light radii of ∼ 10 5−7 M stellar mass dwarfs. Table 2 summarises the final density profiles for the SIDM subhaloes presented in this section. From left to right, it shows the hosted dSph galaxy name, virial mass, concentration parameter, core size and the preferred range of cross section values. Note that the virial mass is calculated by integrating the density profile up to the virial radius, which in turn is estimated as the radius within which the mean density is 200 times ρ crit (z = 0). The core radius, r c , is calculated by fitting an isothermal profile (ρ(r) = ρ 0 /(r 2 c + r 2 )). Figs. 4 and 5 show that all MW dSphs need to be in gravothermal core collapse in order to explain the observational data. This result, however, strongly depends on the initial virial mass of the systems, M 200,init , which is not chosen at random, it is tuned so that the systems, in their final state, have a virial mass that reproduces the observational estimations. If we disregard this and increase the initial mass of the systems, lower values of σ/m χ will be needed to reproduce the observed central DM densities. But again, they would not be a good theoretical representation of the dSphs  because, even considering that the model's rate of mass loss is a lower limit, the systems would too massive. In Appendix C we show that changing M 200,init in 20% can change the final central DM densities of subhaloes in up to 50%. Although the final DM central density is quite sensitive to the choice of initial mass, it also largely depends on σ/m χ . Fig. 4 shows that for a constant M 200,init , changes in σ/m χ of up to 25% lead to changes of 80% in the final central DM densities.

Velocity-dependent cross section
In this section we analyse the range of cross sections that match the observed central DM densities shown in Figs. 4 and 5. To understand its dependence with the particles velocity we calculate the average collision velocities, v , of DM particles within each system's core. We define v = (4/ √ π)v (for a Maxwellian distribution), with v the average 1-D velocity dispersion of each system's core. We find that for the case of Carina for example (Figs. 3 and 4), observations favour a cross section range between 32 and 40 cm 2 g −1 , as its core reaches a stable collision velocity of ≈ 48 km/s. We find that there is a strong correlation between σ/m χ and v . This is shown in Fig. 6, that highlights the range of values for each dSph. From the figure it can be seen that for systems such as LeoII, characterised by v ≈ 21 km/s, observations favour a large cross section, whereas for Draco, which has a v ≈58 km/s, observations favour lower cross sections. We determine the range of cross sections for each dSph by analysing the models with fixed σ/m χ that are able to reproduce the DM central densities (Fig. 5), and also by considering that DM densities from the gravothermal model may differ from N-body simulations by up to a factor of 2 (Essig et al. 2019). We used this uncertainty to further extend the range of σ/m χ . Fig. 6 indicates that the range of cross sections needed to reproduce the observed DM densities are not random, instead they point towards an intrinsic velocity-dependent relation. We investigate this relation in the context of particle physics models for SIDM, where a DM particle χ of mass m χ interacts under the exchange of a light mediator φ, with the self-interactions being described by a Yukawa potential with r the separation between DM particles, α χ the analog of the fine-structure constant in the dark sector, and m φ the mediator mass. There is no analytical form for the differential scattering cross section due to a Yukawa potential, but by using the Born-approximation (valid when the scattering potential can be treated as a small perturbation), the analytical form that approximates the true differential cross section results (see where v is the relative velocity between interacting DM particles, w ≡ m φ /m χ is a characteristic velocity, σ 0 = 4πα 2 χ m 2 χ /m 4 φ is the amplitude of the cross section, and θ is the scattering angle in the frame of the centre of mass. From eq. (14) we calculate the total cross section by integrating over the solid angle and obtain with w = 30 m φ 10 MeV 10 GeV m χ km s −1 .
We fit eq. (15) to the dSphs values in order to determine the values of DM mass and mediator mass that reproduce the relation. We assume α χ = 0.01 and find that the relation is best fitted by a DM particle mass of m χ = 0.648 ± 0.154 GeV and a mediator mass of m φ = 0.636 ± 0.055 MeV. This best-fit relation is shown in Fig. 6 in solid line, the grey region highlights the uncertainty by propagating the errors of m χ and m φ . Fig. 7 extends the cross section-velocity plane to include MW-(v ∼ 150 − 300 km/s) and cluster-size (v ∼ 1000 − 5000 km/s) haloes. In the figure, the values for the MW dSph galaxies are shown in blue symbols and the best-fit relation in solid line. The figure shows that the extrapolation of the best-fit cross section-velocity relation to large scales lies in very good agreement with current observational constraints. For characteristic DM velocities of MW-size haloes, the upper limit of σ/m χ < 10 cm 2 g −1 is set by subhalo evaporation, where collisions between DM particles within subhaloes and in the host are frequent enough to unbind material from the halo. In this case, energy transfer is determined by the relative velocity of the colliding particles, which is of the order of the orbital velocity, therefore the mass loss in subhaloes is enhanced and the subhalo abundance is depleted relative to the CDM case, particularly in the central regions (e.g. Vogelsberger et al. 2012;Rocha et al. 2013;Zavala et al. 2013). Note that the large values of σ/m χ on dwarf-scales does not translate into evaporation of substructure, this is because the relative velocity between the subhalo and its host halo is set by the velocity dispersion of the latter, for which the cross section is suppressed.
The lower limit of σ/m χ > 1 cm 2 g −1 has been imposed to solve the cusp-core and TBTF problem, otherwise dwarfs and low surface brightness galaxies in SIDM models end up too dense and do not produce large enough core radii (e.g. Davé et al. 2001;Kaplinghat et al. 2016). For characteristic DM velocities of cluster-size haloes, the upper limits in the figure correspond to strong lensing measurements of the ellipticity and central density of clusters , measurements of the offset between the DM and galaxy centre (e.g. Kahlhoefer et al. 2015;Harvey et al. 2015;Wittman et al. 2018;Harvey et al. 2019) and measurements of the mass-to-light ratio of the Bullet Cluster (Randall et al. 2008).

Comparison with previous works
Recent works have also investigated the evolution of MW dSph galaxies aiming to constrain the SIDM cross section. Read et al. (2018) and Valli & Yu (2018) considered a subhalo's inner region (limited by a radius r χ ), within which the average scattering rate per particle times the halo age (t age ) is equal to unity. At r = r χ , therefore, the cross section can be constrained from the relation where ρ(r) is the density profile and v(r) the velocity dispersion.
Valli & Yu (2018) analysed the stellar kinematic dataset for the MW dSphs, applying the standard Jeans analysis, to infer the stellar and DM density and velocity dispersion. They assumed t age to be flatly distributed in the range 8-12 Gyr and found that UM, Draco, LeoI and LeoII probed cross sections ∼ 0.1 − 1 cm 2 g −1 , whereas Sextans and Fornax had cross section values that peaked around ∼ 20 − 40 cm 2 g −1 .
Sculptor and Carina probed intermediate cross section values of ∼ 2 − 7 cm 2 g −1 . Read et al. (2018) focused on Draco and first calibrated the parameters r χ and t age using Vogelsberger et al. (2012) SIDM cosmological zoom simulations of MW-mass haloes. They found that Draco's high central density gives an upper bound on the SIDM cross section of σ/m χ < 0.57 cm 2 g −1 .
The model given by eq. (16) does not include the effects of core collapse, we test it by applying eq. (16) to our simulated subhaloes assuming t age in the range 8-12 Gyr, and never recover a σ/m χ larger than 1 cm 2 g −1 . In addition, the model is only valid if the SIDM subhalo density profile follows a cored profile. Observational studies, however, have reported that only Fornax exhibits a prominent core (e.g. Jardel & Gebhardt 2012;Pascale et al. 2018), Draco (e.g. Read et al. 2018) and the remaining MW dSphs are better described by a cuspy profile (e.g. Read et al. 2019 and references therein). Interestingly, for Fornax, Valli & Yu (2018) reported a σ/m χ ≈ 40 cm 2 g −1 , in agreement with our results

Challenges
An important caveat of the gravothermal collapse model is that the parameter C cannot be derived from first principles, instead it needs to be calibrated using N-body simulations. Previous works have done it (Balberg et al. 2002;Koda & Shapiro 2011;Essig et al. 2019;Nishikawa et al. 2020), but have not reach to a consensus of its value, other than it ranges between 0.45 and 0.75. In this section we investigate how changing C from 0.75 (assumed so far) to 0.45 or 0.6 (and b = 0.003 as suggested by Nishikawa et al. 2020) impact on our results. For all models, a cross section of σ/m χ = 40 cm 2 g −1 yields close agreement with the estimations from Kaplinghat et al. (2019) (shown as grey and black symbols). For a fixed cross section and b = 1.38, a larger C accelerates core collapse. In this manner the figure indicates that the range of cross sections (presented in Section 3), derived with a model that assumes C = 0.75, should be taken as a lower limit. A model with b = 1.38 and C = 0.45 requires a factor of ≈1.5 larger cross sections to reproduce the observed central DM densities and virial masses.
The model with C = 0.75 and b = 0.003 largely differs from the rest. It does not lower the central DM density as much during the core expansion phase, and the rate at which the central DM density increases during the core collapse phase seems to agree with the model with b = 1.38 and C = 0.45. For a fixed cross section and C = 0.75, a lower b yields lower central DM densities and a 'slower' core collapse in comparison to the C = 0.75 and b = 1.38 model. Therefore, according to this model, the range of cross sections presented in Section 3 should also be taken as lower limits, since assuming b = 0.003 would result in larger values of σ/m χ being able to reproduce the observations.
Another important caveat to consider is the assumption of mass conservation in the gravothermal collapse model (eq. 1), that is inconsistent with the fact that the subhalo loses mass every 250 Myrs. To solve this we impose that after the subhalo loses mass, it returns to hydrostatic equilibrium, readjusting its density and pressure within each shell. After that the gravothermal model is called back, and the new evolution of the truncated density profile begins. Further details of the numerical implementation of the gravothermal model are presented in Appendix A.

Impact of initial conditions: NFW profile
The SIDM halo model presented in this work evolves the subhaloes hosting the most massive MW dSphs for 10.2 Gyrs, starting when the Universe is 3.5 Gyr old (redshift 1.87) at a point when the subhaloes' initial density (ρ init ) is assumed to follow the NFW profile. This, however, may not be a good assumption. Harvey et al. (2018) analysed the evolution of 19 low-mass dwarf spheroidal galaxies using a SIDM numerical simulation with σ/m χ = 10 cm 2 g −1 , finding that the dwarf galaxies were already forming a core within the first 2 Gyrs of cosmic time.
It is possible that initializing subhaloes with a cored profile will induce an earlier gravothermal collapse, that will yield lower estimates for the cross sections with respect to the ones reported in Section 3.3. In this section we investigate if this occurs by analysing the evolution of Carina and Leo II, that were modelled with σ/m χ = 40 cm 2 g −1 and σ/m χ = 120 cm 2 g −1 , respectively, and initialized with three different density profiles.
The middle panels of Fig. 9 show the initial density profile for the models of Carina (top) and Leo II (bottom). In  . DM density at 150 pc, ρ 150 , as a function of lookback time (left), initial density profile (middle) and initial 1-D velocity dispersion profile (right) for Carina (top panels) and Leo II (bottom panels). Carina was initialized with a cross section of σ/m χ = 40 cm 2 g −1 , whereas Leo II has a cross section of σ/m χ = 120 cm 2 g −1 . The coloured lines correspond to the subhalo model initialised with a NFW density profile (red dashed lines), cored profile with a hot extended halo (yellow solid lines) and cored profile with a cold extended halo (blue solid lines). The symbols show the values of ρ 150 taken from Kaplinghat et al. (2019), who assumed an isothermal cored profile (grey symbol) as well as NFW (black symbol). The figure shows that initializing the subhaloes' density profile with either an NFW profile or cored profile (with hot extended halo) does not impact the results presented in Section 3, however initializing the subhaloes' density with a cored profile surrounded by a hot extended halo slightly alters the evolution of subhaloes.
the panels the red dashed lines correspond to the NFW density profile, characterized by a cuspy and cold inner region surrounded by a hot extended halo. The solid lines correspond to two different core profiles, one where the inner region is hot and it is surrounded by a hot extended halo (yellow solid lines), and the other where the extended halo is cold (blue solid lines). This can be seen in the right panels that show the 1-D velocity dispersion as a function of radius. The difference between these profiles is that they correspond to different evolutionary stages of SIDM haloes. The hot core + hot extended halo corresponds to a SIDM halo that has a hot inner region due to DM-DM particle interactions, but has a hotter periphery due to dynamical heating, induced by mergers and large DM accretion (Colín et al. 2002). On the contrary, the hot core + cold extended halo corresponds to a SIDM halo that has been in isolation.
The left panels of Fig. 9 show the evolution of the DM density at 150 pc, ρ 150 , for Carina (top) and Leo II (bottom). It can be seen from the top-left panel that initializing Carina with either an NFW profile or cored profile does not impact the results presented in Section 3. Differently, the bottom-right panel shows that initializing Leo II with a core profile of hot inner region and cold extended halo, changes the evolution ρ 150 , reaching lower values at present time. We find, however, that the 'hot core + hot extended halo' profile better represents the initial stages of subhaloes hosting the MW dSphs. When the Universe is 3.5 Gyr old, is very likely that the dSph subhaloes have undergone recent mergers or had large rates of mass accretion, since at that point none of them have yet crossed the MW's virial radius.
Another important assumption of the SIDM halo model is the NFW profile for the MW halo. We find, however, that this is a good approximation for our models. At MW halo scale, σ/m χ ∼ 1 − 5 cm 2 g −1 , according to the σ/m χ -velocity relation shown in Fig. 7. For these cross sections, Robles et al. (2019) has shown that the SIDM MW halo embedded in a baryonic potential not only exhibits a remarkably similar density profile to that of a CDM MW-like halo, but it also has no discernible core.

Impact of initial conditions: halo concentration
It has previously been shown that the core collapse timescale, t c , is Essig et al. 2019;Nishikawa et al. 2020) where C is the free parameter described in Subsection 4.2. Eq. (17) indicates that for fixed σ/m χ , c 200,init and M 200,init , the larger parameter C accelerates core collapse as shown in Fig. 8. It also indicates that a larger cross section and/or virial mass accelerates core collapse. Section 3.2 comments that in our model, the initial virial mass of the systems is constrained by the tidal evolution model and observations. This leaves us wondering about the impact of the concentration parameter, c 200,init , on the SIDM subhalo evolution. Sameie et al. (2020) explored the tidal evolution of SIDM subhaloes in the MW's tides. They produced N-body simulations of dwarf spheroidal galaxies orbiting around the MW, modelled with a static potential that included both the disk and bulge components. They found that a constant cross section of σ/m χ = 3 cm 2 g −1 can reproduce the observed DM density profiles of the MW dSphs Draco and Fornax. However, this was only possible if subhaloes were initialised with a large concentration parameter, such as 29.5 for Draco. They showed that if the concentration parameter was lowered to 22.9, not even the model of σ/m χ = 10 cm 2 g −1 was able to reproduce the large DM densities of Draco, and a model with higher cross section was needed (see also Kahlhoefer et al. 2019).
We test the impact of c 200,init by running the models of Draco and Fornax initialized with a concentration of c 200,init = 15 and σ/m χ = 3 cm 2 g −1 . We find that both models reproduce the observed central DM densities, in agree-ment with Sameie et al. (2020). This result is presented in Appendix B.
We believe, however, that setting such large initial concentrations is not well justified. In the starting point of our model, subhaloes represent typical z = 1.87 low-mass subhaloes in the field, whose density profiles follow the median z = 1.87 concentration-mass relation. Correa et al. (2015c) showed that at high redshift (z > 1), the halo has large rates of accretion, with its mass history mainly characterized by an exponential growth. During this time, the scale radius increases simultaneously with the virial radius, hence the concentration hardly grows. At low redshift (z < 1), there is a drop in the accretion and merger rates of small haloes, and the halo mass increases due to the evolution of the reference density used in the spherical overdensity definition of the halo (ρ crit in this case, also referred as pseudo-evolution phase). This leads subhaloes to have roughly constant scale radius but growing virial radius, which produces the rapid growth of concentrations. In our initial conditions, subhaloes have not reached the pseudo-evolution phase, therefore their concentrations should be set to low values.

Uncertainty in orbital parameters
This work analyses whether the anti-correlation between the central DM density of MW dSphs, ρ 150pc , and their pericenter passages, r P , is the result of SIDM effects. The errors in the orbital parameters reported by the Gaia collaboration (Fritz et al. 2018;Helmi et al. 2018;Brown et al. 2018), however, can weaken the anticorrelation. The parameter that produces the largest uncertainties in the pericenter distances is the tangential velocity, v T , whose errors are around 20% for Carina, UM, Draco, Fornax, Sculptor and Sextans, but increase to 60 − 100% for CVnI, LeoI and LeoII.
Like the ρ 150pc − r P anti-correlation, the cross sectionvelocity relation obtained in Section 3.3 can be affected by the uncertainties of v T . In this section we investigate this further by analysing how the central DM densities of the dSphs Carina, Fornax and LeoII depend on 20% (50% for LeoII) changes in the tangential velocity. Fig. 10 shows ρ 150 (left panels), virial mass, M 200 (middle panels) and galactocentric distance (right panels) as a function of lookback time, for Carina (top panels), Fornax (middle panels) and LeoII (bottom panels). The coloured lines correspond to the subhalo models initialised with the same cross section, initial mass, galactocentric distance and radial velocity, but different tangential velocities. The default values of v T for Carina, Fornax and LeoII are 163, 132 and 74 km s −1 , respectively, these are shown in orange solid lines. Blue solid lines indicate the evolution of the models with 20% (50%) larger v T than the default models, whereas green solid lines show the evolution of the models with 20% (50%) lower v T than the default models.
The top panel shows that lowering v T in 20%, decreases the galactocentric distance, increasing the rate of mass loss and accelerating the gravothermal collapse. As a result, Carina reaches a ∼ 50% higher ρ 150 at present time with respect to the default model. Similarly, increasing v T in 20%, increases the galactocentric distance, decreases the rate of mass loss, and Carina reaches a ∼ 60% lower ρ 150 .
For the evolution of Fornax, lowering v T in 20% yields an earlier infall onto the MW gravitational potential, further increasing the rate of mass loss and accelerating the gravothermal collapse. For this case Fornax reaches over an order of magnitude higher ρ 150 at present time with respect to the default model. Increasing v T , on the other hand, results in a final central DM density in close agreement with the default model. The changes in v T for LeoII are of the order of 50%. The bottom panels of the figure show that decreasing v T does not yield a large disagreement between the final DM central densities. This is because both models experience similar rates of mass loss, despite the difference in their orbits. Differently, the model with 50% larger v T , orbits around the MW only once and it does not lose as much mass from tidal interactions, therefore the gravothermal collapse is delayed, reaching an order of magnitude lower ρ 150 at present time, with respect to the default model.
We conclude that the uncertainties in the galaxies tangential velocities can change the range of cross sections that reproduce the central DM densities, altering the cross section-velocity dependence. In a coming study we will analyse the evolution of truncated SIDM density profiles from gravitational tidal stripping and ram pressure, to further improve the modelling of subhaloes hosting the local dwarf galaxies and provide robust constraints of σ/m χ on dwarf galaxy scales. Such constraints will be adjusted by the uncertainties of the orbital parameters.

CONCLUSIONS
Self-Interacting Dark Matter (SIDM) offers a promising solution to the small-scale challenges faced by the otherwiseremarkably successful Cold Dark Matter (CDM) model. However, robust constraints of the SIDM scattering cross section per unit mass, σ/m χ , on dwarf galaxy scales seem to be missing. The anti-correlation between the central DM densities of the bright Milky Way (MW) dwarf spheroidal galaxies (dSph) and their orbital pericenter distances , poses a potential signature of SIDM. In this work, we have investigated such possibility and found that there is a cross section-velocity relation that is able to explain the diverse DM profiles of MW dSph satellites, and is consistent with observational constraints on larger scales.
To model the evolution of SIDM subhaloes hosting the MW dSphs, we have applied the gravothermal fluid formalism for isolated, spherically symmetric self-gravitating SIDM haloes, and extended it to include the orbital evolution around the MW gravitational potential, along with a consistent characterization of gravitational tidal stripping. We have adopted the proper motions from the Gaia mission and used the code galpy (Bovy 2015) to integrate the orbits of the dSphs. We have also used the fitting functions from  to model the truncation of subhaloes' density from tidal mass loss. Our model has the advantage of tracking the subhalo evolution within scales smaller than 100 pc, largely expensive to resolve with N-body simulations, while easily covering a wide range of parameter space for the SIDM scattering cross section per unit mass, σ/m χ .
We have applied the model to the classical dSph galaxies, namely, Ursa Minor, Draco, Sculptor, Sextans. Fornax, Carina, LeoII, LeoI and Canes Venatici I, using the orbital parameters estimated by Fritz et al. (2018) (Fig.1). The subhaloes hosting the dSphs are modelled from an initial NFW profile, virial mass of (0.1 − 4) × 10 9 M and concentration parameter of 6 − 7 (Table 1). Their evolution begins when the Universe is 3.5 Gyr old (redshift z = 1.87), at a point when none the dSphs have crossed the MW's virial radius, and continues for 10.2 Gyrs until present time.
Along with the virial mass and concentration, subhaloes are initialised with a constant σ/m χ . We have shown that the evolution of the subhaloes density profile, as well as the rate of mass loss, largely depends on σ/m χ , so that a large σ/m χ leads to a larger rate of DM-DM collisions that produce a shallower and lower density core. We have exemplified the core expansion phase and gravothermal collapse phase with the subhalo hosting a Carina-like dSph (Fig. 3), and also showed the dependence of the central DM density of Carina with σ/m χ (Fig. 4). As expected, not only for Carina, but for the remaining dSphs, mass loss from tidal interactions leads to a faster gravothermal core collapse (Fig. 5).
We have investigated the range of σ/m χ that produces subhaloes with central DM densities in agreement with the observational estimations. There is not single range of σ/m χ able to reproduce the observed data, instead each subhalo is characterized by a specific range. Draco, for instance, prefers lower values of σ/m χ , ranging between σ/m χ = 23 cm 2 g −1 and σ/m χ = 25 cm 2 g −1 , whereas LeoII requires σ/m χ ranging between σ/m χ = 120 cm 2 g −1 and σ/m χ = 90 cm 2 g −1 .
Interestingly, the range of σ/m χ needed to reproduce the observed DM densities is not random, instead it correlates with the average collision velocity of DM particles within each subhalo's core. This result points towards an intrinsic velocity-dependent relation (Fig. 6), that can be fitted by a particle physics model for SIDM, where DM particles of mass m χ interact under the exchange of a light mediator φ, with the self-interactions being described by a Yukawa potential. By assuming α χ = 0.01 (analog of the finestructure constant in the dark sector), the relation is best fitted by a DM particle mass of m χ = 0.648 ± 0.154 GeV and a mediator mass of m φ = 0.636 ± 0.055 MeV. We show that the σ/m χ -velocity relation is a feasible velocity-dependent model for SIDM that lies in perfect agreement with current observational constraints on larger scales (Fig. 7).
We have extensively analysed the caveats of the model, such as the uncertainty of the free-parameters b and C from the gravothermal modelling, the density truncation time, the MW mass, errors in the orbital parameters, impact of initial subhalo mass, M 200,init , and concentration, c 200,init , as well as different initializations of the density profile. We have shown that our key results are robust to different initializations of the density profile (Fig. 9), that can either be a NFW profile or core profile. Changes of the parameters C and b, the density truncation time and MW mass, modify the ranges of σ/m χ that reproduce the observational estimates, by either increasing or decreasing the overall normalization in up to a factor of 1.5 for C and b (Fig. 8), in up to a factor of 1.13 for the truncation time (Fig. E1), or in up to a factor of 1.2 for the MW mass (Fig. D1). These changes, nevertheless, do not alter the shape of the cross section-velocity relation, only its normalization.
The initial subhalo mass and concentration can also modify the ranges of σ/m χ that reproduce the observations. Models with lower M 200,init reach lower central DM densities at present time (Fig. C1), therefore higher cross sections would be needed to reproduce the observations. We have argued, however, that our choices of M 200,init are not random, they are tuned so that the systems, in their final state, have a virial mass that reproduces the observational estimations. The concentration parameter controls the corecollapse time-scale, which is shorter for low concentration systems. We have shown that initializing the models with higher concentrations require much lower σ/m χ to reproduce the observed central DM densities (Fig. B1). However, we believe that such large initial concentrations are not well justified. At high redshift (z > 1), the halo has large rates of accretion, with its mass history mainly characterized by an exponential growth. During this time, the scale radius increases simultaneously with the virial radius, hence the concentration hardly grows (Correa et al. 2015b).
Finally, we have addressed the impact in the uncertainty of the orbital parameters estimated by the Gaia Collaboration (Fritz et al. 2018;Brown et al. 2018;Helmi et al. 2018). While our results are robust to the small uncertainties in the galactocentric distances and radial velocities, they are not to the larger uncertainties of the galaxies tangential velocities (Fig. 10), that can potentially weaken the cross-section velocity correlation, as well as the anti-correlation of central DM densities with pericenter passage reported by Kaplinghat et al. (2019).
In this paper we have made a first assessment of the viability of a velocity-dependent SIDM model able to explain a specific observable, the anti-correlation between the central DM densities of the bright MW dSph and their orbital pericenter distances. We have found that there is such model, that explains the diverse DM profiles of MW dSph satellites, is consistent with observational constraints on larger scales and predicts that the dSphs are undergoing gravothermal collapse. However more evidence will be gathered in a coming study, to further support or exclude such scenario. We will also assess the impact of baryons, as well as the evolution of truncated SIDM density profiles from gravitational tidal stripping and ram pressure, to further improve the modelling of subhaloes hosting the local dwarf galaxies and provide robust constraints of σ/m χ on dwarf galaxy scales. Such constraints will be adjusted by the uncertainties in the orbital parameters. . The coloured lines correspond to the subhaloes initialised with the same cross section (σ/m χ = 120 cm 2 g −1 for LeoII, 24 cm 2 g −1 for Draco and 40 cm 2 g −1 for Carina), but different initial virial masses. The symbols show the values of ρ 150 (and M 200 ) taken from , who assumed an isothermal cored profile (grey symbol) as well as NFW (black symbol). The figure shows how the initial mass of the models impact on the final DM density profile.
In this section we show that changing the initial virial masses of systems such as LeoII, Draco and Carina in 20%, changes the final central DM densities of these subhaloes in up to 50%. Fig. C1 shows the DM density at 150 pc, ρ 150 , as a function of lookback time for the models of LeoII (topleft panel), Draco (top-right panel) and Carina (bottom-left panel). The coloured lines correspond to the subhaloes initialised with the same cross section (σ/m χ = 120 cm 2 g −1 for LeoII, 24 cm 2 g −1 for Draco and 40 cm 2 g −1 for Carina), but different initial virial masses.
The top panels show that the models of LeoII and Draco lower their final DM densities from 3 × 10 8 M kpc −3 to 1.5 × 10 8 M kpc −3 , and from 4×10 8 M kpc −3 to 2×10 8 M kpc −3 , respectively. The bottom panels shows that for Carina, an initial virial mass of 10 9.3 M leads to a final DM density of 2 × 10 8 M kpc −3 , lowering the initial mass by 20% results in a DM density of 9 × 10 7 M kpc −3 , and further lowering the initial mass by 50% results in a DM density of 4 × 10 7 M kpc −3 .

APPENDIX D: IMPACT OF MILKY WAY MASS
The mass of the Milky Way (MW) can impact the results presented in Section 3 by altering the orbital evolution of the dwarfs (shown in Fig. 1), as well as the rate of mass loss by increasing/decreasing the tidal radius above which mass is tidally stripped. In this Section we investigate how changing the MW mass from 10 12 M (default value adopted through-  Kaplinghat et al. (2019), who assumed an isothermal cored profile (grey symbol) as well as NFW (black symbol). The figure shows that assuming a MW mass of 1.6 × 10 12 M (0.8 × 10 12 M ) requires a factor of ≈ 1.2 lower (higher) σ/m χ to reproduce the observed central DM densities. out the work) to 0.8 × 10 12 M or 1.6 × 10 12 M , changes the final estimates of the DM density of dwarf subhaloes undergoing gravothermal collapse. Fig. D1 shows Carina's DM density at 150 pc, ρ 150 , as a function of lookback time. The various subhalo models shown in the figure were all initialised with a cross section of σ/m χ = 34 cm 2 g −1 and a MW mass of: 0.8 × 10 12 M (shown in blue solid line), 10 12 M (red dashed line) and 1.6 × 10 12 M (yellow solid line). The figure shows that a high MW mass, induces a higher rate of mass loss from gravitational tidal stripping, and therefore further accelerates the subhaloes' gravothermal collapse, resulting in subhaloes with higher central DM densities. On the contrary, a lower MW mass produces subhaloes with lower central DM densities.
For Carina, but this also applies to the other dwarf models, increasing (decreasing) the MW mass in a factor of 1.6 (1.2), increases (decreases) the central DM density in a factor of 2 (1.2). This implies that assuming a MW mass of 1.6 × 10 12 M , requires a factor of up to ≈ 1.2 lower cross sections to reproduce the observed central DM densities, or alternatively assuming a MW mass of 0.8×10 12 M , requires a factor of ≈ 1.2 larger cross sections.

APPENDIX E: IMPACT OF TRUNCATION TIME
The truncation time, t trunc , is a free parameter that regulates the frequency over which the subhaloes' density profile is truncated. In Section 2.3 we showed that the largest changes of the density profile occurs in the outskirt of the halo, beyond the virial radius, such change in the profile is, nevertheless, important because it lowers the pressure of the extended halo, it gets therefore colder as it readjusts to hydrostatic equilibrium. Section 3.2 shows that this favours the conditions for gravothermal core collapse.
In this section we investigate how changing t trunc from 250 Myr (default value) to 350 or 150 Myr, alters the final central DM density of subhaloes. Fig. E1 shows the Carina model initialised with a cross section of σ/m χ = 34 cm 2 g −1 .
The top panel shows the DM density at 150 pc, ρ 150 , as a function of lookback time, whereas the bottom panel shows the evolution of the virial mass, M 200 . The coloured lines correspond to the subhalo model initialised with the same cross section and initial mass, but different values for the truncation time as indicated in the legends. It can be seen from the figure that a more frequent truncation accelerates gravothermal core collapse and subhaloes reach higher ρ 150 at present time. Conversely, a less frequent truncation of the density profile, decelerates gravothermal core collapse and subhaloes reach lower ρ 150 .
For the specific cross section of σ/m χ = 34 cm 2 g −1 , the models with different t trunc reach ρ 150 and M 200 at present time in agreement with the observational estimations (within the uncertainties). Lowering t trunc to 150 Myr, yields 50% higher central DM densities, while increasing t trunc to 350 Myr, yields 30% lower central DM densities. This implies that changing t trunc from 250 to 150 Myr (or 350 Myr), increases (decreases) the cross sections range that reproduces the observed central DM densities in up to a factor of 1.12 (1.13). Therefore the main effect of t trunc in our results is to increase or decrease the normalization of the cross section-velocity relation, but it does not alters the shape of the relation. The coloured lines correspond to the subhalo model initialised with the same cross section (of σ/m χ = 34 cm 2 g −1 ) and initial mass, but different values for the truncation time, t trunc , that determines the frequency over which subhaloes' density profile is truncated due to mass loss. The symbols show the values of ρ 150 (and M 200 ) taken from Kaplinghat et al. (2019), who assumed an isothermal cored profile (grey symbol) as well as NFW (black symbol). The figure shows that a more frequent truncation, accelerates gravothermal core collapse and subhaloes reach higher ρ 150 at present time. Conversely, a less frequent truncation of the density profile, decelerates gravothermal core collapse and subhaloes reach lower ρ 150 .