Search for U1Lμ−Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_{\mu }-{L}_{\tau }} $$\end{document} charged dark matter with neutrino telescope

We study a simple Dirac fermion dark matter model in U1Lμ−Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_{\mu }-{L}_{\tau }} $$\end{document} theory. The new light gauge boson X plays important roles in both dark matter physics and the explanation for the muon g− 2 anomaly. The observed dark matter relic density is realized by a large U1Lμ−Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_{\mu }-{L}_{\tau }} $$\end{document} charge without introducing a resonance effect of the X boson. As a by-product of the model, characteristic neutrino signatures from sub-GeV dark matter ψ are predicted depending on the mass spectrum. We formulate the analysis of ψψ¯→νν¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \psi \overline{\psi}\to \nu \overline{\nu} $$\end{document}, and of ψψ¯→XX\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \psi \overline{\psi}\to XX $$\end{document} followed by X→νν¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ X\to \nu \overline{\nu} $$\end{document} in a model independent way. The energy spectrum of neutrinos in the former process is monochromatic while in the latter process is bowl-shape. We also evaluate sensitivity at Super-Kamiokande and future Hyper-Kamiokande detectors. The analysis is finally applied to the U1Lμ−Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{U}{(1)}_{L_{\mu }-{L}_{\tau }} $$\end{document} dark matter model.


Introduction
High energy collider experiments have confirmed the Standard Model (SM) in particle physics with high precision up to the TeV scale. In spite of the great success, there is some compelling empirical evidence for new physics, one of the most important ones being dark matter (DM). In the last decades, significant attention has been paid to Weakly Interacting Massive Particles (WIMPs) DM, with considerable efforts to directly and indirectly test such DM candidates. It is unfortunate, however, that all of these experimental efforts have obtained no affirmative signals thus far. The resulting stringent constraints often put pressure on the traditional WIMP DM with the electroweak scale mass. This tendency encourages theorists to embark on an unexplored sub-GeV mass region. A simple realization of the sub-GeV DM includes a secluded scenario [1] in which DM couples to the visible sector only via new light mediator particles, such as dark gauge bosons or dark Higgs bosons [1][2][3][4][5][6][7][8][9].
Another attractive hint for new physics is the so-called muon g − 2 anomaly. There is over 3 σ discrepancy between the theoretical prediction and experimental value of the muon anomalous magnetic moment in light of the result measured at the BNL E821 [10,11]. This tension suggests that there exist new interactions in the muon sector, leading to growing interests in lepton flavored new physics scenarios. One well-studied class of extensions is JHEP03(2021)047 gauged U(1) Lµ−Lτ models [12][13][14][15]. The associated gauge boson with the sub-GeV mass has a sizable contribution to accommodate the discrepancy, while avoiding existing constraints from collider experiments and cosmological and astrophysical observations because of the absence of direct coupling with electron and quarks [16][17][18][19]. Next generation low-energy experiments will have the sensitivity to most of the parameter space that addresses the discrepancy.
Given the potential of the muon g − 2 explanation and the increased interest in the secluded DM scenario, it is natural to ask if the U(1) Lµ−Lτ gauge boson can be a mediator between DM and the SM sector. In the present paper, we therefore consider a SM singlet vectorlike fermion DM ψ, charged under the U(1) Lµ−Lτ group with the associated gauge boson X, and scrutinize the phenomenology. Because DM is introduced vectorlike, gauge anomalies relative to U(1) Lµ−Lτ group, which come from DM contributions, are cancelled automatically. Similar DM candidates have extensively been studied in the literature, mostly in the traditional (heavy) WIMP regime [20][21][22][23] and by introducing additional new degrees of freedom for the successful DM production [24][25][26]. The light mass window is studied in [27], under the assumption with no charge hierarchy between DM and muon. It is shown that the resonant DM annihilation via the X boson physical pole is required to correctly produce the DM relic density in the muon g−2 favored region. In contrast to these works, we do not restrict the DM gauge charge to unity and not introduce any additional fields except for DM and the gauge boson. We then scan the whole parameter space. As a result, we find that allowing a charge hierarchy opens a new possibility that both the muon g − 2 and DM relic abundance are explained without the resonance enhancement. It also turns out that the DM mass is restricted to be 10 MeV-1 GeV, in light of the direct detection experiments and the Cosmic Microwave Background (CMB) observations.
We also examine the indirect detection signals of the DM candidate at the Super-Kamiokande (SK) and Hyper-Kamiokande (HK). It is remarkable that the model can be tested with neutrino telescopes by searching for an extra neutrino flux from galactic DM annihilation. The main annihilation mode of DM is the one into neutrinos or X bosons, depending on the mass spectrum of DM and X. In the former annihilation, the produced neutrino flux is monochromatic at E ν = m ψ , providing a sharp DM signal. In the latter one, the neutrinos are produced via the 1-step cascade process, ψψ → XX → 2ν2ν. The shape of the produced neutrino flux depends on the polarizations of the X boson and takes the bowl-like form as a function of neutrino energy [28]. In the present paper, we analyze both the monochromatic and bowl-like neutrino spectra, and derive upper limits on the annihilation cross sections by reinterpreting the results of the supernova relic neutrino (SRN) searches at the SK [29,30]. We would like to emphasize that while there are many similar studies with a special focus on the monochromatic spectrum [31][32][33][34][35][36][37][38][39][40][41][42][43], the analysis for the bowl neutrino spectrum is given for the first time in this paper. Although we analyze the neutrino flux originated from DM annihilation under the U(1) Lµ−Lτ DM model in the present paper, such analyses can be applied to non-flavored U(1) group and other kinds of DM. We further study the future sensitivity of the HK experiment and discuss the impacts on the U(1) Lµ−Lτ DM model. It is possible that the HK with an updated option of the doped Gd will probe the canonical thermal relic cross section in the O(10) MeV region. This paper is organized as follows. In section 2, we introduce the model Lagrangian involving two new particles, a Dirac DM and a massive gauge boson. The constraints on the new gauge boson are also discussed there. We study the DM phenomenology in section 3, and examine the indirect detection signals of the DM annihilation at the SK and HK experiments in section 4. We summarize our results in section 5.

Model
The model is built on an SU(3) c ×SU(2) L ×U(1) Y ×U(1) Lµ−Lτ gauge theory. We introduce a SM singlet Dirac DM ψ that carries an L µ −L τ charge q ψ in unit of the muon (tau) charge being +1 (−1). The charge assignment for DM and leptons is summarized in table 1. The other SM fields have no L µ − L τ charge. The DM stability is guaranteed by an accidental global U(1) symmetry. In this paper, we assume that neutrinos can approximate massless because their masses are much lighter than the other particles. In fact, neutrino oscillation has been observed by various experiments, and then neutrinos have undoubtedly non-zero masses. There are some previous works which discuss neutrino masses in the U(1) Lµ−Lτ model, for example refs. [44][45][46][47][48][49][50]. However, neutrino masses are irrelevant in the following discussion about DM, and thus we neglect them.
The model Lagrangian is given by where X µ denotes the U(1) Lµ−Lτ gauge boson with the field strength X µν . The mass of X is generated by the Higgs mechanism or Stückelberg mechanism. Since the origin of the mass is irrelevant to our discussion, we do not specify it. Now, we have five free parameters, m X , g X , m ψ , q ψ and in the model. In this paper, we will consider q ψ as a free parameter that takes an arbitrary value unless it violates unitarity. As will see in section 3, a large q ψ is indeed required to explain the observed DM abundance. Such a large q ψ is allowed phenomenologically, but it may be unnatural from the theoretical perspective. We do not discuss in this paper any UV origin of the large charge hierarchy between DM and muon, but one of the possible UV extensions includes the addition of another dark U(1) gauge symmetry that is coupled only to DM. If DM is not charged under the U(1) Lµ−Lτ symmetry and there is a small mixing between two U(1)'s, we will be able to apparently realize the large charge hierarchy.

JHEP03(2021)047
We further assume the kinetic mixing of X and the hypercharge gauge field B is vanishing at some high scale. Nonetheless, finite mixings of X with photon A and with the Z boson arise at the one-loop level at low energy: where F µν and Z µν denote the field strength for photon and Z boson, respectively. We will take into account this one-loop level mixing in our analysis. The mixing with photon is calculated below the muon mass scale as To get the canonically normalized gauge fields, we shift the photon field as This shift induces the interaction of X to the electromagnetic current A eX µ J µ em , which is crucial for DM direct detection. Similarly, the mixing with the Z boson is given by with the Weinberg angle θ W . Shifting Z µ → Z µ − Z m 2 X /m 2 Z X µ and X µ → X µ + Z Z µ to get the canonically normalized gauge fields, we find the Z boson couples to L µ − L τ current and the X boson to the neutral current, Since the effect of the mixing Z is suppressed by m 2 X /m 2 Z and is negligible in the phenomenological study below, we ignore it in this paper.
Let us next discuss impacts of the X boson on the muon g − 2 anomaly and other experimental results. It can contribute to the anomalous magnetic moment of muon a µ = (g µ − 2)/2. The one-loop contribution is evaluated as The latest result announced by E821 [10,11] and the theoretical calculation find a 3.3σ discrepancy [51], Since ∆a X µ g 2 X /(8π 2 ) for m X m µ , the discrepancy is resolved with g X 5 × 10 −4 . Thus, a finite but small g X is favored with respect to the muon g − 2. In figure 1, we show the parameter region, in which the muon g − 2 is explained within 1σ (2σ), with the (light) red band. We also illustrate the experimental upper limits on g X from a BABAR search for e + e − → µμX with a subsequent decay X → µμ [52], CHARM-II and CCFR measurements JHEP03(2021)047 CCFR CHARM muon g-2 of the neutrino trident production νN → νN µμ [53][54][55], a Borexino measurement of the interaction rate of the mono-energetic 862 keV 7 Be solar neutrino [56][57][58], and the white dwarf (WD) cooling induced by the plasmon decay via the off-shell X boson [59,60]. Note that the BABAR limit depends on the branching fraction of the X → ψψ decay, because the BABAR experiment searches for a dark photon decaying into a muon pair. As the fraction of the X → ψψ decay is larger, the limit is weaker. We show in the figure the BABAR limit for q ψ = 5 with the dashed blue line, assuming DM is much lighter than X. The neutrino trident, Borexino and WD bounds are not influenced by the existence of DM, since the off-shell DM can only contribute to them. Hereafter, we focus on the parameter space where the muon g − 2 is explained within 2σ while avoiding the other constraints from the BABAR, CHARM, Borexino experiments and WD cooling.
There is a cosmological bound on the X boson in addition to the above experimental bounds. The X boson mainly decays into a pair of neutrinos for m X < 2m ψ . The lifetime is much shorter than the time scale (τ BBN ∼ 1 sec) of Big Bang Nucleosynthesis (BBN) in the parameter space that we are interested in. If the X boson is light enough, the X boson can be in equilibrium with the neutrinos after the neutrino decoupling (T ∼ 1 MeV). Then, it is possible to increase the effective neutrino number N eff . Since the X boson is already nonrelativistic at the BBN, it does not contribute directly to N eff . It can decay into neutrinos, however. The decay releases the energy into neutrinos, reheating the neutrino temperature. Following [61], we estimate the contribution to N eff and obtain the lower mass bound m X 6 MeV for A = 0. When we include the effects of the non-vanishing kinetic mixing, the lower bound becomes 10 MeV for A 7.2 × 10 −6 corresponding to g X = 5 × 10 −4 .

JHEP03(2021)047
For further details of the contribution to N eff and the cosmological implication of the light U(1) Lµ−Lτ gauge boson, see e.g. ref. [62]. We conservatively use m X 6 MeV as the lower mass bound in this paper.

Dark matter physics
In this section, we study the DM thermal production and constraints from direct detection and CMB observations.

Relic density
The DM thermal relic abundance is determined by solving the Boltzmann equation, where H denotes the Hubble parameter, n the total (equilibrium) number density of DM, σ ann v the DM total annihilation cross section averaged over the thermal distributions of the initial state. Study on precise calculation of the thermal relic density [63,64] indicates the canonical value of the cross section to explain the observed DM abundance, for DM mass ranging from MeV to 100 TeV. In this paper, we employ micromegas 4 3 5 [65] to calculate the DM thermal abundance. There are two relevant annihilation processes, ψψ → XX, ff (f = µ, τ, ν µ , ν τ ) in this model. The dominant annihilation mode depends on the mass spectrum of ψ and X. If ψ is lighter than X, the only possible annihilation is ψψ → ff through the s-channel X boson exchanging. The cross section is given by where Γ X denotes the total decay width of X and we only keep the partial s-wave. Since the gauge coupling g X is as small as O(10 −4 ), the resulting cross section is too small to achieve the correct relic abundance for q ψ = O(1). The sufficiently large cross section is realized for q ψ 1 or in the resonant region with m ψ m X /2. A comprehensive study of the resonant production has been made on a benchmark point m ψ = 0.45 m X [27], including various experimental constraints. We study another option of exploiting the large DM charge as well as the resonant case.
If DM is heavier than X, the other annihilation channel ψψ → XX opens. The cross section is (3.4)

JHEP03(2021)047
This is proportional to q 4 ψ , while the cross section for ψψ → ff to q 2 ψ . Since there is no resonance enhancement in this mass regime, q ψ has to be large for the successful thermal production. Therefore, the DM abundance is produced mostly by the ψψ → XX process.

Direct detection constraints
Elastic scattering of DM with a nucleus and an electron is severely constrained by direct detection experiments. The scattering occurs via one-loop kinetic mixing with photon A . Since the mixing with the Z boson is suppressed by a factor of m 2 X /m 2 Z , compared with the leading contribution, we neglect it here.
The cross section of the spin-independent DM-nucleus scattering is given by is the reduced mass of DM and a nucleus with m N being the nucleus mass, and Z the atomic number of the nucleus. This scattering is very similar to the one via the photon exchanging induced by the DM charge radius, b ψ = A q ψ g X /m 2 X . Thus, the cross section is proportional to Z 2 rather than A 2 with A being the atomic mass. The most stringent limit on the cross section is set by XENON1T [69] and a TEA-LAB simulated experiment inspired by the DarkSide50 result (0.05 GeV m ψ 0.1 GeV) [70].
Similarly, the cross section of the DM-electron scattering is given by where µ e denotes the DM-electron reduced mass. The constraints strongly depend on the scenario if the scattering is dependent on the momentum transfer q or not. If the X boson mass is much lighter than keV, the scattering is enhanced at a small momentum transfer. Then, we have to take into account a DM form factor, F DM = (αm e /q) 2 . Since we are interested in m X 6 MeV, however, the scattering is momentum-transfer independent, so that we can use the limit with a DM form factor, F DM = 1, in the literature. The most stringent limit is set by XENON1T (30 MeV m ψ ) [67], DarkSide50 (20 MeV m ψ 30 MeV) [71], XENON10 (10 MeV m ψ 20 MeV) [72], and SENSEI (m ψ 10 MeV) [73].

Indirect bounds
DM pair-annihilates into charged leptons and neutrinos. These annihilations during the cosmic dark ages and the BBN era have impacts on the ionization history and N eff , bringing us cosmological bounds. We will give brief comments on these bounds in the following.
A. CMB. Annihilation of DM into charged particles and photons increases the ionization fraction in the post-recombination era, modifying the CMB anisotropies and in turn providing the strong limit on the cross section. According to refs. [74,75], we JHEP03(2021)047 adopt a conservative limit on the thermal averaged cross section into charged final states, (σv) charged /(2m ψ ) ≤ 5.1 × 10 −27 cm 3 s −1 GeV −1 . 1 It suggests that there is a lower mass bound m ψ 5 GeV if the annihilation into charged particles is only responsible for the thermal production.
Let us see the impact of the CMB limit on the model. For m ψ ≤ m X , DM can only annihilate to muon, tau and neutrinos. Since the cross section for these three processes is almost the same, (σv) µμ,ττ 10 −26 cm 3 /s is suggested if kinematically possible, excluding the DM mass up to several GeV. Thus, the CMB observation totally excludes the region where m µ < m ψ .
For m ψ > m X , the annihilation is dominated by ψψ → XX followed by X → νν. The cross section for the annihilation into charged states is smaller by a factor of (1/q 2 ψ ) than ψψ → XX. In this case, there is less significant energy release from DM annihilation into electron and photon. The CMB bound is much weaker than that for m ψ ≤ m X . Nonetheless, the CMB observations partly limit the parameter space. We shall estimate it by considering m X m ψ for simplicity. The cross section of ψψ → µμ is related to that of ψψ → XX as The canonical cross section is (σv) XX /2 5 × 10 −26 cm 3 /s for Dirac DM [63,64] in the sub-GeV region where the direct detection bound will be avoided. Plugging the canonical value in eq. (3.7), the CMB bound reads It follows from this equation that q ψ 6.9 is excluded for m ψ = 200 MeV. We illustrate the CMB bound expressed by eq. (3.8) in figure 2 (right).
B. Effective neutrino number N eff . The annihilation into neutrinos reheats the SM plasma when DM becomes non-relativistic. If it occurs after the neutrino decoupling, only the neutrino is reheated and the neutrino-to-photon temperature ratio increases, resulting in the higher expansion rate at the BBN than the standard scenario. The expansion rate at the BBN, together with the baryon-to-photon ratio, affects the primordial abundances of helium and deuteron.
The change of the expansion rate is rendered in the effective neutrino number N eff . A significant deviation of N eff from the standard value N eff,SM = 3.046 [76] is faced with the precise Planck observations. In [77][78][79][80], the authors calculate the increase of N eff by a relic particle coupled to neutrinos, assuming that the particle is only in the equilibrium with the neutrinos during the BBN. They obtain the lower mass bound of Dirac DM m ψ 10 MeV, which is drawn in figure 2 with light blue.  Figure 2 shows the allowed region in the (m ψ , q ψ ) plane. The left panel corresponds to m ψ < m X , while the right panel to m ψ > m X . In both panels, we scan two parameters, m X and g X , such that the experimental and cosmological limits given in section 2 are avoided. The DM observed abundance, Ω CDM h 2 = 0.12, can be explained in all the regions in the plots, except for the dark gray region where Ω ψ h 2 > 0.12 is predicted. We also highlight the 1σ (2σ) muon g − 2 favored region with the (light) red band. In figure 2 (left) with m ψ < m X , the muon g − 2 explanation and DM production are simultaneously achieved in a wide DM mass range. In particular, we find that the non-resonant DM production is realized for q ψ O(10), otherwise the large resonant enhancement is needed. We also show the excluded regions from the CMB (green) and N eff (light blue). These restrict the DM mass to be 10 MeV m ψ 100 MeV. Note that the BABAR limit (blue) depends on q ψ for m ψ < m X /2. We estimate this limit by rescaling the announced BABAR limit [52] with the branching ratio: g X → g X · {Br X→µμ | q ψ /Br X→µμ | q ψ =0 } 1/2 . Naively speaking, the limit is weakened by a factor of q ψ , compared to the announced one. The BABAR limit excludes the region of q ψ 5, although it is overlapped with the CMB limit.

Allowed region
It is obvious from figure 2 (right) that allowing the large q ψ opens a new mass regime, m ψ > m X , that has not been pointed out in the previous study. The heavy DM mass, even O(100 GeV), is allowed in this regime. The leading constraints in the heavy mass region are set by the direct detection experiments. Given the current results, the DM mass should be lighter than 2 GeV to resolve the muon g − 2 anomaly. In contrast, if we do not require the muon g−2 explanation, the direct detection limit can be very weak, since we can JHEP03(2021)047 consider the heavy X boson outside the light red band in figure 1. The gain of the X boson mass does not affect the DM production much, while it considerably reduces the elastic scattering with nuclei and electrons, thereby weakening the direct detection bound. 2 As far as the indirect bounds are concerned, the N eff bound disfavors m ψ 10 MeV, while the CMB gives a complementary constraint in the O(100 MeV) region, but only a small part is covered yet. Altogether, it is 10 MeV m ψ 2 GeV where the muon g − 2 anomaly and DM are both addressed. One may wonder if the large DM charge will generate a sizable Sommerfeld enhancement at the freeze-out or in indirect searches. We have confirmed, however, that there is no sizable Sommerfeld effect in the parameters space we consider.
It has turned out so far that allowing the large DM charge opens two new DM production regimes; (i) the non-resonant production for m ψ < m X and (ii) the secluded production for m ψ > m X . These are overlooked possibilities in the previous study. Given the new phenomenological possibilities, it is natural to explore the signals of DM in this new parameter space. In the next section, we will study the indirect signals of this DM candidate at neutrino telescope experiments.

Signature with neutrino telescope
We have seen that the large DM charge revives the successful thermal production without the help of the resonance enhancement. As a by-product of the revival, the model predicts an excess of the neutrino flux from the galactic space. In this section, we discuss the sensitivity of the SK to the neutrino flux excess and the future prospect at the HK.
There are some dedicated studies of neutrino signatures of annihilating DM in the framework of simplified models [36][37][38][39][40][41][42][43]. In these works, it is assumed that DM directly annihilates into a pair of neutrinos. The resulting monochromatic neutrino fluxes are compared with the results of searches for unknown extraterrestrial neutrino sources [82,83], SRN searches [29,30,84] and atmospheric neutrino measurements [85], to derive upper limits on the annihilation cross sections. In the following, we perform an independent analysis for the DM signals measured at the SK, assuming that the extra neutrino flux originates from the direct annihilation ψψ → νν or the secluded annihilation ψψ → XX → 2ν2ν in the Milky Way. We also note that the analysis in this section is applicable as far as the event topologies are the same only by changing the model dependent flux normalization.

Neutrino flux from DM annihilation
We here assume that neutrinos are Dirac particles for definiteness, but the result will be unchanged even if these are Majorana particles. The expected electron (anti-)neutrino flux at the detector from DM pair annihilation in the galactic halo is expressed by

JHEP03(2021)047
where σv i denotes the annihilation cross section into a final state i, κ is a model dependent constant which characterizes electron-neutrino flavor fraction, dN i /dE ν the neutrino spectral function for the final state i, J ∆Ω the astrophysical J-factor which is given by the line-of-sight (l.o.s) integral of the DM density square. In the (b, l) coordinate, the J-factor is expressed by where ρ(r) is the DM density profile in the galactic halo. The galactocentric distance r is expressed in terms of the l.o.s distance s, with r = 8.5 kpc being the distance between the galactic center (GC) and the solar system. We define the direction of the GC as b = 0 • and l = 0 • . The size of the l.o.s distance is restricted by the DM halo size, whose maximum is given by We take the DM halo size R halo = 40 kpc. Note that the J-factor is insensitive to the halo size parameter as far as R halo ≥ 30 kpc. Since the DM profile is spherically symmetric, we define the range of (b, l) as 0 • ≤ b ≤ 90 • and 0 • ≤ l ≤ 180 • and then multiply the result by four. We consider the Navarro-Frenk-White (NFW) profile [86,87], ρ(r) = ρ s r s r where r s = 20 kpc is the scale radius and the density profile is normalized with ρ(r = r ) = 0.4 GeV/cm 3 . The resulting all-sky J-factor is J ∆Ω 1.5 × 10 23 GeV 2 /cm 5 . As an example for calculating the model dependent constant κ, if the neutrino flavor ratio at source is ν e : ν µ : ν τ = 0 : 1 : 1, then the flavor ratio at the detector is 1 : 2 : 2. Thus, κ = 1/5 is obtained for the U(1) Lµ−Lτ DM model.
The neutrino flux exhibits a distinctive feature. There are two important annihilation modes, which are ψψ → νν, XX. The former is dominant for m ψ ≤ m X , while the latter for m ψ > m X . In the former annihilation, the neutrino spectrum is monochromatic at E ν = m ψ in the center-of-mass frame of the annihilation. Since DM is almost at rest in the galactic halo, the monochromatic feature is maintained in the lab frame that is the rest frame of the DM halo. When the annihilation occurs in the Milky Way, the redshift is negligible. Therefore, the neutrino spectrum at the detector keeps the monochromatic form, dN/dE ν ∝ δ(E ν − m ψ ). This spiky spectrum is easy to be disentangled from the smooth background spectrum, providing a good sensitivity.
In the secluded annihilation ψψ → XX → 2ν2ν, the neutrino spectrum takes another characteristic form. The spectral feature from such a cascade annihilation has been studied in [28] for an intermediate state with an arbitrary spin, and the authors have shown that the spectral shape in general depends on the polarization of the intermediate state. In

JHEP03(2021)047
what follows, we restrict ourselves to an annihilation process ψψ → V V → 2ν2ν with a vector boson V (not necessarily be the X boson), and briefly review the spectral form. In the rest frame of the decaying V , each of emitted neutrinos has a monochromatic spectrum at E ν = m V . In the lab frame, the neutrino energy is boosted and reads where r V = m V /m ψ and θ is the angle between momenta of the parent V boson and the emitted neutrino in the lab frame. There is a sharp kinematical cut of the neutrino energy. The maximum (minimum) neutrino energy corresponds to θ = 0 • (180 • ). Since the emission angle θ determines the neutrino energy in the lab frame, we can obtain the neutrino energy spectrum once we know the angular distribution of the neutrino emitted by the V boson decay. If the V boson produced in ψψ → V V is unpolarized, the neutrino emission is isotropic and has no angular distribution. In this case, the neutrino spectrum has no energy dependence between the kinematical endpoints, i.e. it takes the so-called box shape [88,89], where y = E ν /m ψ and y ± = 1 ± 1 − r 2 V /2 and Θ(y) is the Heaviside theta function. On the other hand, if the produced V boson is polarized, the neutrino emission has an angular distribution and hence the resulting neutrino spectrum has an energy dependence.
In general, the neutrino spectrum in the ψψ → V V → 2ν2ν process is expressed by where the summation runs over V boson helicity m, n = +1, 0, −1 and Br m,n denotes the production branching fraction for each helicity state, . (4.9) The spectral functions f m (y) are in the form of where C m are model-dependent coefficients that characterize how the intermediate state with the different polarizations couples to the decay products. For example, when the

JHEP03(2021)047
vector boson couples to the final state fermions as f (g R γ µ P R + g L γ µ P L ) f V µ , we find in the leading order in m f /m V or m f /m V [28]. Here, we would like to mention the relation, It follows from this relation that if the branching fraction is helicity-independent, i.e. the V boson is unpolarized, then we exactly obtain the box-shape spectrum eq. (4.7). Now let us return to our model. The X boson mostly couples to the left-handed neutrinos which can be regarded as massless, so that C 0 C +1 0 and C −1 1. Further, the polarization of the X boson produced in the annihilation ψψ → X m X n is purely transverse, Br = diag(1/2, 0, 1/2). (4.16) As a result, the neutrino spectrum is in the form of a bowl-shape, In figure 3, we show the neutrino spectrum eq. (4.17), which is symmetric with respect to E ν = m ψ /2. The width of the spectrum is determined solely by kinematics and is given by 1 − r 2 V . As ψ and X are more degenerate, the width is narrower and the height is taller. We would like to remark on the DM model dependence of the neutrino spectrum from the secluded annihilation. For example, if DM is a complex scalar charged under the U(1) Lµ−Lτ gauge group as considered in [90], the annihilation into the longitudinally polarized X boson is nonvanishing and the branching fraction depends on the mass degeneracy of DM and X, With a small r X , the annihilation into the longitudinal polarization is sub-leading and the spectral form is similar to the bowl-shape, while all polarizations contribute equally with r X 1 and the spectrum gets close to the box-shape. Moreover, the annihilation into the longitudinal polarization can be even dominant if DM is sizably coupled to the Nambu-Goldstone mode of the U(1) Lµ−Lτ Higgs field. In this case, the spectrum is like an upside-down bowl. With the current and even future-planned precision, however, it will be difficult to distinguish the difference and these spectra will look practically the same shape in experiments. It might be interesting if any effort in future would realize good enough capability to discern the difference of the spectra.

Analysis
We explain our analytical method. We try to reinterpret the result of the 2,853 days SK search for SRN [30], to derive the limits on the DM annihilation cross section. The signal in this experiment is an electron and a positron produced in the reaction of the inverse beta decay (IBD) (ν e + p → e + + n) and the neutrino absorption by Oxygen in the charged current (CC) interactions (ν e (ν e ) + 16 O → e − (e + ) + X A ). The energy of the produced electron/positron is related to the initial neutrino energy, E e E ν − 1.3 MeV (ν e p), E e E ν − 15.4 MeV (ν e O) and E e E ν − 11.4 MeV (ν e O). If no excess of the recoil events is found, one can derive the limit on the neutrino flux, which is in turn translated into the bound on the annihilation cross section.
The signal region is set to E e = 16-88 MeV which is divided into 18 bins with a 4 MeV width. The number of the signal events measured in the i-th bin is expressed by where N SK = 1.5 × 10 33 denotes the number of free protons in the SK detector, T SK the exposure time of SK. We combine SK-I (1,497 days), SK-II (794 days) and SK-III (562 days) data. To take into account the detector efficiency and the finite energy resolution, we introduce the efficiency function, (E vis ) (figure 10 of [30]), and the Gaussian-like resolution function,  with the width σ(E e ) = 0.4 MeV E e /MeV + 0.03E e [37]. Here, E e is the actual electron/positron energy, E vis the measured energy at detector and dΦ νe /dE ν the neutrino flux originated from the galactic DM annihilation in eq. (4.1). We use the NLO analytical expression for the IBD cross section, dσ νep /dE e , in [91] and extract dσ νeO /dE e and dσ νeO /dE e from [92,93]. For the latter, we assume the electron energy is mono-energetic, e.g. dσ νeO /dE ν (E ν , E e ) = σ νeO (E ν ) δ(E e − E ν − 15.4 MeV), and determine σ νeO (E ν ) from figure 2 of [93].
There are four main backgrounds in the signal region E e = [16,88] MeV. The first one is an electron from the Michel decay of a muon that is produced by the CC interaction of the atmospheric muon neutrino in the water of the detector. If the momentum of the produced muon is lower than the Cherenkov threshold, this muon is invisible and its decay electron cannot be removed. This is called the invisible muon background. This has a kinematical edge at E e 60 MeV. The second one comes from the CC interaction of the atmospheric electron neutrino, producing the monotonically increasing events with the increasing E e . This dominates the background above the kinematical endpoint of the invisible muon. The third one is an electron produced via the neutral current (NC) interaction of the atmospheric neutrinos. This is increased at the low energy bins. The last grouping of the background is due to heavy charged particles, pions and muons, created in the NC reactions. Some of them survive the pion and Cherenkov angle cuts, and enter the signal region. These backgrounds are shown in figure 4.

JHEP03(2021)047
To derive the 90% confidence level (C.L.) limit, we first introduce a likelihood function, for each bin and each SK phase. With this likelihood, we define test statistic (TS) as where the summation runs over 18 energy bins and from SK-I to SK-III. The 90% C.L. limit is obtained by solving TS ≤ 2.71. Since TS is a function of the DM mass and annihilation cross section TS(m ψ , σv), we obtain an exclusion curve in the (m ψ , σv) plane. To estimate the HK sensitivity, we consider the 374 kton fiducial volume and 10 yrs livetime with the same efficiency and energy resolution as the SK. The HK detector has an option of doping Gd to reduce the backgrounds. In our analysis, it is assumed that the doped Gd reduces the invisible muon background by 50% or 80%. The sensitivity curves are obtained by solving TS = 2.71 with N obs = N bkg .

SK bounds and HK sensitivity
We compare the DM signal events with the observed events at the SK-I (1,497 days) in figure 4. The left panel corresponds to ψψ → νν and the right panel to ψψ → XX with m X /m ψ = 0.5. We assume 100 % branching fraction for the corresponding annihilation and fix the cross section such that TS = 2.71. The red and blue curves represent the number of events in presence of the DM annihilation. The dashed one is the signal and the solid one is the sum of the signal and backgrounds. The monochromatic flux predicts the distinct event shape, which is peaked at E e m ψ for 50 MeV DM. The shape is broader for 70 MeV DM, on the other hand. This is because there are three reactions that produce the electron/positron with different energy: the IBD with E e E ν − 1.3 MeV and ν e (ν e ) + 16 O scattering with E e E ν − 15.4 (11.4) MeV. The cross sections for the latter two scatterings are still smaller at 50 MeV than the former, while these become comparable at 70 MeV. Thus, another peak appears at an energy 10 MeV below the DM mass. The resulting event shape takes the flatter form. The bowlshape flux induced by the secluded annihilation produces the broader event shape. The event excess in each bin is not sizable, but it equally contributes to multiple bins, not losing the sensitivity. The center of the neutrino spectrum is at half of DM mass. While, the peak appears at slightly higher energy because the cross section grows as the energy.
We show the 90% C.L. upper limits on the annihilation cross section and the future HK sensitivity in figure 5. For the ψψ → νν annihilation, we illustrate the similar limits derived in the previous works [40,41] for comparison. Our SK limit is consistent with their results below m ψ = 50 MeV, while it is a little aggressive for higher mass. We guess the disagreement originates from the difference in modeling the ν e O cross section. We just assume the mono-energetic electron in the reaction, but in general the electron has a broad energy distribution. This will weaken the signal strength and decrease the sensitivity. (σv) [cm 3 /s] Figure 5. The 90% C.L. upper limits on and the future HK sensitivity to the annihilation cross section for ψψ → νν (left) and ψψ → XX (right) cases. Black solid (dashed) line corresponds to the limit on the annihilation cross section from SK (HK). The results given in [40,41] are shown for reference, and are shown in turquoise and brown lines.
For the secluded annihilation, the limits depend on the mass degeneracy. When ψ and X are highly degenerate (m X /m ψ = 0.99), the spectrum has a very narrow width and looks the line shape within the detector resolution, so that the structure of the exclusion curve resembles the monochromatic one. On the other hand, the curve is stretched, which reflects the spectral shape that has the center at E ν = m ψ /2. As ψ and X are less degenerate, the curve is shifting to the lighter mass. It is interesting that the exclusion curves extend to the high mass region above 300 MeV due to the wide spectrum. The sensitivity in this mass region is rapidly decreasing as the lower energy cut of the neutrino flux, E − = y − m ψ , exceeds about 100 MeV. For instance, E − 0.28 m ψ with m X /m ψ = 0.9, so that m ψ ∼ 350 MeV is the boundary. In the high mass region, the sensitivity may be increased by combining other observations, such as measurements of the atmospheric neutrino flux by the SK [85]. In both annihilation modes, the current SK limits are much above the canonical thermal relic cross section, but the future HK sensitivity may reach down to it. We would like to add that ref. [41] also studies the reach of DUNE [94] and JUNO [95] experiments, that is more sensitive than the HK in the same mass range. These sensitivity plots are shown in figure 6 (left). We also note that the difference in modeling the ν e O cross section will affect the exclusion limits in the secluded annihilation as well. Nevertheless, we expect that the modeling error is not so large compared with the monochromatic case unless DM and X are highly degenerate. This is mainly because the neutrino flux is basically broad in the secluded annihilation, in contrast to the monochromatic one, and hence the electron event shape is also broad without the electron energy distribution in the ν e O reaction. Thus, the exclusion limits will not largely be changed even if another modeling generates a broad electron distribution in the reaction. The detailed study of modeling the ν e O cross section and its potential error is beyond the scope of this paper and will be followed in a future work. (σv) [cm 3

JHEP03(2021)047
/s] (σv) [cm 3 /s] Figure 6. The predicted annihilation cross section for ψψ → νν (left) and ψψ → XX (right) cases. The magenta and orange points correspond to m ψ < m X and m ψ > m X , respectively. All these points can explain the observed DM abundance and the muon g − 2 as well as being consistent with the experimental constraints in figure 2. The various sensitivity curves presented in the literature are shown, in addition to our results (black lines with the ♦ mark). The lines with the same mark (♣, ♥) are extracted from the same paper: ♣ Ando et al. [41], ♥ Arguelles et al. [42] and SK20 [34].
To avoid being messy, we omit the similar lines shown in figure 5.
It may be useful to comment on how large impact the polarizations of the intermediate state have on the cross section limits in figure 5 (right). If the intermediate state in the secluded annihilation is an unpolarized vector or a scalar, the neutrino flux is box-shape. We have analyzed the neutrino signature and calculated the upper limit on the cross section for such a neutrino spectrum. Comparing the results with those of the bowl-shape spectrum, we have found the difference is at most 10-20%. Since in general the neutrino spectrum is a medium of the box and bowl-shape, we expect the cross section limit to be comparable with the ones in figure 5 (right), independently of the annihilation fraction into each polarization state.

Implication for U(1) Lµ−Lτ DM
In the above subsections, we have formulated the analysis of neutrino signature and evaluated the experimental reach to the extra neutrino flux from the DM annihilation with neutrino telescopes. In this subsection, we discuss the impact on the U(1) Lµ−Lτ DM.
In figure 6, we plot the predictions of the annihilation cross section for the parameter points that can explain the observed DM abundance and the muon g − 2 as well as being consistent with the experimental constraints in figure 2. The magenta and orange points give the cross sections for m ψ < m X and m ψ > m X , respectively. The velocity of DM is assumed to be a typical virial velocity in the galaxy, v rel ∼ 10 −3 . We also summarize various limits on the cross section derived in the literature, in addition to our results (black lines with the ♦ mark). To avoid the messy figure, however, we omit the similar limits that JHEP03(2021)047 are already shown in figure 5. The lines with the same mark (♣, ♥) are extracted from the same paper: ♣ Ando et al. [41], ♥ Arguelles et al. [42] and SK20 [34]. In referring to these limits, we appropriately adjust the neutrino flavor fraction to κ = 1/5. We see in figure 6 (left) that for m ψ < m X , the model predicts so large cross section that the future experiments can reach. For some parameter points, the cross section is significantly boosted to have even σv = 10 −22 cm 3 /s. Such parameter points correspond to the resonant mass region m ψ m X /2. If the DM mass is fine-tuned to such values at or below % level, the thermal history of the DM production is modified to delay the freeze-out. Such a delay requires the larger cross section than the standard case to deplete the number density down to the observed value. The mechanism for boosting the cross section in this way is discussed in [96,97]. We review the detail of this mechanism in the appendix. It is interesting that some of such resonant parameter sets have already been excluded by the SK data. When the mass tuning to m ψ m X /2 is moderate, on the other hand, suppressed cross sections can also be obtained. These are distributed over σv O(10 −26 ) cm 3 /s in the figure (magenta points). The suppression is caused if the cross section is highly enhanced at the typical freeze-out time due to the physical X boson resonance. In this case, the small DM charge is enough to thermally produce the DM abundance. As a result, the late-time annihilation, e.g. in the galaxy, is suppressed because the DM velocity in the galaxy is too small to produce the X boson resonance. This region corresponds mainly to q ψ = O(1). In the non-resonant region, the DM thermal production works in the standard way, so that we have an almost constant canonical cross section, σv 10 −25 cm 3 /s. None of the current experiments can reach the canonical value, but the future experiments including DUNE and JUNO will be able to cover the DM mass of 20-100 MeV.
In the secluded region (m ψ > m X ), the direct annihilation into a neutrino pair is suppressed and the size of the cross section is at most ∼ 10 −26 cm 3 /s. This tendency is more pronounced as DM is heavier. In this case, however, the secluded annihilation can be large (see figure 6 (right)). We see the model predictions distributed slightly below the canonical value, σv = 10 −25 cm 3 /s. There is no enhancement in this mass region. The cross section can be kinematically suppressed in the mass degenerate case, because the phase-space of the produced X boson becomes small for the low DM velocity. Thus, a moderate mass splitting is favored to observe the signal, although the high degeneracy produces the sharp neutrino spectrum and leads the strong limits. Indeed, we see some parameter points lying much below the canonical value, where DM and the X boson have close mass. The current SK limit is over one order of magnitude larger than the predicted cross section. The estimate of the future HK reach is above the prediction by a factor of 3 or more. We hope that some future updates or improvements of the analysis will grow the experimental sensitivity and reach the thermal relic cross section. Compared the HK sensitivity with the DUNE and JUNO ones in figure 6 (left), the latter ones have better sensitivity. If this is the case for the secluded annihilation, it will be important to analyze the similar signals at these neutrino observatories, that will be pursued in future.
It is worth mentioning other possibilities that predict significant neutrino flux from sub-GeV DM. Indeed, a great deal of effort has been devoted to non-renormalizable DM-JHEP03(2021)047 neutrino interactions. The prime difficulty in embedding them in a renormalizable model is that the model includes interactions with charged leptons, which are isospin partners of neutrinos. Such interactions make DM or a mediator particle visible in terrestrial experiments and cosmological observations and, hence, will strongly limit the categories that achieve the sizable DM-neutrino interactions in a renormalizable manner. The model considered here is one of the simplest renormalizable models, that evades the experimental constraints and does not suffer from theoretical requirements, such as gauge anomalies. The other type of possibility is realized by employing a t-channel mediator. Examples of feasible renormalizable models incorporating a t-channel mediator include [98][99][100][101].

Summary and discussions
We have examined a simple Dirac DM model based on the U(1) Lµ−Lτ gauge theory. In this model, DM ψ interacts with the SM particles only through the X boson. DM is produced thermally in the early Universe and annihilates into the SM particles through two kinds of processes. One is the s-channel ψψ → ff which is important for m X > m ψ . The other is the t-channel ψψ → XX. This annihilation process is allowed if m X < m ψ , and becomes significant for a large U(1) Lµ−Lτ charge q ψ of DM. In both cases, the annihilation cross sections can be large enough even at off-resonance of X boson mass pole. We have shown that the observed relic abundance of DM and the discrepancy in the muon anomalous magnetic moment are explained simultaneously by the U(1) Lµ−Lτ gauge boson X without conflicting the severe experimental bounds thanks to no interaction with SM quarks and electron. As a by-product of a large q ψ scenario, DM annihilations provide characteristic neutrino signatures at SK and HK.
In this paper, we have formulated the analysis of the indirect neutrino signal of DM in a model independent way, and applied it to U(1) Lµ−Lτ DM model. From the ψψ → νν process, neutrinos with monochromatic energy are predicted, while from the ψψ → XX process followed by X → νν, the energy spectrum of neutrinos is bowl-shape. Neutrinos produced in the DM halo can be detected by neutrino telescopes. We have calculated the number of the signal events expected at SK and HK detectors, and estimated the upper limits from the future HK sensitivity to the DM annihilation cross section. As we have shown in figure 6, the future sensitivity to the annihilation cross section at HK almost reaches the canonical thermal relic cross section in the ψψ → νν case. On the other hand, in the ψψ → XX → 2ν2ν case, the obtained future sensitivity is several times larger than the canonical one. Further improvement of the experimental sensitivity is necessary in order to cover the wide range of the model predictions. We hope that the analysis with directional information may help the background subtraction, which is beyond the scope of this paper. It will be shown elsewhere.

JHEP03(2021)047
Indirect detection experiments observe the neutrinos from DM annihilation in our Galaxy. Since the DM velocity is v rel ∼ 10 −3 at most, we can approximate (σv) νν σv T =0 for e.g. δ = 10 −3 . In this case, the boost factor is expected to have Before closing, we derive eq. (A.4). In the standard way, the approximate solution of the Boltzmann equation is obtained by integrating the equation, , where x f is determined by appropriately matching the approximate solution with the actual one. Above, we assumed the s-wave DM annihilation and introduced Y DM = n DM /s with s being the entropy density and Taking into account the resonant behavior of the cross section, the equation is modified as One can readily get eq. (A.4) by solving this equation.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.