title: Elastic and Inelastic Scattering of Cosmic-Rays on Sub-GeV Dark Matter
authors: Gang Guo Yue-Lin Sming Tsai Meng-Ru Wu Qiang Yuan
text: 
  
  
  
  gangg23@gmail.com
  Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan
  
  
  smingtsai@gate.sinica.edu.tw
  Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan
  Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210033, China
  Department of Physics, National Tsing Hua University,
  Hsinchu 300, Taiwan
  
  
  mwu@gate.sinica.edu.tw
  Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan
  Institute of Astronomy and Astrophysics, Academia Sinica, Taipei, 10617, Taiwan
  National Center for Theoretical Sciences, Physics Division, Hsinchu, 30013, Taiwan
  
  
  yuanq@pmo.ac.cn
  Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210033, China
  School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, China
  Center for High Energy Physics, Peking University, Beijing 100871, China
  
  
  
  
  
  
  
  
  
  
  
  We revisit the signatures from collisions of cosmic-rays on sub-GeV dark matter (DM) in the Milky Way. In addition to the upscattered DM component that can be probed by existing DM and neutrino experiments widely discussed, we examine the associated signals in γ-rays and neutrinos that span a wide energy range due to the inelastic scatterings. Assuming a simple vector portal DM model for illustration, 
  we compute both the upscattered DM flux by cosmic-ray protons, and the
  resulting emission of secondary γ-rays and high-energy neutrinos
  from proton excitation, hadronization, and the subsequent meson decay. We derive limits on coupling constants in the vector portal model using data from the γ-ray and high-energy 
  neutrino telescopes including Fermi, H.E.S.S. and IceCube. 
  These limits are
  compared to those obtained by considering the upscattered DM signals at the low-energy DM/neutrino detectors XENON1T/MiniBooNE and the IceCube.
  For this particular model, the limits are set predominantly by non-detection of the upscattered DM events in XENON1T, for most of the DM mass range due to the large scattering cross section at low energies.
  Nevertheless, our study demonstrates that the γ-ray and neutrino signals, traditionally considered as indirect probes for DM annihilation and decay, can also be directly used to constrain the DM–nucleon interaction in complementary to the direct search experiments.
  
   
   
  
  
  
  Elastic and Inelastic Scattering of Cosmic-Rays on sub-GeV Dark Matter
      Qiang Yuan
      December 30, 2023
  ======================================================================
  
  
  
  
  
  § INTRODUCTION
   
  
  
  The extensive searches for
  the interaction between particle dark matter (DM) and the standard model (SM) sectors are conventionally classified into three different categories: the direct detection via the scattering between the DM and SM particles (see, e.g., <cit.> for the most sensitive searches for massive DM particles, and
  for recent reviews see, e.g., <cit.>), the indirect search for signals from DM annihilation or decay (see, e.g., <cit.>), and the search of DM production in colliders (see, e.g., <cit.>).
  Although they can, in principle, be used together to test a given model describing how DM interact with the SM sectors, the involved physical processes in generating the signals are clearly distinct in different categories. Despite great efforts to hunt for DM in all these respects during the past few decades, null signal of DM has been clearly confirmed to date. 
  Usually, the search for DM in each category could be subject to different limitations. 
  For example, the indirect searches focusing on signals from DM annihilation or decay 
  cannot be used to probe asymmetric or stable DM. Moreover, the deep underground direct detection experiments lose sensitivity rapidly for sub-GeV DM with recoil energies below the detection threshold (see, e.g., <cit.>).
  
  
  New detection techniques or possibilities have been recently investigated to detect dark matter lighter than GeV (see, e.g., <cit.>). Among them, the idea of detecting a subdominant fraction of DM particles that are unavoidably upscattered by cosmic-rays (CRs) to high enough energies was recently proposed in Refs. <cit.> and studied subsequently in Ref. <cit.>. Assuming an energy-independent cross section, the subsequent scattering of the accelerated DM inside low-energy DM and neutrino detectors can be used to derive new constraints on DM-nucleon or DM-electron cross sections for light DM. Ref. <cit.> considered the case that DM can be upscattered to high-energies of above ∼PeV that can be detected by IceCube, providing a probe of DM-nucleon cross section at such high-energy (HE) scales. Relevant studies also investigated the diurnal effect of the boosted DM signals <cit.> 
  and those at germanium detectors <cit.>.
  In addition to all these model-independent studies above, Refs. <cit.> explored the CR-boosted DM signals based on simple microscopic DM models and placed limits on the coupling constants and masses related to the dark sectors. Very recently,
  an excess of O(keV) electronic recoil events has been observed in XENON1T <cit.>. The scenarios of CR-boosted DM have also been employed to explain such an anomaly <cit.>. In relevance to the role of CRs,
  similar constraints based on the distortion of the CR spectrum at ∼TeV due to scattering with DM <cit.> or the
  DM signals produced from the collision of CRs with the atmosphere <cit.> were also discussed.
  
  
  We note that most of these studies on CR-boosted DM signals in the existing literature are limited to only elastic χ–p scattering. However, for CRs colliding with DM with a center-of-mass energy ≳ GeV, the nucleons can be excited and, consequently, produce secondary γ-rays <cit.> and neutrinos in the decay sequence.
  Moreover, at higher energy scales, DM can scatter off individual quarks inside nucleons and lead to deep inelastic scatterings (DIS) that will also generate secondary neutrinos and γ-rays from the hadronization and the subsequent meson decay. These inelastic scatterings are just like the pp, pγ, or pν collisions in the SM and are bound to occur regardless of the particle nature of DM. Consequently, the co-produced “indirect” signals testable by HE γ-ray and neutrino telescopes can simultaneously probe the DM–nucleon interaction together with the direct detection of the upscattered DM by CRs without invoking any other assumptions.
   
   
  In this work, we study consistently upscattered DM by CRs, together with the co-produced secondary HE γ-ray and neutrino emissions for the first time.
  For the purpose of including the inelastic scatterings, it is necessary to adopt a specific model describing the interaction between the DM and SM particles.
  Especially, the excitation of nucleons to resonances with different spins or isospins depends explicitly on the interaction type.
  More importantly,
  as we aim to cover a wide energy scale over several orders of magnitudes in this study, relevant to low-energy elastic χ-p scattering as well as γ-ray and HE neutrino productions from inelastic scattering, a consistent treatment based on a specific model taking fully into account the energy dependence as done in Refs. <cit.> is required. This is different from the studies in Refs. <cit.> which adopted either a constant inelastic cross section for π^0 production <cit.> or a rescaled 
  cross section following the energy dependence of pp collision <cit.>. 
  In addition, we also include γ-rays from other hadrons and consider constraints due to the secondary neutrino emission.
  For illustration, we simply consider a fermionic DM which couples to baryons via a new dark vector boson. The corresponding Lagrangian is
  
      ℒ⊃χ(i∂_μγ^μ-m_χ)χ + 
      g_χχγ^μχ V_μ + ∑_f=u,d,s,... g_q_fq_f γ^μ q_f V_μ 
      +1/2 m^2_V V_μ V^μ ,
  
  where χ, V_μ, and q_f refer to the DM, the vector mediator, and quarks, respectively. 
  For simplicity, we consider only that the vector mediator couples to quarks, as only scatterings between DM and baryons are relevant for our study. 
  Furthermore, we consider a universal coupling constant between the vector mediator and quarks of different flavors, i.e., g_q_f = g_q = g_B/3. 
  This can be realized when the vector boson is a U(1)_B gauge boson with B the baryon number. We should point out that the simple U(1)_B model considered in our work suffers from gauge anomalies. To realize a UV completion of the model, new heavy particles charged under the new gauge should be introduced to cancel the anomalies. However, searches for these heavy particles with collider experiments <cit.> as well as the decay of Z and mesons <cit.> can place tight constraints on the viable parameter space. Alternative options are to consider the anomaly free models such as 
  a gauged B-L model (<cit.>) or a kinetically-mixed dark photon model <cit.>. These models with light gauged mediators are also strongly constrained by terrestrial experiments as well as astrophysical and cosmological observations (see, e.g., <cit.> and a recent review on dark photon <cit.>). Even the gauged U(1)_B model without considering the anomalies are tightly constrained by terrestrial experiments (see, e.g., <cit.>), BBN observation <cit.> as well as supernova cooling <cit.>.
  As mentioned above, in this work we simply take the U(1)_B model as an example and demonstrate that all the relevant signatures of elastic/inelastic collisions between HE CRs and DM can be consistently studied, and can be used to constrain the DM models, even though the associated constraints could be weaker than those derived from other methods. To be more specific, we derive exclusion limits on
  the couplings for sub-GeV DM by considering the upscattered DM component with XENON1T <cit.>, MiniBooNE <cit.>, IceCube <cit.>, the secondary γ-rays with data from Fermi <cit.> and H.E.S.S. <cit.>, as well as the secondary neutrinos with IceCube.
  Although we reply on this specific model to explore the inelastic effects, the same method described
  in this work can also be applied
  to other DM models.    
  
  
  
  This paper is organized as follows. In Sec. <ref>, 
  we explicitly compute the cross sections of elastic and inelastic scatterings 
  between χ and p based on the Lagrangian given in Eq. (<ref>).
  The upscattered DM flux by CRs in the Milky Way (MW)
  as well as the fluxes of the secondary HE γ-rays and neutrinos are derived.
  For demonstration, we show
  our results based on two particular choices of vector boson mass:
  (i) m_V=3 m_χ that is widely adopted in the literature (see, e.g., <cit.>) and 
  for which dark bosons decaying to DM pair is kinematically allowed; (ii) a light mediator case with m_V=1 eV. [We point out that our results in this work are actually insensitive to the value of m_V if m_V^2 is much smaller than the typical value of Q^2 involved in DM-nucleon scatterings, where Q^2 is the four-momentum transfer squared. For example, Q^2 ≳ (0.1  GeV)^2 for γ-ray production from resonance production. Here we take a representative value of m_V= 1 eV for the light mediator case so that the low-energy DM-nucleon elastic scattering is not divergent; see also relevant discussions in Sec. <ref>.] In Sec. <ref>, we evaluate the constraints from terrestrial detectors and telescopes mentioned above. We finally conclude in Sec. <ref>.
  
  
  
  
  § ELASTIC AND INELASTIC SCATTERINGS BETWEEN DM AND PROTON
     
  
  
  In this section we study the signatures of collisions between HE CRs and DM inside our Galaxy. Specifically, we consider the elastic scattering that has been commonly studied 
  in previous literature  
  
      χ + p →χ + p,
  
  and, more importantly, the inelastic scattering
  
      χ + p →χ + X →χ +  hadronic showers + γ-rays +  neutrinos.
   
  
  Since CRs span a very wide range in energy, 
  the energies of the resulting DM as well as the secondary γ-rays and neutrinos can be very high. 
  Therefore, we expect that terrestrial neutrino/DM experiments or γ-ray telescopes sensitive to different energy regions can all be used to probe/constrain the DM model. 
  
        
  
  
   §.§ Basic formalisms for DM–proton scatterings
  
   
  
  For convenience, we
  choose to describe DM–proton scattering in the rest frame of the initial proton.
  We introduce the four-momenta, p_χ=(E_χ,  p_χ), p'_χ=(E'_χ,  p'_χ), q =  p_χ-p'_χ, p_p=(m_N,0), and p_X as the incoming momentum of DM, the outgoing momentum of DM, the momentum transfer, the initial momenta of the proton, and the momentum of the outgoing hadronic final state X, respectively. 
  For the case of elastic scattering, p_X is simply the momentum of the outgoing proton. Note that quantities such as energies and three-momenta with a tilde always refer to the rest frame of the initial proton.         
  
  The differential cross section for elastic scattering is given by <cit.>
  
      dσ_ el/dQ^2 =  g_B^2g_χ^2/4πβ^21/(Q^2+m_V^2)^2[ 1 - β^2 Q^2/Q^2_ max + Q^4/8 m_N^2 E_χ^2] G^2_p(Q^2) ,
  
  where Q^2 ≡ -q^2 is the positive four-momentum transfer squared, β=| p_χ|/E_χ is the incoming velocity of DM, and G_p is the form factor of proton. 
  The maximal momentum transfer squared, Q^2_ max, is 
  
      Q^2_ max = 4(Γ^2-1)m_χ^2 m_N^2/m_χ^2 + m_N^2 + 2 m_N E_χ ,
  
  with Γ= E_χ/m_χ=(1-β^2)^-1/2. As in Refs. <cit.>, we assume that G_p takes a dipole form,
  
      G_p(Q^2) = 1/(1+Q^2/Λ_p^2)^2 ,
  
  with Λ_p ≈ 770 MeV <cit.>.
  
  
  
  Similarly to neutrino–nucleon scattering <cit.>, DM–proton inelastic scatterings can be divided into two main categories: resonance excitation of the nucleon (RES) and deep inelastic scattering off the individual quark constituents of the nucleon (DIS). With a lower threshold, resonance production contributes more significantly at E_χ≲ a few GeV, and then DIS starts to dominate at higher energies. To describe the inelastic scattering, one can introduce the following quantities (see, e.g., <cit.>):        
  
  
      ν = q · p_p/m_N = E_χ - E_χ' ,
  
  
      W^2 = p_X^2 = (p_χ+p_p-p'_χ)^2 = m_N^2+2 m_N ν-Q^2 ,
  
  
      x = Q^2/2q · p_p = Q^2/2 m_N ν ,
  
  
      y = q· p_p/p_χ· p_p=ν/E_χ=1-E'_χ/E_χ = Q^2/2m_N E_χ x ,
  
  
  where ν is the energy transfer of DM in the proton rest frame, W is the invariant mass of the hadronic final state X, and x and y are the Bjorken variables introduced for DIS. All these four variables are Lorentz invariant quantities.               
  
   
  For the excitation of a spin-1/2 resonance state, the hadronic vector current is given by <cit.>
  
      H^μ_R,1/2 = ⟨ R(p_X)| J^μ_R,1/2(Q^2) |p(p_p) ⟩  = u̅_R(p_X) Γ^μ_R,1/2(Q^2) u_p(p_p) ,
   
  where u_p and u_R are the Dirac spinors for the proton and the resonance state R, respectively, and the vertex function Γ^μ_R,1/2=V^μ_R,1/2 and V^μ_R,1/2γ^5 for resonance states with positive and negative parities, respectively. The vector part V^μ_R,1/2 can be further parametrized in terms of transition form factors as  
  
      V^μ_R,1/2(Q^2) =  F_1^R(Q^2)/4 m_N^2(Q^2γ^μ+q q^μ) + 
       F_2^R(Q^2)/2m_N i σ^μν q_ν .
  
  
  The formalism for the excitation of a spin-3/2 resonance state (dominated by the Δ-resonance) is more involved. The hadronic current is <cit.> 
  
      H^μ_R,3/2 = ⟨ R(p_X)| J^μ_R,3/2(Q^2) |p(p_p) ⟩  = ψ̅_R,α(p_X) Γ^αμ_R,3/2(Q^2) u(p_p) ,
  
  with ψ_R,α the Rarita-Schwinger spinor for a spin-3/2 state. Introducing Γ^αμ_R,3/2=V^αμ_R,3/2γ^5 and V^αμ_R,3/2 for positive and negative parities, we have 
  
      V^αμ_R,3/2 =   C_3^R(Q^2)/m_N(g^αμ q-q^αγ^μ) + C_4^R(Q^2)/m_N^2(g^αμ q · p_X - q^α p_X^μ) 
          + C_5^R(Q^2)/m_N^2(g^αμ q · p_p - q^α p_p^μ) + g^αμ C_6^R(Q^2) ,
  
  with C_3,4,5^R the transition form factors from the proton to the spin-3/2 states. 
  
  Once the transition form factors to all relevant resonance states are known, the transition matrix elements and the cross sections for resonance productions can be obtained. Taking into account the propagator of the new vector mediator, the differential cross section can be expressed as 
  
      d^2σ_ RES/dν dQ^2 = 9g_χ^2 g_q^2/(Q^2+m_V^2)^2∑_R | T_R|^2  A_R/32π m_N | p_χ|^2 ,
  
  where | T_R|^2 is obtained by contracting the hadronic tensor with that associated with DM, and A_R is the spectral function incorporating the widths of the resonance state R (see more details in Ref. <cit.>). In this work, we simply use the neutrino generator encoded in  <cit.> to calculate the differential cross section d^2σ_ RES/dν dQ^2 given in Eq. (<ref>).
  Specifically, 
  we consider 30 resonance states with masses between 1 and 2 GeV in  and use the electromagnetic (EM) form factors taken from the MAID analysis <cit.>. Note that here we have assumed that the vector form factors related to the new vector mediator are well approximated by the EM form factors.        
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  For DIS, the differential cross section is given by 
  
      d^2σ_ DIS/dν dQ^2 = 1/2m_N E_χ^2 yd^2σ_ DIS/dxdy ,
  
  where 
  
      d^2σ_ DIS/dxdy = g_χ^2 g_q^2/4π(Q^2+m_V^2)^2(1-y+y^2/2 - xy m_N/2E_χ)F_2(x, Q^2) ,
  
  and the corresponding structure function is    
  
      F_2(x,Q^2) = x ∑_f=u,d,c,s,b[f(x,Q^2) + f̅(x,Q^2)] .
  
  with f (f̅) the parton distribution functions (PDFs) of quarks (anti-quarks). In our study, we use the MSTW 2008 NNLO PDFs from <cit.>.  
  
  
  To compare the differential cross sections of the elastic and inelastic scatterings, one can further integrate out the energy transfer ν in the double differential cross sections [see Eq. (<ref>)] and obtain 
  
      dσ_ RES,DIS/dQ^2 = ∫ dσ_ RES,DIS/dν dQ^2 dν .
  
  Fig. <ref> shows the differential cross sections dσ/dQ^2 of 
  elastic scattering, resonance excitation, and DIS, respectively for cases with m_χ=10^-3 GeV, g_χ g_q = 0.03, and a light vector mediator with m_V=1 eV. The differential cross sections for the case with m_V=3m_χ are very close to the light mediator case, as Q^2 explored here is much larger than m_χ^2. For illustration, we have used E_χ=1 and 10^8 GeV as representative low and high energies, respectively. As shown in the figures, all the differential cross sections drop with increasing Q^2 at large Q^2 mainly due to the factor Q^-4 from the propagator of the mediator. Besides that, elastic scattering and resonance excitations are further suppressed by a factor of ∼Q^-8 from the hadronic form factors at Q^2 ≳ 1  GeV^2 [see Eq. (<ref>)], making DIS more important at higher Q^2.
  Related to the values of Q^2 probed at different E_χ, dσ_ DIS/dQ^2 decreases as ∼Q^-2 at Q^2 ≲ 1  GeV^2 and as ∼Q^-4 at higher Q^2. As the typical values of Q^2 increase with E_χ, DIS gives a dominant contribution at high E_χ. This can also be clearly seen in the total cross sections shown in Fig. <ref>, with σ_ RES,DIS=∫ dQ^2 dσ_ RES,DIS/dQ^2. With a finite mediator mass m_V=3 m_χ [see Fig. <ref>(b)], the total cross section is reduced when compared to that with a light mediator with m_V=1 eV [see Fig. <ref>(a)]. The reductions of the total cross sections of both resonance excitations and DIS are mainly due to the suppressed contributions of dσ_ RES,DIS/dQ^2 at Q^2 ≲ m_V^2, with Q^2_ min the order of m_χ^2. We do not show the total cross section for elastic scatterings as it requires a low-Q^2 cutoff to avoid the divergence in the case of m_V=0, and is irrelevant to our studies below.           
  
  
  
  
  
  
   §.§ DM accelerated by cosmic-rays
   
  
  
  DM can be accelerated to high velocities/energies through collisions with HE CRs. In principle, both elastic scatterings and inelastic scatterings contribute to the accelerated DM flux. The scenario involving only the 
  elastic scattering has been well studied in previous literature <cit.>. 
  In this subsection we include the contributions from the inelastic channels 
  and discuss their possible implications.           
  
  
  The differential DM flux per solid angle Ω (in units of GeV^-1 cm^-2 s^-1 sr^-1) upscattered by CRs inside the Milky Way and arriving at the Earth is given by a line-of-sight (l.o.s.) integral <cit.>
  
      ϕ_χ^ MW (T_χ,Ω) ≡d^2N_χ/dT_χ dΩ (T_χ,Ω) = ∫_ l.o.s. dℓ∫ dE_p ρ_χ(r)/m_χϕ_p(E_p) [ dσ_ el/dT_χ + dσ_ inel/dT_χ],
  
  where T_χ is the kinetic energy of the accelerated DM, ϕ_p ≡ d^2N_p/(dE_p dΩ) is the differential flux of CR protons (in units of GeV^-1 cm^-2 s^-1 sr^-1), and
  ρ_χ is the mass density of the DM halo. Following Refs. <cit.>, we have assumed that the galactic HECRs are uniformly and isotropically distributed in a cylinder with radius R= 10 kpc and half-height h=1 kpc. This is a relatively good approximation when compared to the CR distribution simulated by  <cit.> and results in only small differences for the accelerated DM flux (see later discussion and Fig. <ref>).  
  To cover a wide energy range of CRs, we adopt the same flux for CR proton as in Ref. <cit.> for our numerical study.
  The CR flux can be approximately described by a broken power law, ϕ_p ∝ E_p^-γ, with γ≈ 2.7 for E_p ≲ 10^6 GeV (below the knee), γ≈ 3 for 10^6 ≲ E_p ≲ 2× 10^8 GeV (below the 2nd knee), γ≈ 3.3 for 2× 10^8 ≲ E_p ≲ 5× 10^9 GeV (below the ankle), which then followed by a flattening of the CR flux and finally a rapid drop above 3× 10^10 GeV due to Greisen-Zatsepin-Kuzmin (GZK) cutoff <cit.>. As we aim to explore the inelastic effects, we neglect the contributions from helium and other heavy isotopes in the CRs. For ρ_χ, we take the Navarro-Frenk-White (NFW) profile <cit.> with a scale radius r_s=20 kpc, 
  normalized to the local DM density ρ_0=0.3  cm^-3 where the Sun is located at r=8.2  kpc from the galactic center (GC). 
  We assume a single-component DM scenario in this work so that the DM local density is entirely made of χ. With DM initially at rest and, thus, T_χ = Q^2/(2m_χ) in the laboratory frame, the differential cross sections in Eq. (<ref>) are given by 
  
      dσ_ el/dT_χ = 2 m_χdσ_ el/dQ^2 ,  
         dσ_ inel/dT_χ = 2m_χdσ_ inel/dQ^2 = 2m_χ( dσ_ RES/dQ^2 + dσ_ DIS/dQ^2) ,
  
  where dσ_ el/dQ^2 and dσ_ RES,DIS/dQ^2 are given in Eqs. (<ref>) and (<ref>).  
  
  
  
  
  
  
  
  Fig. <ref> shows the boosted DM fluxes from elastic and inelastic scatterings between CRs and DM inside our Galaxy, where we have introduced the solid-angle integrated flux Φ_χ^ MW (in units of GeV^-1 cm^-2 s^-1) as 
  
      Φ_χ^ MW(T_χ) = ∫ dΩ ϕ_χ^ MW(T_χ, Ω) .
  
  Note that, to better present the fluxes, we choose to show T_χ^3 Φ_χ^ MW instead of Φ_χ^ MW. We find that at T_χ≲ 0.1  GeV(m_χ/ GeV)^-1, the fluxes are dominated by elastic scattering. 
  For elastic scattering, dσ_ el/dT_χ∝  m_χ dσ_ el/dQ^2 ∝ m_χ Q^-4∝ m_χ^-1 T_χ^-2 at low Q^2.
  The minimal value of E_p contributing to the flux at T_χ for m_χ < m_N^2/E_p is E_p, min∝ (T_χ/m_χ)^1/2 [see Eq. (<ref>)]. 
  Combining these with a CR proton spectrum of E_p^-2.7, Eq. (<ref>) gives rise to a DM flux at low T_χ being proportional to m_χ^-1(E_p, min)^-1.7 dσ_ el/dT_χ∝ m_χ^-1.15 T_χ^-2.85.
  At higher T_χ or Q^2, dσ_ el/dQ^2 ∝ Q^-12∝ T_χ^-6. The
  corresponding DM fluxes from elastic scattering decrease as ∼T_χ^-6.85 as shown in Fig. <ref>. At T_χ≳ 0.1  GeV (m_χ/ GeV)^-1, resonance excitations and DIS start to contribute dominantly. Based on similar arguments, the fluxes above ∼1 GeV (m_χ/ GeV)^-1 vary as m_χ^-5.15 T_χ^-6.85 from resonance excitations similar to that from elastic scattering, and as m_χ^-2 T_χ^-3.7 from DIS considering that E_p,  min roughly scales as T_χ.                            
  
       
  
  
  
  
  The upscattered DM flux exhibits a strong angular dependence on both the right ascension () and the declination () in the equatorial coordinate system <cit.>.
  For detectors not located at the North or South poles, this can lead to daily modulations due to its dependence on  <cit.>. 
  However, when considering a long exposure time, one can integrate over the -dependence as done in Ref. <cit.>. Fig. <ref> shows the relative DM flux integrated over , ϕ^ MW_χ, = ∫_0^2π d ϕ_χ^ MW(T_χ, Ω), as a function of sin() for our fiducial scenario adopting uniformly distributed CRs and the NFW halo profile.
  
  Note that for this case, we have normalized ϕ^ MW_χ, (T_χ, ) to unity at sin()=0. 
  Since the angular dependence is determined by the spatial distribution of the DM halo and the location of the Earth, the relative angular distribution is independent of T_χ. 
  Also shown in Fig. <ref> is the relative flux using the Burkert halo profile <cit.> with the uniform CR flux. Note that we take the same normalization factor as for the fiducial case. For the Burkert profile, we use a scale radius of 9 kpc and the same local density as the NFW profile. In reality, the CR flux
  varies at different positions inside the Milky Way. 
  To estimate the associated effects, we assume that the spatial distribution of the CR proton flux is independent of the CR energy. We then use the distribution at E_p=1 TeV obtained from  (see Ref. <cit.>) for all values of E_p, and calculate the relative DM flux with the spatial-dependent CR flux and the NFW profile (see the blue dashed line). Obviously, the resulting angular distribution of the DM flux is also independent of T_χ as for the case with a constant CR flux. A comparison of these curves in Fig. <ref> shows that the uncertainty due to the chosen CR model and/or the halo profile is within a factor of ∼2. Hereafter, we stick to the fiducial case, i.e, assuming a constant CR flux and the NFW halo profile.
  
  
  The same angular distribution shown in Fig. <ref> also applies to the secondary γ-rays and neutrinos from the inelastic scatterings. 
  Although the flux peaks at the direction of the GC for cases with the NFW profile and at the direction of sin≈ -0.84 with the Burkert profile,
  the relative differences at different sin() are
  within a factor of ∼3. 
  Hereafter, we approximate the flux of the upscattered DM and those of the secondary neutrinos as isotropic for the sake of simplicity. With this approximation, we assume that the total flux Φ_χ^ MW in Eq. (<ref>) is uniformly distributed in all directions <cit.>. 
  This, on the one hand, simplifies the derivation of 
  constraints from the DM flux. 
  On the other hand, it allows us to
  directly use the existing upper limits on astrophysical or cosmogenic HE neutrinos, which are always assumed to be isotropic.
  Considering the full angular dependence could, in principle, enhance the sensitivity of our study. 
  
    
  
  
  
   §.§ Secondary γ-rays and neutrinos
    
  
  
  The fluxes of the secondary γ-rays and neutrinos summing over all flavors in units of GeV^-1 cm^-2 s^-1 sr^-1 can be expressed as 
  
      ϕ_γ,ν(E_γ,ν, Ω)= 
      ∫_ l.o.s. dℓ∫ dE_pρ_χ(r)/m_χϕ_p(E_p) 
      ∫ dν dQ^2 ×d^2σ_ inel/dν dQ^2×d𝒩_γ,ν(E_p, E_γ,ν, ν, Q^2)/dE_γ,ν ,
  
  where d𝒩_γ,ν/dE_γ,ν is the spectrum of produced γ-rays or neutrinos per pχ collision for given values of E_p, ν, and Q^2. We first use  <cit.> and  <cit.> to generate the γ-ray or neutrino yields in the rest frame of the hadronic final state X, 
  which depend only on the invariant mass W [see Eq. (<ref>)] and are isotropically distributed in three-momentum. We then boost the yields from the rest frame of X to the laboratory frame to obtain d𝒩_γ,ν/dE_γ,ν. The corresponding boost factor is simply Γ_ boost = E_X/W with E_X=E_p-T_χ=E_p-Q^2/(2m_χ) the energy of X in the laboratory frame. 
  
  
  To constrain the DM model using the secondary γ-rays, we consider the observed data from Fermi <cit.> and H.E.S.S. <cit.> with E_γ ranging from ∼GeV to ∼100 TeV. Both the RES and DIS contribute in 
  this energy range. We first rely on  to generate the π^0 yields from resonance excitations, and then calculate the resulting γ-ray yields from π^0 →γ + γ. Note that the γ-rays in the rest frame of π^0 carry an energy of m_π/2 and are isotropically distributed in three-momentum. The resulting γ-ray yields from π^0 decay can be obtained by a Lorentz boost, leading to <cit.>
  
      dn_γ(E_γ)/dE_γ = ∫_E_π, min(E_γ) dE_πdn_π(E_π)/dE_π1/√(E_π^2-m_π^2) ,
  
  where E_π, min(E_γ) satisfies the equation E_π, min(1+β_π, min)/2=E_γ with β_π, min = √(E_π, min^2-m_π^2)/E_π, min for E_γ > m_π/2. For relativistic π^0, E_π, min(E_γ) ≈ E_γ, and the resulting γ-ray takes a same spectral index as the π^0 flux.
  For HE γ-ray production in the DIS region that may go beyond the scope of , we use the  code to obtain the γ-ray yields directly. On the other hand, we consider HE neutrinos above ∼PeV that can be probed by IceCube. Since neutrino production involves three-body decays or decay sequences, we directly use the  code to calculate the total neutrino yields of all flavors from both resonance production and DIS. We further assume that neutrino fluxes are equally distributed in all flavors after oscillation.   
  
  The γ-ray fluxes in Eq. (<ref>) are further integrated over certain solid angles around the GC to be compared with the observed data at Fermi and H.E.S.S. (see  Fig. <ref>). For the secondary neutrinos, we take the isotopic approximation as mentioned in Sec. <ref> and use the averaged neutrino flux per solid angle to derive the bounds (see Fig. <ref>). 
      
  
    
  
  
  § CONSTRAINTS FROM TERRESTRIAL EXPERIMENTS AND TELESCOPE
   
  
  In this section, we derive constraints on the coupling constants of the vector portal model for sub-GeV DM based on the HE γ-ray observation (Sec. <ref>),
  the HE neutrino detection (Sec. <ref>), and 
  the expected DM signals at the low-energy DM/neutrino experiments
  
  (Sec. <ref>).
  
  
  
   §.§ Constraints from Fermi and H.E.S.S.
   
  
  To constrain the coupling constants of the vector portal model with the secondary γ-ray emission, we consider two datasets from Fermi and H.E.S.S..
  For the Fermi data, 
  we consider the Fermi GC excess at 0.5<E_γ<500 <cit.>  extracted from a circular region within 10^∘ from the GC excluding the inner region of radius 2^∘ <cit.>. 
  For the HE H.E.S.S. data, we take those with 0.2<E_γ<60, integrated over an annulus centered at Sgr A^* with inner and outer radii of 0.15^∘ and 0.45^∘, respectively, and a section of 66^∘ excluded <cit.>.
  
  
  In order to compare the produced secondary γ-ray fluxes from the inelastic collisions between the CRs and DM, we introduce the 
  J-factors for Fermi and H.E.S.S. as   
  
      J_ Fermi, H.E.S.S. = ∫_ΔΩ_ Fermi, H.E.S.S. dΩ∫_ l.o.s. dℓ ρ_χ(r)/m_χ
  
  to account for different sky coverages in the two data sets mentioned above. The corresponding J-factors for Fermi and H.E.S.S. are J_ Fermi≈ 2.16  kpc cm^-3 and J_ H.E.S.S.≈ 7.47× 10^-3  kpc cm^-3, respectively.
  Using these J-factors, the expected secondary γ-ray fluxes in units GeV^-1 cm^-2 s^-1 are  
  
      Φ_ Fermi, H.E.S.S.(E_γ)    = ∫_ΔΩ_ Fermi, H.E.S.S.dΩ ϕ_γ(E_γ, Ω)   
        
         =J_ Fermi, H.E.S.S.×∫ dE_pϕ_p(E_p) 
      ∫ dν dQ^2 ×d^2σ_ inel/dν dQ^2×d𝒩_γ(E_p, E_γ, ν, Q^2)/dE_γ .
  
  We point out that these J-factors can be taken out, because we
  assume that the CR proton flux is uniformly distributed in the Galaxy. 
  In reality, the CR flux near the GC could be a factor of a few higher than the local flux, leading to J-factors larger than those based on a constant flux by a factor of ≲2 (see also Fig. <ref>). 
  Therefore, the limits obtained here with γ-rays can be considered as conservative ones.                   
  
  
  We show the resulting γ-ray fluxes Φ_γ(E_γ)≡Φ_ Fermi(E_γ) in Fig. <ref> for different DM masses with m_V=0 [panel (a)] and m_V=3 m_χ [panel (b)]. 
  For comparison, the allowed band of the Fermi GC excess taking into account all possible systematic uncertainties is 
  shown in Fig. <ref> (see Fig. 15 of Ref. <cit.>).
  The flux Φ_ H.E.S.S.(E_γ) to be confronted with the H.E.S.S. data can be obtained from Φ_ Fermi(E_γ) by simply multiplying a scaling factor J_ H.E.S.S./J_ Fermi≈ 3.4× 10^-3.
  Thus, we show instead the rescaled
  H.E.S.S. data at 100  GeV≲ E_γ≲ 100 TeV (see Fig. 3 of Ref. <cit.>), multiplied by a factor of J_ Fermi/J_ H.E.S.S.≈ 290 so that they can be directly compared with the calculated Φ_ Fermi. 
  
  For the range of E_γ considered, resonance productions contribute more significantly than DIS. 
  As the produced γ-ray typically carries a certain fraction of the energy of the primary CR proton, the flux given in Eq. (<ref>) roughly scales as ϕ(E_γ) ∝ϕ_p(E_p) σ_ RES(E_χ) with E_χ = E_p m_χ/m_N; see also discussions below Eq. (<ref>). Considering that σ_ RES firstly increases rapidly above the threshold and then saturates at E^ sat_χ≈ 1 GeV (see Fig. <ref>),
  the γ-ray spectra above the peak shown in Fig. <ref>, can be understood quantitatively. The flux Φ_γ(E_γ) first increases rapidly and then decreases as E_γ^-2.7 above E^ peak_γ≈ 0.1  GeV (m_χ/ GeV)^-1.
  Note that E^ peak_γ is simply determined by E^ sat_χ, as the resulting E_γ from π^0 decay, on average, is ⟨ E_γ⟩= 0.5 ⟨ E_π⟩≈ 0.1E_p = 0.1 E_χ m_N/m_χ, i.e., E_γ^ peak≈ 0.1 E^ sat_χ m_N/m_χ≈ 0.1  GeV(m_χ/ GeV)^-1. The simple argument does not apply to the spectra below the peak, where the γ-ray flux decreases more slowly with decreasing E_γ than that would be indicated from the sharply declined ϕ_p σ_ RES (see Fig. <ref>).
  The reason is that the decay of π^0 with energy E_π generates a flat distribution of γ-rays with E_γ^ max,min = E_π(1±β_π)/2, where β_π is the velocity of π^0 in units of the speed of light. Therefore, the γ-ray flux below the peak receives contributions dominantly 
  from π^0 with E_π≳ E_π^0≈ 2E_γ^ peak produced with large cross sections.
  To understand the γ-ray yields below the peak, one needs to fully consider the energy spread of π^0-decay using Eq. (<ref>). In this case, the lower bound of Eq. (<ref>) is ∼E_π^0 which is independent of E_γ and is larger than E_π, min(E_γ). With the π^0 flux ∝ E_π^-2.7 above E_π^0, the resulting γ-ray flux below the peak would be independent of E_γ, or equivalently, E_γ^2Φ_γ∝ E_γ^2, as clearly shown in Fig. <ref>.
  As m_χ increases, the peak position of E_γ^2Φ_γ shifts to low values of E_γ as E^ peak_γ∝ m_χ^-1.
  Correspondingly, the peak magnitude of E_γ^2Φ_γ scales approximately as m_χ^-0.7. To derive the upper limits at the 90% confidence level (C.L.) from the Fermi data, we require that, at any given energy, Φ_ Fermi(=Φ_γ) is
  smaller than the upper edge of the GC excess band 
  multiplied by a factor of 1.28.[Note that to obtain the 90% C.L. upper limits, we have simply assumed that the Fermi excess data follow a Gaussian distribution with a mean value of zero and a variance given by the upper edge of the band. For the H.E.S.S. data, we also assume a Gaussian distribution. The upper limits derived in our work are insensitive to these assumptions.]
  Similarly, for the H.E.S.S. data, we exclude the regions of large coupling constants if the predicted flux Φ_ H.E.S.S.(≈Φ_γ/290) exceeds the mean values of the measured diffuse flux plus 1.28 times the deviations at any energy. 
  The resulting 90% C.L. upper limits on g_χ g_q from the Fermi and the H.E.S.S. data as functions of
  the DM mass
  are presented in Fig. <ref>. As m_χ increases, the peak of E_γ^2 Φ_γ moves closer to the energies probed by the γ-ray telescopes, resulting in tighter bounds.
  For m_χ≳ 2× 10^-6 and ≳ 2 × 10^-4 GeV, the derived limits using the Fermi and the H.E.S.S. data become weaker with m_χ, due to the suppressed γ-ray fluxes.
          
          
  The constraints using the γ-ray fluxes are roughly consistent with the earlier work <cit.>,[Ref. <cit.> focused on dark matter heavier than 1 GeV.] although a constant inelastic scattering cross section was assumed. The reason is that the relevant energy window for E_γ is relatively narrow, and besides, the U(1)_B model adopted in our work exhibits a relatively flat inelastic cross section above the threshold. To set a scale, both our studies give an upper bound on the inelastic cross section of a few
   times 10^-29  cm^2 for m_χ around 0.1 GeV using the GC γ-ray data. Ref. <cit.> also checked that using the γ-ray data off the GC would only give rise to a relaxed bound.         
          
  
  
  
  
  
  
  
  
  
  
  
   §.§ Constraints from IceCube
  
  
  Both the accelerated DM and the secondary neutrinos can be detected by 
  IceCube.
  As the DM signals are almost indistinguishable from the neutral-current neutrino events <cit.>,
  a combined analysis of these signals should, in principle, be carried out. 
  In this study, however, we choose to do separate analyses for HE DM and neutrino signals. As shown below, the secondary neutrinos can set stronger limits than the DM signals for most of the DM masses considered.        
  
  We first consider the secondary neutrino signals. Fig. <ref> shows the resulting neutrino fluxes summing over all flavors from resonance excitations and DIS for m_χ=10^-6, 10^-3, and 10^-1 GeV, respectively. As discussed in Sec. <ref>, we take the isotropic approximation and show in the figure the averaged flux per solid angle, Φ_ν(E_ν)/(4π) = 1/(4π) ×∫ dΩϕ_ν(E_ν, Ω), with ϕ_ν(E_ν, Ω) given in Eq. (<ref>). For E_ν≳ PeV, the fluxes are mainly from DIS. With E_ν∝ E_p, ϕ_p(E_p) ∝ E_p^-γ, and σ_ DIS∝E_χ^1/3, the neutrino spectra scale as ϕ_p(E_p)σ_ DIS(E_χ) ∝  E_ν^-γ+1/3. With γ≈ 3 and 3.3 below and above the ankle followed by a sharp drop of the CR flux at E_p ≳ 3× 10^10 GeV due to the GZK cutoff, the neutrino spectra can be quantitatively understood.
  Meanwhile, the neutrino flux decreases as m_χ^-1. To constrain the DM model, the existing limits on the diffuse fluxes of HE neutrinos from IceCube <cit.> and Auger <cit.> as well as ANITA <cit.> can be directly used. In our study, we use the relevant limits at the 90% C.L. from Ref. <cit.> for 5× 10^6≲ E_ν≲ 3× 10^7 GeV and Ref. <cit.> for 3× 10^7 ≲ E_ν≲ 8× 10^10 GeV as shown in Fig. <ref>.
  Since E_ν^2 Φ_ν decreases with E_ν while the IceCube upper limit increases,
  the bounds on the couplings, g_χ g_q, are, in fact, determined by comparing the secondary neutrino flux and the upper limit from the IceCube publications at E_ν=5 PeV.
  Fig. <ref> shows the exclusion limits for g_χ g_q as functions of m_χ.
  Since the flux decreases with m_χ approximately as m_χ^-1, the upper limit on g_χ g_q increases as ∼m_χ^1/2. We notice that IceCube has accumulated around 100 events with deposited energies above ∼10 TeV, corresponding to an observed all-flavor flux of E_ν^2 ϕ_ ob≈ 2 × 10^-8  GeV cm^-2 s^-1 sr^-1 for 0.3 ≲ E_ν≲ 3 PeV (see Fig. 2 of Ref. <cit.>). Taking into account the uncertainties for the observed data,
  extending the energy window to E_ν≲ 1 PeV
  does not help to strengthen the bounds on the coupling constants.    
  
  
  The situation for DM signals is more complicated, since the related detector responses are not clear. Considering that no event with deposited energy above 5 PeV has been observed, we simply require the expected number n_χ of DM-induced events at IceCube with deposited energies above 5 PeV to be smaller than 2.4, which leads to the Feldman-Cousins upper limit at the 90% C.L. for the case of negligible background <cit.>. Considering the potential contributions from the cosmogenic HE neutrinos, the limit derived this way is a conservative one. The expected number n_χ can be estimated as 
  
      n_χ≈∫ dΩ∫_5  PeV dν  M_TN_AT Φ_χ,  att^ MW(T_χ, Ω) dσ_ DIS/dν ,
  
  where M_T ≈ 400 Mton is the IceCube effective mass at HE regions <cit.>, N_A is the Avogadro constant, T≈ 3142.5 days is the exposure time taken from Ref. <cit.>, and Φ_χ,  att^ MW is the attenuated DM flux arriving at the detector located 1450 meters below the surface. For 
  detecting the HE DM at IceCube, we need only to consider the scattering cross section from DIS   
  
      dσ_ DIS/dν = ∫ dQ^2 d^2σ_ DIS/dν dQ^2 .
  
  To obtain Φ_χ,  att^ MW, we take into account the energy loss of HE DM when passing through Earth/ice-shell from different declinations and   
  follow Ref. <cit.> to treat the attenuation effects. Again, DIS is the dominant channel to be considered for the energy losses of HE DM. The energy loss rate per unit distance is given by 
  
      dT_χ/dx = - n_N ∫dσ_ DIS/dνν dν ,
  
  where n_N is the number density of nucleons. For DM arriving from different declinations, the lengths traversed by DM inside the Earth are different. We use the Earth density profile from Ref. <cit.> for our study.      
  
  
  The excluded region on g_χ g_q by considering the DM signals at IceCube has been shown in Fig. <ref> (see the gray region). Unlike the limits from the secondary γ-rays and neutrinos, the excluded region from DM signals has an upper bound due to the Earth attenuation effects. For m_χ≳ 10^-7 GeV, the bounds from the upscattered DM are weaker than those from the secondary neutrinos. The reason is that the DM flux decreases rapidly with T_χ and is much smaller than that of the secondary neutrinos at high energies above ∼PeV.
  Without considering the attenuation effects, the expected number of DM signals is proportional to Φ_χσ_ DIS∝ m_χ^-2 (g_χ g_q)^4, which is 
  different from the secondary flux that scales as (g_χ g_q)^2. Consequently, the lower limit of g_χ g_q would grow as m_χ^1/2. However, this simple argument is no longer valid when the attenuation effect is taken into account. Relevant to the lower limits, the expected number of DM events grows more slowly than (g_χ g_q)^4. This finally results in a more rapid increase of the lower bound of g_χ g_q approximately as m_χ and a limited exclusion region. As discussed above, due to a rapidly decreasing differential cross section dσ_ DIS/dQ^2, the upscattered DM fluxes at high energies are suppressed. Therefore, compared to Ref. <cit.>, which assumed a constant cross section and differential cross section, the exclusion region derived in this work shrinks largely.
  
  
  
  
   §.§ Constraints from low-energy neutrino and DM experiments
   
  
  For the vector portal model considered in this work, the differential cross section of χ-p scattering at low Q^2 is dominated by the elastic channel and is orders of magnitude higher than that at high Q^2 (see Fig. <ref>). Therefore, as shown in Fig. <ref>, we expect a much higher DM flux at lower energies that can be detected at low-energy neutrino or DM detectors, resulting in a strong bound on the tested model. Using the subroutines implemented in DarkSUSY <cit.>, we show in Fig. <ref> the constraints on g_χ g_q from XENON1T (green regions) and MiniBooNE (pink regions), respectively. Only the elastic scattering has been considered for calculating the boosted DM flux and for estimating the event rates at detectors. Note that, for the case with a massless mediator, the energy loss rate of DM passing through the Earth, dE_χ/dx ∝∫_0^T_A^ max dT_A  dσ_ el/dQ^2 with Q^2=2m_AT_A, m_A=Am_N the mass, and T_A the recoil energy of nuclei with atomic mass number A, is logarithmically divergent.
  This is the reason why we choose m_V=1 eV instead of a massless mediator, see also footnote <ref>. 
  We also point out that in the original DarkSUSY subroutines a low-energy cutoff T_A^ min=1 eV has been taken, corresponding to a lower cutoff on Q^2 of (2× 10^9 A)  eV^2, which is much larger than m_V^2 for our light mediator case. Thus, when calculating the DM attenuation, we lower T_A^ min so that the related Q^2_ min is far below m_V^2. We have checked explicitly that our results are insensitive to the cutoff.               
  
  
  The searches for CR-boosted DM in XENON1T set the most stringent bounds for most of the DM mass range explored in this work due to the significantly enhanced cross section and flux at low energies. With a shallow location and, thus, a small Earth attenuation effect, the large couplings which cannot be probed by deep underground
  experiments can be further excluded by the MiniBooNE data (see the pink regions). For the case of m_V = 3m_χ,
  the limits on g_χ g_q are quantitatively similar to those for the light mediator case with m_V=1 eV at m_χ≲ 10^-2 GeV. This is simply because the typical values of Q^2 encountered in upscattering the DM by CRs or detecting the DM at detectors are higher than m_V^2 for both cases.
  However, as m_χ increases, m_V^2
  becomes comparable to the typical values of Q^2 for the case of m_V=3m_χ, leading to suppressed cross sections and enhanced upper/lower limits on the coupling constants compared to the light mediator case.       
  
  For the vector portal model considered, the bounds set by the secondary signals are typically weaker than those from XENON1T. However, unlike the DM signals at terrestrial detectors, the secondary signals do not suffer from the attenuation effect and
  can be used to exclude 
  all the strong-coupling regions above the upper limits. Especially for m_χ≳ 10^-2 GeV, 
  the Fermi data exclude the regions of g_χ g_q ≳ 0.1 for the case of a light mediator which cannot be explored 
  via the underground DM
  experiments. For the regions of coupling constant explored, the associated cross sections between nucleons and DM are high enough so that light DM are thermally populated before BBN <cit.>. Therefore, the observed D/H abundance ratio excludes m_χ≲ 7.8 MeV for the case of a Dirac fermion DM <cit.> (see the vertical lines in Fig. <ref>). See also other tighter constraints from beam-dump experiments <cit.>, supernova cooling <cit.>, low energy neutron scattering <cit.>, and rare meson decays <cit.>. 
      
  
  For the sake of comparison to studies which set bounds on χ-p scattering cross section, we also provide some benchmark values of the elastic cross section relevant to the direct detection for given m_χ, m_V, and g_χ g_q. Since the scattering cross sections  
  depend on the energy scales involved, the non-relativistic limit of the scattering cross section σ_χ p^0 ≈ g_B^2g_χ^2μ_χ p^2/(π m_V^4) ∝∫ dQ^2 dσ_ el(Q^2=0)/dQ^2 with μ_χ p the reduced mass is often used as a reference for cross section comparison <cit.>. For m_χ≲ 1 GeV, σ_χ p^0 ∼ g_q^2g_χ^2 (m_χ/ MeV)^2 ( MeV/m_V)^4 × 10^-25  cm^2, where we have used the relation g_B=3g_q. 
  However, we note that 
  σ_χ p^0 is originally introduced for non-relativistic heavy DM and cannot  
  be readily applied to DM direct detection for light DM/mediator. 
  For this reason, we introduce an effective cross section for χ-p scattering as 
  
      σ_χ p^ eff≡∫_Q^2_a^Q^2_b dQ^2 dσ_ el/dQ^2 ∼ g_B^2 g_χ^2 Q̅^2/[π(m_V^2+Q̅^2)^2],
  
  where Q^2_a,b are the lower/upper bounds of Q^2 determined by the detection energy window for nuclear recoils and Q̅^2 ∼ Q^2_a,b is the typical (or averaged) value. The effective cross section should depend on E_χ. Here, we simply assume that E_χ is large enough to produce a detectable recoil energy in the detectors. Taking xenon as the target and a typical value of 10 keV for nuclear recoil, Q̅^2 ∼ (0.05  GeV)^2. For m_V^2 ≪Q̅^2, σ_χ p^ eff∼ g_q^2 g_χ^2 [(0.05  GeV)^2/Q̅^2] × 10^-29  cm^2, and for m_V^2 ≳Q̅^2, σ_χ p^ eff∼ g_q^2g_χ^2 [(0.05  GeV)^2/Q̅^2] ( GeV/m_V)^4 × 10^-34  cm^2.  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSION AND PROSPECT
  
  
  
  As an extension to previous studies, we have investigated the inelastic scattering of CRs on sub-GeV light DM within the Milky Way. At energies high enough to excite the nucleonic resonance states (resonance excitations) or to probe the quark constituent (DIS), the secondary γ-rays or neutrinos can be produced from hadronization and the subsequent meson decay, and act as another indirect probe of DM. To demonstrate these signatures, we consider a simple vector portal DM model in this work.
  We have calculated the expected fluxes of these secondary signals and used the γ-ray and neutrino telescopes such as Fermi, H.E.S.S., and IceCube to derive limits for
  this specific model. On the other hand, the upscattered DM from both elastic and inelastic scatterings are also consistently considered and can be directly detected at low- and high-energy detectors, which can also be used to constrain the same DM model. Because of a much larger cross section of elastic scattering at low energies for this specific model, the low-energy experiments such as XENON1T set stronger bounds than those obtained from the secondary signals and the HE component of the accelerated DM.                                    
  
  The bounds derived are weakly dependent on
  the spatial distribution of the CRs and the DM halo profiles (see Fig. <ref>). 
  We have considered only CR protons in our study and assumed that a fraction of 30% are protons in the HE CRs above ∼PeV. We do not expect our results to be affected too much by this assumption
  unless the proton abundance in HE CRs is extremely small. To calculate the secondary neutrinos from DIS, we have relied on the SOPHIA code using the Lund Monte Carlo generator JETSET to treat hadronization. Since all hadronization models are relatively well constrained by data and can give rise to quantitatively similar yields even at very HE regions,
  the associated uncertainties for the secondary HE neutrino fluxes should not largely affect our conclusion either.                       
  
  Although considering the secondary signals does not 
  improve the bounds for the vector portal model for m_χ≲ 10^-2 GeV, they could be interesting for other DM models with relatively suppressed/enhanced cross sections at low/high energies. One possible scenario is that the DM candidate carries no charge but only couples to the vector mediator via a dipole form; i.e., DM particles have dipole moments. Other possibilities include composite DM scenario or the scenarios that scalar/vector DM particles scatter with nucleons in the s-channel. In these cases, the HE signals could possibly lead to a better sensitivity than the low-energy ones. In this work, we choose to take the simple vector portal DM model as a starting point and show the possible relevance of inelastic channels in HE regions. We plan to extend these studies to other interesting models in the future.      
  
  
  
  
  
  § ACKNOWLEDGMENTS
  
  We thank Hsiang-nan Li for helpful discussions about this project. G. G. and M.-R. W. acknowledge support from the Academia Sinica by Grant No. AS-CDA-109-M11. Y.-L. S. T.
  was funded in part by the Taiwan Young Talent Programme of Chinese Academy of Sciences 
  under the Grant No. 2018TW2JA0005 and 
  the Ministry of Science and Technology, Taiwan under the Grant No. 109-2112-M-007-022-MY3.
  M.-R. W. acknowledges support from the Ministry of Science and Technology, Taiwan under Grants No. 108-2112-M-001-010 and No. 109-2112-M-001-004, and the Physics Division, National Center of Theoretical Science of Taiwan.
  Q. Y. is supported by the National Natural Science Foundation of China (No. 11722328 and No. 11851305), 
  the 100 Talents Program of Chinese Academy of Sciences and the Program for Innovative Talents and
  Entrepreneur in Jiangsu.
  
  
  99
  
  
  
  Akerib:2016vxi 
    D. S. Akerib et al. [LUX Collaboration],
    
    Phys. Rev. Lett.  118, 021303 (2017)
    doi:10.1103/PhysRevLett.118.021303
    [arXiv:1608.07648 [astro-ph.CO]].
    
    
  
  
  
  Akerib:2017kat 
    D. S. Akerib et al. [LUX Collaboration],
    
    Phys. Rev. Lett.  118, 251302 (2017)
    doi:10.1103/PhysRevLett.118.251302
    [arXiv:1705.03380 [astro-ph.CO]].
    
    
  
  
  
  Cui:2017nnn 
    X. Cui et al. [PandaX-II Collaboration],
    
    Phys. Rev. Lett.  119, 181302 (2017)
    doi:10.1103/PhysRevLett.119.181302
    [arXiv:1708.06917 [astro-ph.CO]].
    
    
  
  
  
  Fu:2016ega 
    C. Fu et al. [PandaX-II Collaboration],
    
    Phys. Rev. Lett.  118, 071301 (2017)
    Erratum: [Phys. Rev. Lett.  120, 049902 (2018)]
    doi:10.1103/PhysRevLett.120.049902, 10.1103/PhysRevLett.118.071301
    [arXiv:1611.06553 [hep-ex]].
    
    
  
  
  
  Aprile:2018dbl 
    E. Aprile et al. [XENON Collaboration],
    
    Phys. Rev. Lett.  121, 111302 (2018)
    doi:10.1103/PhysRevLett.121.111302
    [arXiv:1805.12562 [astro-ph.CO]].
    
    
  
  
  
  Aprile:2019dbj 
    E. Aprile et al. [XENON Collaboration],
    
    Phys. Rev. Lett.  122, 141301 (2019)
    doi:10.1103/PhysRevLett.122.141301
    [arXiv:1902.03234 [astro-ph.CO]].
    
    
  
  
  
  Liu:2017drf 
    J. Liu, X. Chen and X. Ji,
    
    Nature Phys.  13, 212 (2017)
    doi:10.1038/nphys4039
    [arXiv:1709.00688 [astro-ph.CO]].
    
    
  
  
  
  Schumann:2019eaa 
    M. Schumann,
    
    J. Phys. G 46, 103003 (2019)
    doi:10.1088/1361-6471/ab2ea5
    [arXiv:1903.03026 [astro-ph.CO]].
    
    
  
  
  
  Gaskins:2016cha 
    J. M. Gaskins,
    
    Contemp. Phys.  57, 496 (2016)
    doi:10.1080/00107514.2016.1175160
    [arXiv:1604.00014 [astro-ph.HE]].
    
    
  
  
  
  Leane:2020liq 
    R. K. Leane,
    
    arXiv:2006.00513 [hep-ph].
    
    
  
  
  
  PerezdelosHeros:2020qyt
  C. Pérez de los Heros,
  
  Symmetry 12, 1648 (2020)
  doi:10.3390/sym12101648
  [arXiv:2008.11561 [astro-ph.HE]].
  
  
  
  
  Buchmueller:2017qhf 
    O. Buchmueller, C. Doglioni and L. T. Wang,
    
    Nature Phys.  13, 217 (2017)
    doi:10.1038/nphys4054
    [arXiv:1912.12739 [hep-ex]].
    
    
  
  
  
  Kahlhoefer:2017dnp 
    F. Kahlhoefer,
    
    Int. J. Mod. Phys. A 32, 1730006 (2017)
    doi:10.1142/S0217751X1730006X
    [arXiv:1702.02430 [hep-ph]].
    
    
  
  
  
  Boveia:2018yeb 
    A. Boveia and C. Doglioni,
    
    Ann. Rev. Nucl. Part. Sci.  68, 429 (2018)
    doi:10.1146/annurev-nucl-101917-021008
    [arXiv:1810.12238 [hep-ex]].
    
    
  
  
  
  Knapen:2016cue 
    S. Knapen, T. Lin and K. M. Zurek,
    
    Phys. Rev. D 95, 056019 (2017)
    doi:10.1103/PhysRevD.95.056019
    [arXiv:1611.06228 [hep-ph]].
    
    
  
  
  
  An:2017ojc 
    H. An, M. Pospelov, J. Pradler and A. Ritz,
    
    Phys. Rev. Lett.  120, 141801 (2018)
    Erratum: [Phys. Rev. Lett.  121, 259903 (2018)]
    doi:10.1103/PhysRevLett.120.141801, 10.1103/PhysRevLett.121.259903
    [arXiv:1708.03642 [hep-ph]].
    
    
  
  
  
  Ibe:2017yqa
  M. Ibe, W. Nakano, Y. Shoji and K. Suzuki,
  
  JHEP 03, 194 (2018)
  doi:10.1007/JHEP03(2018)194
  [arXiv:1707.07258 [hep-ph]].
  
  
  
  
  Berlin:2018sjs 
    A. Berlin, D. Hooper, G. Krnjaic and S. D. McDermott,
    
    Phys. Rev. Lett.  121, 011102 (2018)
    doi:10.1103/PhysRevLett.121.011102
    [arXiv:1803.02804 [hep-ph]].
    
    
  
  
  
  Akesson:2018vlm 
    T. Åkesson et al. [LDMX Collaboration],
    
    arXiv:1808.05219 [hep-ex].
    
    
  
  
  
  Matsumoto:2018acr 
    S. Matsumoto, Y. L. S. Tsai and P. Y. Tseng,
    
    JHEP 07, 050 (2019)
    doi:10.1007/JHEP07(2019)050
    [arXiv:1811.03292 [hep-ph]].
    
    
  
  
  
  Berlin:2019uco 
    A. Berlin, R. T. D'Agnolo, S. A. R. Ellis, P. Schuster and N. Toro,
    
    Phys. Rev. Lett.  124, 011801 (2020)
    doi:10.1103/PhysRevLett.124.011801
    [arXiv:1908.06982 [hep-ph]].
    
    
  
  
  
  Depta:2019lbe 
    P. F. Depta, M. Hufnagel, K. Schmidt-Hoberg and S. Wild,
    
    JCAP 04, 029 (2019)
    doi:10.1088/1475-7516/2019/04/029
    [arXiv:1901.06944 [hep-ph]].
    
    
  
  
  
  Hertel:2019thc 
    S. A. Hertel, A. Biekert, J. Lin, V. Velan and D. N. McKinsey,
    
    Phys. Rev. D 100, 092007 (2019)
    doi:10.1103/PhysRevD.100.092007
    [arXiv:1810.06283 [physics.ins-det]].
    
    
  
  
  
  Dror:2019onn 
    J. A. Dror, G. Elor and R. Mcgehee,
    
    Phys. Rev. Lett.  124, 18 (2020)
    doi:10.1103/PhysRevLett.124.181301
    [arXiv:1905.12635 [hep-ph]].
    
    
  
  
  
  Dror:2019dib 
    J. A. Dror, G. Elor and R. Mcgehee,
    
    JHEP 02, 134 (2020)
    doi:10.1007/JHEP02(2020)134
    [arXiv:1908.10861 [hep-ph]].
    
    
  
  
  
  Bringmann:2018cvk 
    T. Bringmann and M. Pospelov,
    
    Phys. Rev. Lett.  122, 171801 (2019)
    doi:10.1103/PhysRevLett.122.171801
    [arXiv:1810.10543 [hep-ph]].
    
    
  
  
  
  Ema:2018bih 
    Y. Ema, F. Sala and R. Sato,
    
    Phys. Rev. Lett.  122, 181802 (2019)
    doi:10.1103/PhysRevLett.122.181802
    [arXiv:1811.00520 [hep-ph]].
    
    
  
  
  
  Cappiello:2019qsw 
    C. Cappiello and J. F. Beacom,
    
    Phys. Rev. D 100, 103011 (2019)
    doi:10.1103/PhysRevD.100.103011
    [arXiv:1906.11283 [hep-ph]].
    
    
  
  
  
  Guo:2020drq
  G. Guo, Y. L. S. Tsai and M. R. Wu,
  
  JCAP 10, 049 (2020)
  doi:10.1088/1475-7516/2020/10/049
  [arXiv:2004.03161 [astro-ph.HE]].
  
  
  
  
  Ge:2020yuf 
    S. F. Ge, J. L. Liu, Q. Yuan and N. Zhou,
    
    arXiv:2005.09480 [hep-ph].
    
    
  
  
  
  Zhang:2020htl 
    B. L. Zhang, Z. H. Lei and J. Tang,
    
    arXiv:2008.07116 [hep-ph].
    
    
  
  
  
  Bondarenko:2019vrb 
    K. Bondarenko, A. Boyarsky, T. Bringmann, M. Hufnagel, K. Schmidt-Hoberg and A. Sokolenko,
    
    JHEP 03, 118 (2020)
    doi:10.1007/JHEP03(2020)118
    [arXiv:1909.08632 [hep-ph]].
    
    
  
  
  
  Dent:2019krz 
    J. B. Dent, B. Dutta, J. L. Newstead and I. M. Shoemaker,
    
    Phys. Rev. D 101, 116007 (2020)
    doi:10.1103/PhysRevD.101.116007
    [arXiv:1907.03782 [hep-ph]].
    
    
  
  
  
  Wang:2019jtk 
    W. Wang, L. Wu, J. M. Yang, H. Zhou and B. Zhu,
    
    arXiv:1912.09904 [hep-ph].
    
    
  
  
  
  Cho:2020mnc 
    W. Cho, K. Y. Choi and S. M. Yoo,
    
    arXiv:2007.04555 [hep-ph].
    
    
  
  
  
  Cao:2020bwd 
    Q. H. Cao, R. Ding and Q. F. Xiang,
    
    arXiv:2006.12767 [hep-ph].
    
    
  
  
  
  Aprile:2020tmw
  E. Aprile et al. [XENON],
  
  Phys. Rev. D 102, 072004 (2020)
  doi:10.1103/PhysRevD.102.072004
  [arXiv:2006.09721 [hep-ex]].
  
  
  
  
  Jho:2020sku
  Y. Jho, J. C. Park, S. C. Park and P. Y. Tseng,
  
  Phys. Lett. B 811, 135863 (2020)
  doi:10.1016/j.physletb.2020.135863
  [arXiv:2006.13910 [hep-ph]].
  
  
  
  
  
  Bloch:2020uzh 
    I. M. Bloch, A. Caputo, R. Essig, D. Redigolo, M. Sholapurkar and T. Volansky,
    
    arXiv:2006.14521 [hep-ph].
    
    
  
  
  
  Cappiello:2018hsu 
    C. V. Cappiello, K. C. Y. Ng and J. F. Beacom,
    
    Phys. Rev. D 99, 063004 (2019)
    doi:10.1103/PhysRevD.99.063004
    [arXiv:1810.07705 [hep-ph]].
    
    
  
  
  
  Alvey:2019zaa 
    J. Alvey, M. Campos, M. Fairbairn and T. You,
    
    Phys. Rev. Lett.  123, 261802 (2019)
    doi:10.1103/PhysRevLett.123.261802
    [arXiv:1905.05776 [hep-ph]].
    
    
  
  
  
  Plestid:2020kdm 
    R. Plestid, V. Takhistov, Y. D. Tsai, T. Bringmann, A. Kusenko and M. Pospelov,
    
    arXiv:2002.11732 [hep-ph].
    
    
  
  
  
  Su:2020zny 
    L. Su, W. Wang, L. Wu, J. M. Yang and B. Zhu,
    
    arXiv:2006.11837 [hep-ph].
    
    
  
  
  
  Cyburt:2002uw 
    R. H. Cyburt, B. D. Fields, V. Pavlidou and B. D. Wandelt,
    
    Phys. Rev. D 65, 123503 (2002)
    doi:10.1103/PhysRevD.65.123503
    [astro-ph/0203240].
    
    
  
  
  
  Hooper:2018bfw 
    D. Hooper and S. D. McDermott,
    
    Phys. Rev. D 97, 115006 (2018)
    doi:10.1103/PhysRevD.97.115006
    [arXiv:1802.03025 [hep-ph]].
    
    
  
  
  Dobrescu:2014fca
  B. A. Dobrescu and C. Frugiuele,
  
  Phys. Rev. Lett. 113,
  061801 (2014)
  doi:10.1103/PhysRevLett.113.061801
  [arXiv:1404.3947 [hep-ph]].
  
  
  
  
  Dror:2017ehi 
    J. A. Dror, R. Lasenby and M. Pospelov,
    
    Phys. Rev. Lett.  119, 141803 (2017)
    doi:10.1103/PhysRevLett.119.141803
    [arXiv:1705.06726 [hep-ph]].
    
    
  
  
  
  Foot:1990mn 
    R. Foot,
    
    Mod. Phys. Lett. A 06, 527 (1991).
    doi:10.1142/S0217732391000543
    
    
  
  
  
  Okun:1982xi
  L. B. Okun,
  
  Sov. Phys. JETP 56, 502 (1982).
  
    
    
  
  Holdom:1985ag 
    B. Holdom,
    
    Phys. Lett.  166B, 196 (1986).
    doi:10.1016/0370-2693(86)91377-8
    
    
  
  
  An:2013yfc
  H. An, M. Pospelov and J. Pradler,
  
  Phys. Lett. B 725, 190-195 (2013)
  doi:10.1016/j.physletb.2013.07.008
  [arXiv:1302.3884 [hep-ph]].
  
  
  
  Batell:2014mga
  B. Batell, R. Essig and Z. Surujon,
  
  Phys. Rev. Lett. 113, 171802 (2014)
  doi:10.1103/PhysRevLett.113.171802
  [arXiv:1406.2698 [hep-ph]].
  
  
  
  Heeck:2014zfa 
    J. Heeck,
    
    Phys. Lett. B 739, 256 (2014)
    doi:10.1016/j.physletb.2014.10.067
    [arXiv:1408.6845 [hep-ph]].
    
    
  
  
  
  Bilmis:2015lja 
    S. Bilmis, I. Turan, T. M. Aliev, M. Deniz, L. Singh and H. T. Wong,
    
    Phys. Rev. D 92, 033009 (2015)
    doi:10.1103/PhysRevD.92.033009
    [arXiv:1502.07763 [hep-ph]].
    
    
  
  
  Alexander:2016aln 
    J. Alexander et al.,
    
    arXiv:1608.08632 [hep-ph].
    
    
  
  
  Fayet:2016nyc 
    P. Fayet,
    
    Eur. Phys. J. C 77, 53 (2017)
    doi:10.1140/epjc/s10052-016-4568-9
    [arXiv:1611.05357 [hep-ph]].
    
    
  
  
  Knapen:2017xzo 
    S. Knapen, T. Lin and K. M. Zurek,
    
    Phys. Rev. D 96, 115021 (2017)
    doi:10.1103/PhysRevD.96.115021
    [arXiv:1709.07882 [hep-ph]].
    
    
  
  
  Chang:2018rso
  J. H. Chang, R. Essig and S. D. McDermott,
  
  JHEP 09, 051 (2018)
  doi:10.1007/JHEP09(2018)051
  [arXiv:1803.00993 [hep-ph]].
  
  
  
  Bauer:2018onh 
    M. Bauer, P. Foldenauer and J. Jaeckel,
    
    JHEP 07, 094 (2018)
    doi:10.1007/JHEP07(2018)094
    [arXiv:1803.05466 [hep-ph]].
    
    
  
  
  Sung:2019xie
  A. Sung, H. Tu and M. R. Wu,
  
  Phys. Rev. D 99, 121305 (2019)
  doi:10.1103/PhysRevD.99.121305
  [arXiv:1903.07923 [hep-ph]].
  
  
  
  DeRocco:2019njg
  W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares, and S. Rajendran,
  
    JHEP 02, 171 (2019)
    doi:10.1007/JHEP02(2019)171
    [arXiv:1901.08596 [hep-ph]].
  
  Li:2020roy 
    J.-T.. Li, G. Fuller and E. Grohs,
    arXiv:2009.14325[astro-ph.CO].
  
  
  Fabbrichesi:2020wbt 
    M. Fabbrichesi, E. Gabrielli and G. Lanfranchi,
    
    arXiv:2005.01515 [hep-ph].
    
    
  
  
  
  Barbieri:1975xy
  R. Barbieri and T. E. O. Ericson,
  
  Phys. Lett. B 57, 270-272 (1975)
  doi:10.1016/0370-2693(75)90073-8
  
  
  
  Batell:2014yra
  B. Batell, P. deNiverville, D. McKeen, M. Pospelov and A. Ritz,
  
  Phys. Rev. D 90, 115014 (2014)
  doi:10.1103/PhysRevD.90.115014
  [arXiv:1405.7049 [hep-ph]].
  
  
  
  
  Aguilar-Arevalo:2017mqx
  A. A. Aguilar-Arevalo et al. [MiniBooNE],
  
  Phys. Rev. Lett. 118, 221803 (2017)
  doi:10.1103/PhysRevLett.118.221803
  [arXiv:1702.02688 [hep-ex]].
  
  
  
  Dror:2017nsg
  J. A. Dror, R. Lasenby and M. Pospelov,
  
  Phys. Rev. D 96, 075036 (2017)
  doi:10.1103/PhysRevD.96.075036
  [arXiv:1707.01503 [hep-ph]].
  
  
  
  
  Aguilar-Arevalo:2018wea
  A. A. Aguilar-Arevalo et al. [MiniBooNE DM],
  
  Phys. Rev. D 98, 112004 (2018)
  doi:10.1103/PhysRevD.98.112004
  [arXiv:1807.06137 [hep-ex]].
  
  
  
  Krnjaic:2019dzc
  G. Krnjaic and S. D. McDermott,
  
  Phys. Rev. D 101, 123022 (2020)
  doi:10.1103/PhysRevD.101.123022
  [arXiv:1908.00007 [hep-ph]].
  
  
  
  Rrapaj:2015wgs
  E. Rrapaj and S. Reddy,
  
  Phys. Rev. C 94, 045805 (2016)
  doi:10.1103/PhysRevC.94.045805
  [arXiv:1511.09136 [nucl-th]].
  
  
  
  Aprile:2015uzo 
    E. Aprile et al. [XENON Collaboration],
    
    JCAP 04, 027 (2016)
    doi:10.1088/1475-7516/2016/04/027
    [arXiv:1512.07501 [physics.ins-det]].
    
    
  
  
  
  Aprile:2017aty 
    E. Aprile et al. [XENON Collaboration],
    
    Eur. Phys. J. C 77, 881 (2017)
    doi:10.1140/epjc/s10052-017-5326-3
    [arXiv:1708.07051 [astro-ph.IM]].
    
    
  
  
  
  AguilarArevalo:2008qa 
    A. A. Aguilar-Arevalo et al. [MiniBooNE Collaboration],
    
    Nucl. Instrum. Meth. A 599, 28 (2009)
    doi:10.1016/j.nima.2008.10.028
    [arXiv:0806.4201 [hep-ex]].
    
    
  
  
  
  Aartsen:2016nxy 
    M. G. Aartsen et al. [IceCube Collaboration],
    
    JINST 12, P03012 (2017)
    doi:10.1088/1748-0221/12/03/P03012
    [arXiv:1612.05093 [astro-ph.IM]].
    
    
  
  
  
  TheFermi-LAT:2017vmf 
    M. Ackermann et al. [Fermi-LAT Collaboration],
    
    Astrophys. J.  840, 43 (2017)
    doi:10.3847/1538-4357/aa6cab
    [arXiv:1704.03910 [astro-ph.HE]].
    
    
  
  
  
  Abramowski:2016mir 
    A. Abramowski et al. [H.E.S.S. Collaboration],
    
    Nature 531, 476 (2016)
    doi:10.1038/nature17147
    [arXiv:1603.07730 [astro-ph.HE]].
    
    
  
  
  Uehling:1954wp
  E. A. Uehling,
  
  Ann. Rev. Nucl. Part. Sci. 4, 315-350 (1954)
  doi:10.1146/annurev.ns.04.120154.001531.
  
  
  
  Angeli:2004kvy 
    I. Angeli,
    
    Atom. Data Nucl. Data Tabl.  87, 185 (2004).
    doi:10.1016/j.adt.2004.04.002
    
    
  
  
  
  Formaggio:2013kya 
    J. A. Formaggio and G. P. Zeller,
    
    Rev. Mod. Phys.  84, 1307 (2012)
    doi:10.1103/RevModPhys.84.1307
    [arXiv:1305.7513 [hep-ex]].
    
    
  
  
  Zuber
  Kai Zuber, 
  Neutrino Physics (Institute of Physics Publishing, Ltd., London, 2004).
  
  
  Leitner:2009zz 
   T. J. Leitner, Ph.D. thesis, UniversitätGießen, 2009. 
    
    
    
  
  
  
  Buss:2011mx 
    O. Buss et al.,
    
    Phys. Rept.  512, 1 (2012)
    doi:10.1016/j.physrep.2011.12.001
    [arXiv:1106.1344 [hep-ph]].
    
    
  
  
  
  Tiator:2006dq
  L. Tiator and S. Kamalov,
  
  arXiv:nucl-th/0603012.
  
  
  
  
  Drechsel:2007if 
    D. Drechsel, S. S. Kamalov and L. Tiator,
    
    Eur. Phys. J. A 34, 69 (2007)
    doi:10.1140/epja/i2007-10490-6
    [arXiv:0710.0306 [nucl-th]].
    
    
  
  
  MSTW
  .
  
  
  Martin:2009iq 
    A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt,
    
    Eur. Phys. J. C 63, 189 (2009)
    doi:10.1140/epjc/s10052-009-1072-5
    [arXiv:0901.0002 [hep-ph]].
    
    
  
  
  
  Strong:1998pw 
    A. W. Strong and I. V. Moskalenko,
    
    Astrophys. J.  509, 212 (1998)
    doi:10.1086/306470
    [astro-ph/9807150].
    
    
  
  
  Zyla:2020zbs
  P. A. Zyla et al. [Particle Data Group],
  
  PTEP 2020, 083C01 (2020)
  doi:10.1093/ptep/ptaa104.
  
  
  
  Navarro:1995iw 
    J. F. Navarro, C. S. Frenk and S. D. M. White,
    
    Astrophys. J.  462, 563 (1996)
    doi:10.1086/177173
    [astro-ph/9508025].
    
    
  
  
  
  Navarro:1996gj 
    J. F. Navarro, C. S. Frenk and S. D. M. White,
    
    Astrophys. J.  490, 493 (1997)
    doi:10.1086/304888
    [astro-ph/9611107].
    
    
  
  
  
  Burkert:1995yz 
    A. Burkert,
    
    IAU Symp.  171, 175 (1996)
    [Astrophys. J. Lett.  447, L25 (1995)]
    doi:10.1086/309560
    [astro-ph/9504041].
    
    
  
  
  
  Mucke:1999yb 
    A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe and T. Stanev,
    
    Comput. Phys. Commun.  124, 290 (2000)
    doi:10.1016/S0010-4655(99)00446-4
    [astro-ph/9903478].
    
    
  
  
  
  Aartsen:2016ngq 
    M. G. Aartsen et al. [IceCube Collaboration],
    
    Phys. Rev. Lett.  117, 241101 (2016)
    Erratum: [Phys. Rev. Lett.  119, 259902 (2017)]
    doi:10.1103/PhysRevLett.117.241101, 10.1103/PhysRevLett.119.259902
    [arXiv:1607.05886 [astro-ph.HE]].
    
    
  
  
  
  Kopper:2017zzm 
    C. Kopper [IceCube Collaboration],
    
    PoS ICRC 2017, 981 (2018).
    doi:10.22323/1.301.0981
    
    
  
  
  
  Aartsen:2018vtx 
    M. G. Aartsen et al. [IceCube Collaboration],
    
    Phys. Rev. D 98, 062003 (2018)
    doi:10.1103/PhysRevD.98.062003
    [arXiv:1807.01820 [astro-ph.HE]].
    
    
  
  
  
  Aab:2015kma 
    A. Aab et al. [Pierre Auger Collaboration],
    
    Phys. Rev. D 91, 092008 (2015)
    doi:10.1103/PhysRevD.91.092008
    [arXiv:1504.05397 [astro-ph.HE]].
    
    
  
  
  
  Allison:2018cxu 
    P. W. Gorham et al. [ANITA Collaboration],
    
    Phys. Rev. D 98, 022001 (2018)
    doi:10.1103/PhysRevD.98.022001
    [arXiv:1803.02719 [astro-ph.HE]].
    
    
  
  
  
  Feldman:1997qc 
    G. J. Feldman and R. D. Cousins,
    
    Phys. Rev. D 57, 3873 (1998)
    doi:10.1103/PhysRevD.57.3873
    [physics/9711021 [physics.data-an]].
    
    
  
  
  
  Aartsen:2013jdh 
    M. G. Aartsen et al. [IceCube Collaboration],
    
    Science 342, 1242856 (2013)
    doi:10.1126/science.1242856
    [arXiv:1311.5238 [astro-ph.HE]].
    
    
  
  
  McDonough:2003
  W. F. McDonough,
  Treatise Geochem. 2, 547 (2003).
  
  
  Bringmann:2018lay
  T. Bringmann, J. Edsjö, P. Gondolo, P. Ullio and L. Bergström,
  
  JCAP 07, 033 (2018)
  doi:10.1088/1475-7516/2018/07/033
  [arXiv:1802.03399 [hep-ph]].
  
  
  
  
  
  