Right-handed neutrino dark matter under the B − L gauge interaction

We study the right-handed neutrino (RHN) dark matter candidate in the minimal U(1)B−L gauge extension of the standard model. The U(1)B−L gauge symmetry offers three RHNs which can address the origin of the neutrino mass, the relic dark matter, and the matter-antimatter asymmetry of the universe. The lightest among the three is taken as the dark matter candidate, which is under the B − L gauge interaction. We investigate various scenarios for this dark matter candidate with the correct relic density by means of the freeze-out or freeze-in mechanism. A viable RHN dark matter mass lies in a wide range including keV to TeV scale. We emphasize the sub-electroweak scale light B − L gauge boson case, and identify the parameter region motivated from the dark matter physics, which can be tested with the planned experiments including the CERN SHiP experiment.


Introduction
In understanding nature, the gauge symmetry and its spontaneous breaking play a core role. The standard model (SM) of particle physics is an extremely successful model so far in explaining the data. Needless to say, its beauty is ascribed to the gauge principle which not only regulates the interactions among particles but also organizes the content of particles by means of the anomaly cancellation conditions.
There are, however, various issues that the SM fails to address. For instance, although the existence of the dark matter (DM) is quite certain to explain many independent astrophysical observations, it is convinced that the dark matter does not belong to the SM, and its identity has been still unknown. Due to the fact that the neutrinos are massive, there plausibly exist their chiral partners, the right-handed neutrinos (RHNs), which are also not a part of the SM. Unlike the other SM fermions, they can be Majorana particles which can exploit the seesaw mechanism to explain their small masses [1][2][3][4]. Through the seesaw mechanism, the RHN can stay effectively as a sterile neutrino, decoupled from the active neutrinos.
The RHN, which is neutral under the SM gauge symmetries, has a potential to be a viable dark matter candidate. This has been realized in the ν minimal standard model (νMSM) [5,6], 1 which sets the lightest RHN (N 1 ) mass around keV scale such that it can be naturally long-lived against its decay, N 1 → νγ, induced through the mixing between JHEP02(2017)031 N 1 and active neutrinos, where the mixing angle is conventionally denoted by θ 1 . The framework of the νMSM can also address the baryon asymmetry of the universe (BAU) [6] through the GeV scale RHNs and active neutrino oscillations [10].
In the νMSM, the sterile neutrino DM can be produced through the mixing between the N 1 and the active neutrinos, which is known as Dodelson-Widrow mechanism [11] (see also refs. [12,13]). However, the non-observation of the X-ray signal from the N 1 decay (N 1 → νγ) [14] and the phase space density constraint on the N 1 mass [15] excluded this simple approach (for the Lyman-α forest constraint, see, e.g., ref. [16]) except for turning to the resonant effect which requires an anomalously large lepton asymmetry [17]. As another way out, introducing extra interactions can provide a viable dark matter production mechanism that is independent of the mixing angle θ 1 . For instance, the N 1 can be produced by the decay of a scalar particle [18][19][20][21][22][23][24][25][26][27][28] through the freeze-in mechanism [29,30]. (For a discussion on the freeze-in scenario for the hidden sector dark matter that communicates with our sector through the kinetic mixing and/or the scalar mixing, see refs. [31,32].) In this paper, we present the minimal U(1) B−L gauge extension of the SM with the RHN dark matter candidate, which we call the U(1) B−L extended νMSM or the UνMSM. We also explore a comprehensive picture of the sterile neutrino DM candidate in this model. In the light of the success of the gauge principle in the SM, the U(1) B−L gauge symmetry is expected to play a similar role for the DM. 2 In fact, due to the anomaly cancellation conditions, the U(1) B−L regulates the number of the RHNs to be three. The lightest RHN can be a DM candidate with its mass scale from keV to TeV, or even higher. The other two RHNs may be responsible for the BAU, which will be studied elsewhere. The interaction can be mediated by both a B − L gauge boson Z and an associated scalar S that is associated with the spontaneous symmetry breaking. To gain the control in the number of free parameters, we will consider only the Z interaction in this work, unless specifically stated, which is valid in the limit the S is heavy enough and/or inefficiently communicate with the SM sector so that its contribution to the DM production is greatly suppressed. The new gauge interaction can play an important role in the sterile neutrino production especially via the Z mediated freeze-in mechanism, and provide distinguishable implications that can be tested experimentally.
There are some related works such as refs. [38][39][40]. They impose a Z 2 protective symmetry on some sterile neutrino while requiring two others to accommodate realistic neutrino phenomenology. In this scenario, the sterile neutrino can be an ordinary cold DM candidate around the weak scale, i.e., it has a weak interaction (say, the U(1) B−L gauge interaction) and gains a correct relic density via the conventional freeze-out mechanism. In our study, the most interesting case (also the main case) actually is a very light RHN which does not necessarily call for a Z 2 protective symmetry, although for the sake of a global picture we also include the heavy RHN dark matter case, which then may require a flavor symmetry as in refs. [38][39][40]. We also exploit the freeze-in mechanism to account for correct relic density for the RHN DM. A scalar DM candidate in a similar framework was studied JHEP02(2017)031 in ref. [41]. We also note a larger gauge group SU(3) C × SU(2) L × SU(2) R × U(1) B−L based on the Left-Right gauge symmetry was considered before [42,43]. Heavy gauge bosons (W ± R and Z ) and usual freeze-out production method with a dilution was used, which is a different approach from ours. The rest of this paper is organized as follows. In section 2, we describe our framework, the UνMSM. In section 3, we discuss possible DM production scenarios and the relevant constraints on the model. In section 4, we discuss implications for various phenomena including the SHiP experiment. In section 5, we summarize our study.

The framework of the UνMSM
Following the success of the gauge principle in the SM, we consider a model with the U(1) B−L gauge symmetry as a minimal choice in terms of the matter contents, which offers three RHNs N i , a U(1) B−L gauge boson Z , and a single scalar Φ S being responsible for spontaneous breakdown of U(1) B−L . The Lagrangian of the UνMSM is given by where α = e, µ, τ , i = 1, 2, 3, and D µ = ∂ µ − ig B−L Q Z µ with g B−L and Q being the B − L gauge coupling and B − L charge (Q = −1 for the SM leptons and N 's, Q = 1/3 for the SM quarks, Q = 2 for Φ S ). Z µν is the field strength of the Z . We take fourcomponent fermion notations, by which N i represents a four-component fermion having only the right-handed part. The gauge kinetic mixing of ( /2)Z µν B µν is highly constrained, and for the simplicity we take it zero in this paper. The gauge kinetic mixing [44] has been a great source of research interests in the past decade [45] and also branched out some variant forms such as the one in ref. [46]. See ref. [47] for the details of the physics related to this term in the gauged B − L model.
The Higgs potential is given by where Φ H and Φ S develop the vacuum expectation values (VEVs), Φ H = v H and Φ S = v S , so that the electroweak and the B − L gauge symmetries are spontaneously broken. After diagonalizing the mass matrix, we obtain the masses

JHEP02(2017)031
for the physical states H and S, respectively, where the mixing angle is given by tan 2θ = 2λ The VEV of Φ S gives the mass of Z and N i as follows: The coupling κ i is in general a complex value. Our following discussion is, however, independent from its CP phases, and thus we take κ i as a real parameter in what follows. The N 2 and N 3 are not directly related to the DM production and their masses are not bounded by the DM relic density as in the νMSM, as the resonant production through the large lepton asymmetry is not necessary in this model. The dominant decay mode is N 1 → 3ν given by [48,49] Requiring the N 1 lifetime is longer than the universe age (τ U ∼ 13.7 × 10 9 years), we get the following constraint.
where we have used θ 2 1 α |y α1 | 2 v 2 H /M 2 N 1 from the the see-saw mechanism. Thus, the low mass of the N 1 (not too larger than the eV scale) can satisfy the DM lifetime constraint easily, but the heavier N 1 would require α |y α1 | 2 1 to be sufficiently stable. Although the heavier the N 1 DM may mean the less natural setup, we will include the heavier N 1 in our study that expands the relevant phenomenology significantly (e.g., see section 4). As a matter of fact, the N 1 will be stable as long as the α |y α1 | 2 0. In this limit, which might invoke a flavor symmetry like refs. [38][39][40], the lightest neutrino would be massless (m ν 1 = 0) or almost massless, which is still consistent with the experimental constraints [50]. Throughout the rest of this paper, we will discuss in the zero N 1 mixing angle (θ 1 = 0) limit, which also allows us to leave out of account the constraints from the X-ray observations with the N 1 → γ + ν process.

Dark matter production and constraints
We now turn to discussing how the B − L gauge boson Z makes an impact on the N 1 dark matter production. The dark matter scenario drastically changes, depending on whether the Z can decay into the dark matters (M Z > 2M N 1 ) or not (M Z < 2M N 1 ).
In the rest of this section, we will approach the dark matter issues from very general points. First, we will discuss how and where the N 1 and Z can be thermalized (section 3.1). Then, we will discuss various constraints including the Big Bang nucleosynthesis (BBN), lab experiments, and astrophysical bounds (section 3.2) before we discuss the dark matter relic density. Although some of the discussions and constraints may not be directly relevant to the parameter region that gives the right relic density for the N 1 , it might be still instructive JHEP02(2017)031 to have them as they might be relevant when we consider somewhat altered scenario such as the late time entropy injection. In section 3.3, we briefly go over the issues for the keV scale N 1 dark matter for the thermal production. We discuss mainly the non-thermal N 1 dark matter production for the M Z > 2M N 1 (M Z < 2M N 1 ) case in section 3.4 (section 3.5).
3.1 Thermalization of the N 1 and Z Before heading towards the production of the correct relic density of the N 1 , we describe how the dark sector, the dark matter as well as its portal Z , is thermalized. For the thermalization of the N 1 and Z , the relevant reactions among the N 1 , Z and the SM particles are (a) of which the reaction rates are denoted by r a , r b , and r c , respectively. The relevant formulae are given in appendix A. If r i (i = a, b, c) is larger than the Hubble expansion parameter, H = (g * π 2 /90) 1/2 (T 2 /M Pl ) with M Pl 2.4 × 10 18 GeV being the reduced Planck mass, at some time, the N 1 and/or Z enter the thermal bath (reaching the relative chemical equilibrium of the SM sector and/or dark sector). In the following discussion, we take the numbers of degrees of freedom for the energy density and the entropy density to be the same value g * since they are very close, and g * is evaluated as a function of the temperature according to refs. [51,52].
It should be noted that N 1 N 1 ↔ Z Z mediated by s-channel S also exists. As we will discuss later, however, S can be always heavier than the N 1 and Z in the parameter regions of our interest, and this process will be suppressed as we will take a very heavy S. For other possible processes, N 2 N 2 , N 3 N 3 ↔ N 1 N 1 mediated by S may become significant when κ i is strong. In such a case, SS ↔ N 1 N 1 may also be relevant for the thermalization. On the other hand, as we will see, we can take M N 1 , M Z < M N 2 , M N 3 , M S in the parameter region of our interest, and these processes can be omitted by taking a specific reheating temperature T R as max{M In what we follow we take this case for the sake of simplicity. In order to focus on the Z interaction, we turn off the other possible reactions involving scalars, such as HH, SS, SH ↔ N 1 N 1 , by taking S very heavy and λ HS vanishingly small in a similar way to refs. [18][19][20][21]. Figure 1a shows whether the N 1 is thermalized or not depending on the Z mass and coupling, where we take M N 1 = 10 keV for an illustration purpose. In the deep blue regions, the N 1 never enters the thermal bath; in the other regions, the N 1 becomes thermal at some time. For the thermal N 1 , there are two distinct regions depending on if the N 1 is relativistic (hot or warm dark mater case) or non-relativistic (cold dark matter case) at its decoupling temperature T dec In the light yellow region, the N 1 satisfies M N 1 /T dec N 1 < 1, namely, it is a relativistic particle, while in the light green region, the N 1 is a non-relativistic particle.
In thermalization of the N 1 , the reaction rate r a is the dominant contribution to take the N 1 into the thermal equilibrium with the SM particles. 3 As mentioned in the beginning of this section, the DM production is sensitive to the critical line thermalized even when g B−L is very small; there is also a region in the bottom right corner of the parameter space where the N 1 is not thermalized because the mediator Z is too heavy and suppresses the reaction rate. On the other hand, r a gets suppressed for M Z < 2M N 1 since the process (a) becomes off-resonance, and the required g B−L for thermalization becomes larger. Figure 1a would not change even if there is a late time entropy injection, and clearly illustrates the distinction between M Z < 2M N 1 region and M Z > 2M N 1 region.

Collecting relevant constraints
In figure 1b, we collect relevant constraints in the g B−L and M Z parameter space, for a choice of M N 1 = 10 keV (which is considered to be a conservative value for the lowest M N 1 [15]). The constraint from Big Bang nucleosynthesis (BBN) can eliminate a large part of the parameter space shown in dark blue. The existence of additional relativistic degrees of freedom can speed up the expansion of the universe, which leads to the earlier decoupling of the active neutrinos, and hence a higher yield of 4 He and so on. The extra radiation density is included in the conventional parametrization of where ρ γ is the photon energy density, and N eff counts 3 for three active neutrinos. 4 4 Here we ignore the flavor dependence of the neutrino decoupling temperature, and take T dec ν ∼ 1 MeV. In reality, νµ and ντ might decouple before νe, which would induce a small correction to N eff .

JHEP02(2017)031
In our case, the deviation from 3 contains the contributions from the N 1 and Z (if it is relativistic at T dec ν ∼ 1 MeV), which is given by where g * (T dec ν ) = 10.75. By demanding ∆N eff < 1 [53], we obtain the exclusion region shaded in dark blue for the range of 1 MeV M Z 10 MeV in figure 1b. For masses 2M N 1 M Z 1 MeV, we impose that the N 1 enters the thermal bath after T ∼ 1 MeV so that the N 1 does not affect the SM neutrino decoupling [54,55], which leads to the bound for the coupling, g B−L 3 × 10 −9 -10 −10 . 5 For M Z 2M N 1 and g B−L < several × 10 −6 , only the thermal Z contributes to ∆N eff because the N 1 is non-thermal.
The other individual constraints shown in figure 1b are following. 6 1. LEP experiments. The high mass regions are sensitive to the LEP experiments which give the exclusion limit depicted by the brown region [56]. The constraint for the contact interactions [57] is valid only for the M Z much larger than the collision energy at LEP, 209 GeV, while the initial state photon radiation process, e + e − → γνν [58], can be used for the M Z lower than the collision energy.
3. Beam dump (BD) experiments. The orange regions are excluded by the electron and proton BD experiments, where the regions from top to bottom correspond to E774 [60], E141 [61], Orsay [62], ν-Cal I (proton bremsstrahlung) [63], E137 [64], respectively. The black solid curve shows the expected reach of the SHiP experiment based only on the proton bremsstrahlung [8,65], which we will discuss in section 4. We have followed the method in ref. [66] to calculate the bounds from the electron beam dump experiments. For the proton beam dump experiments, the relevant calculation is shown in refs. [63,65]. 4. ν − e scattering at Borexino. The Borexino experiment has reported the interaction rate of neutrino-electron scattering from 867 keV 7 Be solar neutrinos [67]. The observed value is consistent with the SM prediction, which gives the bound denoted by the dark gray region, by imposing that the ratio between the cross section involving Z and the SM contributions should not exceed the maximum error [68]. This constraint is very powerful as it applies to a wide region of M Z . See also ref. [69] 5 When the thermalization temperature of the N1 is lower than the temperature at which the BBN is completed, observations of the light elements can not give any constraints. On the other hand, when the thermalization of the N1 occurs after the recombination (T ∼ 0.1 eV), the thermalized N1 may leave an imprint on the cosmic microwave background. This temperature range is beyond the scope of this paper though. 6 We did not take into account the Z → N1N 1 branching ratio for the BABAR, BD, SHiP, LEP bounds, which depend on it, and these bounds are taken as the same as figures 2-3. The change will be small nevertheless.

JHEP02(2017)031
for the similar level of the constraint from theν − e scattering based on the reactor experiments.
5. ν − q scattering at NuTeV. The mass range of Z above 10 GeV is constrained by the neutrino-nucleon scatterings. The NuTeV experiment measured ν µ (ν µ ) − q scattering, where ν µ andν µ were provided by the beamline at the Fermilab [70]. Since there is relatively large systematic errors, we take a conservative limit: M Z /g B−L > 0.4 TeV [55,71], which is depicted by the light yellow region.
6. Horizontal-branch (HB) stars. For the lighter Z , the energy loss rate of the stars in the globular clusters gives the more restrictive constraints, where the larger energy loss shortens the lifetime of the stars, and hence the observed population of the stars would be changed [72]. Here, we show the constraint from HB stars represented by the red region [73].
7. Supernova 1987A (SN1987A). The green region shows the constraint from the supernova explosion. The extra light particle taking energies from the center of the supernova can affect the signal duration of the neutrinos [72], in which the energy loss argument puts the bound [74]. An updated constraint [75], although not taken in our paper, is similar to the one in ref. [74] for the parameter regions we plot. (Cf. for a discussion on the potential way out, called the Chameleon effect, see ref. [76].) The latest LHC bound on the Z through the Drell-Yan process comes into the region above TeV scale [38], which is beyond the region of our interest, and we omit it in the figure.
3.3 Thermal production of the keV scale N 1 As a warm-up, we first consider a well-known case that the N 1 is around the keV scale, specifically 10 keV, which can be a candidate for a warm dark matter. 7 As one can see from figure 1a, the N 1 can be thermalized in the bulk space of the M Z − g B−L plane, and we concentrate on this case. The N 1 that once entered the plasma can be a warm or cold relic, depending on its mass and the decoupling temperature. The light yellow region in figure 1a indicates that the N 1 is relativistic (M N 1 /T dec N 1 < 1) at T dec N 1 , while it is non-relativistic (M N 1 /T dec N 1 > 1) in the light green region, where T dec N 1 is the decoupling temperature of the N 1 . When the N 1 is non-relativistic, T dec N 1 is lower than T dec ν , and the BBN constraint excludes this region. (The HB and SN1987A bounds also ruled out some part of this region independently.) When the N 1 is relativistic at T dec N 1 , the relic abundance of the N 1 is given by

JHEP02(2017)031
where n N 1 is the number density of the relativistic N 1 , n N 1 = (3/2)(ζ(3)/π 2 )T 3 , and ρ c = 1.05368 × 10 −5 h 2 GeV cm −3 is the critical density of the universe. s = (2π/45)g * T 3 and s 0 = 2889.2 cm −3 are the entropy density and its present day value. In this case, the abundance of the N 1 exceeds the observed value of the dark matter abundance Ω DM h 2 0.12 [53], and the universe is overclosed. This parameter space is depicted by the gray region above the dashed curve in figure 1b, excepting the non-relativistic region. Such a large abundance could be diluted if we take into account the late time entropy production by, e.g, the decay of N 2,3 as studied in refs. [42,43] although they employed a different gauge extension. 8 We note large parameter regions including a new window much below the weak scale can be viable in the case of the dilution, which has low energy laboratory implications. This can be compared to the refs. [42,43] where only the weak scale or above was considered. This is manifest in figure 1b, which shows that BBN, BD and BABAR already excluded a large portion of the parameter space, and the SHiP experiment is able to cover more space.

M Z > 2M N 1 case
We here discuss the case of M Z > 2M N 1 . It is well known that when the N 1 is around the electroweak scale while the Z is at TeV scale, the N 1 can be a thermal relic dark matter. This scenario was addressed in the context of the classically conformal models [77,78], and collider signatures of such a heavy Z were studied in, e.g., refs. [38,79]. We do not purse to study the thermal N 1 DM with a heavy Z in this paper.
On the other hand, there is another possibility that the N 1 is produced by the freezein mechanism [29,30], where the Z is produced as an on-shell state, and subsequently decays into a pair of the N 1 . In this scenario, the N 1 never enters the thermal bath, and is produced by the annihilations of a pair of the SM particles and also a pair of the Z if it is thermalized. This implies that the Z gauge coupling is very small compared to the thermal dark matter scenario.
We also require that the N 1 does not exist at the time when the universe is reheated up to the temperature T R after the inflation, namely n N 1 (T R ) n N 1 (∞) = 0, and thus the Boltzmann equation for n N 1 is given by where √ s is the center of mass energy. (K 1 is the modified Bessel function of the first kind.) For the annihilation cross section σv, the process (a) is the dominant contribution, which is given by

JHEP02(2017)031
where we have utilized the narrow width approximation. 9 Substituting eq. (3.5) to the right hand side of eq. (3.4), we obtain where we have used the yield Y N 1 ≡ n N 1 /s and d/dt = −HT d/dT , and take g * as a constant in the following. By replacing T with x ≡ M N 1 /T and integrating x from 0 to ∞ in eq. (3.6), we end up with the non-thermal abundance where f (τ ) = τ (1 − τ 2 ) 3/2 with τ = 2M N 1 /M Z taking 0 < τ < 1, and f (τ ) takes the maximal value f (τ ) 0.19 at τ = 2/5. We also approximate the total decay width as Γ Z ∼ C f g 2 B−L /(12π)M Z where C f is a coefficient in taking massless limit for final state particles. If Z decays into all the SM fermions (and N 1 ), C f becomes 7. We will approximate our results using C f = 7, and the parameter region where the right DM relic density is satisfied will be slightly changed if we use the exact values.
In figure 1b, we also depict the region of Ω nt N 1 h 2 > 0.12 as the gray region below the dashed curve. Therefore, the gauge coupling should be extremely small in order to obtain the observed dark matter abundance in this case, and it is quite challenging to test such a feebly interacting Z experimentally.

M Z < 2M N 1 case
Next, let us further focus on a possible dark matter scenario for M Z < 2M N 1 , where the BBN bound gets relaxed significantly because the reaction (a) is suppressed. This can be seen in the region M Z 20 keV in figure 1b, where the BBN bound on the gauge coupling becomes weak since the N 1 is hardly thermalized. In our setup, there are two scenarios for the dark matter depending on whether the relic abundance is produced in a thermal or non-thermal way. Figure 2 shows the N 1 relic density for a couple of examples of the M Z < 2M N 1 case. In the region above dashed curves in figure 2, the N 1 comes into the thermal bath at some time. In this parameter region, we find numerically the N 1 is always non-relativistic at the decoupling temperature T dec N 1 , and thus, we can evaluate the dark matter abundance in the same way as the usual cold dark matter case, which is given by 9 We here consider the case that TR is sufficiently large compared to the masses of the N1 and Z . As another possibility, the scenario with TR < M Z was discussed in refs. [80][81][82][83]. where the thermally averaged annihilation cross section, σv , includes the processes (a) and (c). The gray regions above the dashed curves in figure 2 show the parameter space of Ω th N 1 h 2 > 0.12, where we have given two benchmark cases, figure 2b]. In both cases, however, the thermal dark matter scenario is ruled out by various experiments such as the Borexino. 10 As a viable dark matter scenario, let us consider the non-thermal case where the N 1 is produced by the freeze-in mechanism discussed earlier. By demanding the condition n N 1 (T R ) n N 1 (∞) = 0, we obtain the abundance given by An important feature of this case is that the abundance is almost independent from the N 1 mass. To see this, let us approximately derive the analytical expression of the relic abundance. Since we consider the off-resonance reactions here and only the reaction (a) is sufficient in most of the cases, we can take σv ∼ g 4 B−L /(3πs). Substituting σv to the right hand side of eq. (3.11), we obtain dY nt

JHEP02(2017)031
It should be noted that the right hand side of eq. (3.12) takes the maximum value around T ∼ M N 1 , and thus, the produced number density is not sensitive to higher temperatures. Because of this, it turns out to be Y N 1 ∝ 1/M N 1 after integrating over the temperature, and hence the abundance is independent of M N 1 . By replacing T by x ≡ M N 1 /T and integrating over x from 0 to ∞ in eq. (3.12), we end up with the non-thermal abundance (3.13) This estimate well coincides with our numerical calculation shown as the gray regions below the dashed curves in figure 2, where the small fluctuations are caused by the temperature dependence of g * whose value is roughly given by g * (T ∼ max{M Z , M N 1 }).
We briefly comment on the BBN bound in figure 2, which is depicted by the dark blue regions. Since the BBN bound is sensitive only for the relativistic spices at around the neutrino decoupling temperature, it can eliminate up to M N 1 , M Z T dec ν . Below g B−L ∼ 10 −5 , the thermalization temperature of the N 1 and Z is below T dec ν or they never come into thermal bath, and thus the BBN can not constrain this region.
Before closing this section, we note perturbative unitarity on the coupling κ i for i = 2, 3, which can be expressed as , and κ 1 < κ 2 , κ 3 should hold as the N 1 is the lightest among the three in our setup. By demanding κ i < 4π, the masses of N 2 and N 3 are bounded from above as

Implications
In the non-thermal scenario for 2M N 1 > M Z , the dark matter abundance given by eq. It turns out that, e.g., for the mass regions 500 keV M Z 1 MeV and M Z 0.1 GeV, the scalar mass is at most 200 GeV M S 400 GeV and M S 4 TeV, respectively, when we take the perturbatively allowed maximum value λ S = 4π. Although scrutinizing the effect of the S is beyond the scope of this paper, our analysis is valid when we take λ HS vanishingly small so that the S does not come into the thermal bath and the non-thermal production of the N 1 through the decay of the S is sufficiently small [18][19][20][21].
For direct searches of the dark matter, the scattering between the N 1 and a nucleon can be induced by the Z and S mediated processes. However, since an effective operator of the Z mediation is given by (N 1 γ 5 γ µ N 1 )(qγ µ q), the scattering cross section is suppressed by velocity or momentum in the non-relativistic limit [84], which makes the measurement of this process difficult. We do not consider the S mediated process [39,40] as this interaction is turned off by taking λ HS 0 in this paper. Next, let us discuss possible experiments to test the freeze-in region with the right relic density (roughly, g B−L ∼ 10 −6 region) in figure 2. Beam dump experiments are a powerful tool to look for the small coupling regions. We estimate the expected reach of the SHiP experiment [8] using only the proton bremsstrahlung mode. The SHiP experiment utilizes the CERN SPS 400 GeV proton beam, where the Z can be produced via bremsstrahlung in proton scattering off the target. The SHiP is designed to have a 60 m muon shield and a 50 m detector area, and the Z decaying into the dileptons inside the detector may be observed. To estimate the signal events, we take the same kinematic parameters shown in ref. [65]. Following a similar approach to ref. [65], we take no background, in all figures 1-3 we show the expected region with the signal events more than three, which is depicted by the black solid curves. 11,12 Figure 3 shows the region around the BD constraints. In the blue area in figure 3a, the non-thermally produced N 1 can explain the correct relic density with the N 1 with M N 1 < M Z /2. The SHiP might barely have a chance to test this case when the M N 1 is near 10 keV. The blue band in figure 3b shows the case that the N 1 with M N 1 > M Z /2 can explain the whole amount of the observed DM abundance. The bottom side of this band is determined by taking M N 1 M Z /2, and the top side by taking M N 1 M Z . 11 The actual SHiP experiment sensitivity curves will be somewhat different from the ones given in our figures for the higher Z mass region as they should include additional production channels and the parton level analysis. 12 A study on how the decaying N2,3 signals with the B − L gauge boson can appear in the experimental searches at the LHC and the SHiP can be found in ref. [85].

JHEP02(2017)031
As a result, for the freeze-in region, the SHiP is expected to cover the mass range of 1 MeV M Z 200 MeV. The two cases in figure 3a and figure 3b are distinguishable in the sense that the Z decaying into a pair of the light DM is applicable in the former but not in the latter. While the blue area in the former case is not easily accessible with the planned experiments, the blue band in the latter case is quite well accessible partly because its coupling is larger.
There are some other forthcoming experiments that might be sensitive to our scenario. The NA64, one of the beam dump experiments at the CERN SPS looks for a missing energy carried away by a light gauge boson [86], and it may be sensitive to the M Z < 2M e region too as the Z can decay into the neutrinos in the B −L model. Also the Belle II experiments using the mono-photon trigger has a sensitivity that can cover 10 times smaller than the BABAR results in terms of the gauge coupling [87]. Detailed analysis for these experiments and developing methods to cover the whole blue regions in figure 3 are called for.

Summary and outlook
Success of the SM has been astonishing and it is amazingly consistent with high precision experiments. Yet, there are reasons to believe the complete description of nature requires the SM to be extended.
Following the success of the gauge principle in the SM, we have investigated the minimal gauge U(1) B−L extension of the SM, where three RHNs, a U(1) B−L gauge boson, and a singlet scalar are introduced. In particular, we have discussed the possibility that the lightest RHN N 1 is a dark matter candidate. Due to the presence of the Z interaction, the production mechanism of the dark matter does not need to rely on the mixing between active and sterile neutrinos, i.e., Dodelson-Widrow mechanism, and thus the X-ray bounds can be evaded.
For the keV scale dark matter, the U(1) B−L gauge interaction can bring the N 1 into the thermal bath, and thus the dark matter abundance is determined by the freeze-out mechanism. The produced N 1 is, however, relativistic at its decoupling in most parameter regions, which requires extra entropy production to dilute the overproduced number density. Note that even if the N 1 is never thermalized, non-thermal production such as the freeze-in mechanism can work. However, the produced number density is fairly small in this case.
As another viable possibility, we have considered heavier mass scales for the N 1 DM candidate based on two different relative mass spectrum: For the 2M N 1 < M Z case, we have discussed the freeze-in production of the N 1 , and found that extremely small g B−L is required for the correct number density for the DM candidate, which makes it difficult to be covered by the planned experiments except for a tiny region in the parameter space.
For the 2M N 1 > M Z case, the N 1 can be produced either in a thermal or nonthermal way depending on the parameter region. In the parameter region where the N 1 is thermalized, it is always non-relativistic at its decoupling, and thus becomes thermally produced cold dark matter in a typical way. The thermal abundance, however, requires a rather large gauge coupling, and such regions are already excluded by various experiments JHEP02(2017)031 (Borexino, etc.). On the other hand, a non-thermal production is still allowed for a smaller gauge coupling. We found that the appropriate value of the DM abundance can be obtained for g B−L ∼ 10 −6 largely independent of the N 1 mass. Interestingly, this parameter region (indicated as the blue band in figure 3b) can be sensitive to the planned beam dump experiments, and we found this freeze-in scenario can be tested by the light gauge boson searches at the SHiP experiment up to M Z ∼ 200 MeV.
We recall that the muon g − 2 favored parameter region (of the mass and coupling) in the dark photon scenario [88][89][90] has been a target of the active experimental searches in the past decade [45]. The parameter space was completely excluded by 2015 through the collaborative efforts of many different experiments [91]. The blue band in our study is specifically determined parameter region in our scenario (for the case the Z does not decay into a pair of the light DM), and some part of it is testable with the planned experiments. It would be worth investigating the possible ways to completely cover this parameter region of the B−L gauge boson, motivated by the the relic dark matter, the neutrino mass, and the BAU.
We have not scrutinized the interaction through an additional Higgs singlet assuming that λ HS is vanishingly small so that it does not contribute to the dark matter production. We also have not discussed the effect of the θ 1 , taking it negligibly small. The effect of these additional contributions will be discussed elsewhere.

A Reaction rates
We summarize the formulae used to calculate the relic abundance. The relevant processes are The squared amplitudes of these processes are given by

JHEP02(2017)031
where M f represents the SM particle masses, s, t, u are the Mandelstam variables, and N C is the color factor (N C = 3 for quarks, otherwise N C = 1). All the squared amplitudes are summed over spins. For the left-handed neutrinos, we take the massless limit in our numerical analysis. In particular, under this limit, the squared amplitudes of N 1 N 1 ↔ νν and Z Z ↔ νν become a half of eq. (A.4) and eq. (A.5), respectively. It should be mentioned that the full expression of the |M c | 2 would contain both the longitudinal component contribution which diverges in the high energy region, and the S contribution which cancels the divergence. Since the dark matter production discussed in this paper is not sensitive to the high energy behavior of M c , we did not include them in eq. (A.6). They will be presented and discussed in the subsequent work when we discuss the S contribution in detail.
The total decay width of Z is written by Γ Z of which the hadronic decay channel is obtained by Γ(Z → hadrons) = Γ(Z → µ + µ − )R(s = M 2 Z ) with R(s) being the usual R ratio (at a collision energy √ s) defined by R(s) = σ(e + e − → hadrons)/σ(e + e − → µ + µ − ) [92,93]. The partial decay widths are given by For the decay of Z into three massless neutrinos, its partial decay width becomes Γ(Z → νν) = 3g 2 B−L M Z /(24π). The reaction rates can be defined by using thermally averaged cross sections. For instance, the reaction rate of the process N 1 N 1 → ff is given by r a = σv(N 1 N 1 → ff ) ×n eq N 1 where n eq i = g i (2π 2 ) −1 M 2 i T K 2 (M i /T ) is the number density of particle i (having the mass M i and the degrees of freedom g i , e.g., g N = 2 and g Z = 3) in the equilibrium state. (K 2 is the modified Bessel function of the second kind.)