The Dispirited Case of Gauged $U(1)_{B-L}$ Dark Matter

We explore the constraints and phenomenology of possibly the simplest scenario that could account at the same time for the active neutrino masses and the dark matter in the Universe within a gauged $U(1)_{B-L}$ symmetry, namely right-handed neutrino dark matter. We find that null searches from lepton and hadron colliders require dark matter with a mass below 900 GeV to annihilate through a resonance. Additionally, the very strong constraints from high-energy dilepton searches fully exclude the model for $ 150 \, \text{GeV}<m_{Z'}<3 \, \text{TeV}$. We further explore the phenomenology in the high mass region (i.e. masses $\gtrsim \mathcal{O}(1) \, \text{TeV}$) and highlight theoretical arguments, related to the appearance of a Landau pole or an instability of the scalar potential, disfavoring large portions of this parameter space. Collectively, these considerations illustrate that a minimal extension of the Standard Model via a local $U(1)_{B-L}$ symmetry with a viable thermal dark matter candidate is difficult to achieve without fine-tuning. We conclude by discussing possible extensions of the model that relieve tension with collider constraints by reducing the gauge coupling required to produce the correct relic abundance.


Introduction
The nature of the dark matter in the Universe and the origin of the tiny neutrino masses are two unsolved fundamental questions in particle physics. Small neutrino masses are naturally generated by extending the Standard Model (SM) with massive Majorana neutrinos, singlets under the SM gauge group, via the well-known seesaw mechanism. While the scale of such sterile neutrino masses is largely unconstrained (it can be the GUT scale for Yukawa couplings Y of O(1), a few GeV for Y ∼ 10 −6 , etc.), the mixing among active and sterile neutrinos is generically very small, making it difficult to detect sterile neutrinos at current colliders, even if kinematically accessible. Interestingly, this simple extension of the SM can also provide a dark matter candidate if the lightest of the sterile neutrinos has a keV-scale mass [1][2][3][4]. Although a very appealing possibility, this extension involving only the addition of right-handed neutrinos is disfavored by the absence of an X-ray signal arising from the sterile neutrino decay N → νγ [5][6][7] -a process necessarily induced by the required mixing with the active neutrinos [8].
A slightly enlarged scenario that also links sterile neutrinos and dark matter consists of charging both under U (1) B−L and coupling them to a SM singlet scalar, also charged under this global U (1) B−L symmetry, which generates Majorana masses for all the SM singlet fermions. The spontaneous breaking of the U (1) B−L leads to the appearance of a Goldstone boson, the Majoron [9,10]. In this scenario there are two potential dark matter candidates: 1) the Majoron, assuming it acquires a mass e.g. via non-perturbative gravitational effects which break the global symmetry [11,12] (see [13] for a recent review on the subject, and for a fresh perspective on global symmetries see [14]), and 2) one (or several) of the massive SM singlet fermions. As before, the simplest solution involves taking the dark matter to be a keV sterile neutrino. However, unless the Majoron is more massive than the sterile neutrino, the sterile neutrino is unstable 1 . Thanks to the interactions with the scalar singlet, a heavier sterile neutrino stabilized by an additional symmetry (e.g. a Z 2 symmetry which forbids a Yukawa coupling to SM leptons) can also play the role of the dark matter. This possibility has been extensively studied by two of the authors in [16], where it was shown that the observed dark matter relic abundance can be obtained by the annihilation freeze-out of a Majorana fermion with mass in the range of a WIMP (∼1 GeV -2 TeV). Within this scenario, the common origin of the dark matter mass and the masses of the other sterile neutrinos (i.e. those that mix with the active neutrinos) naturally links the unknown seesaw scale to the scale of a thermal WIMP.
Despite containing an electroweak scale dark matter candidate, these scenarios are often quite elusive, a consequence of the fact that the production mechanism in terrestrial experiments (other than the possible mixing of active and sterile neutrinos) arises from the mixing between the SM Higgs and the CP-even component of the singlet scalar. Given current limits on the aforementioned mixing from e.g. precision Higgs measurements at the LHC (see e.g. [17][18][19][20][21]), significant production of the sterile neutrinos at colliders is unattainable. Furthermore, the presence of the massless (or light) Majoron allows for invisible decays of both the sterile neutrinos and the CP-even scalar, implying that even should these particles be produced, their dominant decay modes make this model extremely hard to test.
As a consequence, considerably more attention has been given to models in which the SM is extended by a U (1) B−L gauge symmetry, which is interesting from both a phenomenological perspective (in the sense that it is easier to probe experimentally) and a theoretical point of view, since canceling the anomaly arising from thew new gauge symmetry naturally requires the existence of three righthanded neutrinos. The absence of the Majoron makes a keV sterile neutrino again cosmologically stable, and its new gauge interaction allows for the production of the observed dark matter relic abundance by the freeze-in mechanism, without conflicting with the non-observation of an X-ray signal from its decay [22,23]. Alternatively, one can consider a more massive sterile neutrino dark matter candidate, once again stabilized by an additional symmetry, whose relic abundance is now determined via the more conventional freeze-out mechanism. Historically, studies focused on this scenario have only studied a small part of the viable parameter space, analyzing just one particular mediator and/or mass range (e.g. the Z portal has been explored in [22,[24][25][26], the Higgs portal in [27,28], and the classically conformal minimal model in [29]).
We find that a thorough analysis of the minimal gauged U (1) B−L sterile neutrino 2 dark matter model is missing from the literature, and thus we present here a comprehensive analysis of the viability of a sterile neutrino WIMP dark matter candidate within the context of a gauged U (1) B−L extension of the SM. We focus here explicitly on the case where the dark matter fermion has a lepton charge L = 1, but comment in the text on alternative possibilities [30][31][32][33]. Moreover, we also study generic phenomenological features appearing in extensions of the minimal scenario via the addition of singlet scalars.
The paper is organized as follows. In Section 2 we briefly review the minimal U (1) B−L dark matter framework considered here. After summarizing in Section 3 all current constraints on the U (1) B−L gauge boson Z , we focus in Section 4 on the dark matter phenomenology and identify, via a comprehensive scan, the viable parameter space within the minimal scenario for which a dark matter candidate can still be accommodated. In Section 5 we summarize the results of the scan and discuss various theoretical concerns regarding the allowed parameter space, e.g. the potential appearance of a Landau pole or induced instability of the scalar sector occurring at low energies. Section 6 is devoted to the analysis of generic extensions of the minimal framework. We conclude in Section 7.

Model Outline
In this work we consider possibly the simplest gauged U (1) B−L extension of the Standard Model with sterile neutrinos and a dark matter candidate. We construct the model based on the minimal particle content required for U (1) B−L to be non-anomalous. Namely, we introduce three right handed leptons with charge L = 1. Two (N 1 , N 2 ) of these three states will be the ones responsible for the masses and oscillatory pattern of the active neutrinos via the type-I seesaw mechanism. The other (χ) will be our dark matter particle, which is stabilized via a Z 2 symmetry. In order to give masses to these fermions we also introduce a scalar field (φ) with charge L = −2, which upon symmetry breaking will render Majorana masses for these right-handed leptons. The relevant part of the Lagrangian then reads: Additionally, it should be understood that the inclusion of the new gauge symmetry necessitates a modification to the covariant derivatives which subsequently produce interactions with the new Z . Dirac dark matter candidates within the context of this (or a similar) model have been studied in e.g. [25,[34][35][36][37]. We choose here to focus on a Majorana candidate because this maximally suppresses constraints from direct dark matter detection; collider constraints, on the other hand are largely insensitive to the difference between Dirac and Majorana dark matter, thus making our analysis here maximally conservative. We also note that one could consider a slight non-minimal extension of this model, with a scalar dark matter candidate [25,38,39]; however, in addition to being non-minimal, this model suffers from strong direct detection constraints similar to that of the Dirac dark matter case.

Constraints on the Z
In the following subsections we enumerate the constraints that arise from having a Z with nonnegligible couplings to both baryons and leptons. We emphasize that the collider constraints described below implicitly assume that the invisible branching fraction of the Z is negligible, however we have explicitly checked that for any choice of m χ and m N1,N2 this invisible channel does not exceed ∼ 25%, and thus the constraints are never significantly altered. Nevertheless, we verify in the following sections the extent to which reducing collider bounds, such that invisible decays are saturated, ameliorates existing tension on the model. A summary of the collider constraints on the B − L gauge coupling as a function of the Z mass is presented in Fig. 1.

LHC
A number of searches from the LHC produce strong constraints on Z models, in particular for models in which the Z couples to charged leptons. We briefly review the most stringent constraints of our model below.
In the high mass region, i.e. 150 GeV m Z 4 TeV, bounds from ATLAS on the dilepton final state using 36fb −1 of data at 13 TeV lead to strong limits on the B − L coupling (shown in blue in Fig. 1) [40] (see [41] for similar CMS constraints on this channel). Constraints on the B − L model are produced by the collaboration, and thus presented bounds are inferred using Fig. 5 of [40].
Constraints on lower Z masses (30 GeV m Z 300 GeV) can be derived by analyzing the modification that would be produced to the Drell-Yan (DY) differential cross section. The triple differential cross section for DY production has been measured to remarkably high accuracy; here, we focus on the modifications to the invariant mass spectrum with muon final states (in the central rapidity channel), using measurements from ATLAS at √ s = 8 TeV with 20.2fb −1 [42]. Note that a similar analysis was performed in the context of a different model in [43]. Explicitly, we use MadGraph5 aMC@NLO [44,45] to compute the modified DY differential cross section for invariant dimuon masses 46 GeV < m µµ < 200 GeV. Cuts on the pseudo-rapidity and transverse momentum of the Z are imposed using MadAnalysis5 [46], with 10 equally spaced bins in the dimuon spectrum and assuming uncorrelated errors between bins, following the analysis of the collaboration [42]. The derived constraint is shown in Fig. 1 in red. At yet lower masses (10.6 GeV < m Z < 70 GeV), recent constraints from LHCb on the identical channel, i.e. pp → µ + µ − , provide the dominant constraint on the Z [47]. Re-deriving such constraints using MadGraph5 aMC@NLO is non-trivial, as the analysis performed by the collaboration is complicated. However, since the results are interpreted in terms of a dark photon model, such bounds can be effectively translated into the B − L model, as is done e.g. in [48]. This can be accomplished by noting that the production cross section from uū and cc are suppressed by a factor of 4 relative to that of the dark photon model, and production cross sections from the dd, ss, and bb channels are identical in both models [48]. Since the branching fractions to muons in both models are approximately equal within the mass range of interest, conservative limits can be set on the B − L coupling by diving the cross section by a factor of 4 (note that this is conservative because only some, not all, of the production channels are actually suppressed). The bounds shown in cyan in Fig. 1 accurately reproduce those derived in [48].
The above list certainly does not account for the entirety of the constraints on this model; less stringent constraints from other channels are not shown here to maintain clarity (see [49] for some other current and future constraints).
It is also worth noting that in generic models with multiple U (1) gauge symmetries, mixing between the gauge bosons may occur. Such a mixing typically leads to even stronger constraints on the B − L gauge coupling, see e.g. [50] for derived limits on the mixing in a U (1) B−L model and [51] for general light Z constraints from kinetic mixing. Throughout this paper we have implicitly taken this mixing to be zero, however in principle the running of the couplings may induce non-zero mixing at scales probed by colliders. We have verified that for our purposes, the evolution of this coupling for the parameter space of interest is small over the relevant scales.
Before continuing we would like to briefly comment on the potential reach of the future highluminosity LHC (HL-LHC) in constraining the existence of such a Z . From the discussion above, it should be clear that dilepton searches provide the most stringent constraints; this is simply a consequence of the fact that the branching ratio to leptons is large and electron/muon searches are clean. In order to understand the extent to which the high mass region of parameter space may be probed (within the authors' lifetimes), we preform the following simple first-order analysis: using MadGraph5 aMC@NLO, we determine the distribution of dimuon events as a function of invariant dimuon mass arising from the Z , taking the B − L coupling to be at the perturbative limit (i.e. g B−L = √ 4π), setting the center-of-mass energy √ s = 14 TeV, and assuming a luminosity of L = 3 ab −1 . We find that roughly O(5) events can be observed for m Z 7 TeV; thus we will consider m Z 7 TeV to be untestable in near-future colliders.

BaBar
The BaBar collaboration has published strong constraints on both visible decays [52], with e and µ final states, and invisible decays [53] of a dark photon in the mass range ∼ 2m µ − 10 GeV. Since the Z in the U (1) B−L model decays predominantly to visible particles, we use this search to set bounds. The published limits on the visible decay modes can be straightforwardly translated into constraints on the B − L gauge coupling by equating the dark photon mixing with the gauge coupling, and multiplying the constraint by the relative branching fractions in each model (for each decay channel). Since the branching fractions of the dark photon decay to charged leptons are approximately equal to those of the Z , the limits are nearly identical in strength. Furthermore, limits from BaBar on invisible decays are of comparable strength, implying these bounds cannot differ by more than a factor of 2.

LEP
The bound from LEP originates from the agreement of the measured e + e − → + − cross sections at various energies with the SM prediction [54]. In particular, for the B − L gauge boson, it was shown in [55,56] that LEP observations require (m Z /TeV) > 6 g B−L .

Dark matter phenomenology
In this section we comment on the relevant phenomenology of our Majorana dark matter candidate. In particular we discuss the computation of the relic density (Section 4.1), direct detection constraints (Section 4.2), the algorithm used to determine the viable parameter space (Section 4.3), and relevant annihilation phenomenology (Section 4.4).

Dark Matter Relic Density
The dark matter relic density has been computed by implementing the model in FeynRules [57] and using micrOMEGAs 4.3.5 [58]. Since the parameter space is quite large, it is not possible to show the diverse range of possible couplings that give rise to the correct dark matter abundance. For illustrative purposes, we show in Figure 2 the value of the gauge coupling providing the correct relic abundance as a function of the Z mass for various choices of the dark matter, scalar and right-handed neutrino masses, plotted against the constraints discussed in Section 3. It is worth pointing out that the values of the gauge coupling producing the correct relic abundance are fairly independent of the masses of the CP-even scalar and the sterile neutrinos, provided they are not on-resonance with χ. Thus, the examples provided in Fig. 2 are broadly representative of the model.

Direct Detection
The direct detection cross section for the dark matter candidate contains a spin-independent component that is suppressed by the square of the perpendicular component of the velocity of the recoiling where q is the momentum transfer and µ the dark matter-nuclei reduced mass) and a spin-dependent component suppressed by the momentum transfer squared (see e.g. [59] for dark matter-nuclei differential cross section); thus one naturally expects direct detection constraints to be significantly reduced relative to the constraints from colliders. To verify this, we have computed the constraints from XENON1T and PandaX-II explicitly [60, 61] 3 . The cross section for this interaction has been calculated using the tools provided in [63]. The bounds have been computed for each experiment by taking the efficiency function and exposure from [60] and [61], and by applying Poisson statistics. This analysis reproduces to an accurate degree the published limits on the spin-independent interaction.
The constraints on the coupling as a function of dark matter and Z mass from the Xenon1T experiment are shown in Fig. 3. We note that PandaX-II produces nearly identical limits. As expected, they are reduced with respect to collider searches in the parameter space of interest.
In order to illustrate the sensitivity of future direct detection experiments, we also derived a 90% CL sensitivity limit for the DARWIN experiment [64], assuming no events are observed (see right panel of Fig. 3). This sensitivity limit is constructed considering perfect detection efficiency for nuclear recoils between 5 keV and 40 keV, a 40 ton fiducial volume and a 5 year runtime.
In our analysis we take the scalar-Higgs mixing to be zero; note that this is valid from a phenomenological perspective as this mixing is significantly constrained from direct detection searches and LHC Higgs measurements [17][18][19][20][21]. We have included for completeness the direct detection cross section arising from Higgs-mixing in the Appendix A.

Full Parameter Space Scan
Here, we describe the sampling method that allows for a comprehensive understanding of the viable parameter space without relying on simplifying assumptions that may bias the results.
This parameter space scan begins by discretely sampling the Z mass at 200 logarithmically spaced points between 1 GeV and 30 TeV. Note that we have truncated our analysis at 1 GeV as strong constraints exist for MeV scale Z s that couple to both baryons and leptons (see e.g. Ref. [48]) 4 . For each fixed value of m Z , randomly generated points are sampled in the 3-dimensional parameter space defined by (m χ , m N1 = m N2 , m ρ ), with all masses restricted to reside between 0.1 GeV and 30 TeV. The B − L gauge coupling is chosen so as to produce the correct relic abundance. A summary of the range of sampled parameters is contained in Table 1. Since the interaction rate scales ∝ g 4 B−L , this scan strategy is guaranteed to find the unique coupling giving rise to the correct relic abundance. The chosen values of the parameters are only retained should they pass the following tests: (i) the gauge coupling needed to produce the correct relic density is allowed by current constraints (see Fig. 1) and (ii) the values of the couplings remain perturbative (the perturbative limit is taken here to be √ 8π for the Yukawas and √ 4π for the B − L gauge coupling). This sampling procedure is run continuously for 24 hours, resulting in a total of ∼ 7.5 × 10 5 points that are capable of providing a viable dark matter candidate. We emphasize that for a large region of parameter space, namely 150 GeV m Z 3 TeV, there did not exist a single point passing the scan; the strength of the LHC pp → + − bound accounting for this exclusion has recently been pointed out in [26].
One might wonder whether a slight reduction of the constraints arising from a partial branching fraction to invisible final states might significantly alter the results of this initial parameter scan. In order to investigate this possibility we repeat the scan, but only retaining values of g B−L that reside between the current bound g lim B−L and the reduced bound, where the reduced bound is conservatively defined as g lim B−L /0.7. This factor has been chosen so as to reflect the strength of the bound that would be derived should the maximally allowed invisible branching ratio be realized. In this case, an additional constraint is imposed that requires either m χ or m N1,N2 < m Z /2, since these are the decays responsible for the reduced bounds.
The most important results of the parameter space scan can be summarized in the following: • if the dark matter mass m χ 900 GeV, the model only remains viable if annihilations are on-resonance; • strong constraints from high energy dilepton searches rule out the model for 150 GeV m Z 3 TeV, even for annihilations proceeding through a resonance; • as will be discussed shortly in Section 5, the viable parameter space passing resonance cuts for m χ 5 TeV and m Z 4 TeV often have deep theoretical flaws arising from the appearance of low-energy Landau poles or an instability of the potential.
We defer a more comprehensive discussion of the parameter scan results to Section 5.

Parameter
Description Table 1. Ranges of the parameters explored in the parameter scan. As discussed in Section 4.3, values of m Z were sampled using 200 logarithmically spaced points. For each point random values of mχ, mN 1 = mN 2 , mρ were drawn from the range detailed in the table above (assuming a log-flat distribution), and gB−L was set to be consistent with the measured relic abundance at 2σ, i.e. ΩDM h 2 = 0.120 ± 0.003 [66]. We further require all Yukawa couplings to be ≤ √ 8π and gB−L ≤ √ 4π in order for the model to remain perturbative. The total number of points in parameter space that meet such requirements are 7.5 × 10 5 .

Annihilation Phenomenology
By analyzing the viable parameter space identified with the scan described in Section 4.3, we can determine the relevant annihilation phenomenology in various regimes. We begin by presenting some of the relevant annihilation channels in Fig. 4 -this is not intended to be a comprehensive list but rather it serves as a representative sample. In Fig. 4, we have also listed the annihilation cross section's leading order behavior in velocity. Assuming the dark matter abundance does not proceed through a resonance, the annihilation phenomenology can generically be summarized with two points: (i) if kinematically accessible, then the annihilation proceeds to Z ρ, and (ii) if this channel is not accessible, the annihilation predominantly to N N . Both channels are s-wave and lead to SM final states that are detectable by indirect searches. In order to highlight some of the most important features of the annihilation phenomenology for the viable parameter space, we identify three distinct regimes that are differentiated by the mass of the Z and motivated by the weakness of collider constraints in those regimes; namely, we consider 1 GeV < m Z < 10 GeV, m Z ∼ m Z and m Z > 3 TeV. (a) m Z ∼ m φ(1020) . We have found that dark matter of m χ > 50 GeV with m χ 0.002 g B−L 50 GeV annihilates at a rate of σv ∼ 2 × 10 −26 cm 3 /s to either Z ρ or N N in the contemporary Universe. Thus, limits from dark matter annihilation in dwarf galaxies [68] rule out m χ 100 GeV, or equivalently g B−L 9 × 10 −4 .
2. m Z ∼ m Z . Constraints derived in this region in parameter space are strongly suppressed due to the large backgrounds from the SM Z (as shown in Section 3.1). However, this region of parameter space can potentially be tested with indirect searches. For 86 GeV m Z 97 GeV and m χ 500 GeV, the annihilation cross section is σv ∼ 2 × 10 −26 cm 3 /s, and proceeds either to Z ρ (if 2m χ > m Z + m ρ ), or to N N is the former is kinematically forbidden. Although this region of parameter space is currently unconstrained by indirect searches, future searches could potentially probe this region. It is also worth mentioning that a dedicated study of the the dilepton search channel probing invariant mass distributions close to the mass of the SM Z could conceivably test this region, since relatively large gauge couplings are required to produce the relic abundance (i.e. g B−L 0.034 500 GeV mχ ).
(a) 2m χ > m ρ + m Z . Since the Z ρ channel is accessible, it will always dominate, and if 2m χ > 1.5 (m ρ + m Z ) the annihilation rate as relevant for indirect searches is σv > 10 −26 cm 3 /s, and thus potentially testable with indirect observations from future groundbased observatories.
(b) 2m χ < m ρ + m Z . If the Z ρ channel is closed, the annihilation final state can be ρρ, Z Z , f f and N N , depending on the masses of these particles. In particular if g B−L < 1, then the annihilation proceeds predominantly to N N via ρ exchange, orf f via resonant Z exchange; the Z Z channel, on the other hand, only occurs in the resonant region m χ ∼ m ρ /2. We note that the annihilation to SM fermions is p-wave suppressed, while the annihilation to N N is s-wave, but when mediated by the ρ is Yukawa suppressed. Thus, only if m ρ m χ and m N m χ the annihilation rate will be comparable to σv ∼ 10 −26 cm 3 /s, which represents a very small portion of parameter space. This can explicitly be seen in Figure 7.

Summary & Theoretical Considerations
Thus far we have identified the parameter space for which collider constraints on the U (1) B−L gauge boson do not exclude the possibility that one of the right-handed neutrinos constitutes the dark matter of our Universe. The important question that remains is whether or not the viable parameter space is theoretically well-motivated. In this section we attempt to address exactly this question. Here, we show that much of viable parameter space exists only on-resonance, and most of the viable parameter space not on-resonance exhibits deeper theoretical issues related to the appearance of a Landau pole or an instability of the scalar sector at remarkably low energies, which become apparent only through an analysis of the RGEs.
We begin by considering how many of the viable points identified in the parameter scan require production through a resonance. Models annihilating through the Z or ρ resonance present the most likely candidates for surviving this analysis, and thus we begin here. For m χ ≤ 900 GeV, the requirement that m χ ≥ 0.6m Z or m χ ≤ 0.4m Z removes more than 70% of the originally identified viable points. Further imposing that m χ ≥ 0.6m ρ or m χ ≤ 0.4m ρ removes an additional ∼ 23%. A quick look at the remaining parameter space reveals that the mass of the Z must be near the mass of the Z (i.e. m Z ∼ m Z ), a consequence of the fact that collider constraints are weakened near the mass of the Z. Imposing the requirement that m Z > 100 GeV or m Z < 80 GeV removes all remaining points. Thus, one can concretely say that for dark matter masses m χ ≤ 900 GeV the only viable parameter space is produced via a resonance or is nearly degenerate with the SM Z.
One can perform a similar analysis as above, systemically removing resonance points, but for larger values of the dark matter mass. We have performed this analysis for two independent regions, 900 GeV ≤ m χ ≤ 5 TeV and 5 TeV ≤ m χ ≤ 6 TeV; removing these resonances eliminates ∼ 55% and ∼ 20% of the originally identified viable points, for each region respectively (a summary of these cuts is provided in the top half of Table 5).
The two more conceptional issues we address in this section deal with the appearance of a Landau pole in the running of the B − L gauge coupling, and the stability of the scalar sector, both of which require analyzing the RGEs of this model.
The two-loop RGEs are calculated using SARAH [69,70]. For simplicity, only the one-loop RGEs for the gauge coupling, the quartic coupling of the new scalar, and the couplings of the right-handed neutrinos to the new scalar are provided below 5 : where Yukawa contributions Y αβ have been set to zero as these are negligible. A more comprehensive discussion of the RGEs in the U (1) B−L model can be found in e.g. [76,77]. We define the couplings at the freeze-out scale, i.e. m χ /20. If m Z is above this scale, the B − L gauge coupling is not evolved below the mass of the Z . After preforming the aforementioned resonance cuts, we use the RGEs to determine at what energy either (i) a Landau pole appears, (ii) the scalar sector becomes unstable (i.e. when either 4λ H λ φ − λ 2 Hφ becomes negative or when λ H , λ φ < 0), or (iii) one of the couplings becomes nonperturbative (i.e. exceeds √ 8π). If this energy is below Max(m χ , m Z , m ρ , m N ), we deem this model to have an inconsistency and reject this point in parameter space. The results for these cuts, in the order they have been presented above, is listed in Table 5 and shown graphically in Fig. 6. In Fig. 5 we show two examples of the RGEs that highlight the instability of the scalar sector (left) and the appearance of a Landau pole (right).
For dark matter masses in the range 900 GeV ≤ m χ ≤ 5 TeV, the perturbativity cut, and in particular the cut on scalar stability, drastically reduces the viable parameter space for m Z 5 TeV; in particular, after all aforementioned cuts, only ∼ 7% of the originally identified points remain. The same cut applied for the 5 TeV ≤ m χ ≤ 6 TeV region similarly removes nearly all points, except for those in the range 5 TeV m Z 10 TeV (in total, roughly 30% of the originally identified points remain). 5 Note that in principle there exist a total of three gauge couplings that enter a model with two U (1) symmetries, the extra coupling arising from the potential kinetic mixing term. We have implicitly assumed in the above discussion that this mixing is effectively zero, as non-zero Z − Z constraints can be stringent (see e.g. [71][72][73][74][75]). However, in principle, the RGEs may evolve so as to induce a kinetic mixing at higher energies; e.g. defining the mixing to vanish at a scale µ( = 0) induces a non-zero mixing of magnitude (µ) ∼ 0.015 × g B−L log µ( =0) µ . We have verified that for this model, and for the energies of interest, this does not happen. Thus for simplicity we have set this kinetic mixing term in the RGEs shown in Eq. 5.1 to zero, however its evolution has been included in the calculations. In the reduced bound analysis, the low m χ regime remains effectively unchanged; we note that a small amount of parameter space has become viable in the vicinity of m Z ∼ m Z , however a slightly wider (i.e. ∼ 5 GeV) resonance cut would have eliminated all of these points. In the intermediate dark matter mass region, a majority of the viable points lie at large values m Z (i.e. m Z 4 TeV), although the number of viable points at m Z 100 GeV has also marginally increased. Finally, for m χ > 5 TeV there exists a large amount of parameter space for both m Z ∼ O(GeV) and m Z 4 TeV (also, as before, there exists a region near the mass of the Z only slightly evading the resonance cuts). An intriguing feature that appears in the high dark matter mass reduced bound analysis (5 TeV ≤ m χ ≤ 6 TeV) is the non-existence of viable points for m Z 15 TeV 6 . This is a feature that arises from the fact that the only constraint on the Z in this region is from LEP, and those points producing the correct relic abundance are not naturally occurring in the vicinity of this bound (this is shown in Fig. 7 for the high-mass parameter space scan, to be discussed below).
Finally, before concluding this section, it is worth commenting on constraints coming from unitarity. The conditions for maintaining unitarity in dark matter models with additional U (1) symmetries have been previously derived in e.g. [78,79]. These require that the dark matter mass and scalar mass satisfy and the gauge and Yukawa couplings to satisfy While we do not find that enforcing unitarity constrains any of the viable parameter space after applying the aforementioned cuts, we emphasize that such constraints may strongly limit the dark matter mass from above, further constraining this model. More specifically, we find that for values of m Z ∼ 30 TeV, typical couplings required to meet the relic density constraints are of order g B−L 2.5, implying m χ , m ρ 20 TeV. In order to make more conclusive statements about the validity of the high-mass region, we perform a specialized scan in which we randomly sample in a linear scale m χ , m Z , m ρ , m N ∈ 0.5−100 TeV, and find the g B−L coupling that provides the correct relic abundance (without violating constraints from LEP and the LHC). In Fig. 7, we show that off-resonance, the mass of the Z , the dark matter, and the scalar must be 40 TeV. If one considers resonances, values as high as m Z ∼ 80 TeV are attainable, while both m χ , m ρ 40 TeV. The color of the points identified in Fig. 7 also illustrates the scale Λ at which either the Landau pole appears, or the scale at which the scalar potential becomes unstable. In order to maintain the validity of this model to moderately large energy scales (Λ ∼ 10 4 TeV), the Z mass must be 20 TeV. Comparing this number with the maximal Z mass that can be probed by near-future colliders (i.e. m Z 7 TeV), we see that future HL-LHC constraints are unlikely to provide significant insight into this potentially viable window of parameter space.

Non-Minimal Extensions
We have shown that in the minimal model, null collider searches for the B − L gauge boson negate the possibility of thermal dark matter if m χ < 900 GeV, except for on-resonance regions (i.e. m χ m Z /2, m ρ /2) or near the mass of the SM Z (i.e. m Z ∈ (m Z ± 5 × Γ Z )). However, such strong constraints only probe the coupling and Z mass. In the minimal model presented in Section 2, the vev of the scalar is directly related to the mass and coupling of the Z : m Z = 2g B−L v B−L . However, upon introduction of an additional scalar it is possible to alter this relation, such that the relation between m Z and g B−L is given by 7 Here, v 2 corresponds to the vev of the scalar with charge −2 under B − L, which has been introduced in order to provide the Majorana mass terms to the sterile neutrinos, and v q corresponds to the vev of  Figure 6. Histograms showing the fraction of sampled points removed by each of the cuts described in Section 5 for mχ < 900 GeV (top), 900 GeV < mχ < 5 TeV (middle), and 5 TeV < mχ < 6 TeV (bottom) as a function of m Z . Left (right) panels denote the results for the full (reduced) bound analysis. Viable points surviving all cuts are shown in pink. When no histogram is shown (e.g. 150 GeV < mχ < 3 TeV) it is understood that no points in the sampling routine were capable of producing a viable dark matter candidate. Numerical values detailed the percentage of points surviving each cut are shown in Table 5. some new scalar. The pertinent question at hand is: how tuned must the ratio v 2 /v q be such that for m χ < 900 GeV there are sizable regions of parameter space where the dark matter can be thermally produced without requiring m χ m Z /2, m ρ /2 or m Z ∼ m Z .
If the charge q is such that the Lagrangian respects the symmetry φ q → φ q , then a physical Goldstone boson, the Majoron [9], will appear in the spectrum. On the other hand, if this symmetry of the scalar potential is absent, then a massive pseudoscalar will form upon symmetry breaking. In order to understand the impact of each scenario we will consider two models. One with q = 1/2, which therefore includes a physical Goldstone boson in the spectrum, and other with q = 1 in which a massive pseudoscalar will instead be present.
It is worth stressing that the modification of the g B−L to m Z relation through the presence of new scalar states is present in a variety of models in the literature, and thus deserves a dedicated study. Note that in both cases we assume that dark matter is one of the three sterile neutrinos (i.e. has L = 1), so that it couples to φ 2 and thus the annihilation channel χχ → N N mediated by the scalar and pseudoscalar fields is not further suppressed by the singlet scalar's mixing angle, tan β = 2v 2 /(qv q ), as in models where the dark matter fermion has L = 1 [35]. In this sense, the scenarios analyzed below can be considered as the most favorable ones within the gauged U (1) B−L symmetric paradigm. It is also worth pointing out that the inclusion of additional particle content in a U (1) B−L model leads to an increase in the β function for the B − L gauge coupling, and thus generically many extensions may induce a Landau pole at rather low scales.
6.1 Additional Sterile Scalar with q = 1/2 Upon the introduction of a new scalar with charge q = 1/2, one can write the potential as: where H represents the Higgs doublet. For such charge assignment, the potential has two global symmetries U (1) D × U (1) B−L , where the first corresponds to φ q → e q D iθ φ q and the second to φ 2 → e 2iθ φ 2 , φ q → e qiθ φ q [80]. When the scalars develop a vev φ 2 = 1 √ 2 (v 2 + ρ 2 + iη 2 ) and φ q = 1 √ 2 (v q + ρ q + iη q ), both symmetries will be spontaneously broken and the minimization of the potential yields two CP-odd massless scalar fields η 2 and η q . The mass matrix for the CP-even scalars will be given by in the h, ρ 2 , ρ q basis. In order to simplify the phenomenology we will assume that the mixing with the Higgs is negligibly small (i.e. λ H2 = λ Hq = 0), thus negating various measurements at the LHC.
To further simplify the discussion we consider the case where λ 2q = 0. We note that this last step is not required by any phenomenological reasons, but reduces the number of parameters to be explored. This therefore leads to the following relations: Turning to the phenomenology of the Goldstone boson, after the breaking of both the U (1) D and U (1) B−L symmetries, the two massless CP-odd states originally contained in φ 2 and φ q mix. One linear combination of these scalars forms the longitudinal mode of the Z gauge boson and the orthogonal combination renders a massless Goldstone boson. This can be explicitly expressed as where η G is the physical Goldstone boson. Such a relation can be expressed as a rotation of the fields: This formulation is convenient as it allows for a linear realization of the scalar fields. Namely, we can write the scalar fields as In the Unitary Gauge, one rotates away the longitudinal modes of the gauge fields, and thus from the scalars one can fully express the interacting Lagrangian for the fermion fields. Since only χ and N have Yukawa interactions with the φ 2 scalar, they are the only fields that interact directly with such a boson. The interacting Lagrangian that concerns the fermions and scalars in the model therefore reads In the above expression we haven't explicitly written the scalar-scalar and gauge-scalar interactions that will originate from both the scalar potential and the kinetic interaction. We will assume here that m η G = 0, as we expect the gravity breaking contributions to be non-perturbative [81], and thus negligible. In the case that m η = 0, the presence of such a massless particle in the spectrum only leads to a measurable impact in CMB Stage IV experiments [82] if it is still in thermal contact with the SM bath for T < 0.5 GeV, however such Goldston boson interacts weakly and has momentum suppressed couplings to fermions and thus likely decouples considerably earlier (see e.g. [80] for cosmological implications in the case that m η G = 0).

Additional Sterile Scalar with q = 1
Here, we consider the case of an additional sterile scalar with charge q = 1. After taking λ H2 = λ Hq = λ 2q = 0 for simplicity, the scalar potential reads where we take the trilinear coupling µ 21 to be real and positive. For such a charge assignment, the mass matrix in the basis {h, ρ 2 , ρ q , η 2 , η q } is given by where the minimization condition requires Upon diagonalization, this leads to the following physical CP-even scalar mixing where the fields ρ 2 and ρ q are the CP-even scalars in the interaction basis, and ρ 2 and ρ q are mass eigenstates. In terms of the couplings, the mass of the scalars read On the other hand, for the CP-odd states: In this model we denote the physical pseudoscalar with a in order to remind the reader that in this case it is massive. Again, as for the model described in Section 6.1, one can express the scalars as In the limit m a → 0, it follows that µ 21 → 0, and thus λ 2 → m 2 ρ2 /(2v 2 2 ) and λ 1 → m 2 ρ1 /(2v 2 1 ). It is worth mentioning that in order to have a well-defined mass matrix, namely, that the masses of the CP-even scalars are positive, one must require that det[M] CP-even > 0. This results in the following condition for the mass of the pseudoscalar: (6.20) Since the masses of the scalar fields are related to the couplings, we will scan λ 1,2 ∈ 10 −5 − √ 8π rather than the masses directly. The mass of the pseudoscalar will be varied from 0.1 − 30 × 10 3 GeV subject to the consistency requirement given in Eq. 6.20.

Results
The inclusion of an additional scalar field changes the dark matter phenomenology due to the appearance of two new annihilation channels. The first, χχ → a a is p-wave, and if m a = 0, it is always open. The second, χχ → a ρ 2 is s-wave, and if open, dominates over the χχ → a a channel. In addition to the new aforementioned channels, the annihilation χχ → N N mediated by a and ρ could be enhanced by a factor v B−L /v 2 , a consequence of the fact that v 2 may be lower than in the minimal scenario. In order to draw quantitative conclusions about the required values of v 2 /v q needed to obtain the correct relic abundance without requiring m χ m Z /2, m ρ /2, we performed a parameter scan for both the case m a = 0 and m a > 0 case assuming that the dark matter is relatively light m χ < 900 GeV and the Z mass is ≤ 3 TeV; we have chosen this region as it is directly testable in terrestrial experiments. The precise details of the parameter scan are presented in Table 3.
We found that if the scalar masses are such that 2m χ > m ρ2 + m a , then the χχ → a ρ 2 channel will always dominate the annihilation in both the early and the contemporary Universe. If that condition is not met, then the χχ → N N channel will dominate the annihilation, accompanied by the χχ → a a annihilation should it be kinematically allowed. Therefore, models with a non-minimal scalar sector will allow the dark matter candidate to be formed in the early Universe through the annihilation to the aforementioned channels in a manner similar to what happens in the global B − L scenario (studied at depth in [16]). The required level of tuning of the vevs, in order not to require m χ m Z /2, m ρ /2, depends on both the Z mass and the dark matter mass, and typically ranges between v 2 /v q ∈ 0.001 − 0.1 as shown in Fig. 8 (note that in Fig. 8 we have chosen to remove resonant points, defined by 0.4 < m χ /m ρ2 , m χ /m ρq , m χ /m a , m χ /m Z < 0.6). s-wave: Z ′ ρ 2 NN ρ 2 a Z ′ Z ′ p-wave: aa f f Figure 8. Viable values of v2/vq identified for the non-minimal models described in Section 6.1 (top) and Section 6.2 (bottom) as a function of the Z mass. The colors represent the dominant annihilation channel in the contemporary Universe. Resonant points, defined by 0.4 < mχ/mρ 2 , mχ/mρ q , mχ/ma, mχ/m Z < 0.6 have been removed.
As shown in [16], for the case m a = 0 unless m χ m N > 100 GeV the final state of s-wave annihilations will lead to invisible particles. The reason is that the main decay of light sterile neutrinos (m N m W ) will be N → ν a which does not render gamma-rays. Similarly, the decay of the ρ 2 scalar will be ρ 2 → a a, and is thus also invisible. If the pseudoscalar is not massless the scenario may differ. In this case, the massive pseudoscalar will always decay to light active neutrinos unless m N < m a /2. Therefore, if the annihilation is dominated by χχ → N N the sterile neutrinos will decay to SM particles, and m χ 100 GeV are excluded from dwarf galaxies observations [83][84][85]. Dark matter masses m χ 100 GeV are also excluded from gamma-ray observations in the case that the annihilation is dominated by χχ → a ρ 2 , provided that the sterile neutrinos are light enough, i.e. m N < m ρ /2, m a /2. In the case that χχ → a ρ 2 dominates and m N > m ρ /2, m a /2, then the final state will be entirely composed of light active neutrinos. Current bounds form neutrino telescopes on such annihilation channels are about 3 orders of magnitude above σv = 3 × 10 −26 cm 3 /s [86,87] and thus this portion of parameter space will remain elusive.

Parameter
Description Table 3. Ranges of the parameters explored for the non-minimal models. Values for m Z were sampled using 154 logarithmically spaced points. The value of the gB−L coupling was not taken as a free parameter, but instead fixed to the 2σ exclusion value for each m Z value. For each point random values of mχ, mN 1 = mN 2 , mρ 2 , mρ q , ma were drawn from the range detailed in the table above (assuming a log-flat distribution, except for ma in the q = 1/2 case where ma = 0), and tan β or equivalently v2 was set to be consistent with the measured relic abundance at 2σ, i.e. ΩDM h 2 = 0.120 ± 0.003 [66]. We further require all coupling to be < √ 8π in order for the model to remain perturbative. The total number of points in parameter space that meet such requirements are 1.54 × 10 5 for both the q = 1/2 and q = 1 cases.

Conclusions
In this work we have presented a thorough investigation of the simplest extension of the Standard Model that generates active neutrino masses and contains a viable dark matter candidate within the context of a gauged U (1) B−L symmetry. A number of groups have addressed aspects of similar models, e.g. by considering a Dirac dark matter candidate (for which direct dark matter detection experiments are strongly constraining) [34][35][36][37] or a simplified subset of parameter space [22,[24][25][26][27]29], however until now a comprehensive understanding of the viability of this model has been absent from the literature. Specifically, we take one of the right-handed neutrinos required to achieve anomaly cancellation to be a Majorana particle and stable under a new symmetry, thus providing a dark matter candidate with suppressed direct detection constraints. We investigate the unavoidable limits from various colliders on the existence of the new gauge boson, and consider deeper theoretical concerns such as the appearance of a Landau pole or instability of the scalar sector arising at energies as low as the TeV scale.
After performing a full parameter space scan, our results reveal a number of important points: (i) dark matter with mass 900 GeV must annihilate on-resonance to remain viable, (ii) high energy dilepton searches rule out 150 GeV < m Z < 3 TeV (even when annihilations proceed on-resonance), and (iii) the theoretical considerations, such as annihilations that proceed on-resonance, the appearance of a Landau pole, or the instability of the scalar sector at very low energies strongly disfavors models in which m χ , m Z 5 TeV. While unitarity and perturbativity constraints require m χ , m Z 40 TeV, we have also shown that avoiding low-energy Landau poles and instability in the high-mass region likely require m Z 20 TeV. Unfortunately, constraints from the HL-LHC cannot close this window (we estimate the maximal sensitivity to be approximately m Z 7 TeV), implying the small remaining viable and theoretically well-motivated parameter space is largely untestable. Thus, we believe it is safe to say that the minimal gauged U (1) B−L extension of the Standard Model with right-handed neutrino dark matter is under siege, and should at this point be considered either a disfavored or untestable extension of the Standard Model.
The stringent constraints on the minimal model raise the question of whether minor modifications could relieve existing tension with experimental bounds. Perhaps the simplest, however entirely representative, extension for which this can be achieved involves introducing a new scalar particle; this modification allows one to decouple the relation between the vev of the scalar and the ratio of m Z to g B−L . We have considered two such extensions, one which gives rise to a physical Goldstone boson and the other which gives rise to massive pseudoscalar. We investigate the viable parameter space in these models, and show that typically the tuning between the ratio of the vevs of the B − L charged scalars must be ∼ 10 −2 . The addition of more complicated particle content to resolve the aforementioned problems inherently increases the instability of the β functions, potentially yielding problems at rather low scales.
In summary, constraints on the minimal gauged U (1) B−L extension of the Standard Model with thermal right-handed neutrino dark matter significantly reduce its current appeal. Minor modifications via the introduction of a new scalar can revive the model, however only at the cost of a potential finetuning of the vevs of the scalar fields. Given the theoretical motivations for this model, one might wonder whether non-thermal dark matter production mechanisms could revitalize interest in righthanded dark matter within a gauged U (1) B−L extension.

A Appendix: Cross Sections
For completeness we list the various annihilation cross sections, in expansions of velocity v, for annihilations to SM fermions, right-handed neutrinos, two Z s, two scalars ρ, and one Z and one scalar ρ: