Origin of Sterile Neutrino Dark Matter via Vector Secret Neutrino Interactions

Secret neutrino interactions can play an essential role in the origin of dark matter. We present an anatomy of production mechanisms for sterile neutrino dark matter, a keV-scale gauge-singlet fermion that mixes with active neutrinos, in the presence of a new vector boson mediating secret interactions among active neutrinos. We identify three regimes of the vector boson's mass and coupling where it makes distinct impact on dark matter production through the dispersion relations and/or scattering rates. We also analyze models with gauged $L_\mu-L_\tau$ and $B-L$ numbers which have a similar dark matter cosmology but different vector boson phenomenology. We derive the parameter space in these models where the observed relic abundance is produced for sterile neutrino dark matter. They serve as well-motivated target for the upcoming experimental searches.


I. INTRODUCTION
The nature of dark matter is one of the biggest puzzles of our universe. Many experimental searches for dark matter are ongoing but with no success yet. New theoretical targets are needed to guide the path forward. Neutrinos, in many aspects, are the least well-known particles within the Standard Model owing to their weak interacting nature. In particular, neutrino self-interactions have never been directly measured in laboratories. Only indirect bounds have been inferred [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. There is room for neutrinos to interact with themselves much more strongly than the ordinary weak interaction. This also accommodates a potential connection to dark matter. We investigate such a possibility in the present work.
Sterile neutrino, a gauge singlet fermion that mixes with the Standard Model neutrinos, is one of the simplest dark matter candidates. It is requires very little effort of model building to realize such a candidate. Its cosmological longevity does not even require a symmetry but is simply attributed to the small mass and mixing. As fermionic dark matter, a sterile neutrino cannot be arbitrarily light [17,18], which grants the opportunity to observe monochromatic photons from its decay as an indirect detection signal [19,20]. The attractiveness of sterile neutrino dark matter is further enhanced by a novel finding of how its relic abundance could be produced. Dodelson and Widrow [21] pointed out that the same mixing angle allowing sterile neutrino to decay also enables it to be efficiently produced in the early universe. The corresponding parameter space has been scrutinized by various astrophysical probes and is virtually ruled out by the existing constraints [22][23][24][25][26]. A minimal solution to ameliorate this tension is to allow for a non-zero lepton-number asymmetry [27]. However, the required asymmetry does not have to be large and is challenging to observe elsewhere [28]. Alternative proposals include inflaton decays [29], freeze-in via additional particles [30,31], new interactions in the sterile sector [32,33], and gauge extensions of the Standard Model [34][35][36]. These mechanisms often make dramatic changes to the ultraviolet side of the story, i.e., the initial conditions of the early universe.
The focus of this work are theories that keep the most salient feature of the Dodelson-Widrow mechanism, the infrared dominance of dark matter production, and make testable predictions experimentally. The intuition is quite simple. Dark matter relic abundance is set by the product of the two ingredients in this mechanism, active-sterile neutrino mixing and the active neutrino reaction rate. While the former is strongly constrained by indirect searches, the latter is allowed to be much stronger than the weak interaction, especially when neutrinos are interacting with themselves. This idea was recently entertained within a model where neutrinos self-interact via a lepton-number charged scalar particle [37]. In the present work, we explore the impact of a new vector boson with coupling to neutrinos on the production of sterile neutrino dark matter. During the dark matter production epoch, such a vector boson could manifest as a heavy mediator, or a light and thermalized degree of freedom. These possibilities help us to discover new variations beyond the original Dodelson-Widrow mechanism. We consider three incarnations of the vector boson, in a neutrinophilic model, as well as in models with gauged L µ − L τ and B − L symmetries. We carry out detailed calculations of the relic density of sterile neutrino dark matter in each model, and derive the parameter space for the observed value to be reproduced. Furthermore, we confront our results with the existing and future experimental searches for each type of the new vector boson. This paper is organized as follows. In the next section, we give an overview of sterile neutrino dark matter, and use the tension between Dodelson-Widrow and astrophysical constraints as the motivation for introducing new neutrino forces. In section III, we first analyze a model with a neutrinophilic vector boson, deriving the parameter space for dark matter relic abundance and comparing it with a number of experimental probes. In sections IV and V, we do the same but in the gauged U (1) Lµ−Lτ and U (1) B−L models, respectively. We draw conclusions in section VI.

II. MOTIVATION FOR SECRET NEUTRINO INTERACTIONS
We start with the following Lagrangian which introduces a sterile neutrino field, where ν s is a left-handed fermion and a Standard Model gauge singlet, C = −iγ 2 γ 0 , and L T = (ν, − ) is the Standard Model lepton doublet with hypercharge −1/2. The parameter y ε is a tiny Yukawa coupling between ν s and L that makes a negligible contribution to the observed neutrino mass. The main role of this term is to generate a mixing between the active and sterile neutrinos, after electroweak symmetry breaking, H T = (0, v/ √ 2). In the case where both active and sterile neutrinos are Majorana particles, their mass terms take the form, where m d = y ε v/ √ 2 is generated by the Yukawa coupling in Eq. (2.1), and m ν and m s are the Majorana masses generated from lepton-number-violating sources beyond the above Lagrangian. Assuming the mass hierarchy m 2 d /m s m ν m s , the active-sterile neutrino mixing angle is approximately θ m d /m s . The two mass eigenvalues are m 1 m ν and m 4 m s , respectively. The heavier mass eigenstate and dark matter candidate is the following linear combination, Alternatively, if the neutrino masses are Dirac (both active and sterile), we must introduce partner fields for ν and ν s for writing down the corresponding Dirac mass terms, where m d = y ε v/ √ 2, and m ν and m s are Dirac masses. We set the (12) element to zero for simplicity. In this case, ν mixes with the ν s field, the resulting dark matter state is ν 4 = ν s cos θ + ν sin θ, where θ m d /m s . Its mass is approximately m 4 m s , assuming that m d m s . In both cases, when θ 1, the mass eigenstate ν 4 has a large sterile neutrino component (i.e., ν 4 ν ( ) s ) and a very small active component. Hereafter, we refer to ν 4 as the sterile neutrino. The orthogonal linear combination, denoted by ν 1 , is a mostly-active neutrino mass eigenstate.
With a mass between keV to MeV scales, ν 4 is a viable dark matter candidate. The angle θ, parametrizing the presence of its active component, is assumed to be small enough so that ν 4 never reaches thermal equilibrium with the Standard Model sector in the early universe. Assuming the universe was born without a population of ν 4 , its relic density could be produced through neutrino oscillation effects. In the early universe, at temperatures before the decoupling of weak interaction, active neutrinos are constantly produced and destroyed by Standard Model weak interactions, as the flavor eigenstate ν. Since the produced ν is a linear combination of mass eigenstates ν 1,4 , neutrino oscillations will occur. In the early universe, the active neutrino encounters a refractive potential due to interaction with the thermal bath, which alters its dispersion relation from the zero temperature case, thereby changing the oscillation probabilities. However, as the neutrino oscillates on a timescale given by its energy and the mass-squared difference, it can undergo weak interactions, that reset only the active component. By this time, the original ν state has already developed a ν s component, which survives such a "measurement", and eventually contributes to the relic density of ν 4 dark matter. The above oscillation process can repeat for many times until weak interaction decouples. This is the crux of the Dodelson-Widrow mechanism [21]. It was demonstrated that a proper choice of the mixing parameter θ can produce the correct dark matter relic density.
The probability for each active neutrino ν to oscillate into a ν s (the latter mostly ends up as dark matter ν 4 ) is dictated by the effective mixing angle participate only in the Standard Model weak interaction, as assumed in the original work by Dodelson and Widrow, the effective thermal potential V T results from the self-energy of the neutrino, as depicted by the diagrams in Fig. 1. The resultant potential is given by [20,38] where G F is Fermi's constant and T is the temperature of the universe. The label "SM" indicates that this expression is only valid in the absence of any new neutrino interaction beyond the Standard Model (BSM). This is the result in the absence of any lepton number asymmetry.
The interaction or collision rate Γ, due to exchange of W and Z bosons, is given by [20,38] We will derive new expressions of these quantities in BSM frameworks in the upcoming sections.
Assuming the initial dark matter abundance to be negligible, the amount dark matter produced is given by the integration of active neutrino production rate times the oscillation probability over the cosmological time. Because the oscillation probability is neutrino energy dependent, we write down the equation that governs the phase-space density evolution for the sterile neutrino with fixed energy E [21,32,38] where H is the Hubble parameter and the parameter z ≡ µ/T is introduced to label the cosmological time. Clearly, the choice of µ does not affect the result, and for convenience, we choose µ ≡ 1 MeV throughout this work. f ν is the Fermi-Dirac distribution for an active neutrino or antineutrino, of the form f ν (x) = 1/(1 + e x ), where x ≡ E/T and it is z-independent. Unlike the proposal of [27], we will not consider the presence of a lepton asymmetry in the universe, thus the chemical potential for active neutrino is set to zero. This is an elegant mechanism, given its simplicity and predictability. Indeed, sterile neutrino dark matter has been explored extensively, from the model building to its various aspects in detection (see e.g., [39,40]), and strong tensions emerge from the indirect detection. In particular, a nonzero ν − ν s mixing θ makes the sterile neutrino dark matter ν 4 unstable.
It could decay into either ν 1 plus a photon, or three ν 1 . The former decay channel could lead to extra X-ray radiation from regions of accumulating dark matter. The non-observation by X-ray telescopes disfavors the mixing angle θ needed for the Dodelson-Widrow mechanism to work in regions where the sterile neutrino is heavier than a few keV. Meanwhile, as a fermionic dark matter, it is found difficult to successfully fill lighter sterile neutrino into the known dwarf galaxies [17,18,41]. As a result, the entire parameter space that produces the correct relic density for sterile neutrino dark matter via the Dodelson-Widrow mechanism is almost closed [22][23][24][25][26].
It is also worth noting that there have been claims of an unresolved ∼ 3.55 keV X-ray line from various nearby galaxy clusters [42,43], which is under thorough scrutinization nowadays [25,[44][45][46]. While the dust has not settled, the relevant takeaway for our work is that if sterile neutrino dark matter decays into this line, the corresponding mixing angle θ is too small to account for its relic density via the Dodelson-Widrow mechanism.
A recent letter [37] has suggested that a new, secret interaction among the active neutrinos is effective for alleviating the above tensions. There, it was postulated that the secret neutrino interaction is mediated by a lepton-number charged scalar and much stronger than the ordinary weak interaction, thus allowing the sterile neutrinos to be produced more efficiently in the early universe. It was shown that with the new interactions, one can explain the relic density without violating all the existing constraints, for dark matter mass between a few keV to MeV scale. This includes the point favored by the observation of the 3.55 keV X-ray line. As the most intriguing aspect of this model, it addresses the relic density of dark matter but with the new physics being introduced in the active neutrino sector. There are predictions in low energy experiments, such as precision measurement of pion, kaon decays, and accelerator neutrino facilities with a near detector (e.g., DUNE). In turn, these probes of secret neutrino interaction can test the fate of dark matter.
In this work, we take this enticing idea further by analyzing models for secret neutrino interactions mediated by mediated by a new vector boson V . Such a new vector boson could naturally arise from U (1) gauge extensions of the Standard Model. In the upcoming sections, we explore three of its incarnations: i) a model with a neutrinophilic V ; ii) gauged U (1) Lµ−Lτ model; and iii) gauged U (1) B−L model. In each case, we confront the parameter space favored by relic density with existing and future experimental constraints.
Strictly speaking, in the last two models, the neutrino interactions are not so secret because other Standard Model particles also see them, and naively they are already tightly constrained. Interestingly, we find that there still exists a viable parameter space to accommodate the correct dark matter relic density in these models. Future experiments will fully probe the remaining parameter space, to either discover or falsify our proposal. In contrast, the neutrinophilic vector boson model with a genuine secret neutrino self interaction is less constrained and calls for new experiments for it to be fully covered.

III. MODEL WITH A NEUTRINOPHILIC VECTOR BOSON
In the first model, we consider a new vector boson V which couples only to the active neutrinos. Because neutrinos exist in an SU (2) L doublet, such a neutrinophilic nature of V is achieved via a higher dimensional operator, added to the Lagrangian in Eq. (2.1), where the cutoff scale Λ αβ characterizes the interaction strength of V with different combinations of lepton flavors. After electroweak symmetry breaking, the vacuum expectation value of the Higgs boson projects out the neutrino field from the lepton doublets. At low energies, the V couplings are neutrinophilic.
where the couplings are λ αβ = v 2 /(2Λ 2 αβ ). Such a new interaction, if strong enough, could keep the neutrinos in thermal equilibrium with themselves longer than the weak interaction, thus facilitating the production rate of sterile neutrino dark matter in the early universe. The task of this section is to quantify this statement and find the favored model parameter space for which the sterile neutrino has a relic density that matches today's observed amount. Meanwhile, the neutrinophilic vector boson could also lead to observational effects in various processes in the laboratories where neutrinos interact. In the following subsections, we derive a list of existing and near-future experimental coverage on the model parameter space which have interesting interplay with the relic density favored region. Up to subsection III G, we consider a representative case where V couplings only to the muon neutrino ν µ , and the sterile neutrino also mixes with ν µ . The flavor dependence of our analysis will be commented afterwards, where it is pointed out that the relic density results remain similar for other choices of flavors. Subsection III H comments on the phenomenological implications of allowing V to couple to different flavors of neutrinos. In subsection III I, we present a possible ultraviolet (UV) completion for the effective operator introduced in Eq. (3.1).

A. The anatomy of sterile neutrino dark matter production
In the original Dodelson-Widrow mechanism, the relic density of sterile neutrino dark matter depends on the neutrino weak interaction rate Γ, and the effective active-sterile neutrino mixing angle that is also controlled by the weak interaction, through the thermal potential V T . In the presence of the new neutrino self interaction mediated by V , both V T and Γ are modified, although the relic density can still be calculated using Eq. (2.8).
The new contribution to the active-neutrino thermal potential is shown in Fig. 2, and takes the following form [47,48]: In the very heavy or very light mediator V limit, the above potential simplifies into It is worth noting that this potential changes sign as T passes m V . These asymptotic expressions are useful for us to infer the parametric dependence in the final relic density. In the numerical integration of the Boltzmann equation, we keep the most general form of the thermal potential, Eq. (3.3). The new interaction mediated by the V boson provides new scattering channels among the neutrinos, including ν µ ν µ → ν µ ν µ and ν µνµ → ν µνµ , as shown by the Feynman diagrams in Fig. 3. The corresponding cross sections are where √ s is the center of mass energy. If the initial-state four-momenta are denoted as where the tiny active neutrino masses are neglected, and ϑ is the relative angle between p and p . These cross sections are used for calculating the thermal averaged reaction rates for a neutrino with fixed energy E, by averaging over the phase space of the other neutrino (or anti-neutrino) it scatters with, is the Fermi-Dirac distribution for active neutrino or antineutrino, σ tot is the sum of the two cross sections in Eq. (3.5), and the relative velocity v rel is equal to 2(1 − cos ϑ) for ultra-relativistic neutrinos.
Here, Γ V also includes the reaction rate of antineutrino, which occurs via the chargeconjugate channels of those reported in Eq. (3.5). At high temperatures when helicity is approximately equivalent to chirality, neutrinos are left-handed and antineutrinos are righthanded. Both of them can contribute to the dark matter relic density via the mixing and oscillation (in a helicity preserving way). Thus, the prefactor of 2 captures both helicity states of the active neutrinos and antineutrinos.
It is useful to show the asymptotic forms of the cross sections and thermal rates in the heavy or light V limits. For m V T , we find and In contrast, for T m V , the typical center-of-mass energy of the scattering is much higher than m V . In this limit, the cross sections in Eq. (3.5) takes the forms These cross section blow up in the m V → 0 limit, similar to the case of Bhabha scattering in QED. Here the scattering range is regularized by a non-zero V mass. Adding the two, the corresponding thermal rate is (3.10) In addition, for T m V , the process ν µνµ → ν µνµ could also occur with an on-shell V exchange. In the early universe, this can also be interpreted as the decay and inverse decay of fully thermalized V particles. In the narrow width approximation, where γ V is the decay rate of V , γ V = |λ µµ | 2 m V /(12π). The corresponding thermal reaction rate is approximately (here we use the Boltzmann distribution instead of Fermi-Dirac in order to obtain an analytic expression) , the on-shell rate has an extra factor (m V /T ) 4 , which is more suppressed at high T . This corresponds to a phase space (collinear) suppression for two energetic neutrinos to scatter at center-of-mass energy equal to m V . However, this factor is not important when the temperature drops to T ∼ m V . In addition, the on-shell rate only involves only two powers of the |λ µµ | coupling, in contrast to the other rates, making it easier to stand out in the small coupling regime.
Numerically, we find that at T ≥ m V , Γ V is approximately given by Γ on-shell These asymptotic rates in the two regimes capture the dominant contributions, and they match each other well at T = m V because the Γ on-shell V term dominates, as long as |λ µµ | < 1.
To summarize, in the presence of the V boson mediated neutrino self interaction, the total thermal potential and reaction rate are given by where V T, SM and Γ SM are the Standard Model contribution, given in Eqs. (2.6) and (2.7), respectively. The new physics contribution V T, V is given by Eq. (3.3), and Γ V given by Eq. (3.6). With the above potential and rate, we numerically solve the Boltzmann equation for the phase space distribution f ν 4 , Eq. (2.8), up to a sufficiently low temperature T f such that the value for f ν 4 saturates. This temperature should be much lower compared to the mediator mass m V but still higher than the dark matter mass m 4 . In practice, we take T f = 100 keV. The dark matter relic density today corresponds to where ρ 0 = 1.05 × 10 −5 h −2 GeV/cm 3 is the critical density, s is the entropy density as a function of temperature of the universe, and s 0 = 2891.2 cm −3 is the entropy density today. The result of the numerical calculation described above is shown in Fig. 4, in the parameter space of λ µµ versus m V , where the orange contours stand for constant values of the sterile neutrino dark matter relic density Ωh 2 . We choose a set of input parameters, m 4 = 7.1 keV and sin 2 (2θ) = 7 × 10 −11 . This point corresponds to the potential dark matter explanation of the 3.55 keV X-ray line. Interestingly, all the curves exhibit an "S-shape", with three regimes having distinct slopes of the curve. This diverse parametric dependence on λ µµ and m V suggests different detailed production processes on which we shall elaborate on below. Along the curve where the correct relic density is reproduced (Ωh 2 = 0.12), we mark three points A, B, C, one from each regime.
The different parametric dependence mentioned above can be understood by examining Eq. (2.8), which in the very small θ limit takes the form,  We plot different factors on the right-hand side of the above equation in Fig. 5. In the figure, each column corresponds to case A, B, or C, respectively. The first row depicts Γ/H as a function of time (labeled by z ≡ 1 MeV/T ) in each case. This ratio dictates the epoch when the active neutrino self interaction drops out of thermal equilibrium (when Γ/H < 1). In all cases the rate Γ features a bump which corresponds to the exchange of an on-shell V in the neutrino self interaction. This contribution is most important when T ∼ m V but becomes Boltzmann suppressed at T m V . At low temperatures, the off-shell V exchange dominates.
The second row of Fig. 5 depicts the effective mixing angle θ ef f [defined in Eq. (2.5)] divided by the vacuum mixing angle θ. The value of θ 2 ef f dictates the probability for each active neutrino in the thermal plasma to oscillation into a sterile neutrino. At very early time, θ ef f /θ is highly suppressed because the oscillation frequency ∆ is much smaller compared to Γ and V T . The production of the sterile neutrino state become more efficient at later time when θ ef f catches up with θ.
The above temperature dependence lead to interesting interplay in the third row of Fig. 5, which depicts the actual dark matter production rate d f ν 4 (x, z)/d ln z and is proportional to the product of Γ/H and θ 2 ef f . Here we focus on the energy E = T (or x = 1) which has the highest population of active neutrinos as the source term. We introduce three useful time scales: i) z 0 : when Γ/H = 1. At z z 0 , the active neutrino self interaction falls out of thermal equilibrium and dark matter production ceases; ii) z 1 : when ∆ Max{|V T |, Γ a }. production; and iii) z 2 ≡ µ/m V . At z z 2 , the vector boson V is too heavy to be produced in the universe.
First of all, if z 0 is smaller than both z 1 and z 2 , the new neutrino interaction via V decouples too early and is irrelevant for dark matter production. We do not consider this case. The three cases A, B, C mentioned above correspond to the following different hierarchies of z 0 , z 1 , z 2 and thus different sterile neutrino dark matter production scenarios.
• Case A corresponds to heavy V and large coupling λ µµ . In this case, the hierarchy of time scales satisfy, z 2 < z 1 < z 0 . As a result, when the temperature is high enough to thermalize V , the mixing angle θ ef f is so suppressed that dark matter produced around this time (z ∼ z 2 ) is negligible. The universe has to wait until a later time (z ∼ z 1 ) for more efficient dark matter production. Dark matter is mainly produced from active neutrino scattering via heavy off-shell V exchange. We find df νs /dlnz ∝ z 9 for z z 1 , and df νs /dlnz ∝ z −3 for z z 1 . They lead to the triangle peak in the third row of case A in Fig. 5 and the following parametric dependence in the final relic density, (3.16) Case A also represents the standard Dodelson-Widrow mechanism with |λ µµ | and m V replaced by the SU (2) L gauge coupling g 2 and the W -boson mass M W respectively.
• Case B has the hierarchy of time scales: z 2 < z 0 < z 1 . The production of dark matter is restricted to the time window z 2 z z 0 , where the the effective mixing angle θ ef f is still suppressed compared to the vacuum angle θ, but the universe remains hot enough compared to V mass. Therefore, dark matter can be produced through active neutrino scattering via on-shell V exchange, or equivalently, the decay (and inverse decay) of V particles in the plasma. In this case we find df νs /dlnz ∝ z 7 and the key parametric dependence in the final relic density is, This is a new production regime beyond the Dodelson-Widrow mechanism.
• Case C corresponds to a very light V and tiny coupling λ µµ . In this case the time scales satisfy, z 1 < z 2 < z 0 . As a result, when the effective mixing angle θ ef f relaxes to the vacuum angle θ, the universe is still hot enough compared to V mass. Thus, dark matter is mainly produced from active neutrino scattering via on-shell V exchange. In the dominant production window, z 1 z z 2 , we find df νs /dlnz ∝ z 3 . This implies the z integral is dominated by z z 2 , where T ∼ m V . The key parametric dependence in the final relic density is, Like case B, this is another new production regime beyond Dodelson-Widrow.
In practice, we find that cases A, B, C cover all the possible production scenarios by varying the values of m ν 4 and θ. The S-shape of the relic density contours is an interpolation through the three cases. In case A, the dominant production temperature is T m V , whereas in cases B and C, the dominant production temperature is T ∼ m V . Thanks to this, thermal corrections to the V mass has little impact on the result.
Mathematically, another hierarchy z 1 < z 0 < z 2 is also possible. However, for the phenomenologically allowed values of λ µµ and m V considered in this work, there is no room for z 0 < z 2 to occur, i.e., the V boson always becomes heavy before the new interaction decouples.
For varying values of dark matter mass m 4 and mixing angle θ, the S-shape relic density curve sweeps across the λ µµ versus m V plane. In Fig. 6, we plot two choices of parameters: • Gray curve: m 4 = 7.1 keV, sin 2 (2θ) = 7 × 10 −11 .
For heavier sterile neutrino dark matter (black curve), we choose a correspondingly smaller mixing angle in order to satisfy the X-ray constraints. We find that heavier dark matter generically requires larger |λ µµ | and/or smaller m V for the relic density, making it more constrained experimentally. On the other hand, due to the dwarf galaxy and Lyman-α constraints, there is limited room to make the dark matter ν 4 much lighter than 7.1 keV. As a result, the gray curve in Fig. 6 roughly indicates the lower margin of parameter space where ν 4 comprises the total dark matter relic density in the universe. I.e., all the target parameter space lies above the gray curve.

B. Constraint from the Higgs boson invisible decay
In this and the next few subsections, we present experimental constraints on the neutrinophilic V boson. To obtain the correct dark matter relic density, the mass of V is required to be smaller than the weak scale. It mainly decays into neutrinos and thus appears invisible after production in laboratories. This leads to a number of constraints by making precision measurements of processes that involve an active neutrino.
The first process we consider is the Higgs boson decay. The operator introduced in Eq. (3.1) opens up a new Higgs invisible decay channel, h → ν µνµ V . The corresponding partial decay rate is where h V ≡ m 2 V /m 2 h . We see that in the small m V limit, this partial width is enhanced by a factor of 1/m 2 V , corresponding to the Higgs decaying into the longitudinal component of the V boson. The operator responsible for Higgs decay, λ(h/v)νγ µ νV µ , does not admit any conserved U (1) gauge symmetry for V µ to be promoted as the corresponding gauge boson. For example, the U (1) of neutrino number is explicitly broken by the presence of the Higgs field.
To set a limit using this channel, we note that the invisible decay branching ratio of the Higgs boson in this model is calculated as, Γ h→νµνµV /Γ total , where the Higgs total width Γ total is the sum of the Standard Model Higgs width (∼ 4 MeV) and Γ h→νµνµV . The present LHC upper bound on Higgs invisible branching ratio, Br(H → invisible) < 24% [49], leads to a constraint in the λ µµ versus m V parameter space, corresponding to the purple shaded region in Fig. 6. The future running of high-luminosity (HL) LHC is expected to further improve the above limit to Br(H → invisible) < 2.5% [50], which corresponds to the purple dashed curve in Fig. 6. We find this is the leading constraint on the model parameter space for m V > 5 MeV. In particular, for the case m 4 = 50 keV and sin 2 (2θ) = 10 −15 , the present Higgs invisible limit still allows a region with m φ ∼ 100 MeV, but that will be covered by the HL-LHC. It is also worth pointing out that the Higgs invisible decay constraint equally applies to V coupling to other neutrino flavors.

C. Constraint from W boson decay width
The next channel we examine is the W − → µ −ν µ V where V is radiated from the final stateν µ via its neutrinophilic coupling. Ignoring the muon mass, this decay rate is where w V ≡ m 2 V /M 2 W . As with the Higgs decay, we see that in the small m V limit, this decay width is enhanced by 1/m 2 V . To derive a conservative constraint, we simply require this decay rate to be smaller than the uncertainties in the W total width measurement, ∼ 42 MeV [49]. The resulting limit is found to be much weaker than the one from Higgs invisible decay, and we do not show it in Fig. 6. Including the muon mass in this calculation would drive the constraint to be slightly weaker than this calculation.
We note that the W decay limit might be further improved by studying the kinematic edge observables of W → µ + MET decay at the LHC. In the Standard Model, when the W -boson is singly produced, the final charged lepton transverse momentum distribution features a Jacobian peak [51]. This feature is absent when V is present in the decay product. We leave a careful study of this for a future work.

D. Constraint from Z boson invisible width
We also examine the Z → ν µνµ V decay where V is radiated from either ν µ orν µ in the the final state. Because V decays into neutrino and appears invisible at colliders, this decay channel contributes to additional invisible Z-boson decay width. This decay rate is (3.21) Unlike the Higgs and W cases, the new Z boson decay width does not feature the 1/m 2 V enhancement at small m V . This could be understood as follows. In this decay process, the vector boson V couples toν µ γ µ ν µ which can be defined as the neutrino current. The U (1) symmetry for ν µ neutrino number associated with this current is preserved by the Standard Model Zν µνµ coupling. This symmetry is only broken by the other Standard Model couplings (e.g., the W µν coupling considered earlier) which do not take part in this decay rate at leading order. Due to the lack of the 1/m 2 V factor, numerically, we also find the constraint from Z invisible decay is much weaker than that from Higgs invisible decay derived above. It is shown as the red shaded region in Fig. 6.

E. Constraints from exotic meson decays
For the vector boson V µ lighter than GeV scale, the V νν coupling could also lead to new decay channels of charged mesons will exist, m + → µ + ν µ V , where m = π, K, B, etc. Because V will invisibly decay back to neutrinos, such processes would contribute to the branching ratios of the m + → µ + ν µ νν and get constrained. The decay width m + → + νV is where G F is the Fermi constant, F m is the decay constant of m, and |V qq | is the relevant CKM matrix element for this decay. In the small m V limit, this integral may be performed analytically, In practice, we have considered the charged kaon and pion decays, and apply the experimental constraints, Br(K + → µ + ν µ νν) < 2.4 × 10 −6 and Br(π + → µ + ν µ νν) < 5 × 10 −6 [49]. The resulting limit is shown by the blue shaded region in Fig. 6. It is weaker than that from Higgs invisible decay.

F. Constraints from BBN and Supernova 1987A
If the new vector boson V is thermalized in the early universe and light enough to remain relativistic by the time of big bang nucleosynthesis, it will contribute to ∆N eff. and affect the primordial element abundances. A conservative constraint rules out m V 5 MeV [13]. For a detailed computation of the effects of non-standard neutrino self-interactions on big bang nucleosynthesis, see [52].
If the vector V is lighter than ∼ 100 MeV, it can be produced from neutrino scatterings in the explosion of Supernova 1987A. It could carry away significant fraction of energy and modify the observed time scale of neutrino emission. Ref. [53] has estimated this effect in the gauged U (1) Lµ−Lτ model (see also Section IV). We rescale and apply their bound by restricting V to couple to only one neutrino flavor.

G. Mono-neutrino probes at DUNE
Neutrinophilic particles like V may be emitted in beam neutrino experiments, if kinematically allowed. The processes ν µ n → V µ − p + occurs via initial state radiation. In contrast to standard neutrino quasi-elastic scattering, ν µ n → µ − p + , the radiation of V carries away both energy and transverse momentum. As a result, the final state µ − carries lower energy than expected in the quasi-elastic scattering case. The resulting µ − p + system appears to have a non-zero p T with respect to the beam direction. These are the mono-neutrino signals introduced in [9,11] in a different context. In contrast to the lepton-number charged scalar radiation studied carefully in [11], here the radiation of V does not carry away leptonnumber. The final state muon in the above signal process has the same electric charge as the quasi-elastic background, thus it will not benefit form the muon charge identification capability of the neutrino detector. In contrast, the emission cross section for a light V features a E 2 /m 2 V longitudinal enhancement, making the signal stronger at low m V . We derive an expected 95% CL constraint assuming five years of DUNE [54] data collection in the neutrino mode, as shown by the dashed green line in Fig. 6. An additional five years of DUNE data collection in antineutrino mode does not improve this limit significantly, as the background processes are more difficult to resolve when an antineutrino scatters. Our obtained sensitivity is comparable to the current Higgs invisible limit but will be surpassed by the high-luminosity LHC search.

H. On V coupling to other neutrino flavors
Throughout Section III we have focused on λ µµ , the V coupling to muon-flavored neutrinos. In principle, the couplings λ αβ comprise a 3 × 3 matrix including couplings to other flavors and even flavor-violating couplings. Since all of the early-universe reactions that lead to the production of ν 4 are driven by the self interaction of neutrinos, the relic density results shown in Figs. 4 and 6 are independent of the choice of flavor, as long as one only element of λ is turned on each time. Experimentally, the strongest constraint on λ µµ comes from the Higgs boson invisible decay, as shown in Fig. 6. It equally applies to all other λ αβ couplings.
On the other hand, we note that if the V coupling is flavor off diagonal, there will be strong constraints from loop induced µ → eV , τ → eV , τ → µV decays. Given order-ofmagnitude estimates of this loop process, we find that these constraints would be stronger than all of those discussed above and would rule out the desired parameter space for relic sterile neutrino dark matter.

I. A possible UV completion
In this subsection, we present a UV completion for the effective operator introduced in Eq. (3.1). We extend the Standard Model with a pair of chiral fermions N L , N R and a complex scalar φ. All are Standard Model gauge singlets. The N L and φ fields are oppositely charged under a new U (1) gauge symmetry whereas N R is neutral. With such a particle content the U (1) still possess a gauge anomaly that could be canceled by resorting to including additional heavy fermions without direct couplings to the lighter fields. We will keep those heavy particles implicit.
With the N L , N R and φ fields, we could write down the following renormalizable Lagrangian governing their interactions, (3.24) where L, H are the Standard Model lepton and Higgs doublets, and D µ = ∂ µ ± ig V µ . We assume φ has a potential V (φ) such that φ gets a vacuum expectation value v φ and breaks the U (1) symmetry. This gives a mass to the new gauge boson V , and also give a Dirac mass for the N L , N R fermions.
Next, we go to energy scales below M N and integrate out N L , N R . From the square bracket in Eq. (1) the equation of motion ofN R yields (neglecting its kinetic term) As a result, the gauge interaction of N L with the V in the first line of Eq. (1) leads to the effective interaction This is same as the higher dimensional operator introduced in Eq. (3.1), with Λ 2 = M 2 N /(g y 2 ). The presence of the new Dirac fermion N in this UV completion could be subject to further experimental constraints.
The second model we explore is an extension of the Standard Model with gauged U (1) Lµ−Lτ symmetry. It is one of the simplest anomaly-free U (1) models known to provide an explanation to the discrepancy between the Standard Model prediction of muon g −2 and the experimental measurement [55]. The phenomenology of this model has been explored extensively [56][57][58][59][60][61][62]. Here, we point out that the model also contains a target parameter space for sterile neutrino dark matter relic density.
The couplings of the U (1) Lµ−Lτ gauge boson V with Standard Model fermions are where g µτ is the gauge coupling. Like the previous section, we assume that the sterile neutrino dark matter ν 4 has a small ν µ component, characterized by a mixing angle θ. To generate such a mixing, one cannot rely on the Yukawa coupling in Eq. (2.1). Instead, the following higher-dimensional operator must be introduced, where φ is a complex scalar field that Higgses the U (1) Lµ−Lτ symmetry and carries an opposite charge to that of L = (ν µ , µ) T . This is because the ν s field ( ν 4 ) must be a singlet under U (1) Lµ−Lτ in order to prevent dark matter from being fully thermalized in the early universe and violating the initial condition for the production mechanism discussed in this work.
The relic density calculation in this model proceeds in similar fashion as the neutrinophilic one in Sec. III. The contribution to the thermal potential V T is the same as Fig. 2 with the gauge boson V and ν running in the loop. It is calculated using Eq. (3.3) with λ µµ replaced by the gauge coupling g µτ . In contrast, the thermal reaction rate of neutrino receives contributions from the scattering with ν µ ,ν µ , ν τ ,ν τ , µ ± , τ ± particles, via the new L µ − L τ gauge interaction.
To simplify the thermal averaging in the numerical calculation, we make the instantaneous decoupling approximation. For temperatures above the mass of a fermion participating in a scattering process, we calculate the thermal rates assuming the fermion to be massless. In contrast, when the temperature drops below any of the fermion mass, we assume the process vanishes immediately. Because sterile neutrino dark matter is dominantly produced at temperatures below a few hundred MeV (see, e.g., Fig. 5 in subsection III A), the τ ± leptons are already decoupled from the universe, and the presence of muons is marginal. This is a good approximation to simplify the numerical workload.
Under this approximation, we report the cross sections and thermal rates relevant to the relic density calculation in this model: .

(4.3)
From the experience gained in the previous section (see Fig. 5), the most relevant neutrino reactions for dark matter production occur at T m V . We present the thermal rate results for the heavy V and on-shell V cases. In the heavy-V limit (m V T ), T < m µ .

(4.4)
For on-shell V contribution (m V T ), the corresponding rate is where w ≡ m 2 V /(4ET ). γ V is the decay width of the L µ − L τ gauge boson, where r α = m 2 α /m 2 V and Θ is the Heaviside theta function. The above thermal potential and reaction rates are added to their counterparts in the Standard Model and then inserted to the Boltzmann equation (2.8). The dark matter relic density favored parameter space is shown in the g µτ versus m V plane in Fig. 7. Similar to Fig. 6, the gray and black curves corresponds to two choices of parameters: m 4 = 7.1 keV, sin 2 (2θ) = 7 × 10 −11 and m 4 = 50 keV, sin 2 (2θ) = 10 −15 , respectively. The relic density curves exhibit similar S-shapes, corresponding to the three production regimes discussed in subsection III A.
Sterile neutrino dark matter production in this model has been explored in [63] in the light V limit, which corresponds to the case C defined in subsection III A. In this regime, our result is consistent with theirs. Fig. 7 also shows the experimental constraints, mostly adapted from Ref. [64], where the colored shaded regions are already excluded and the regions enclosed by colored dashed curves will be covered by future experiments, correspondingly labeled. We also include the cosmological/astrophysical constraints from BBN and supernova 1987A [53], which mainly covers the region with m V below a few MeV scale. Compared to the neutrinophilic case, the U (1) Lµ−Lτ model has a stronger interplay between the dark matter relic density favored regions and experimental reaches. The present limits almost rule out the entire parameter space for the m 4 = 50 keV case. The reach of future experiments including SHiP [65,66] and NA64µ [67,68], if they occur, will be able to cover nearly the entire region above the m 4 = 7.1 keV curve, that is, roughly the whole viable parameter space for dark matter relic density. The DUNE experiment, with a liquid argon near detector, can also probe the parameter space overlapping with the (g−2) µ favored region, by measuring the neutrino trident production [69,70]. DUNE, if equipped with a gaseous multi-purpose near detector, is also able to probe the displaced decay of V into e + e − via a loop-generated kinetic mixing [71]. Finally, a proposed muon-missing-momentum experiment [72] can probe a parameter space similar to NA64µ, for m V < 2m µ . The third model we consider is an extension of the Standard Model with U (1) B−L gauge symmetry [73]. This is another popular framework where many phenomenological works have been done. The interactions of the B − L gauge boson V take the form universal for all quark and lepton flavors. The new gauge interaction, under which active neutrinos are charged, could potentially facilitate the production of sterile neutrino dark matter. Like the U (1) Lµ−Lτ model, the active-sterile neutrino mixing has to be generated by a higher dimensional operator similar to Eq. (4.2). The relic density calculation proceeds the same as before. The thermal potential is calculated using Eq. (3.3) with λ µµ replaced by the gauge coupling g BL . To calculate the thermal reaction rates of neutrinos, we make the same instantaneous decoupling approximation as described in the previous section. The thermal reaction rate in the heavy mediator V limit takes the form T < m e .

(5.2)
For on-shell V contribution, the rate is 2m e ≤ T < 2m µ 3, T < 2m e (5.3) where w ≡ m 2 V /(4ET ), and γ V is the decay width of the B − L gauge boson [74]. The above rates are calculated based on the interactions of neutrinos with themselves and the charged leptons. For temperatures above GeV scale, the processes where neutrino scatters with/annihilates into baryons also contribute. 1 However, they are unimportant because the sterile neutrino dark matter is dominantly produced at temperatures well below a few hundred MeV. See discussion in subsection III A and Fig. 5. At higher temperatures, the effective active-sterile neutrino mixing angle is suppressed.
The result of our numerical calculation is presented in Fig. 8, where the gray and black curves corresponds to the two same choices of parameters as before: m 4 = 7.1 keV, sin 2 (2θ) = 7 × 10 −11 and m 4 = 50 keV, sin 2 (2θ) = 10 −15 . The relic density curves exhibit similar Sshape, with three production regimes having distinctive parametric dependence, as discussed in subsection III A. In the same figure, we also show the experimental constraints on the U (1) B−L model parameter space Ref. [64,75], where the colored shaded regions are already excluded. Because the gauge boson V couples to the electron, constraints are the strongest among the three models we have studied. We find out that the entire curve in the m 4 = 50 keV case has already been firmly excluded. The ongoing Belle-II [76], the upcoming FASER [77] at LHC, and the proposed LDMX experiment [72] will cover the region enclosed by the blue and purple dashed curves. Combined, they will test all the relic density favored parameter space. 2 To summarize, the available parameter space for sterile neutrino dark matter production is already quite small in the U (1) B−L model. It can be fully covered by experiments in the near future.

VI. CONCLUSION
In this work we considered a sterile neutrino as the dark matter candidate, and studied the origin of its relic abundance. The Dodelson-Widrow mechanism fails to work, given the severe constraints from X-ray indirect detection and small scale structure observations. We explored simple extensions of the Standard Model where a new vector boson V mediates new interactions for active neutrinos that facilitates efficient production of sterile neutrino dark matter. We calculated the dark matter production via neutrino scattering and oscillation in the early universe and identify three scenarios (A, B, C as defined section III A). Depending whether V plays the role of a heavy or light mediator, it leads to different parametric dependence in the final dark matter relic density. The viable mass of V lies between MeV to GeV and the corresponding neutrino coupling ranges from 10 −6 to 10 −2 . For comparison, we also carried out similar calculations in the gauged U (1) Lµ−Lτ and U (1) B−L models. The main results are presented in Figs. 6, 7 and 8.
The parameter space for viable dark matter relic density found in this work is a well motivated target for the next experimental probes. Among the three models studied, the U (1) B−L model is the most strongly constrained. Its parameter space for sterile neutrino dark matter is already narrow and will be fully covered by future high intensity experiments, including BELLE-II, FASER, and LDMX. The U (1) Lµ−Lτ case is relatively less constrained but could also be covered by the proposed experiments, including SHiP, NA64-µ, M 3 , and DUNE. In contrast, the model with a neutrinophilic vector boson is the least constrained.
The strongest present limit comes from the Higgs boson invisible decay search, which will be improved by HL-LHC. A portion of its parameter space has not been targeted by any experiment, to our knowledge.
On the theory side, the same dark matter production mechanisms could also work in, e.g., gauged U (1) Le−Lµ or U (1) Le−Lτ models, or even gauged lepton-number that has other motivations [78,79]. The presence of electron coupling for the vector boson V in these models implies strong experimental probes similar to U (1) B−L .