Sterile neutrino portals to Majorana dark matter: effective operators and UV completions

Stringent constraints on the interactions of dark matter with the Standard Model suggest that dark matter does not take part in gauge interactions. In this regard, the possibility of communicating between the visible and dark sectors via gauge singlets seems rather natural. We consider a framework where the dark matter talks to the Standard Model through its coupling to sterile neutrinos, which generate active neutrino masses. We focus on the case of Majorana dark matter, with its relic abundance set by thermal freeze-out through annihilations into sterile neutrinos. We use an effective field theory approach to study the possible sterile neutrino portals to dark matter. We find that both lepton-number-conserving and lepton-number-violating operators are possible, yielding an interesting connection with the Dirac/Majorana character of active neutrinos. In a second step, we open the different operators and outline the possible renormalisable models. We analyse the phenomenology of the most promising ones, including a particular case in which the Majorana mass of the sterile neutrinos is generated radiatively.


Introduction
Non-zero neutrino masses and the existence of dark matter (DM) in the Universe constitute the two most compelling pieces of evidence for physics beyond the Standard Model (SM). Many extensions of the SM able to account for non-zero neutrino masses involve SM gauge singlet fermions N R , referred to as sterile or right-handed (RH) neutrinos. Apart from the Yukawa coupling y ν to the SM singlet operator LH, with L and H being the SU (2) L lepton and Higgs doublets, respectively, N R could also interact with a dark sector containing, in particular, a DM candidate. In this case, N R would serve as a portal between the visible and dark sectors. This is the idea of the sterile neutrino portal to DM [1]. Due to the stringent limits on Weakly Interacting Massive Particles (WIMPs) from direct detection (DD), indirect detection (ID) and collider searches, it is natural to consider the case in which the dominant coupling of the dark sector to the SM is via sterile neutrinos.
In what follows, we will assume that DM is produced through the thermal freeze-out mechanism. Then, depending on the mass of sterile neutrino, m N , one can identify two distinct regimes. If m N is smaller than the mass of DM, m DM , the DM relic abundance is set by its annihilations into sterile neutrinos. The latter subsequently decay into SM particles leading to ID signatures in photon, charged particle and neutrino spectra. This is the so-called secluded regime [1]. The phenomenology of this regime has been studied in Refs. [2][3][4][5] considering a simple renormalisable model for the dark sector containing a fermion and a real/complex scalar, both charged under a Z 2 /U (1) "dark" symmetry. In the opposite regime, when m N > m DM , DM annihilation into a pair of sterile neutrinos is kinematically forbidden, unlike the annihilation into SM particles, in particular, active neutrinos, which proceeds via active-heavy neutrino mixing. In this case, this mixing has to be sizeable to provide the annihilation cross section required to explain the observed relic abundance. The same annihilation process could lead to ID signatures at neutrino experiments, see e.g. [6][7][8][9][10][11]. The phenomenology of this regime has been investigated in Refs. [12][13][14] assuming that the dark sector comprises a fermion and either a scalar or a vector boson. Furthermore, the freeze-in mechanism of DM production in neutrino portal scenarios has been studied in Refs. [15][16][17][18]. Further examples of studies investigating DM-neutrino connections can be found in Refs. [19][20][21][22].
In the present work, we will concentrate on the secluded regime, assuming that (i) the DM candidate is a Majorana fermion χ charged under a Z 2 symmetry responsible for its stability, and (ii) DM interactions with N R are given by effective four-fermion interactions generated at the new physics scale Λ. The effective field theory (EFT) approach for interactions of DM with the SM extended with RH neutrinos has been developed in Ref. [23], 1 and the four-fermion operators we will focus on in this work form part of the basis of dimension-six operators derived in Ref. [23]. As we will see, restricting ourselves to the interaction of DM with the lightest of RH neutrinos, there are three different four-fermion operators. After studying them in detail, we will discuss simple UV completions generating these operators at tree level. Some of the UV models will lead to other operators as well, resulting in either DM and sterile neutrino self-interactions, and/or direct interactions of DM with the SM. We will identify the most interesting/promising models and study their DM phenomenology, making also a connection to the mechanism of neutrino mass generation.
It is worth noting that what is often meant in the literature as neutrino portal is the operator (LH) × O dark , where O dark is a singlet of a dark symmetry group [28]. In particular, O dark can be composed of a dark fermion and a dark scalar [29,30]. The operator (LH) × O dark may result from integrating out a heavy RH neutrino, and thus, the latter is not present in the corresponding EFT. This is different from what we will discuss in the present work. See Fig. 1 for a schematic diagram of our set-up.
The paper is structured as follows. In Sec. 2, we describe the EFT approach to the interactions of DM with RH neutrinos. In Sec. 3, we classify the UV models generating the four-fermion operators of interest at tree level. In Sec. 4, we discuss the phenomenology of several selected Figure 1: Schematic diagram of the considered set-up, in which the sterile neutrino N serves as a portal between the SM and the dark sector containing the DM χ. NP stands for new physics generating effective four-fermion interactions between N and χ at the scale Λ. models. In Sec. 5, we summarise our findings and draw our conclusions. Finally, some technical details are given in the appendix.

Effective field theory approach 2.1 Four-fermion sterile neutrino portal operators
We add to the SM particle content two chiral fermions, N R and χ L , transforming as (1, 1) 0 with respect to the SM gauge group (SU (3) C , SU (2) L ) Y . The first one, N R , is the usual RH neutrino, while we consider the case in which the latter, χ L , has a discrete symmetry Z 2 that stabilises it and makes it a potential DM candidate. The most general renormalisable Lagrangian reads where L SM is the SM Lagrangian (with massless neutrinos). We can always re-phase N R and χ L in such a way that the masses m N and m χ are real and positive. If the lepton number U (1) L symmetry is imposed with L(N R ) = 1 and L(χ L ) = 0, then the Majorana mass term for N is forbidden. After electroweak symmetry breaking (EWSB), the Higgs field takes a vacuum expectation value (VEV), H = (0, v h / √ 2) T , with v h = 246 GeV. In order to generate two non-zero light neutrino masses (either of Dirac or Majorana type), at least two RH neutrinos should be included. For simplicity, in the following we assume that only one of the sterile neutrinos is lighter than the DM particle, so that m N < m χ , and there is an open annihilation channel χχ → N N , while the heavier ones decouple from the spectrum. Here χ and N stand for Majorana fields, i.e.
We assume that the lightest of the RH neutrinos contributes to active neutrino masses and is also responsible for the DM phenomenology we are interested in. The results can be easily extended to the case of all of them being lighter than DM. We focus on four-fermion effective operators describing interactions between χ and N , which have masses below the new physics scale Λ, so that m N < m χ < Λ. At dimension D = 6, there are three such operators that connect DM with the SM through the sterile neutrinos. We will refer to them as sterile neutrino portal operators or simply portal operators. The corresponding D = 6 Notation Operator Dimension effective Lagrangian reads

Majoron interactions
with where in the second equalities we have used Fierz identities. 2 The Wilson coefficient c 1 is real, whereas c 2 and c 3 are, in general, complex. In general, many other operators at mass dimensions D = 5 and D = 6 involving χ L exist [23]. In Sec. 3, we will consider renormalisable models that generate the portal operators, generically together with other ones. In Tab. 1, we summarise all the D ≤ 6 operators generated in these models.
The operator O 1 preserves lepton number (we will refer to this operator as lepton-numberconserving, LNC). In particular, it is allowed if light neutrinos are Dirac particles. In this case m N = 0, and light neutrino masses are given by m ν = y ν v h / √ 2. The operators O 2 and O 3 break lepton number in two units (we will refer to these operators as lepton-number-violating, LNV). 3 2 Note that the r.h.s of Eq. (5) would in general involve +1/2(NRσµν N c R )(χ c L σ µν χL). However, such a term vanishes for one generation of N or χ. 3 As we will see in Sec. 3, in some models χL is also charged under U (1)L. Depending on its charge, O2 (O3) If lepton number is broken, then in general, also m N = 0, and in the basis (ν c L , N R ) the neutrino mass matrix reads where m D = y ν v h / √ 2 is the Dirac neutrino mass matrix. Assuming for simplicity one generation of active and sterile neutrinos, the eigenvalues of M tree so the active neutrino mass is In the following, unless stated otherwise, we assume that there is a contribution to active neutrino masses via the standard seesaw mechanism [31][32][33][34][35][36], 4 so that for m N m D , we obtain Notice that the mass eigenstates are approximately equal to the weak eigenstates because the neutrino mixing with the heavy states ∼ m light /m heavy is always very small, namely, it is smaller than 10 −5 for m N 2 GeV. This lower limit on m N stems from the fact that sterile neutrinos should decay before Big Bang nucleosynthesis (BBN) in order not to spoil light-nuclei abundances. More specifically, following Ref. [3], we require the N lifetime τ N 1 s. In the case of Dirac neutrinos (when N R is the RH component of a Dirac neutrino field), the Yukawa coupling y ν ∼ 10 −13 is such that N R are never in thermal equilibrium with the SM bath and their energy density at BBN is negligible, leading to no contribution to the number of relativistic degrees of freedom (N eff ≈ 3), see e.g. Ref. [38].

Thermal equilibrium
The only new state that directly couples to the SM is the sterile neutrino, through the Yukawa coupling y ν . In this scenario, this coupling may be suppressed, as unitarity, EFT validity and DM annihilations require and therefore, imposing m ν 0.05 eV and using the seesaw relation in Eq. (10), we get will preserve (break) lepton number or vice versa. More specifically, if L(χL) = 1, then O2 is LNC and O3 is LNV, whereas if L(χL) = −1, the roles are interchanged. In both cases, the LNV operator will break lepton number in four units. Note that O1 always preserves lepton number. 4 In the context of the SM gauge group, the seesaw mechanism has been discussed in Ref. [37].
One may wonder whether such a sterile neutrino would be in thermal equilibrium in the early Universe and have standard freeze-out. Indeed, such small Yukawa coupling does not suffice to keep N R , and thus, DM in equilibrium with the SM all the way down to the freeze-out of DM. 5 However, as we will see in Sec. 3, the openings of the effective operators include scalars and gauge bosons that always have interactions with the SM, via Higgs portals and/or kinetic mixing. Moreover, other EFT operators, even if irrelevant with respect to DM annihilations compared to those in Eqs. (4)-(6), may be strong enough to keep DM in thermal equilibrium with the SM. Therefore, it is safe to assume that, early on, N R and χ L were in thermal equilibrium with the SM.
Even if early on all particles are in equilibrium, the dark sector may kinetically decouple from the SM before the freeze-out of DM. Thermal evolution of a such decoupled dark sector has been studied in Refs. [43][44][45]. The corrections to the cross section needed to obtain the observed relic abundance in this case depend, in particular, on whether dark sector particles are relativistic or not at the time of DM chemical freeze-out. Particularising to our dark sector, if sterile neutrinos are relativistic at the time of kinetic decoupling and down to the freeze-out (i.e. for m N m χ /20), the temperature of the dark sector, T D , is similar to that of the SM bath, T , leading to a very mild effect on the relic abundance [43]. This is the case realised for the most of the parameter space investigated in our analysis in Sec. 4. On the other hand, if sterile neutrinos and/or dark matter are relativistic at kinetic decoupling but become non-relativistic at freeze-out (i.e. for m N m χ /20), the dark sector will be reheated. In this case, T D /T could reach a factor of a few [43], and the impact on the final relic abundance would be larger. In App. A, we discuss some implications of the departure from the standard assumptions of chemical and kinetic equilibrium. However, a full treatment of these effects on the relic abundance is beyond the scope of the present study.

Dark matter annihilations
As there is thermal equilibrium, the DM relic abundance is set by the standard freeze-out of annihilations χχ → N N . In the non-relativistic limit, which is appropriate for freeze-out temperatures, the cross section only depends on the relative velocity v = | v 1 − v 2 | of the DM particles and can be expanded as: The coefficients a and b are associated to the s-and p-wave contributions to the annihilations. Performing the calculation, we find: with r N ≡ m N /m χ . In the limit of negligible m N , these formulae reduce to These results can be understood from the arguments based on the discrete symmetries of a pair of Majorana fermions and conservation of the total angular momentum, J. They agree with the conclusions of the general analysis performed in Ref. [46]. Below we apply them to our case.
The wave function of the Majorana DM particles in the initial state should be anti-symmetric. Since this is defined by (−1) L i (−1) S i +1 = −1, with S i and L i being the spin and the orbital angular momentum of the initial pair, we conclude that L i + S i must be even. Choosing the z-axis to lie along the direction of motion of the outgoing particles, we have L f z = 0 and J f z = S f z , where the index f refers to the pair in the final state. Now, we can study the final states generated by the portal operators in Eqs. (4)-(6).
• If m N = 0, N R can be described by a Weyl fermion. Then, the operator O 1 produces a pair N R , N R , with opposite helicities, +1/2 for N R and −1/2 for N R . Therefore, the spins are aligned and |S f z | = 1. Conservation of the total angular momentum implies that |J i z | = |J f z | = 1. Since L i + S i must be even, the lowest order combination that can realise |J i z | = 1 is S i = L i = 1. Therefore, we conclude that O 1 generates a p-wave suppressed DM annihilation cross section, cf. Eqs. (17) and (18). If m N = 0, there can be a helicity flip and S f z = 0 can be attained leading to the s-wave proportional to m 2 N . This agrees with Eq. (15). Moreover, when the N R , N R pair is coupled in an s-channel, the mediator should be a vector boson. This can also be easily seen by observing the Fierz transformed version of O 1 in Eq. (4).
• If m N = 0, O 2 (O 3 ) creates a pair N R , N R (N R , N R ) with both states having the same helicity, +1/2 (−1/2), and therefore, the spins are anti-aligned and S f z = 0. Conservation of J implies that J i z = J f z = 0, and at the lowest order we have S i = L i = 0, leading to the s-wave DM annihilation cross section. This is in agreement with Eq. (15). Moreover, when the N R , N R (N R , N R ) pair is coupled in an s-channel, the mediator should be a scalar boson. This can also be easily seen by observing the Fierz transformed version of O 2 in Eq.   6 Here χ and N are the Majorana fields defined in Eq. (2).
For a state formed by a pair of Majorana particles, parity is given by P = (−1) L+1 , where L is the orbital angular momentum of the state. In particular, the bilinear χχ annihilates a pair of DM particles with P = +1 which at the lowest order corresponds to L = 1 (see e.g. Ref. [46]), 7 hence the DM annihilation cross section χχ → N N is p-wave. Conversely, the bilinear iχγ 5 χ annihilates a pair of DM particles with L = 0, so the DM annihilation cross section is s-wave. Finally, the zeroth component of the bilinear χγ µ γ 5 χ has P = −1, whereas the spatial components have P = +1, and thus it contributes to both s-and p-waves. In view of that, having p-wave DM annihilation cross section requires that the terms in the Lagrangian in Eq. (19) involving the bilinear iχγ 5 χ and χγ 0 γ 5 χ vanish. This implies that c 1 = 0 and c 2 = −2c * 3 . For completeness, we also consider the scenario in which N R is the RH counterpart of the LH SM neutrino ν L , i.e. ν is Dirac with mass m ν ∼ 0.05 eV. Lepton number is conserved in this case, and we have only one operator, O 1 . The coefficients a and b in the annihilation cross section, Eq. (14), are now given by with to a very good approximation. As expected, these expressions agree with those obtained in the Majorana case in the limit m N → 0 with c 2 = 0 and c 3 = 0, cf. Eqs. (17) and (18). When the DM annihilation cross section in Eq. (14) is thermally averaged one obtains an expansion in inverse powers of x = m χ /T : Notice that relativistic corrections will only affect higher order terms in this expansion [47]. The typical value for x at freeze-out is around 20-25. The observed relic abundance corresponds to σv ≈ 2.2 × 10 −26 cm 3 /s if a = 0 [48]. If the cross section is p-wave (a = 0 and b = 0), then σv ≈ 4.4 × 10 −26 cm 3 /s is required at freeze-out [49] (see also [9,50]). More precisely, the quoted values of σv apply for m χ 10 GeV, whereas for smaller DM masses larger values of σv are needed [45,48]. Hence, in our numerical analysis we employ the results of Ref. [45], where σv that reproduces the observed relic abundance is given in Figs. 1 and 4 as a function of m χ for the cases of s-and p-wave DM annihilation cross section, respectively.
In the left panel of Fig. 2, we show values of m χ and the new physics scale Λ 1 (associated with the LNC operator O 1 ) leading to the correct thermal relic cross section assuming standard freeze-out taking place at x 20. For m N = 0, the cross section is p-wave, whereas for m N = 0, it is s-wave. As can be seen, the impact of m N on the scale of new physics needed to get the correct relic abundance is rather moderate. For m N = 0 (blue line), the lower limit on m χ 100 MeV is set by requiring that χ is in thermal equilibrium until the freeze-out, which we take at x 20.   Near the threshold m χ m N one can not use the expansion of the DM annihilation cross section in terms of the coefficients a and b, as discussed in Ref. [47]. In this case, we use the general relativistic expression for σv given in Eq. (3.8) of the aforementioned reference (see also [51]). We provided this formula in Eq. (70) in App. A. When we take into account the thermal average for DM masses slightly smaller than m N , the σv is Boltzmann suppressed. Therefore, we need to decrease Λ 1 in order to reproduce the observed relic abundance. In the right panel, we work in the limit m N → 0 and turn on one operator at a time. The LNV operators lead to the annihilation cross sections not suppressed by 1/x, cf. Eq. (17). Thus, for a given DM mass and O 2 (O 3 ), a new physics scale a factor of 3 (4) larger than that for O 1 is needed to reproduce the observed relic abundance. We also display the region where Λ i < m χ (assuming c i = 1), i.e. where the EFT description is not valid.

Tree-level UV completions of the portal operators
In this section, we consider tree-level UV completions of the four-fermion neutrino portal operators. They can be divided into models involving a real/complex heavy scalar mediator in either t-channel (type A) or s-channel (type B), or else a heavy vector mediator (type C), as depicted in Fig. 3. In what follows, we assume that the mass of the mediator is always larger than the mass of DM. This forbids annihilation of DM into the mediators and ensures the neutrino portal regime. In Tab. 2, we summarise the dark sector particles of each model along with their Z 2 and B − L charges. We note that in type-A models, the scalar mediator is charged under the Z 2 symmetry stabilising DM, whereas in Models B and C, the mediators are neutral under this symmetry. In the following, the models that at D ≤ 6 generate only the portal operators O 1 and/or O 2 and/or O 3 will be referred to as genuine. As we will see below, these are Models A. In addition, in Tab. 3 we provide the tree-level matching conditions for the Wilson coefficients of the effective operators generated by Figure 3: Diagrams for DM annihilation into sterile neutrinos in the three types of renormalisable models considered. φ and σ stand for a real and complex singlet scalar, respectively, whereas Z is a vector mediator. In Models A, the mediator, exchanged in the t-channel, is charged under the Z 2 . We termed these models genuine, as they generate only the four-fermion effective operators with sterile neutrinos and DM.

Models C
the UV models. See details in the next subsections.

Model A1: real scalar
The model includes a real scalar φ, charged under the same Z 2 symmetry as χ. The Lagrangian reads where L 4 is given in Eq. (1) and V (φ, H) is the most general scalar potential invariant under Z 2 : For the sake of simplicity, in what follows, we consider one generation of N R and χ L , such that the new Yukawa couplings are numbers. Since we have already re-phased N R and χ L to render the masses m N and m χ real and positive, and φ is a real field, the phase of the Yukawa coupling f cannot be absorbed. Hence, in general, it is complex. The new Yukawa interaction breaks U (1) L . Thus, when integrating out φ, both LNC and LNV operators are expected to be present. Indeed, we find that O 1 and O 2 are generated with the tree-level matching conditions given in Tab. 3. No other operators at D ≤ 6 are induced.
The DM phenomenology of this model has been studied in Ref. [2]; in addition, Ref.
[52] analyses the cases when N and χ have similar masses and Ref. [5] when also φ does. A wide range of masses for N and χ was studied in Ref. [3]. On the other hand, the freeze-in mechanism of DM production in this model was considered in Refs. [15,17,18].
Model • A2c. Majorana neutrinos, but with m N = 0 in Eq. (1) and lepton number being softly broken only in the scalar potential: We can absorb the phase of µ 2 σ in σ, hence making µ 2 σ real and positive. Further, in the absence of m N , we can redefine N R (as well as L and e R ) in such a way that f also becomes real. This model has some interesting features. For example, finite m N is generated at one loop. We will discuss this feature in Sec. 4.2. Integrating out the complex scalar σ, we find that both O 1 and O 2 are generated with the matching conditions given in Tab. 3. As expected, the Wilson coefficient of O 2 is proportional to the LNV parameter µ 2 σ . Alternatively, we can integrate out the real, ρ, and imaginary, θ, parts of σ. This leads to the following matching relations: where For µ 2 σ /m 2 σ 1, these matching conditions reduce to those given in Tab. 3. No other non-renormalisable operators of D ≤ 6 are induced.

Model B1: real scalar
It includes a real scalar φ, this time not charged under the Z 2 symmetry. In this case, where V (φ, H) is the most general scalar potential: 8 In general, the Yukawa couplings f and g are complex. Integrating out φ leads to the LNV portal operators O 2 and O 3 . As can be inferred from the matching conditions given in Tab As can be seen from Tab. 3, the Wilson coefficients of these operators are controlled by the scalar coupling µ φH . Thus, if µ φH is small, these operators are suppressed with respect to the neutrino portal operators. Upon EWSB, the first operator contributes to the (mostly-)sterile neutrino mass and Higgs or N decays, depending on the value of m N . The second operator contributes to the DM Majorana mass and provides the fermionic Higgs portal with associated DM phenomenology, see e.g. [56][57][58]. In addition, at D = 6, we find the four-fermion self-interactions: with the matching conditions provided in Tab. 3. Since we consider one generation of N R and χ L , the operators (N c

Model B2: global U (1) B−L
Instead of a real scalar φ, this model includes a complex scalar σ. A complex scalar calls for an associated U (1) symmetry. In this case, we will consider lepton number, or rather U (1) B−L , since the latter is an anomaly-free global symmetry of the SM. The corresponding lepton charges are L(N R ) = L(χ c L ) = 1 9 and L(σ) = −2. The Lagrangian of this model is given by where V (σ, H) is the LNC potential given in Eq. (27). The couplings f and g can be rendered real. This model has been studied in detail in Ref. [59]. If the complex scalar acquires a VEV, 8 We note that the term linear in φ can always be removed by a shift. 9 Other charge assignments are possible since the new fermions are singlets under the SM gauge group.
v σ , the U (1) B−L symmetry gets broken spontaneously, and Majorana masses m N = √ 2f v σ and m χ = √ 2gv σ are generated. We can parameterise the complex scalar as Then J corresponds to the (massless) Goldstone boson, the Majoron, and s is the radial excitation.
In this parameterisation, J is not present in the potential and appears in the Lagrangian only through the kinetic term and the Yukawa interactions in Eq. (37). Further, we can rotate the fields carrying non-zero lepton number, namely, Ψ = N R , χ c L , L and e R , as and remove J from all Yukawa interactions. After this transformation, the kinetic terms for the fermions Ψ will induce In this way, the derivative nature of Goldstone boson's couplings is manifest. It is worth noting that, despite having D = 5, this operator is not related to integrating out a heavy mediator, but is a consequence of the non-linear field redefinition performed in Eq. (40).
In the spirit of the EFT approach we are pursuing in the current study, it is interesting to see which effective operators are generated at low energies if s is heavy (for concreteness, we assume that its mass is larger than the electroweak scale). In what follows, we work in the unbroken phase of the electroweak symmetry. Minimising the potential in Eq.
It is interesting to note that the |H| 6 operator is not generated at tree level due to a peculiar cancellation coming from the s 3 and s 2 |H| 2 terms in the potential upon using the equation of motion for s and the relation between m s and v σ . (This had been previously noted in Ref. [61].) Finally, the parameters of the SM potential, get shifted as

Model C1: massive vector boson
The vector form of the operator O 1 , see Eq. (4), suggests that it can be generated by the exchange of a heavy neutral vector boson, Z µ . The Lagrangian that could (effectively) describe such an exchange has the following form: where Z µν is the corresponding field strength tensor. In general, kinetic mixing among Z and Z, Z µν Z µν , as well as mass mixing, δm 2 Z µ Z µ , are also allowed. However, to ensure the neutrino portal regime (and reduce the number of independent parameters), we will set and δm 2 to zero.

Model C2: gauged U (1) B−L
It is well known that promoting U (1) B−L to a local symmetry requires the addition to the SM particle content of three RH neutrinos to cancel gauge anomalies. 11 In the considered case, one of them is traded by the chiral fermion χ c L odd under Z 2 . This is why it is important to have L(N 1,2 R ) = L(χ c L ) = 1. For concreteness, we assume that one of the two sterile neutrinos is lighter than DM, whereas the second one has a mass around the scale of U (1) B−L symmetry breaking, v σ . The Lagrangian of this model reads 12 where in L B2 the (covariant) derivatives are modified to include the piece associated with the new gauge symmetry: with g and Q B−L being, respectively, the new gauge coupling and the B − L charge of the field D µ acts upon. We provide the B − L charges of this model's fields in Tab. 4. Upon spontaneous 10 For the conditions of applicability of such type of models see e.g. Ref. [62]. 11 This is not the only accidental (global) symmetry of the SM that can be gauged. It is well known that differences of individual lepton flavour numbers, such as Lµ − Lτ , are also anomaly free in the pure SM (with no RH neutrinos) [63]. In SM extensions with additional fermions (like our scenario), Lµ − Lτ and its variants can also be gauge symmetries if the new fermions have the proper charge assignments. However, flavour symmetries have strong implications for neutrino masses and mixings, and thus, deserve further studies. Therefore, although they may have interesting implications, in the following we discuss the flavour-blind symmetry B − L. Let us stress that a distinct feature of this gauge symmetry is that cancellation of gauge anomalies calls for the addition to the SM of three chiral fermions, unlike the cases of gauged differences of individual lepton flavour numbers. 12 In general, the kinetic mixing term among Z and B, Z µν B µν , is allowed. Here B µν is the field strength tensor of U (1)Y . However, to ensure the neutrino portal regime we assume that the physical kinetic mixing is negligible. breaking of U (1) B−L , Z , N R and χ L acquire their masses: . From an EFT point of view, if v σ is larger than the weak scale and the couplings g and λ σ are not too small, we can integrate out both s and Z . It is convenient to go to the unitary gauge, rendering σ from Eq. (38) real in each point of spacetime. In this gauge, the would-be Goldstone boson J is removed from the theory. We list in Tab. 3 the Wilson coefficients of the operators generated in the EFT. As in Model B2, we find c 2 = −2c * 3 , so the contributions of these operators to the s-wave part of the cross section for χχ → N N cancel, see Eq. (15). It is interesting to note that c 4 and c 5 vanish if λ σ = f 2 and λ σ = g 2 , respectively (see Tab. 3).
Apart from the operators summarised in Tab. 3, we find the following four-fermion interactions: where ψ stands for the SM fermions, i.e. ψ = L, e R , Q, u R , d R . The corresponding Wilson coefficients read: Here Q ψ , Q N and Q χ denote the B − L charges of ψ, N R and χ L , respectively, see Tab. 4. The phenomenology of this model has been investigated in detail in Ref. [64], and its parameter space has been shown to be severely constrained. For other studies exploring DM-neutrino connections in gauged U (1) B−L models, see references in [64] as well as [65,66].

Phenomenology of selected renormalisable models
In this section, we study the phenomenology of the UV completions presented in Sec. 3. The main focus of this work is on the neutrino portal regime, where the relic abundance is set by the DM annihilations into sterile neutrinos χχ → N N , with no connection with the SM through the Higgs or vector portals. However, we need some interaction that guarantees the thermal equilibrium of N with the SM particles in the early Universe, as has been mentioned in Sec. 2.2.1. Therefore, we assume a (small) value for the Higgs portal coupling, 10 −6 λ φH 10 −3 , that keeps the dark sector in kinetic equilibrium with the SM up to a certain temperature (see App. A). 13 For genuine Models A, this coupling does not affect the DM relic abundance. However, in Models B2 and C2, when σ and H develop VEVs, the two scalars mix, and the coupling λ σH as small as 10 −4 -10 −3 would contribute significantly to the relic abundance around the resonance m χ = m h /2, where m h is the mass of the physical Higgs boson, see e.g. Ref. [44]. We will not discuss this effect in what follows.
We summarise the main phenomenological features of the models described in the previous section in Tab. 5, while some further characteristics are described below.
• Type-A models. The Higgs portal term could produce DD signals at one loop. For m χ < m h /2, this coupling would also induce invisible Higgs decay, h → χχ. In any case, if λ φH 1, these constraints are evaded. Furthermore, there will always be a one-loop contribution to DD through the exchange of Z boson. If m χ < m Z /2, this will also lead to invisible Z decay, Z → χχ. However, both processes are suppressed by the small neutrino Yukawa coupling. For a detailed analysis of such one-loop contributions see Ref. [67]. Model A2a has interesting features avoiding ID bounds because the DM annihilation cross section is p-wave, and in Model A2c finite m N is generated at one loop. We will discuss the latter in Sec. 4.2.
• Type-B models. For Model B1, the most general scalar potential written in Eq. (32) includes the µ φH term, that generates a VEV for the scalar φ upon EWSB. In that case, there is a mixing between the scalar and the Higgs, and elastic scattering of DM off nuclei occurs at tree level. In addition, if m χ < m h /2 we also have invisible Higgs decay. However, as we are interested in the neutrino portal regime, we take µ φH = 0 in the phenomenological analysis; also the small value for the Higgs portal term does not generate any additional contribution to the relic abundance with processes involving SM particles, except for around the resonance m χ = m h /2, if λ φH 10 −4 .
In Model B1, if the Yukawa coupling g of DM to the scalar mediator is real, the annihilation cross section for χχ → N N is p-wave. This is a reflection of the fact that DM requires a pseudo-scalar coupling (given by the imaginary part of g) to annihilate through s-wave, see e.g. Ref [68]. In Model B2, since both f and g are real, the annihilation cross section is also p-wave. Thus, in these cases, the ID limits are avoided.
Finally, this kind of models has four-fermion self-interactions of N and χ. However, the DM self-interactions χχ ↔ χχ are very suppressed in the parameter space considered in our analysis, i.e. σ χχ→χχ /m χ 10 −6 cm 2 /g for Model B1, well below current limits [69].
• Type-C models. Model C1 is not UV-complete and has to be understood as an effective description of the interaction between N R and χ L via a new massive vector boson. Note that in the presence of kinetic mixing, there would be other processes like Z → χχ if m χ < m Z /2. Even if the tree-level parameter = 0, this kinetic mixing will be induced at one loop but will be further suppressed by the small neutrino Yukawa coupling. Model C2, instead, leads to direct interactions of DM with the SM particles, cf. Eq. (50), and thus, it is severely constrained [64].
Given that the phenomenology of Models A1, B2 and C2 have been studied in detail in Refs. [2,3,5,52], [59] and [64], respectively, and that Model C1 is not UV-complete, we will analyse in detail Models A2 (b and c) and B1 in the next subsection. Moreover, in the last part we will comment on the one-loop generation of the RH neutrino mass in Model A2c.

Feature
Model A1 A2a A2b A2c B1 B2 C1 C2 s-wave σv χχ→N N * DD @ tree level Self-interactions Table 5: Classification of the phenomenological features of the UV completions discussed in Sec. 3. Notice that when the DM annihilation cross section σv χχ→N N is p-wave, one can easily avoid indirect bounds. The asterisk * in Model B1 implies that if the coupling g is real, the thermally averaged DM annihilation cross section is p-wave.

Dark matter phenomenology
In this subsection, we focus on models A2b, A2c and B1. As discussed in Sec • Model A2b, with coupling f real: • Model A2c, with coupling f real: • Model B1, with couplings f and g complex: We consider the following two situations: (i) Real f and g, for which (ii) Purely imaginary f and g, resulting in We see that, in all the cases except Model B1 with real couplings, the DM annihilation cross section is s-wave.
In Fig. 4, we depict the thermally averaged DM annihilation cross section corresponding to the observed value of DM relic abundance Ωh 2 = 0.12 [70] using the results from Ref. [45]. 14 The cross section is computed both (i) in a renormalisable model (red line) and (ii) in the corresponding EFT (blue dotted line), using the matching conditions from Tab. 3. Therefore, red and blue dotted lines correspond to the values of σv that reproduce the observed relic abundance, as explained before. The four panels correspond to Models A2b (top left), A2c (top right) and B1 with real (bottom left) and purely imaginary (bottom right) Yukawa couplings. 15 We observe that the EFT approach to the calculation of DM relic abundance works for m χ m σ,θ /4 in the type-A models and m χ m φ /6 in Model B1. For the latter, the resonance behaviour of the cross section clearly cannot be captured by the EFT. Our analytical results for the relic abundance agree with the results obtained using micrOMEGAs [71,72]. Regions in red stand for values of the relic abundance which would overclose the Universe, i.e. Ωh 2 > 0.12. We show the values of the parameters that are fixed for each model in the upper region of the plots, taking for the RH neutrino mass in models A2b and B1 the minimal value allowed by the BBN constraint, i.e. m N = 2 GeV. As we will discuss in Sec. 4.2, in Model A2c the sterile neutrino mass m N is generated radiatively. In Fig. 4, we take µ σ = 10 4 GeV in order to have m N 2 GeV in the part of the parameter space of interest, and contours of fixed m N (in GeV) are also shown as black dotted lines. In addition, the brown region is excluded by the BBN constraint. We assume f = 1 as an illustrative example for models A2b and A2c without loss of generality. For Model B1, we distinguish between the cases of real and purely imaginary couplings, f r = g r = 1 and f i = g i = 1, respectively.
We also add ID constraints from Ref. [3], which we briefly discuss below. Planck cosmic microwave background (CMB) measurements set bounds on the DM annihilations into SM particles. The related production of particles leads to homogenisation of the CMB power spectra and the modification of the ionisation history of the Universe. In addition, the Fermi analysis of dwarf spheroidal galaxies (dSphs) provides limits on the DM annihilation cross section by non-observation of excess above the astrophysical backgrounds in the gamma-ray flux, with photon energies in the 500 MeV-500 GeV range. In the models discussed here, DM has negligible interactions with the SM particles, in particular with quarks. Therefore, it is not captured in the Sun and no associated ID constraints exist. The bounds from CMB and dSphs are represented in the plots by the blue and orange hatched regions, respectively. Notice that for Model B1 in the case of real g (coupling Finally, for the parameters of the models that we take in Fig. 4, the white regions correspond to points that avoid all the experimental bounds and provide some fraction of the total relic abundance of DM. For Model A2b, a small region of the parameter space with 100 GeV m χ 300 GeV and 200 GeV m σ 300 GeV is open. In Model A2c, DM masses between approximately 100 GeV and 800 GeV, in conjunction with 300 GeV m θ 800 GeV are allowed. In this case, the values for the one-loop generated RH neutrino mass are 2 GeV m N 10 GeV. For Model B1, due to the resonance behaviour of the annihilation cross section, a larger part of the parameter space is open. For real couplings, χ with mass between 2 GeV and 10 TeV, and φ with mass between 2 GeV and 20 TeV can be responsible for the totality of observed DM relic abundance, whereas for purely imaginary couplings, the allowed interval of DM masses is 30 GeV m χ 50 TeV and that of the scalar mediator masses is 1 TeV m φ 100 TeV.
Differences in the results shown in Fig. 4 could, in principle, come from allowing other particles of the dark sector to evolve out of the thermal equilibrium and enabling particles to decay, as it is detailed in App. A for Model A2b. However, the evolution of the full set of Boltzmann equations shows that the deviations are not significant in almost all of the parameter space analysed in this work.

Neutrino masses in Model A2c
The Lagrangian of the model is given in Eq. (28). The U (1) L symmetry, under which the complex scalar σ has charge (+1), is softly broken by a quadratic term in the potential, whereas m N = 0 at tree level. The soft breaking term splits the masses of the real, ρ, and imaginary, θ, components of σ, see Eq. (30). We choose µ 2 σ > 0, such that m ρ > m θ . The lighter of the Z 2 -odd fields, i.e. either χ or θ, yields a DM candidate. However, here we only focus on fermionic DM. We notice that other U (1) L breaking terms are possible, as e.g. λ σH σ 2 |H| 2 or m N N c R N R , but they are harder (higher dimension). Even if these terms are absent at tree level, finite contributions to both of them are generated at one loop. In the case of λ σH , this reads The splitting in the masses of ρ and θ leads to a finite m N being generated at one loop, see Fig. 5. A similar mechanism has been proposed in Ref. [20]. It is analogous to that of the scotogenic model [73] and its generalisations [74][75][76], where light (mostly-active) neutrinos acquire their mass through one-loop diagrams, for a review see Ref. [77]. After EWSB the tree-level Majorana mass term for electrically neutral fermions reads where m D = y ν v h / √ 2. The DM candidate χ is decoupled from the rest of neutral leptons, since it is charged under the Z 2 .
For the computation of m N at one loop, we assume that there are n N generations of N R and n χ generations of χ L . Furthermore, we work in a basis for χ L in which m χ is diagonal with positive and real elements m χ k . Performing the computation (which is identical to the one in the scotogenic model), we find where the loop function F is defined as follows: From this formula, it is evident that in the limit of lepton number conservation (µ 2 σ = 0, and hence, m ρ = m θ ), the mass matrix m N = 0. In fact, in the limit m χ k m ρ , m θ and m ρ m θ m σ one finds This result, which depends linearly on m χ k , can easily be estimated by using dimensional analysis and symmetry arguments.
Depending on the number of generations n N and n χ , some of the RH neutrinos may remain massless. For instance, if n N = 3 and n χ = 1, only one of the eigenvalues of m N is non-zero, since m N ∼ f * f † in this case. n χ = 2 leads to two massive N , a minimal number (within the type I seesaw mechanism) needed to explain low-energy neutrino oscillation data. For n χ = 3, all three RH neutrinos get masses. In the case of only one generation n N = n χ = 1, and assuming m N m D , we can express the coupling f as with m ν ∼ 0.05 eV the mass of the light active neutrino and v h = 246 GeV the Higgs VEV.

Summary and conclusions
Motivated by the lack of WIMPs signals, in the present work we have revisited the possibility of SM singlet DM interacting primarily with sterile neutrinos. The latter can explain the light neutrino masses. We have extended the SM with a Majorana fermion χ and RH neutrinos N R , assuming that their interactions are described by effective four-fermion operators. The stability of χ is ensured by a Z 2 symmetry. Restricting ourselves to the case in which DM interacts with the lightest of sterile neutrinos, we have shown that there are three independent four-fermion operators. One of them, O 1 , always preserves lepton number, whereas the remaining ones, O 2 and O 3 , may either preserve or violate it (depending on the lepton number of χ L ). We refer to O 1,2,3 as sterile neutrino portal operators.
Assuming that the mass of the lightest sterile neutrino, m N , is smaller than that of χ, the observed DM relic abundance can be entirely explained by the freeze-out of χ due to the annihilation process χχ → N N triggered by the portal operators. For O 1 , the s-wave part of the corresponding annihilation cross section is proportional to m 2 N , and thus, is suppressed for small values of m N . For O 2 and O 3 , there is no such a suppression. Turning one operator at a time we have derived the scale of new physics Λ required to reproduce the observed relic abundance.
Further, we have formulated simple UV completions that lead to one or more portal operators when integrating out a heavy mediator at tree level. Depending on the Lorentz nature of the mediator and whether it propagates in t-or s-channel of the χχ → N N process, we classified the UV models into • Model A1 (A2) containing a real (complex) scalar φ (σ) propagating in t-channel; • Model B1 (B2) involving a real (complex) scalar φ (σ) propagating in s-channel; • Model C1 (C2) having a massive vector (gauge) boson Z propagating in s-channel. 16 In Models A, the mediator is charged under the Z 2 symmetry stabilising DM, whereas in Models In Model A2a, where m N = 0 and light neutrinos are Dirac, the annihilation cross section is effectively p-wave, and the ID bounds from annihilation to neutrinos [10] are avoided. This scenario is interesting since it allows for light thermal DM, with masses as small as 100 MeV. However, it is very difficult to probe it. In Models A2b and A2c, m N = 0 at tree and one-loop level, respectively. Model A2c possesses this interesting feature of finite m N being generated at one-loop level, analogously to the generation of light neutrino masses in the scotogenic model. In the limit of DM mass being smaller than the mass of the scalar mediators, m N ∼ m χ /(4π) 2 . In other models considered in the present work, m N is a free parameter. In any case, it should be larger than approximately 2 GeV for N to decay before BBN. Decays of N will modify the spectra of charged particles (in particular, antiprotons and positrons) as well as photons. These modifications can be looked for in ID as discussed in Ref. [3]. We have adopted the constraints from CMB measurements by Planck and dSphs observations by Fermi derived in Ref. [3], showing that in Model A2b (A2c) a Majorana fermion χ with m χ between approximately 100 GeV and 300 GeV (800 GeV) can constitute 100 % of the observed DM relic abundance, respecting the ID constraints.
In Model B1, the annihilation cross section is p-wave if the Yukawa coupling of DM to a scalar mediator, g, is real. In this case, the ID bounds are avoided, and DM masses between 2 GeV and 10 TeV are allowed. On the contrary, for complex g, the annihilation cross section is s-wave, and thus, the ID constraints apply. For purely imaginary f = g = i, we find that the viable range of DM masses is 30 GeV m χ 50 TeV. Larger DM masses in this model are accessible due to the resonance enhancement of the cross section.
In conclusion, we have shown that DM-sterile neutrino interactions described by effective fourfermion operators constitute a viable option. They can be generated in a number of UV-complete models possessing somewhat distinct phenomenology. This scenario provides a possible connection between neutrinos and dark matter, which arguably are among the most feebly interacting sectors of nature.

A.1 Departure from chemical equilibrium within the dark sector
In this appendix we review the comparison of calculating the relic abundance with the full set of Boltzmann equations (BEqs) and the standard approximation (STD). We denote as STD the case considering just the evolution of the DM particle χ and the 2 ↔ 2 processes. For example, in Model A2b described in Sec. 3.1.2, the complete BEqs for m θ = m ρ > m χ , m N , can be expressed as: in terms of the yields Y i = n i /s, where n i is the number density for species i and s is the total entropy density, x = m χ /T , H is the (x-dependent) Hubble rate, and the superscript "eq" denotes equilibrium distributions with zero chemical potential, as in Refs. [47,49].Γ i are the thermal decay rates for a decaying particle i given bỹ and Γ i are the zero-temperature decay rates. When m N < m h (T ), decays into N should be taken into account by means of the following substitutionΓ N →Γ h with, Here we consider the decay with the approximation as in Ref. [78], where all the four states of the Higgs doublet have the Higgs boson mass m h (T ); thermal masses were taken from [52]. The thermally averaged cross section σv is given by [47]: where K 1,2 are modified Bessel functions, and m is the DM mass. In Fig. 6, the relative deviation of the BEqs' result from the STD approach is shown by plotting the quantity (Ω − Ω STD )/Ω STD [%] as a function of m χ /m σ,N for Model A2b, for given values of f = 1, m σ = 10 3 GeV and m N = 2 GeV (left) and 100 GeV (right), in three different cases: (i) assuming that N is in equilibrium, and following the evolution of σ and χ (orange line); (ii) same as before but taking into account 1 ↔ 2 processes (grey line); and finally (iii) the general case solving the BEqs for σ, χ and N , i.e. assuming that N is not in equilibrium, and also allowing for 1 ↔ 2 processes (brown line). Notice that we show the case m N = 100 GeV in order to illustrate the case when masses for χ and N are degenerate, m χ ∼ m N . Departures from Ω STD in the plot can be understood as follows.
From Eqs. more significant when N and/or σ and the DM start to become non-relativistic at nearly the same temperature, so for similar masses. Then it tends to make the freeze-out happen earlier. 18 This is shown in the plot when the three lines depart from zero for m χ ∼ m σ or when the brown line does so for m χ ∼ m N . Indeed, the latter feature is present only in the right panel of Fig. 6, where m N = 100 GeV is relatively close to m χ . Moreover, the addition of decay widths allows for the production/decay of N from/to SM particles, and the decays of σ to χ and N . This would produce, in contrast, the opposite effect by making the particles follow the equilibrium for longer if the decaying particle is not excessively Boltzmann suppressed. This can be noticed by the fact that the grey line is closer to zero than the orange line. In conclusion, the deviation of the full BEqs' result from the STD is below 5 % (10 %) in almost all of the parameter space for m N = 2 (100) GeV, except for m χ ∼ m σ (m χ ∼ m σ , m N ).

A.2 Departure from kinetic equilibrium of the dark sector with the SM
Here we briefly discuss kinetic decoupling of the dark sector from the SM. For concreteness, we focus on Model A2b. In Fig. 7, we show the thermal rates of the most relevant 1 ↔ 2 and 2 ↔ 2 processes, normalised to the Hubble rate. The values of the fixed parameters correspond to a point in Fig. 4 yielding the observed DM relic abundance and avoiding all the experimental constraints. As can be seen, the neutrino Yukawa coupling fixed by the seesaw relation in Eq. (10) to approximately 5 × 10 −8 (for m N = 2 GeV and m ν 0.05 eV) cannot keep the dark sector in equilibrium with the SM. At the same time, the Higgs portal coupling λ σH does ensure kinetic equilibrium between the dark sector and the SM as long as it is larger than ∼ 10 −6 , at least in some range of temperatures. For the example shown in Fig. 7, kinetic decoupling of the dark sector from the SM happens before the DM chemical freeze-out. However, kinetic equilibrium within the dark sector is maintained through χN ↔ χN process. 19 From the moment of kinetic decoupling, the dark sector and the SM bath evolve with two different temperatures, T D and T , respectively. We assume that entropy is conserved independently in both sectors [43]: with s D and s SM being the entropy densities of the dark sector and the SM bath, respectively, and T kd the temperature of kinetic decoupling. The evolution of ξ ≡ T D /T can be obtained using s SM = (2π 2 /45)g * (T )T 3 and s D = (ρ D (T D ) + p D (T D ))/T D , where ρ D and p D are the energy and pressure densities of the dark sector, and g * is the effective number of relativistic degrees of freedom in the visible sector. In Fig. 8, we display the evolution of ξ as function of the dark temperature, T D , for the same point in the parameter space as in Fig. 7. We are interested in the value of ξ at chemical freezeout. For blue (red) line corresponding to the Higgs portal coupling λ σH = 10 −3 (4 × 10 −6 ), kinetic decoupling takes place at T D ≈ 20 (200) GeV, cf. also Fig. 7. As can be seen from Fig. 8, if T kd < m χ , the temperature of the dark sector is very similar to that of the SM bath, whereas if T kd > m χ , the ratio of temperatures reaches approximately 1.2 at the freeze-out of DM. In both cases, N is relativistic at the freeze-out, and according to Ref. [43], the DM relic abundance is modified with respect to the standard solution by a factor ξ g eff * /g * , where g eff * = g * + g D ξ 4 , with g D being the effective number of relativistic degrees of freedom in the dark sector. Since g D g * and 1 ξ 1.2, the correction to the relic abundance can reach up to approximately 20 %. 20 We have checked that this holds in a large part of the parameter space presented in Fig. 4; e.g. if λ σH = 10 −3 , this condition is fulfilled for m χ 10 GeV. After the freeze-out and until N becomes non-relativistic, ξ is constant, whereas after T D drops below m N , the ratio ξ ∝ T [43]. More precisely, ξ can be expressed as ξ = ξ fo T D /m N , where ξ fo is the value of ξ at the freeze-out.
If sterile neutrinos become non-relativistic before the freeze-out, the dark sector may be significantly reheated. In that case, there could be an order one correction to the relic abundance. For a precise determination of the relic abundance in the presence of decoupled dark sectors and the cases where the impact of such a decoupling can be sizeable we refer the reader to Refs. [44,45].