Constraints on dark matter-neutrino scattering from the Milky-Way satellites and subhalo modeling for dark acoustic oscillations

The elastic scattering between dark matter (DM) and radiation can potentially explain small-scale observations that the cold dark matter faces as a challenge, as damping density fluctuations via dark acoustic oscillations in the early universe erases small-scale structure. We study a semi-analytical subhalo model for interacting dark matter with radiation, based on the extended Press-Schechter formalism and subhalos' tidal evolution prescription. We also test the elastic scattering between DM and neutrinos using observations of Milky-Way satellites from the Dark Energy Survey and PanSTARRS1. We conservatively impose strong constraints on the DM-neutrino scattering cross section of σ DM–ν,n ∝ En ν (n = 0,2,4) at 95% confidence level (CL), σ DM–ν,0 < 10-32 cm2 (m DM/ GeV), σ DM–ν,2 < 10-43 cm2 (m DM/ GeV)(Eν /Eν 0)2 and σ DM–ν,4 < 10-54 cm2 (m DM /GeV)(Eν /Eν 0)4, where Eν 0 is the neutrino energy and Eν 0 is the average momentum of relic cosmic neutrinos today, Eν 0 ≃ 6.1 K. By imposing a satellite forming condition, we obtain the strongest upper bounds on the DM-neutrino cross section at 95% CL, σ DM–ν,0 < 4 × 10-34 cm2 (m DM/ GeV), σ DM–ν,2 < 10-46 cm2 (m DM/ GeV)(Eν /Eν 0)2 and σ DM–ν,4 < 7 × 10-59 cm2 (m DM/GeV)(Eν /Eν 0)4.


Introduction
Dark matter (DM) is the most dominant but elusive matter in the universe.The DM mass and interactions beyond gravity still remain a mystery.The standard collisionless cold dark matter (CDM) cosmological model is that DM interacts only gravitationally with other particles.While this model explains the cosmic microwave background (CMB) and large-scale structure (LSS) data [1,2], it has challenges in explaining observations in the small-scale structure [3] such as the missing satellite [4,5], core-cusp [6][7][8][9] and too-big-to-fail [10,11] problems.Galaxy formation physics may solve these problems [12][13][14][15], but observations on small scales still provide a chance to explore DM scenarios beyond the CDM cosmological model.
In this work, we study the constraints on the elastic scattering cross sections of DM and neutrinos using the latest observational data of Milky-Way satellite galaxies from the Dark Energy Survey (DES) and PanSTARRS1 (PS1) [72,73]. 1 To quickly derive theoretical prediction of the number of Milky-Way satellites, we study and use a modified version of a semi-analytical subhalo model based on the extended Press-Schechter (EPS) formalism [74][75][76][77] and subhalos' tidal evolution prescription, as developed in refs.[78,79].The primordial density fluctuations collapse gravitationally, forming DM halos.Halos grow by accreting smaller halos, and the accreted halos are called subhalos.After the accretion of subhalos onto a host halo, they lose mass through gravitational tidal stripping, changing their internal structure [80][81][82].The model developed in refs.[78,79] describes the complementary semi-analytical evolution of a host halo and subhalos, including tidal stripping, with analytical expressions for the hierarchical assembly of DM halos provided by the EPS formalism.In this study, we extend it to interacting dark matter models by modifying the EPS formalism, in particular the window function and the critical overdensity.To obtain the expected number of subhalos, we use a modified version of publicly available SASHIMI codes. 2 3  We compare the computed number of Milky-Way satellite galaxies with the observed number of satellite galaxies to constrain the cross section of DM and neutrinos.We use observational data of 270 Milky-Way satellite galaxies from DES and PS1 after completeness correction [72] as well as a subset of 94 satellite galaxies that contain kinematics data.We consider the mass of the Milky-Way halo in the range of (0.6-2.0) × 10 12 M ⊙ , which is determined in refs.[83][84][85][86] based on the latest Gaia data.For our conservative constraints, we assume that all the subhalos host satellite galaxies.We also derive stronger limits, imposing a galaxy forming condition in subhalos that reduces the expected number of satellite galaxies.
We find the strongest constraints on the DM-neutrino scattering cross section of σ DM-ν,n ∝ E n ν (n = 0, 2, 4) at 95% confidence level (CL) in figures 3, 4 and 5, where E ν is the neutrino energy.These results are independent of specific particle physics model.The constraints from CMB [23,[26][27][28]32] and Lyman-α forests [27,34] are also shown in these figures for comparison.To compare with astrophysical constraints from observations of supernova (SN) 1987A neutrinos [26] and high-energy neutrinos [87][88][89][90][91][92][93][94], we need to consider a more specific model rather than the cross sections as discussed in section 6.2.This paper is organized as follows.In section 2, we review DM-neutrino interactions and their effects on primordial density fluctuations.Section 3 presents our semi-analytical model of subhalo evolution for the DM-radiation interaction scenarios.In section 4, we compare the results of our model with publicly available data of N-body simulations for interacting DM model.In section 5, we show the constraints on the DM-neutrinos cross sections using the presented semi-analytical model and the population data of Milky-Way satellites.In section 6, we discuss uncertainties of our results and compare our constraints with astrophysical ones on the DM-neutrino scattering.Here we also compare our results with the results in the literature.Section 7 concludes the paper.Supplemental materials on our semi-analytical subhalo model are given in appendices A, B and C.

DM-neutrino interactions
We consider three types of the scattering cross section of DM and neutrinos, σ DM-ν .One is constant and the others are proportional to neutrino-energy squared and quartered, σ DM-ν,n ∝ E n ν ∝ a −n (n = 0, 2, 4), where a is the scale factor of the universe normalized to unity today.A constant cross section can be obtained when the DM and mediator masses are strongly degenerate (like Thomson scattering).The cross section would typically depend on neutrino-energy squared or quartered (see appendix C in ref. [87] for the cross sections in specific DM-neutrino interactions).
DM-radiation interactions result in additional collision terms in the perturbation equations for DM and radiation.The modified Euler equations for DM and neutrinos are given by4 [26,27,95] where θ ν and θ DM are the velocity divergence for the neutrinos and DM, respectively, k is the comoving wavenumber, ψ is the curvature potential, δ ν is the neutrino density fluctuation, σ ν is the neutrino anisotropic stress potential and H = ( ȧ/a) is the conformal Hubble parameter.An overdot denotes a derivative with respect to conformal time.The DM-neutrino interactions also modify the Boltzmann hierarchy for neutrinos.For brevity, we do not show the full modified Boltzmann equations for DM and neutrinos [31,32].Γ ν-DM (Γ DM-ν ) is the conformal neutrino (DM) collision rate with DM (neutrinos), where n DM = ρ DM /m DM is the DM number density and m DM is the DM mass.The unknown parameter of DM in eqs.( 3) and ( 4) is σ DM-ν /m DM , which can be written as the dimensionless quantity, Figure 1: The linear matter power spectra for CDM, the DM-neutrino interaction scenario with u 0 DM−ν,0 = 10 −7 , u 0 DM−ν,2 = 10 −20 and u 0 DM−ν,4 = 10 −33 , which are defined in eqs.( 5) and ( 6).
where σ Th = 6.65 × 10 −25 cm 2 is the Thomson cross section and u 0 DM-ν,n is a value of u DM-ν,n normalized by relic neutrino momentum today.
In figure 1, we show the effects of DM-neutrino interactions on the linear matter power spectrum.We obtain the matter power spectrum using publicly available modified versions of the Boltzmann code CLASS [96] described in refs.[31,32,97] for DM-neutrino interactions.For higher n, the decoupling time scale is more instantaneous, inducing the higher frequency of oscillations and the flatter spectrum shape.
For the interaction strength of our interest, DM decouples with neutrinos long before recombination [97].We expect that the DM-neutrino interactions affect only the initial matter power spectrum and have negligible direct effects on the late-time evolution of DM.
Note that we assume the DM energy density behaves as ρ DM = ρ 0 DM a −3 at the DM-neutrino decoupling epoch, where ρ 0 DM is the DM energy density observed today.As a result, for m DM ≲ 10 keV with u DM-ν of our interest, DM may be enough heated by neutrinos to be relativistic at the DMneutrino decoupling [97] (see also eq. ( 44) in appendix A), changing the evolution of ρ DM .We focus on the DM mass of m DM ≳ 10 keV.

Semi-analytical subhalo model for dark acoustic oscillations
To quickly and comprehensively estimate the effects of the suppressed matter power spectrum due to DM-radiation interactions on the number of subhalos via DAO, we consider a modified version of semi-analytical model for subhalo (and host halo) evolution proposed in refs.[78,79].Our semianalytical model is similar to the case of WDM particles [79], where the small scale structure is similarly suppressed due to their free streaming, but we modify the EPS formalism, in particular the window function and the critical overdensity, as described in sections 3.1 and 3.2.

Variance of density fluctuations and window function
An important ingredient in the semi-analytical model based on the EPS formalism is the linear matter power spectrum characterized by the variance, S = σ 2 at z = 0, where z is the redshift.In the EPS method, the variance has to be spherically smoothed on a mass scale M weighted by a window function W (kR), where M = M (R).W (kR) needs to be fixed by comparing with N-body simulations.The top-hat filter is a successful choice in the CDM case, while the sharp-k filter successfully reproduces the prediction with N-body simulations for the truncated power spectra such as WDM [98,99].However, both window functions fail to reproduce the halo abundance in the universe for interacting DM models [65,100] since they cannot properly account for the suppressed power spectrum on small scale.We adopt another window function, called the smooth-k filter proposed in ref. [101], where β is a free parameter that controls the small-scale slope.For β → ∞, the smooth-k filter coincides with the sharp-k filter, , where Θ is the Heaviside step function.For large k modes, the smooth-k filter has non-zero values, unlike the sharp-k filter.It successfully reproduces the halo abundance in some models of interacting DM with dark radiation (DR) [65][66][67], where the matter power spectrum is similar to figure 1.
The mass M associated with the filter scale R is less well defined.We assign the mass M as M = 4π ρ(cR) 3 /3, where ρ is the average matter density of the universe and c is a free parameter.The parameter set (c, β) must be calibrated using simulations.We adopt (c, β) = (3.7,3.5) [65,66].Similar but different choices of (c, β) for interacting DM are also proposed in refs.[67,101].

Critical overdensity for collapse
Another important ingredient in the EPS formalism is the critical overdensity δ c , above which the DM linear density spherically collapses into a halo.In the DM-radiation interaction scenario, radiation may heat DM, producing pressure for DM and delaying the collapse.This may also change the critical overdensity from the CDM case and be similar to the WDM case.However, for m DM ≳ MeV with the interaction strength of our interest, DM would be heavy enough that we adopt the CDM case to δ c (see also the detailed discussion in appendix A), where D(z) is the linear growth factor normalized to unity at z = 0.For m DM ≲ MeV, using eq.( 9) is still a conservative choice.

Evolution of subhalo structure
After the accretion of subhalos onto the host halo, they lose mass by tidal stripping in the host halo.
The average mass-loss rate of subhalos is given by [80] ṁ where m(z) and M (z) are the subhalo and host halo masses, respectively, and τ dyn (z) is the dynamical timescale [80].A and ζ are functions of M (z) and z.The mass evolution of the host halo is the same in the CDM and WDM universes [78,79] because the accretion of larger halos may mainly contribute to the evolution of the host halo.The functions A and ζ are also almost the same in the CDM and WDM universes [78,79] because the mass evolution of subhalos would mainly depend on the host halo mass.We expect that even for interacting DM models, the subhalo mass evolution is the same with the CDM and WDM cases.We adopt the analytical expression of the host halo evolution M (z) in ref. [102], where with M (z) is further generalized to obtain the host halo mass M (z i ) at redshift z i as [103] M We also adopt the fitting functions of A and ζ proposed in ref. [79], We assume the subhalo density profile follows the Navarro-Frenk-White (NFW) profile [104] with a truncation radius r t due to tidal stripping [105], We also assume the NFW density profile for (sub)halos before their accretion onto a host halo.At accretion redshift z a when the halo becomes a subhalo, we determine the virial radius r vir,a for a subhalo with mass m a as where [106] and ρ c (z) is the critical density.The scale radius r s,a at z a is determined as r s,a = r vir,a /c vir,a , given a virial concentration parameter c vir,a .Then the characteristic density ρ s,a at z a is obtained by where The suppression of power spectrum on small scales causes DM halos to form later than in the CDM case.The concentration of halos would be lower than the CDM case due to a lower background density owning to cosmic expansion.We adopt the mean values of the concentration parameter obtained in ref. [107].In ref. [107], the formula for the concentration is applied to CDM and WDM.Recently, ref. [66] confirms that this formula is in a good agreement with simulations in the DM-DR interaction case.The parameters r s and ρ s are related to the maximum circular velocity V max and the corresponding radius r max as The subhalo structure before and after tidal stripping is related to V max and r max at accretion redshift z a and any later redshift z 0 [78,108], where m 0 is the subhalo mass at z 0 after tidal stripping.We obtain r s,0 and ρ s,0 at z 0 through eqs.(20).
Then the truncation radius r t,0 is obtained by solving and eq.( 10).When r t,0 /r s,0 < 0.77, subhalos may be completely disrupted [109], but it has also been pointed out that its complete disruption may not occur [110].In this work, we take into account its complete tidal disruption.

Evolution of the number of subhalos
The number of subhalos as a function of mass can be computed as where δ(x) is the Dirac delta function, c a are the concentration parameter at accretion redshift z a and P (c a |m a , z a ) is the distribution of c a .We can calculate the number of subhalos as a function of other quantities such as r s , ρ s and V max by replacing the argument of the delta function in eq. ( 23) appropriately.For P (c a |m a , z a ), we adopt the log-normal function with the mean value c(m a , z a ) [107] and standard deviation of σ log c = 0.13 [111].
The number of subhalos accreted onto the main branch of host halo with mass m a at accretion redshift z a , d 2 N a /dm a dz a , can be computed based on the EPS formalism as [112] where M a (z a ) is the mean mass of the main branch of host halo at z a given by eq. ( 11).s a ≡ σ 2 (m a ) and S 0 ≡ σ 2 (M 0 ) are the variances of liner matter density on mass scales of m a and M 0 , respectively.δ a ≡ δ c (z a ) and δ 0 ≡ δ c (z 0 ) are the critical overdensities.F(s a , δ a |S 0 , δ 0 ; M a ) is the fraction of accreted subhalos with mass m a at z a in the mass range of the main branch of the host halo [M a − dM a , M a ] defined as where P (M a |S 0 , δ 0 ) is the distribution of the main branch of host halo M a and [...] −1 is the normalization factor.m max ≡ min[M a , M 0 /2] is the maximum mass of subhalos accreting onto the main branch and M = M max ≡ min[M a + m max , M 0 ] is the maximum mass of the main branch of host halo after the mass of the main branch has increased due to the accretion of subhalos with m a .δ M are defined at the redshift when the main branch has mass M max and S M = σ 2 (M max ).

Comparison with N-body simulations
We test our semi-analytical model with numerical simulations.Unfortunately, there is no numerical simulation of subhalo evolution for interacting DM models with neutrinos in the literature.However, there are some results of N-body simulations for subhalo population for interacting DM models with DR [70].We compare the results from our model with these numerical results.Any radiation would also affect only the initial matter power spectrum directly.
The DM-DR interactions are parametrized by the DM and mediator masses, the abundance of DR and the DM-mediator and DR-mediator couplings.In refs.[59,70], the effective theory of structure formulation (ETHOS) including DM-DR interactions and DM self-interaction is proposed 5 .In the following, we focus on ETHOS1-ETHOS3 models, where the above parameters are fixed to concrete values (see ref. [70] for details).In figure 2, we compare the cumulative number of subhalos as a function of maximum circular velocity with the numerical results for the host halo mass of M 200 = 1.6×10 12 M ⊙ .We find our model to be in very good agreement with the N-body simulations within a factor of 1.8.
We should note that there are also some results of N-body simulations for the DM-photon interaction scenarios [17,25].In ref. [17], the authors compute the number of subhalos as a function of the maximum circular velocity V max in the mass bin of (2.3-2.7)× 10 12 M ⊙ but they calculate V max by the observed stellar line-of-sight velocity dispersions, V max = √ 3σ, which may be underestimated.In ref. [25], the authors compute the number of subhalos as a function of subhalo mass in the large mass bin of (0.3-1.4) × 10 12 M ⊙ /h.Some details of precision of the simulations in refs.[17,25] are also lower than in the simulations in the case of ETHOS models [70].Thus, since the simulations in refs.[17,25] might include relatively large uncertainties, we will not compare our results with them.
Our semi-analytical model is in a very good agreement with the ETHOS models with DM-DR interactions within a factor of 1.8.In addition, our models with the smooth-k filter of different parameter values (c, β) [67,101] also predict very similar results of figure 2 as discussed in appendix B. These results indicate our model is a valuable tool to explore DM substructure in various interacting DM models quickly.Note that the parameter set of (c, β) in the smooth-k filter is determined to account for the power spectrum on small scale properly.As in figure 1, the small-scale shape of the power spectrum is flatter for higher n of radiation-energy dependence.This indicates that for better predictions of subhalo properties, the parameter set of (c, β) might slightly change, depending on the nature of interactions (interacting species, the radiation-energy dependence of the cross section, the DR abundance, etc).

Constraints on DM-neutrino scattering 5.1 Constraints on DM-neutrino scattering cross sections
We compare the expected number of subhalos in the Milky Way obtained by integrating eq. ( 23) with the observed number of Milky-Way satellite galaxies from DES and PS1 reported in ref. [72].Ref. [72] corrects the number of detected satellites in DES and PS1 by volumic weights relative to the observable volume.This completeness correction predicts the Milky Way contains 270 satellite galaxies within 300 kpc with absolute V -band magnitude of M V < 0.
We also consider the stellar kinematics data of the Milky-Way satellites.The observed satellite with the smallest line-of-sight velocity dispersion is Leo V with σ = 2.3 km s −1 [113].Assuming that the velocity dispersion is isotropic, the halo circular velocity of Leo V is V circ = √ 3σ = 4 km s −1 .We conservatively regard V circ = 4 km s −1 as the threshold of maximum circular velocity.In ref. [79], the authors found 94 satellite galaxies with V max > 4 km s −1 , which contains 82 satellite galaxies after completeness corrections and 12 satellite galaxies without the corrections.
We compute eq. ( 23) taking 500 logarithmic steps from m a = 10M ⊙ to 0.1M 200 in the subhalo mass and steps of dz = 0.1 from z a = 7 to 0.1 in redshift, where M 200 is the Milky-Way halo mass.M 200 is the mass contained within the radius whose mean density 200 times the critical density.We have confirmed that the results converge enough in these steps.We consider the Milky-Way mass in Figure 2: Cumulative number of subhalos as a function of the maximum circular velocity with the host halo mass of M 200 = 1.6 × 10 12 M ⊙ for DM-DR interaction scenarios (ETHOS1-ETHOS3 models described in ref. [70]).Solid lines denote our results and dashed lines denote the results of N-body simulations from ref. [70].
To constrain the elastic scattering between DM and neutrinos, we define the Poisson probability of obtaining the observed number of satellites N for the expected number of satellites in our models µ as P (N |µ) = µ N exp(−µ)/N !.We rule out the elastic scattering of DM and neutrinos at 95% CL when P (> N obs |µ) = N =∞ N =N obs P (N |µ) < 0.05.For the observational data of satellite galaxies of N obs = 270, µ ≤ 243 is ruled out while for the kinematics data of N obs = 94 with V max > 4 km s −1 , µ ≤ 78 is ruled out.
For our conservative constraints, we assume that all subhalos host satellite galaxies and use the kinematics data of 94 Milky-Way satellites with V max > 4 km s −1 .In figures 3, 4 and 5, we show the upper bounds on the DM-neutrino scattering cross section of σ DM-ν,n ∝ E n ν (n = 0, 2, 4) at 95% CL as a function of the Milky Way mass, respectively.We constrain the DM-neutrino cross section of σ DM-ν,0 < 10 4 , where E 0 ν ≃ 3.15T 0 ν ≃ 6.1 K. E 0 ν and T 0 ν are the average momentum and temperature of relic cosmic neutrinos today, respectively.These constraints are stronger than the constraints from CMB [23,32] and Lyman-α [34], which are also shown in figures 3 and 4.
More realistically, not all subhalos would host satellite galaxies.An ultraviolet background produced by quasars and stars heats gas in small halos, allowing it to escape from their shallow gravitational potential and suppressing star formation in small halos.Some satellite forming conditions reduce the expected number of satellite galaxies, imposing more stringent upper bounds on the scattering of DM and neutrinos.We adopt a threshold halo mass at accretion of m a ≃ 10 8 M ⊙ [12,71,114], below Figure 3: Constraints on the DM-neutrino cross section as σ DM-ν,0 = const at 95 CL as a function of the Milky-Way mass considering the kinematics data of 94 Milky-Way satellites with V max > 4 km/s (blue) as well as the data of 270 Milky-Way satellites imposing the satellite forming condition of m a > 10 8 M ⊙ (yellow).The constraints from CMB (solid line) [23,32] and Lyman-α (dashed line) [34] are shown for comparison.which we assume that subhalos do not host satellite galaxies.In this case we can impose the strongest upper bounds of σ DM-ν,0 < 4 × 10 −34 cm 2 (m DM /GeV), σ DM-ν,2 < 10 −46 cm 2 (m DM /GeV)(E ν /E 0 ν ) 2 and σ DM-ν,4 < 7 × 10 −59 cm 2 (m DM /GeV)(E ν /E 0 ν ) 4 .

Constraints on specific DM-neutrino interactions
Following appendix C in ref. [87], we map our strongest constraints to couplings of specific DMneutrino interactions. 6We denote the DM and mediator masses as m DM and m ϕ , respectively.We consider the case of m ϕ ≳ m DM ≫ E ν .
Scalar DM, scalar mediator: Figure 4: Constraints on the DM-neutrino cross section as σ DM-ν,2 ∝ E 2 ν at 95 CL. σ 0 DM-ν,2 is the cross section normalized by relic neutrino momentum today, The other details are the same as figure 3. Dirac fermion DM, scalar/vector mediator: Scalar DM, fermion mediator: Scalar DM, fermion mediator with m DM = m ϕ : 6 Discussions

Uncertainties of constraints
We have used the observational data of 270 Milky-Way satellites after completeness corrections [72].A companion paper [73] in ref. [72] predicts 220 ± 50 Milky-Way satellites after completeness corrections at 68% confidence level.Adopting 170 Milky-Way satellites, we find that our constraints on the DMneutrino cross section with a satellite forming condition of m a > 10 8 M ⊙ become weaker by a factor of 3, 40, 6 × 10 2 for n = 0, 2, 4, respectively.Theoretical uncertainties of the completeness corrections are marginalized except for the Milky-Way halo mass.The uncertainty of the Milky-Way halo mass is expected to change the number of the total satellites from 170 to 280 [73], based on the Gaia data of the Milky-Way halo mass.
Our subhalo model coincides with N-body simulations for the ETHOS models with DM-DR interactions within a factor of 1.8 as in figure 2, but the result of cumulative number of satellites on small scales of V max ≲ 10 km/s by our model is in better agreement with the N-body simulations than a factor of 1.8.This factor might also induce the same order of uncertainties as the above uncertainties.

Comparison with results from astrophysical neutrinos
Observations of astrophysical neutrinos propagating in DM also constrain the scattering cross sections.Due to the different energy scales of neutrinos, there is no simple mapping to the cross section between cosmological constraints and astrophysical constraints of the DM-neutrino cross sections.
In the following, we consider mainly the DM scattering with high-energy neutrinos [87][88][89][90][91][92][93][94].We will comment on the constraints from observations of SN 1987A neutrinos in the final paragraph of this section.The most stringent constraints from high energy neutrinos are imposed by neutrino events with E ν ∼ 10 TeV from active galaxy NGC 1068.The upper bounds of the scattering cross section from these events are roughly σ DM-ν,n ≲ 10 −30 cm 2 (m DM /GeV)(E ν /10 TeV) n [94]. 7e consider several types of scattering between DM and neutrinos via a mediator particle, following appendix C of ref. [87], and compare with our constraints derived in section 5.2.For scattering of DM with high-energy neutrinos, we consider the parameters of E ν ≫ m ϕ ≳ m DM , where m ϕ and m DM are the mediator and DM masses, respectively.For m ϕ ≫ E ν and/or m DM ≫ E ν , the cross section would be further suppressed by their masses, which is of less interest.We also leave the case of m ϕ ≫ m DM as future work for simplicity.See ref. [87] for a discussion of comparison in the general case.As a result, the constraints on the DM-neutrino scattering cross section from cosmology and high-energy neutrinos are highly complementary.
Scalar DM, scalar mediator: where we assume the log term is ∼ O (10).The constraint from high energy neutrinos is weaker than that from the Milky-Way satellite observations, eq. ( 29).
Dirac fermion DM, scalar mediator: The above constraint is weaker than that from the Milky-Way satellite observations, eq. ( 31).
The constraints are also provided by observations of neutrinos from SN 1987A.The thickness of DM through which SN 1987A neutrinos propagate is 50 kpc 8 kpc ρ(r)dr ∼ 10 25 MeV cm −2 , where ρ(r) is the DM density in the Milky Way.Then the constraints from SN 1987A observations read σ DM-ν,n ≲ 10 −25 cm 2 (m DM /MeV)(E ν /10 MeV) n [26], where the average energy of SN 1987A neutrinos is ∼ 10 MeV.In many of the parameter regions of the above examples, the constraints on the couplings from SN 1987A neutrinos are comparable or weaker than those from the Milky-Way satellite or highenergy neutrino observations.

Comparison with results from the Milky-Way satellites in the literature
Strictly speaking, there is no previous constraint on the DM-neutrino scattering from the Milky-Way satellite observations, but refs.[17,21] provided the constraints on the DM-photon scattering cross section, which might be similar to the limits on the DM-neutrino cross section.We compare our methods and results with refs.[17,21].
In ref. [17], using the N-body simulations of subhalos and the kinematics data of 40 Milky-Way satellites with V max ≥ 8km/s, the authors constrain the DM-photon constant cross section of σ DM-γ < Using this semi-analytical model, we have provided stringent upper bounds on the elastic scattering of DM and neutrinos from the latest population data of the Milky-Way satellite galaxies by DES and PS1 surveys.
The modifications to our semi-analytical model compared with the model for WDM [79] are only using the smooth-k filter and changing the critical overdensity.Our model is in a reasonable agreement with the publicly available results of N-body simulations for three DM-DR interaction scenarios.This agreement indicates that our model is a simple and powerful tool to study subhalo properties in presence of DAO induced by DM-radiation interactions quickly and comprehensively.
By comparing the data of the Milky-Way satellites from DES and PS1 surveys with the results from our model, we have obtained the constraints on the DM-neutrino scattering cross sec- . By imposing a satellite forming condition, we obtain the upper bounds on the DM-neutrino cross section at 95% CL, σ DM-ν,0 < 4 × 10 −34 cm 2 (m DM /GeV), σ DM-ν,2 < 10 −46 cm 2 (m DM /GeV)(E ν /E 0 ν ) 2 and σ DM-ν,4 < 7 × 10 −59 cm 2 (m DM /GeV)(E ν /E 0 ν ) 4 .The results obtained in this analysis are independent of specific particle physics models.In section 5.2, we have mapped their constraints to couplings for specific DM-neutrino interactions.
interacting DM scenario is defined as where γ = 5/3 for DM particles and G is the gravitational constant.The Jeans mass can be written, After the matter-radiation equality of z eq ∼ 3000, the density fluctuations begin to grow as ∝ a.
In ref. [117], the authors discussed the overcritical density for WDM.Using a model for WDM, where WDM is an adiabatic gas with a gas temperature, T ∝ a −2 , they solved one-dimensional hydrodynamical simulations for spherical collapse with pressure (i.e., with the gas temperature T ).The hydrodynamical equations for the WDM model are mathematically identical to those for the DMradiation interaction scenario.Thus, we can apply the results in ref. [117] to our model.The fitting formula of the critical overdensity for the DM-radiation scenario is, based on refs.[117] and [98], where x = log[M/M IDM J (z eq )] and h(x) = (1 + exp[(x + 2.4)/0.1])−1 .For x ≫ 1, we obtain δ c,IDM (M, z) ≃ δ c,CDM (t).In section 5, we consider the number of Milky-Way satellites N sh with V max ≥ 4 km/s or m a ≥ 10 8 M ⊙ .The relation between the subhalo mass m and their maximal circular velocity V max in presence of the DM-radiation interactions is still non-trivial, but in the WDM case, which has a truncated power spectrum, V max ∼ 8 km/s corresponds to m ∼ 10 8 M ⊙ [118].In the interaction strength of our interest with m DM ≳ MeV, we typically obtain M J (z eq ) ≲ 10 6 M ⊙ .To compute N sh with V max ≥ 4 km/s or m a ≥ 10 8 M ⊙ , we expect to be able to approximate δ c,IDM (M, z) ≃ δ c,CDM (t). ( We have confirmed this approximation does not affect N sh with V max ≥ 4 km/s or m a ≥ 10 8 M ⊙ for m DM ≥ MeV and u 0 DM-ν,0 = 10 −6 .To compute the properties of smaller subhalos or lighter DM particles, we should use the fitting formula of the critical overdensity given by eq. ( 49).Eq. ( 50) is a conservative choice, where the number of subhalos is less suppressed, compared with eq. ( 49).
In figure 6, we show the number of subhalos predicted by our models as a function of the maximum circular velocity V max (left panels) and the virial subhalo mass m 0 (right panels) for ETHOS models.

C Distributions of the number of subhalos
In figure 7, we show the number of subhalos predicted by our model as a function of the maximum circular velocity V max (left panels) and the virial subhalo mass m 0 (right panels) for several examples.In this figure, we consider the constant DM-neutrino cross section for simplicity.We take u DM-ν = 10 −6 (blue lines), u DM-ν = 10 −7 (red lines) and u DM-ν = 10 −8 (yellow lines).In each panel, we assume the Milky-Way mass of M 200 = 2.0 × 10 12 M ⊙ (top panels), M 200 = 1.6 × 10 12 M ⊙ (middle panels) and M 200 = 0.8 × 10 12 M ⊙ (bottom panels).