Indirect Detection of Neutrino Portal Dark Matter

We investigate the feasibility of the indirect detection of dark matter in a simple model using the neutrino portal. The model is very economical, with right-handed neutrinos generating neutrino masses through the Type-I seesaw mechanism and simultaneously mediating interactions with dark matter. Given the small neutrino Yukawa couplings expected in a Type-I seesaw, direct detection and accelerator probes of dark matter in this scenario are challenging. However, dark matter can efficiently annihilate to right-handed neutrinos, which then decay via active-sterile mixing through the weak interactions, leading to a variety of indirect astronomical signatures. We derive the existing constraints on this scenario from Planck cosmic microwave background measurements, Fermi dwarf spheroidal galaxies and Galactic Center gamma-rays observations, and Alpha Magnetic Spectrometer - 02 antiprotons observations, and also discuss the future prospects of Fermi and the Cherenkov Telescope Array. Thermal annihilation rates are already being probed for dark matter lighter than about 50 GeV, and this can be extended to dark matter masses of 100 GeV and beyond in the future. This scenario can also provide a dark matter interpretation of the Fermi Galactic Center gamma ray excess, and we confront this interpretation with other indirect constraints. Finally we discuss some of the exciting implications of extensions of the minimal model with large neutrino Yukawa couplings and Higgs portal couplings.

A wide array of gravitational phenomena over a range of cosmological scales strongly supports the hypothesis of dark matter (DM) [1][2][3]. There is, however, no firm evidence that DM couples to ordinary matter other than through gravity, and the search for such non-gravitational DM interactions has become one of the main drivers in particle physics today. Neutrinos (ν) in the Standard Model (SM) may be identified as a component of DM, since they are color-singlet, electrically neutral cosmic relics. However, the smallness of the lightest neutrino mass makes them relativistic at freeze-out in the early universe, and thus incompatible with current observations to account for the majority of the cold DM. One therefore must seek a solution beyond the SM. Since we do not know how DM couples (if at all) to the SM, it is important to explore a variety of models to understand in a comprehensive manner how non-gravitational DM interactions may manifest [4].
Since DM is presumably electrically neutral, it may be either the neutral component of an electroweak multiplet, as in the well motivated weakly interacting massive particle (WIMP) paradigm, or alternatively it may be a Standard Model (SM) gauge singlet state. In the latter case of gauge singlet DM, an economical and predictive mechanism for mediating DM interactions to the SM is provided by the so-called "portals"− renormalizable interactions of DM through gauge singlet SM operators. There are only three such portals in the SM−the Higgs portal [5,6], the vector portal [7,8], and the neutrino portal [9]. As applied to DM, the Higgs portal [5,[10][11][12], and the vector portal [13][14][15] have been extensively investigated, while the neutrino portal option has received comparatively little attention, despite the strong motivation due to its connection to neutrino masses. In this paper we will examine a minimal model of neutrino portal DM in the simplest setup of a Type-I seesaw scenario [9].
The neutrino portal to DM relies on DM interactions being mediated by the right-handed neutrinos (RHNs). Since the RHNs are responsible for generating neutrino masses, one may typically expect the DM interaction strength with the SM to be very small since it is governed by the neutrino Yukawa coupling. In this case it is challenging to probe neutrino portal DM in accelerator experiments or in direct detection experiments. On the other hand, the DM coupling to the RHN can be sizable, thereby facilitating the efficient annihilation of DM to pairs of RHNs. This allows DM to be produced thermally in the early universe with the observed relic abundance and furthermore presents an opportunity to test the scenario through a variety of indirect detection channels. In this work we investigate the indirect detection signatures of neutrino portal DM. The scenario investigated here was first proposed in Ref. [13] and falls into the class of "secluded" DM scenarios. Some aspects of the thermal cosmology were investigated in Ref. [16]. In regards to indirect detection signatures, Ref. [17] explored a possible interpretation of the Fermi Galactic Center gamma ray excess [18][19][20][21][22] in terms of the DM annihilation to RHNs. Recently, Ref. [23] investigated the limits from gamma ray observations on DM annihilation to RHNs, although did not explore the implications for specific particle physics models. Extensions of the simplest scenario, which include additional states and/or interactions have also been discussed in Refs. [24][25][26][27][28][29][30][31][32]. Our work provides a comprehensive and updated analysis of the indirect detection phenomenology of neutrino portal DM. In particular, we present constraints from Planck cosmic microwave background (CMB) measurements, Fermi dwarf spheroidal galaxies and Galactic Center gamma-rays studies, and AMS-02 antiproton observations, and also describe the future prospects for Fermi and the Cherenkov Telescope Array. Thermal relic annihilation rates are already constrained for DM masses below about 50 GeV. This scenario can also provide a DM interpretation of the Fermi Galactic Center gamma ray excess, although we demonstrate that such an interpretation faces some tension from dSphs and antiproton constraints. We also describe extensions of this scenario beyond the minimal model, including scenarios with large Yukawa and Higgs portal couplings, and highlight the potentially rich physics implications in cosmology, direct detection, and collider experiments. The outline of the paper is as follows. In Section II we describe a minimal neutrino portal DM model, outline the expected range of couplings and masses, and discuss the cosmology.
The primary analysis and results concerning the indirect detection limits and prospects are discussed in Section III. In Section V we describe several features and phenomenological opportunities present in non-minimal neutrino portal DM scenarios. Our conclusions are presented in Section VI.

II. NEUTRINO PORTAL DARK MATTER
The simplest construction beyond the Standard Model to account for the neutrino masses is the introduction of right-handed neutrinos (RHN). Beside the normal Dirac mass terms with the Yukawa interactions, the RHN can also have a Majorana mass term since it is a SM gauge singlet. This is the traditional Type-I seesaw mechanism [9]. For the same reason of its singlet nature, N can serve as a mediator to the dark sector via the neutrino portal. A simple model of neutrino portal DM based on the Type-I seesaw contains three new fields, N, χ, φ, where N and χ are two component Weyl fermions and φ is a real scalar field. They are charge-neutral with respect to the SM gauge interactions. The fermion N is identified as a RHN. We will assume that χ is lighter than φ, and they are charged under a Z 2 symmetry, which renders χ stable and a potential DM candidate.
The Lagrangian has the following new mass terms and Yukawa interactions where L and H are the SM SU(2) L lepton and Higgs doublets, respectively. There are two central features of this model. First, the RHN field N serves as a mediator between the dark sector fields χ, φ and the SM fields, due to the couplings λ and y. This mediation allows for non-gravitational signatures of the DM and a thermal DM cosmology. Second, after the Higgs obtains a vacuum expectation value, H = v/ √ 2 with v = 246 GeV, a small mass for the light SM-like neutrinos is generated via the seesaw mechanism: Given the observed neutrino masses 1 , the Yukawa coupling y depends on the RHN mass, m N . For instance, fixing m ν ∼ (∆m ν ) atm ∼ 0.05 eV suggsts a small neutrino Yukawa coupling of order As we will discuss in more detail shortly, the requirement of thermal freeze-out of the DM puts an upper bound on the DM and RHN mass less than 20 TeV. Therefore, the Yukawa couplings that we will be interested in will generally be quite small. It will thus be extremely difficult to produce the DM at accelerators, or directly detect it through its scattering with SM particles. However, there is an opportunity to probe this type of DM via indirect detection, and this will be the primary focus of this paper.
As alluded to already we will be interested in DM that is thermally produced in the early universe. The RHN mediator allows for the dark sector to couple to the SM thermal bath in the early universe. Then, provided that m χ > m N and that all of the particles are sufficiently light, say below O(10 TeV), the DM can efficiently annihilate to RHNs, and achieve the correct relic abundance. The process Eq. (4) is governed by the coupling λ, which is a priori a free parameter. The thermally averaged annihilation cross section is We observe that the annihilation cross section Eq. (5) depends on the coupling λ and the masses m χ , m N , m φ . However, the indirect detection signatures that we will investigate will depend in a detailed way only on the size of the annihilation cross section σv , which determines the rate, as well as the masses m χ and m N , which will affect the energy spectrum of the SM annihilation products. Thus, it will be more convenient to simply work with the three parameters { σv , m χ , m N }. Note that for a given set of masses, one can always obtain the desired cross section by an appropriate choice of the coupling λ through Eq. (5), provided the coupling remains perturbative. We will discuss this point in detail shortly.
We can restrict the parameter space further if we demand that the DM saturates the observed relic density. For Majorana fermion DM the observed relic abundance is obtained for [33] σv thermal = 2.2 × 10 −26 cm 3 s −1 .
Once we fix the annihilation cross section to saturate the observed relic abundance, then all of the physics can be characterized in terms of the two masses m χ and m N . Parameter choices that predict cross sections smaller than (6) overproduce the DM.
We now discuss the expected range of masses and couplings of the new states in the model. A first constraint comes from demanding that the coupling λ be perturbative and thus the theory be predictive. Assuming m N m χ , the partial-wave perturbative unitarity bound for the DM annihilation amplitude requires that λ < √ 4π. The over-closure and perturbative unitarity constraints lead to the bound m χ π 4 σv thermal ≈ 20 TeV, which is in broad agreement with the general analysis of Ref. [34]. Furthermore, there are a variety of limits on the right-handed neutrinos N , which depend on its mass and mixing angle with active neutrinos. In particular, for seesaw motivated mixing angles, the lifetime of N is typically longer than O(1 s) for m N 1 GeV, and is thus constrained by Big Bang Nucleosynthesis [35,36]. Then, considering m χ > m N in order to obtain and efficient DM annihilation cross section we will consider in this paper masses in the range The discussion above assumes a standard thermal history for the DM particle χ, which relies on χ being in equilibrium with the plasma. Since the dark sector particles χ and φ have no direct couplings to the SM, it is the RHN that is ultimately responsible for keeping χ and φ in equilibrium. It is therefore important that N remain in equilibrium with the SM during the freezeout process. The relevant processes to consider are the decay and inverse decays of N to the SM. This question has been investigated recently in Ref. [17] 2 .
For Yukawa couplings dictated by the naive seesaw relation, these process are very efficient when m N m W , since N decays through a two body process. However, if N is light, m N m W , the three body decays of N become inefficient and N can fall out of equilibrium. As a consequence, this fact requires an annihilation cross section that is larger than the canonical thermal relic value by some order one factor in the early universe to efficiently deplete the χ abundance, as explored in detail in Ref. [17]. A detailed investigation of the cosmology is beyond the scope of this paper, but we will take the standard thermal value for the annihilation cross section as a motivated benchmark.
Besides the terms in Eq. (1), an additional Higgs portal coupling, φ 2 |H| 2 is allowed in the model. This interaction provides an alternative means to keep φ, χ, and N in thermal equilibrium with the SM. We will assume for now that this coupling is small so that the phenomenology is dictated by the minimal neutrino portal interaction. However, a large Higgs portal coupling can lead to a variety of interesting effects, and we will discuss this topic in Section V.
We now come to the main subject of this work: the constraints and prospects for indirect detection of neutrino portal DM. We will investigate several indirect signatures of DM annihilation in this scenario, including observations of the CMB, gamma rays, and antiprotons. For each of these indirect probes the relevant underlying reaction is DM annihilation to RHNs as in Eq. (4), followed by the weak decays of the RHNs to SM particles due to mixing. We will thererfore require the energy spectrum dN/dE per DM annihilation in the photon, electron and antiproton channels as an input to our further analysis below. To compute these spectra we first simulate the decay of RHNs to SM particles in the N -rest frame using MadGraph5 aMC@NLO [38] in conjunction with the SM HeavyN NLO model files [39,40].
These partonic events are then passed to Pythia 8 [41] for showering and hadronization, thereby yielding the prediction for the resulting photon, electron, and antiproton spectrum coming from the N decay, dN i /dE for i = γ, e − ,p. These events are then boosted to the DM rest frame (see, e.g., Refs. [42][43][44][45]) thus giving the prediction for the required spectrum in each channel. We note that spin correlations are not accounted for in our simulation, but these are expected to have only a modest effect on the broad continuum spectra of interest to us (see Ref. [44] for an explicit example where this expectation is borne out).
We display in Figure 1 examples of the predicted continuum γ-ray, electron, and antiproton spectra for ( GeV (dashed), 100 GeV (dotted). Here we have assumed that N couples solely to the first generation (electron-type) lepton doublet. In the case of the γ-ray and antiproton spectrum, one observes a broad spectrum that peaks in the O(10 GeV) range. The location of the peak is largely dictated by the DM mass, which controls the total injected energy. There is a mild sensitivity to the RHN mass, with harder spectra resulting from a larger mass gap between the DM and RHN. For the electron case, in addition to the continuum component, there is a hard component resulting from the primary N → W e decay, which is clearly seen in Figure 1.
In this work we will restrict to the case in which N couples to the electron-type lepton doublet, but it is worth commenting on the cases of couplings to muon and/or tau flavor. In these cases, we have checked that the continuum spectra is very similar to the electron-flavor case, as is expected since these particles dominantly originate from decay of the electroweak bosons. The primary difference for muon or tau-flavor couplings will be the absence of the hard electron component from the primary N decay. The electron spectrum will be used below as an input to the CMB bounds, so one may expect a mild difference in the resulting limits in the case of muon or tau flavor couplings.
We now present in turn the constraints on neutrino portal DM from the Planck cosmic microwave background measurements, Fermi observations of gamma rays from the Galactic Center and from dwarf spheroidal galaxies, and AMS-02 observations of antiprotons. A summary of these constraints, as well as a discussion of other indirect searches not considered here, and an analysis of the future prospects, is presented below in Section III E.
The relevant quantity of interest for DM annihilation during recombination is the energy absorbed by the plasma per unit volume per unit time at redshift z, where ρ c is the critical density of the Universe today and Ω χ is the DM density parameter today. Production of neutrinos as daughter particles and free-streaming of electrons and photons after creation until their energy is completely deposited into the intergalactic medium (IGM) (via photoionization, Coulomb scattering, Compton processes, bremsstrahlung and recombination) affect the the efficiency of energy deposition. This is accounted for in Eq. (10) by the efficiency factor, f (z), which gives the fraction of the injected energy that is deposited into the IGM at redshift z and depends on the spectrum of photons and electrons arising from DM annihilations. Furthermore, since the CMB data are sensitive to energy injection over a narrow range of redshift, i.e., 1000−600, f (z) can be well-approximated by a constant parameter f eff .
The additional energy injection from DM annihilation in Eq. (10) alters the free electron fraction (the abundance ratio of free electrons to hydrogen), which in turn affects the ionization history. These effects are quantitatively accounted for with new terms in the Boltzmann equation describing the evolution of the free electron fraction. The additional terms are added to the baseline ΛCDM code and used to derive limits on the energy release from DM annihilation. Planck sets a limit on the particle physics factors in Eq. (10) which is obtained from temperature and polarization data (TT,TE,EE+lowP) [50].
To apply the Planck constraints of Eq. (11) to the neutrino portal DM model, it remains to compute the efficiency factor f eff (m χ ) in our model. We use the results of Ref. [64], which provides f γ(e − ) eff (E) curves for photons and electrons to compute a weighted average with the photon/electron spectrum (dN/dE) γ,e − predicted in our model according to The photon and electron spectra for each DM and RHN mass point are computed with Monte Carlo simulation described at the beginning of this section and are displayed for a few benchmarks in Figure 1. Using these spectra and Eqs. (11) and (12) GeV. A small feature in the limit contour is apparent in the region near m W m N m Z .
This is a consequence of the dominance of the two body decay to N → W in this small mass window.

B. Gamma rays from the galactic center
One of the primary signatures of DM annihilation are high-energy gamma rays. In comparison to other cosmic ray signatures involving electrically charged particles, gamma rays are essentially unperturbed by magnetic fields and the astrophysical environment as they travel to us from their source, yielding information about both the energy and location of the underlying DM reaction. One can search for both gamma ray line signatures as well as a continuum signal. While a line signature is unfortunately not present in the neutrino portal DM model, there can be a distinct continuum gamma ray signal, and this will be the subject of investigation here. Significant advances in our study of the gamma-ray sky have been achieved over the past several years by the Fermi Gamma Ray Space Telescope, and ������ ��� ��% ���� ����� data from the Fermi collaboration can be used to probe DM annihilation over a wide range of models and DM masses. In this section we will consider gamma ray signatures from the center of the Milky Way. The Galactic Center has long been recognized as the brightest source of DM induced gamma rays, a consequence of its proximity and the rising DM density in this region. At the same time extracting a signal from this region is challenging due to significant and not well-understood astrophysical backgrounds. Below we will also investigate gamma ray signals from dwarf spheroidal galaxies, which provide a cleaner, albeit dimmer, source of gamma rays. The quantity of interest for gamma ray signals of DM annihilation is the gamma ray flux per unit energy per unit solid angle in a given direction, Φ γ (E,n), where E is the energy andn is a unit vector along the path of the line of sight. The gamma ray flux can be written The term in square brackets in Eq. (13) above depends only on the underlying particle physics properties of the DM model, including m χ , σv , and the spectrum of photons emitted per DM annihilation dN γ /dE. This spectrum is shown in Figure 1 for the channel χχ → N N for several choices of χ and N masses. The quantity J(n) in Eq. (13), also called the J-factor, depends only on astrophysics and involves an integral over the DM density profile ρ χ (r) that runs along the path of the line of sight defined byn: In practice, the J-factor is averaged over a particular region of interest relevant for the analysis. The J-factor depends sensitively on the DM distribution and can vary by several orders of magnitude depending on this assumption, which translates into a substantial uncertainty in the derived annihilation cross section limit. At present, there is no consensus on the expected DM halo profile. Cuspy profiles such as NFW [66,67] or Einasto [68] find support from N -body simulations [69,70]. These simulations only involve DM, and the inclusion of baryonic processes may significantly impact the shape of the profile, especially towards the inner region of the Milky Way. However, even the qualitative nature of the resulting DM distribution is a matter of debate, and it is possible that the resulting profile is either steepened [71][72][73][74] or flattened [75] due to baryonic effects. Besides the assumption of the DM distribution, a separate, smaller O(1) uncertainty arises from the overall normalization of the profile, which is fixed to match the local DM density ρ 0 [76].
The current situation regarding the observed gamma ray flux from the Galactic Center is somewhat murky. A number of analyses, starting from the works of Goodenough and Hooper [19,20] and culminating most recently in the Fermi analysis [18], have found a broad excess of gamma rays from the Galactic Center, which peaks in the 1 − 3 GeV range.
All analyses conclude that there is a highly statistically significant excess above the currently accepted diffuse background models (see for example Refs. [21,22]). However, the origin of these gamma rays is still not clear. While there has been a significant effort devoted to possible DM interpretations, recently it has been argued that the excess is more likely to be a new population of unresolved point sources, which would disfavor the simplest DM interpretations [77][78][79][80] (see however [81]). It is certainly interesting to speculate on a possible DM origin, and we will carry out this exercise below in Section IV. Here we will instead take a conservative approach and use the Fermi data to place limits on DM annihilation.
To obtain limits on the neutrino portal DM scenario, we use the model independent results of Ref. [82]. In that work, four years of data from the Fermi Large Area Telescope was used to construct maps of the gamma ray flux in the region around the Galactic Center in four energy bins in the range from 300 MeV−100 GeV. Backgrounds templates from known point sources and emission from the Galactic Disk are then subtracted to yield the residual ����� �� ��% ���� ����� flux. Assuming that DM annihilation accounts for the remaining emission, the authors then place limits on DM annihilation for several choices of halo profiles. This procedure yields conservative limits since it is expected that additional background sources, such as the central supermassive black hole, unresolved point sources, and cosmic ray interactions with the gas, also contribute significantly to the residual emission. Limits on the the particle physics factor that governs the gamma ray flux, ( σv /m 2 χ ) dE dN γ /dE, are provided in Ref. [82].
For the neutrino portal DM model, we can use these results to derive a limit on the annihilation cross section for the process χχ → N N as a function of the DM and RHN mass. In Figure 3 we show contours of the 95% C.L. upper limit on the annihilation cross section in the m χ − m N plane labelled by the black curves. These limits are derived under the assumption of an NFW profile and local DM density ρ 0 = 0.3 GeV cm −3 . We see that under these assumptions, the Fermi data probes the thermal relic cross sections of Eq. (6) for m χ 10 GeV (thick red contour). The constraints on the annihilation cross section are again translated to limits on the minimum value of the coupling constant λ as shown by the vertical (blue) lines. The shaded (blue) region indicates the perturbative unitarity bound.
However, we again emphasize that there are significant uncertainties associated with halo profile, and the limits will become stronger (weaker) by a factor of a few to 10 (depending of course on the detailed shape) if one assumes a contracted (cored) DM distribution [82].
We observe a small feature near m W m N m Z where the two body decay N → W dominates.
C. Gamma rays from dwarf spheroidal galaxies for kth dwarf and jth energy bin. For each dwarf and energy bin, F ermi provides the likelihood, L k,j as a function of ϕ k,j . The likelihood function accounts for instrument performance, the observed counts, exposure, and background fluxes. For a given DM annihilation channel, the energy flux depends on m χ , σv , and J k (the J-factor of the dSph -see Eq. (14)) according to Eqs. (13,14,15), i.e., ϕ k,j = ϕ k,j (m χ , σv , J k ). The likelihood for a given dwarf, L k , is where LN accounts for statistical uncertainty in the J-factor determination (from the stellar kinematics in the dSphs), incorporated as a nuisance parameter in the likelihood. The where J k is the true value of the J-factor andJ k is the measured J-factor with error σ k on the quantity log 10J k . The combined likelihood for all the dwarfs is then where {J i } is the set of J-factors.
Given that no significant excess is observed, a delta-log-likelihood method is used to set limits on DM model parameters, treating the J-factors as nuisance parameters. The delta-log-likelihood ∆ ln L is given by where  lines and the associated numbers show the limits on the minimum value of the coupling constant λ. The shaded (blue) region indicates the perturbative unitarity bound. In the region m W m N m Z the two body decay N → W opens up and saturates the branching ratio, which is clearly seen in Figure 4.

D. Antiprotons
Antiprotons (p) have long been recognized as a promising indirect signature of DM.
While DM annihilation typically produces equal numbers of protons and antiprotons, the astrophysical background flux of antiprotons is very small in comparison to that of protons.
On the other hand, describing the production and propagation of these charged hadrons is a challenging task, and any statement regarding DM annihilation rests on our ability to understand the associated astrophysical uncertainties. The Alpha Magnetic Spectrometer (AMS-02) experiment has provided the most precise measurements of the cosmic ray proton and antiproton flux to date [84], and here we will explore the implications of this data on our neutrino portal DM scenario. Since DM annihilates to RHNs, which subsequently decay via W , Z, and Higgs bosons, the resulting cascade decay, showering and hadronization produce a variety of hadronic final states including antiprotons. AMS-02 will therefore provide an important probe of the model.
The propagation of antiprotons through the galaxy to earth is described by a diffusion equation for the distribution of antiprotons in energy and space (see, e.g., Ref. [85] and references therein). The transport is modeled in a diffusive region taken to be a cylindrical disk around the galactic plane and is affected by several physical processes. These include diffusion of the antiprotons through the turbulent magnetic fields, convective winds that impel antiprotons outward, energy loss processes, solar modulation, and a source term describing the production and loss of antiprotons. The source term accounts for astrophysical sources such as secondary and tertiary antiprotons, and antiproton annihilation with the interstellar gas, as well as primary antiprotons produced through DM annihilation. The propagation depends on a number of input parameters, and a set of canonical models, called MIN, MED, MAX are often employed [86]. The diffusion equation is solved assuming the steady state condition to find the flux of antiprotons from DM annihilation at earth, where dNp/dK is the kinetic energy (K) spectrum of antiprotons per DM annihilation, vp is the antiproton velocity, and ρ 0 is the local DM density. The propagation function R(K) accounts for the astrophysics of production and propagation, and we use the parameterization provided in Ref. [87].
AMS-02 has provided precise measurements of the proton flux, Φ p (K) [88], and the antiproton-to-proton flux ratio, r(K) [84], which can be used to place constraints on DM annihilation. To proceed, we require an estimate of the secondary background antiproton flux originating from astrophysical sources. For this purpose we use the best-fit secondary flux, Φp ,bkg (K), from [89], which provides an acceptable fit to the AMS-02 data. With the total antiproton flux, Φp ,tot (K, m χ , σv ) = Φp ,bkg (K) + Φp ,χ (K, m χ , σv ), and the measured proton flux from AMS-02, Φ p (K), in hand, we form the ratio of these two fluxes and fit it to the observed ratio. The test statistic is where i runs over energy bins, and σ i is the reported uncertainty of the flux ratio [84]. Following Ref. [89], we define a limit on σv as a function of m χ , m N according to the condition where χ 2 0 is the best fit chi-squared statistic assuming no primary DM antiproton source from Ref. [89]. The limit is derived under the assumption of a Einasto profile and using the MED propagation scheme. Contours of the limit on the annihilation cross section in the m χ − m N plane are displayed in Figure 5. For DM masses in the range of 20 -80 GeV, AMS-02 is able to probe the thermal cross section Eq. (6), as indicated by the thick (red) line. The vertical (blue) lines show the limits on the minimum value of the coupling constant λ. The shaded (blue) region indicates the perturbative unitarity bound. It is important to note again that there are significant uncertainties associated with the DM halo profile and the propagation scheme, which can lead to a variation in the cross section limits by one order of magnitude or more [89]. Note that for a fixed m χ , the limits in Figure 5 become stronger as m N is increased. This is because for fixed m χ , heavier RHNs tend to produce more low energy antiprotons (see Figure 1). However, the ratio r(K) shows good agreement with the astrophysical background model at low value of kinetic energy K and a slight excess at larger values of K, explaining the behavior seen in Figure 5. There are several other notable indirect DM searches that we wish to comment on here.
AMS-02 has provided detailed measurements of the cosmic ray positron spectrum [90]. Much attention has been paid to these results (and those of its forerunner PAMELA [91]) due to the observation of a striking rise in the fractional positron flux, which potentially points to a new primary source of positrons. While it is true that DM annihilation in our scenario produces a significant positron flux, the cross section limits from Fermi dSphs gamma rays and AMS-02 antiproton observations are expected to be stronger than those from AMS-02 positron measurements by an order of magnitude or more, and thus we have chosen to focus on these stronger tests.
Another well-known indirect DM probe is high energy neutrinos from DM annihilation in the sun, which can be probed with the IceCube experiment [92]. But under the minimal assumption of typical seesaw values for the neutrino Yukawa coupling (see Eq. (3)) the DMnucleon scattering rate will be too small to allow for the efficient capture of DM in the sun, so we do not consider this possibility further.
As we have demonstrated, the data collected so far by Fermi-LAT already leads to stringent limits on DM parameter space, and the sensitivity will improve significantly in the coming years. The projected sensitivities for 10 and 15 years of data taking has been studied in detail by the collaboration in Ref. [93]. The fast discovery of new dSphs is the primary upcoming change in dSph targeted DM searches. The identification of new dSph candidates by the Dark Energy Survey (DES) [94] over the past two years, if confirmed, will double the number of known dSphs. Following on important discoveries of the Sloan Digital Sky Survey (SDSS) [95], which covered 1/3 of the sky and discovered 15 ultra-faint dSphs, surveys like DES and especially the Large Synoptic Survey Telescope (LSST) [96] will cover complementary regions of the sky which are expected to discover potentially O(100) dSphs. Ref. [93] takes 60 total dSphs as an estimate of the number of dSphs that can be used for LAT searches. They find that the sensitivity of searches targeting dwarf galaxies will improve faster than the square root of observing time. Following Ref. [93] we expect an improvement on the cross section limit from Fermi-LAT 15 years dSph observations by a factor of a few, which will probe thermal relic DM with masses m χ 100 GeV in the neutrino portal DM scenario. Due to their large effective areas, ground-based imaging air Cherenkov telescopes (IACTs), such as H.E.S.S. [97], VERITAS [98], and MAGIC [99], and in the future CTA [100] and HAWC [101], are well suited to search for higher energy gamma rays originating from heavy DM annihilation. In particular, H.E.S.S. has presented a search for DM annihilation towards the Galactic Center using 10 years of data [97]. Assuming a cuspy NFW or Einasto profile the search sets the strongest limits on TeV mass DM that annihilates to W W or quarks, and almost reach thermal annihilation rates. Taken at face value, the H.E.S.S. limits are indeed stronger than the Fermi dSphs limits for DM masses above a few hundred GeV, but are however less robust due to the inherent astrophysical uncertainties associated with the central region of the Milky way, both in terms of conventional gamma-ray sources and the DM distribution. The H.E.S.S data is not publicly available, so unfortunately we are not able to properly recast their limit. However, for a fixed DM mass, the continuum photon spectrum produced in our model from χχ → N N is qualitatively similar to the spectrum produced by χχ → W W . We can therefore obtain a rough estimate of the H.E.S.S. sensitivity by translating their limits in the W W channel to our parameter space The H.E.S.S.
limits are approaching the canonical thermal relic annihilation rate for DM masses around 1 TeV.
In the future, the Cherenkov Telescope Array (CTA) will be able to further probe heavy TeV-scale DM annihilation, with the potential to improve by roughly an order of magnitude in cross section sensitivity over current instruments depending on the annihilation mode and DM mass. Here we estimate the sensitivity of future CTA gamma-ray observations of the Galactic Center using a "Ring" method technique [102]. Our projections are based on a simplified version of the analysis carried out in Ref. [103] that we now briefly describe. The analysis begins with the definition of signal (referred to as "ON") and background ("OFF") regions. A binned Poisson likelihood function is constructed in order to compare the DM model µ µ µ to a (mock) data set n: where µ ij is the predicted number of events for a given model µ µ µ in the ith energy bin and the jth region of interest, corresponding to ON (j = 1) and OFF (j = 2) regions.
These model predictions are compared to the corresponding observed counts n ij . We use 15 logarithmically-spaced energy bins, extending from 25 GeV to 10 TeV. The number of gamma-ray events predicted by each model consists of three components: a DM annihilation signal, an isotropic cosmic-ray (CR) background, and the Galactic diffuse emission (GDE) background: The details for the regions of interest that have been used in our analysis, including the corresponding solid angles and J-factors, can be found in Ref. [100]. Furthermore, we have used the effective area produced by MPIK group [104] and fixed the time of observation to be 100 hours.
We account for differential acceptance uncertainties (i.e. acceptance variations across different energy bins and regions-of-interest) by rescaling the predicted signals µ ij by parameters α ij and profiling the likelihood over their values. Following Ref. [103] we assume Gaussian nuisance likelihoods for all α with respective variance σ 2 α independent of i and j. Our limits correspond to differential acceptance uncertainties of 1%. The mock data n we employ includes a fixed isotropic cosmic-ray background component in all bins, and no signal from DM annihilation. We derive 95% CL upper limits (sensitivity) on the annihilation cross-section σv in the usual way by requiring −∆ ln L ≤ 2.71/2. Our projections are shown in Figure 7. We have not included systematic uncertainties for the background components, which can be as large as order one and thus significantly degrade the CTA sensitivity. However, this can be partially overcome through a more sophisticated morphological analysis, which leverages the shape differences between the galactic diffuse emission and DM signal [103]. In the end, we expect that Figure 7 provides a reasonable ballpark estimate of the CTA sensitivity, which can improve over H.E.S.S. by a factor of a few to ten in the 100 GeV -TeV DM mass range. We expect Fermi dSphs observations to provide superior limits for lower mass DM, m χ 100 GeV.

IV. GALACTIC CENTER GAMMA RAY EXCESS INTERPRETATION
As mentioned in Section III B, various analyses of Fermi-LAT data show a spherically symmetric excess of gamma rays coming from the central region of the Milky Way peaking in the 1-3 GeV energy range [18][19][20][21][22]. Since DM annihilation to RHNs abundantly produces gamma rays, it is interesting to explore a possible interpretation of this excess in the context of the neutrino portal DM model. In fact, this possibility was previously investigated in Ref. [17], which found that DM annihilation to RHNs could indeed provide a good fit to the Galactic Center excess. Here we will additionally confront this interpretation with existing constraints from other indirect probes, and notably Fermi gamma-ray observations from dSphs and AMS-02 antiproton observations. We fit the neutrino portal DM model parameters to the Galactic Center excess spectrum given in Ref. [22]. We adopt Navarro-Frenk-White (NFW) profile with γ = 1.2. Following [22] we define the χ 2 as where θ={ σv , m χ , m N }, Φ i ((Φ i ) obs ) is the predicted (observed) γ-ray flux (see Eq. (13)) in the i th energy bin, and Σ is the covariance matrix. We find that the best-fit point is degrees-of-freedom. Figure 8 displays 1σ, 2σ, and 3σ CL regions in the m N − m χ parameter space. We see that neutrino portal DM can provide an acceptable fit over a significant range of mass parameters.
Next, we would like to confront this interpretation with the other constraints derived in Section III. To this end, we perform the Galactic Center excess while fit fixing the annihilation cross section to its thermal value, and overlay the limits derived from Planck CMB, Fermi dSphs, and AMS-02 antiproton observations. The result is displayed in the right panel of Figure 8. We see that this interpretation faces some tension with limits from dSphs and antiprotons. However, it is too early to conclude from this analysis that the DM interpretation of the excess is not viable given the significant astrophysical uncertainties in the local DM density, dSphs DM densities, and the modeling of the antiproton propagation.

V. BEYOND THE MINIMAL SCENARIO
We have explored what is perhaps the simplest scenario of neutrino portal DM. The primary probe of this model comes from indirect detection, and we have presented a comprehensive picture of the current constraints. However, it is possible that the neutrino mass model is more complex than the simplest Type-I seesaw, or that there are additional interactions of the scalar mediator with the Higgs, in which case a much richer phenomenology is possible. In this section we will highlight some of these possibilities.

A. Large neutrino Yukawa coupling
Taking the naive seesaw relation in Eq.
(2) as a guide, one generally expects very small active-sterile mixing angles, θ ∼ m ν /m N 10 −6 × (m N /100 GeV) − Such large Yukawa couplings not only offer increased chances to probe the RHNs directly (see, e.g., Ref. [106,107] for a revew), but will also enhance the detection prospects of the DM sector. For instance, one can induce sizable DM couplings to the Z and Higgs boson at one loop that mediate large scattering rates with nuclei, which is relevant for direct detection experiments and capture of DM in the sun. One can also potentially produce the RHNs directly in accelerator experiments.
This also opens up the possibility for the RHN to be heavier than the dark sector particles, while still having a thermal cosmology. Due to the large mixing angle, it is possible for DM to annihilate efficiently into light active neutrinos, and furthermore the DM may annihilate to other SM particles through the loop-induced Z and h couplings. We refer the reader to Refs. [27,28,32] for recent investigations of these issues.

B. Higgs portal coupling
The scalar particle φ can couple to the Higgs portal at the renormalizable level We have so far assumed that this coupling is small. The reason we have made this assumption is primarily for simplicity, as then the phenomenology and cosmology is solely dictated by the neutrino portal link to the SM. However, this assumption can certainly be questioned.
Restricting to the fields and interactions of our scenario in Eq. (1), we observe that the Higgs portal coupling (26) will be induced at one loop with strength of order λ φH ∼ λ 2 y 2 /16π 2 , which is very small due to the small neutrino Yukawa coupling. Still, one may expect unknown UV physics to generically induce a larger coupling. This is because there is no enhanced symmetry in the limit λ φH → 0, and so even though the operator (26) is marginal, we cannot rely on technical naturalness ensure a small value without further information about the UV physics. That being said, one can certainly imagine completions in which the Higgs portal coupling is suppressed. For example if φ is a composite scalar state of some new strong dynamics, then the Higgs portal operator would fundamentally be a higher dimension operator and could be therefore be naturally suppressed.
Another good reason to consider the Higgs portal operator is that it provides additional opportunities to probe the dark sector in experiment. A one loop coupling of the DM to the Higgs will be induced and this can mediate scattering of DM with nuclei, or invisible decays of the Higgs to DM [27][28][29][30].
pair of light scalars, h → φφ. These scalars, once produced would then cascade decay via φ → N χ. The resulting RHN N , being lighter than the W boson, will have a macroscopic decay length and could leave a striking displaced vertex signal (see, e.g., [108]). The signature would thus be an exotic Higgs decay with two displaced vertices.

VI. SUMMARY AND OUTLOOK
In this paper, we have investigated a simple model of neutrino portal DM, in which the RHNs simultaneously generate light neutrino masses via the Type-I seesaw mechanism and mediate interactions of DM with the SM. The model, presented in Section II, is quite minimal and contains a dark sector composed of a fermion χ (the DM candidate) and scalar φ, along with the RHN N . Given the generic expectation of tiny neutrino Yukawa couplings, testing this model with direct detection or accelerator experiments is likely to be challenging.
However, it is possible in this model that DM efficiently annihilates to RHNs, which allows for a number of indirect probes of this scenario.
We have carried out an extensive characterization of the indirect detection phenomenology of the neutrino portal DM scenario in Section III. Restricting to an experimentally and theoretically viable mass range, 1 GeV m N < m χ 10 TeV, we have derived the constraints on the χχ → N N annihilation cross section from Planck CMB measurements, Fermi gamma-ray observations from the Galactic Center and from dSphs, and AMS-02 antiproton observations. Currently, the dSphs and antiproton measurements constrain DM masses below 50 GeV for thermal annihilation rates. In the future, Fermi dSphs observations will be able probe DM masses above the 100 GeV range for thermal cross sections, while CTA will be able to approach thermal cross section values for DM masses in the 100 GeV -1TeV range.
This model can also provide a DM interpretation of the Fermi Galactic Center gamma ray excess as discussed in Section IV. We have verified that the predicted spectrum of gamma rays is compatible with the observed excess for RHN and DM masses in the 20−60 GeV range and annihilation rates close the the thermal value. However, we have also shown that this interpretation faces some tension with the existing constraints from Fermi dSphs and AMS-02 antiprotons (subject of course to various astrophysical uncertainties). It will be interesting to see how this situation develops as Fermi and AMS-02 collect more data. However, at least in the simplest model explored here, it will be challenging to find complementary probes outside of indirect detection.
It is possible that the neutrino mass generation mechanism is more intricate than the simplest Type-I seesaw, as discussed in Section V. If so, the implications for neutrino portal DM could be dramatic, particularly if the neutrino Yukawa coupling is large, as this could lead to direct detection prospects, accelerator probes, and new annihilation channels. Additionally, it is possible in this scenario for additional Higgs portal couplings to be active, which could yield further phenomenological handles.
Portals provide a simple and predictive theoretical framework to characterize the allowed renormalizable interactions between the SM and DM. Furthermore, the existence of neutrino masses already provides a strong hint that the neutrino portal itself operates in nature. These two observations provide a solid motivation for testing the neutrino portal DM scenario, both through the generic indirect detection signals investigated in this paper, and also the additional signals present in more general models. It is worthwhile to broadly explore these scenarios and their associated phenomenology in detail, and we look forward to further progress in this direction in the future.