Gauged Inverse Seesaw from Dark Matter

We propose an economical model addressing the generation of the Inverse Seesaw mechanism from the spontaneous breaking of a local $U(1)_{B-L}$, with the Majorana masses of the sterile neutrinos radiatively generated from the dark sector. The field content of the Standard Model is extended by neutral scalars and fermionic singlets, and the gauge group is extended with a $U(1)_{B-L}$ and a discrete $\mathbb{Z}_4$ symmetries. Besides dynamically generating the Inverse Seesaw and thus small masses to the active neutrinos, our model offers two possible dark matter candidates, one scalar and one fermionic, stable thanks to a remnant $\mathbb{Z}_2$ symmetry. Our model complies with bounds and constraints form dark matter direct detection, invisible Higgs decays and $Z'$ collider searches for masses of the dark sector at the TeV scale.


Introduction
The origin of neutrino masses is one of the big open questions in Particle Physics and, consequently, a plethora of explanations have been proposed, see e.g. Ref. [1]. The type-I seesaw mechanism [2][3][4][5][6][7][8] is maybe the simplest and most popular one, which requires the extension of the Standard Model (SM) particle field content with new neutral leptons, known as right-handed (RH) or sterile neutrinos. Although the most minimal seesaw realization capable of accounting for neutrino oscillation data [9,10] requires two RH neutrinos, the case of three RH neutrinos is particularly interesting, as it restores the balance between the number of quarks and leptons, canceling the gauge anomalies of U (1) B−L . In this case, U (1) B−L can be promoted from being an accidental global symmetry of the SM to a local gauge group. This kind of models have been extensively studied in the literature, see e.g. Refs. [11][12][13][14][15][16][17][18][19][20].
In order to accommodate sub-eV neutrino masses [21], the canonical type-I seesaw model requires either a high lepton number violating (LNV) scale, or tiny Yukawa couplings if realized at low scale. Either way, the phenomenology of this model is extremely suppressed, which makes it difficult to probe experimentally. An alternative to have a low-scale realization with large Yukawa couplings is to consider additional sterile fermions and an approximate symmetry. This is the case of the the inverse seesaw (ISS) [22][23][24] or the linear seesaw (LSS) [25,26] realizations, where the B − L global symmetry is used to protect active neutrino masses, linking their smallness to small parameters that quantify the breaking of the B −L symmetry. Here, the neutral leptons are introduced in pairs of opposite B − L, so they cannot cancel the anomalies of U (1) B−L . Therefore, gauging low-scale seesaw models requires to extend the particle content with new exotic fermions, as it was done for instance in Refs. [27][28][29][30][31][32][33][34][35].
Additionally, to generate tree-level Majorana masses for the active neutrinos, low-scale seesaw models introduce a new scale that breaks explicitly lepton number symmetry either with a Majorana mass term (as in the ISS) or via a Yukawa interaction for the new fermion singlets (as done in the LSS). In these scenarios, LNV parameters are ad-hoc and are assumed to be small, an hypothesis that, despite being technically natural and thus stable under radiative corrections, lacks of a more fundamental motivation. Possible explanations for the origin of these parameters, in particular the small Majorana mass term in the ISS construction, have been proposed in several models [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54].
In this work we explore the possibility of gauging the ISS model by adding only sterile neutrinos in the fermionic sector and, at the same time, providing a dynamical one-loop generation of the Majorana mass µ for the sterile neutrinos, which justifies its smallness. Furthermore, this model provides viable dark matter (DM) candidates either from the sterile neutrino sector or from the extended scalar one, which are actually the particles behind the generation of the µ mass. As a consequence, the neutrino and dark sectors are connected in this setup.
Our model considers an extension of the SM by a local U (1) B−L symmetry and an additional Z 4 symmetry, which protects active neutrinos from obtaining tree-level masses and ensures the stability of the DM candidates. The scalar sector is extended with three scalar singlets, σ that breaks U (1) B−L , χ that breaks Z 4 to a conserved Z 2 symmetry, and an extra inert one ζ. We will assume that the breakings occur at the TeV scale in order to be consistent with LHC bounds on the new heavy vector boson [55,56]. The fermionic sector contains new singlets: three RH neutrinos N R that are in charge of canceling the gauge anomalies of U (1) B−L , and two additional pairs of sterile neutrinos (ν R , ν S ) with opposite B − L numbers, so that the low-scale seesaw is realized without spoiling the anomaly cancellation. Due to the Z 4 symmetry, the N R fields do not mix with the active neutrinos and the ν R,S neutral leptons are not allowed to acquire Majorana mass terms, consequently the active neutrinos remain massless at the tree level, even after the spontaneous symmetry breaking of U (1) B−L ⊗ Z 4 . Nevertheless, our setup generates Majorana masses for the ν S fields radiatively at the one-loop level, providing the origin of its smallness, and triggering the ISS mechanism. Furthermore, we will see that this loop contribution is proportional to small parameters protected by an accidental symmetry, and thus they could additionally suppress the µ term. Moreover, the model provides two different DM candidates, whose stability is ensured by the residual Z 2 symmetry. Depending on the mass hierarchy, we could have a fermionic DM candidate from the lightest of the three N R , or a scalar candidate from the lightest component of the inert scalar ζ (real ζ R or imaginary ζ I component). Interestingly, the same parameter needed to dynamically generate the µ term of the ISS will break the mass degeneracy of ζ R and ζ I , connecting thus both the DM and neutrino sectors. It is worth stressing that the loop generating the µ term involves the new DM candidates, making our mechanism similar to the scotogenic model for neutrino masses [57], i.e. our model provides a scotogenic origin of the ISS µ term. We find that in order to comply with neutrino data, DM direct detection searches, invisible Higgs decays and Z collider searches, the masses of the DM candidates are larger than few TeV. Furthermore, in the fermionic DM scenario, the U (1) B−L gauge coupling needs to be above 0.7 and the Z gauge boson heavier than about 6-20 TeV.
The paper is organized as follows: in Section 2 we introduce the model, providing a detailed description of the scalar sector. The fermionic sector is discussed in Section 3, with details on the active and sterile neutrino mass generation. The viability of both scalar and fermionic DM candidates is discussed in Section 4. Details on the anomaly cancellations, stability and unitarity conditions are collected in the Appendices. We summarize our findings in Section 5.

The model
We propose an extension of the SM where the gauge symmetry is extended by the inclusion of spontaneously broken U (1) B−L gauge and Z 4 discrete symmetries, and the particle content is enlarged with new scalar and fermionic singlets. The particle spectrum and their charge assignments are shown in Table 1.
The scalar sector contains three new electrically neutral scalar gauge singlets σ, χ and ζ. The spontaneous symmetry breaking (SSB) of the U (1) B−L gauge symmetry is triggered by the vacuum expectation value (VEV) of σ, which we will take above the TeV scale in order to comply with collider bounds on the Z gauge boson [55,56]. On the other hand, the VEV of χ breaks the Z 4 symmetry down to a preserved Z 2 one, which will ensure the stability of the DM particles. Finally, the scalar singlet ζ does not acquire a VEV, providing thus a scalar DM candidate with the lightest of its CP-even and -odd components. As we will see, ζ plays a key role in implementing the active neutrino mass generation mechanism, inducing the ISS µ term at the one-loop level.
The fermionic sector is extended with new sterile fermions ν R , ν S and N R , all of them singlets Starting with the SM Higgs boson doublet φ, the scalar sector is collected in the last four columns, preceded by three columns for the new sterile fermions, while de SM fermions are collected on the four first ones. More details are provided in the main text. The indices run as follows: i = 1, 2, 3 and k = 1, 2.
under the SM group but with different U (1) B−L ⊗ Z 4 charges, thus playing different roles in the model. The N R fields have the same B−L charge as the SM leptons, allowing to cancel the gauge anomalies of U (1) B−L . Indeed, this requirement imposes the number of N R to be three, one per generation. 1 Due to the Z 4 symmetry, they do not mix to the active neutrinos, preventing them from generating tree-level neutrino masses and providing a fermionic DM candidate, the lightest of the three states N R . On the other hand, the ν R and ν S fields have opposite B − L charges and do not contribute to the anomaly cancellation of B − L as long as they are introduced in pairs. Interestingly, this same condition allows us to implement an ISS mechanism for the neutrino mass generation. Notice that in principle we do not have any constrain on the number of (ν R , ν S ) pairs in the model, nevertheless we will consider only two pairs, the minimal ISS scenario able to accommodate neutrino oscillation data [58].
Regarding the neutral fermion masses, given the charge assignments displayed in Table 1, they are all dynamically generated by the SSB of U (1) B−L ⊗Z 4 . The scalar σ provides Majorana masses for the N R fields and the χ generates Dirac-like masses for the (ν R , ν S ) pairs. Notice that the Z 4 symmetry prevents the ν R and ν S fields from acquiring Majorana masses at tree level, and thus keeps the active neutrinos massless. Nevertheless, as it is discussed later, the N R fields together with the real and imaginary parts of scalar singlet ζ generates a Majorana mass for ν S at the one-loop level, triggering an ISS mechanism responsible for the light (active) neutrino masses.

The scalar potential
The most general scalar potential invariant under the symmetries of the model reads Notice that the χ field, given its charge assignment, can be taken to be real. Since we are interested in the case where the fields φ, σ and χ acquire VEVs (v, v σ and v χ , respectively), but ζ does not, we consider in Eq. (1) m 2 , m 2 σ and m 2 χ to be negative, whereas m 2 ζ is taken positive. We focus on the CP-conserving case where the dimensionful trilinear η µ and the quartic couplings are real. Additional constraints on the parameter space come from stability and unitarity conditions, which are discussed and summarized in Appendix B. The scalar fields are given by, where the VEVs v, v σ and v χ break the electroweak symmetry, U (1) B−L and Z 4 , respectively.
Here, h,σ andχ are the three physical scalars remnant from each SSB, while φ ± , φ Z and σ Z are the Goldstone bosons associated to the longitudinal components of the gauge fields W ± , Z and Z , respectively. Minimizing the scalar potential, we find that the VEVs of the scalar fields are solutions of the following equations After SSB, the Z 4 is broken down to a residual Z 2 under which h,σ,χ are even, and ζ R , ζ I are odd, thus decoupling the two sectors. We can write the mass matrix for the Z 2 -even states as The quartic couplings λ φσ and λ φχ control the size of the deviations between h and the SM Higgs boson, which are strongly suppressed from LHC Higgs data [59]. Consequently, we will focus on the decoupling scenario λ φσ = λ φχ = 0, in which h behaves as the SM Higgs. Then, one can analytically find the masses for the eigenstatesσ andχ , defined as with On the other hand, the masses of the real and imaginary part of the ζ field are given by and their squared mass difference is proportional to the dimensionful parameter η µ , Therefore, depending on the sign of η µ , either of ζ R or ζ I could be a viable DM candidate, whose stability is ensured as they are odd under the residual Z 2 symmetry. As we will see in the next section, the same parameter η µ is the key ingredient behind the ISS mechanism, as it will control the dynamical generation of the one-loop Majorana mass for the sterile neutrinos.

Active and Sterile neutrino mass generation
With the symmetries and particle content of our model, the neutrino Yukawa terms in the Lagrangian are given by: Here, the superscript C stands for charge conjugation Ψ C = CΨ T , φ = iτ 2 φ * , and summation over repeated indices must be understood with i, j = 1, 2, 3 and k, r = 1, 2. After the SSB of the local SU (2) L ⊗ U (1) Y ⊗ U (1) B−L ⊗ Z 4 symmetry we obtain a Majorana mass matrix for the neutrinos, which in the basis (ν L , ν C R , ν C S , N C R ) T has the following structure, Notice that it is separated in two sectors: a low-scale seesaw one and a decoupled Majorana mass matrix m N . The dimensions of the submatrices are 3 × 2 for m D , 2 × 2 for M and µ, and 3 × 3 for m N . At tree level, they are given by The structure of the mass matrix in Eq. (15) is a consequence of the Z 4 charge assignment, which also forbids a tree-level mass for ν S even after the SSB. In this case of µ = 0, the active neutrinos remain massless, as the (ν R , ν S ) neutrinos become degenerate, forming a Dirac neutrinos of mass M , and their contribution to the light active neutrino masses vanishes. On the other hand, the Majorana neutrinos N R are decoupled from the active neutrinos, again due to the Z 4 symmetry, they therefore do not contribute to the active neutrino masses at tree level. They are however responsible for generating the Majorana µ term at one-loop level. We show in Fig. 1 the relevant diagram generating the Majorana mass term for the ν S sterile neutrinos, which involves a loop with the N R fields and the inert scalar ζ. Notice that both real and imaginary components of the scalar singlet ζ participate in the diagram, coupled to the scalar singlet χ via the trilinear coupling η µ defined in Eq. (1), which is at the same time responsible for generating the mass splitting between ζ R and ζ I in Eq. (13). From this diagram, we obtain the one-loop contribution given by where we have assumed, without lost of generality, that the submatrix m N is diagonal. This non-zero µ matrix triggers the ISS mechanism in the (ν L , ν C R , ν C S ) sector, generating masses for the active neutrinos of the order 2 m ν ∼ µ m 2 D /M 2 and mass splittings for the pairs of sterile neutrinos of O(µ), thus implying that the sterile neutrinos form pseudo-Dirac pairs.
As usual, a successful ISS realization requires a light µ scale, in particular µ m D , M , which occurs naturally in our model. In order to explore this feature in detail, let us start assuming a small mass splitting between the real and imaginary components of ζ, i.e.
Then, we can write From this equation, we see that the generated µ scale is suppressed by a loop factor, and it is proportional to the y ζ and η µ parameters. Interestingly, in the limit when these parameters are zero, together with the quartic coupling λ ζ , the Lagrangian becomes invariant under a global U (1) ζ symmetry under which only the ζ fields are charged. In this sense, it is technically natural to consider these three parameters, and in particular y ζ and η µ , to be small, a fact that further suppresses the µ scale. Therefore, our model predicts a small scale for the ISS µ term, which is suppressed by a loop factor and three powers of technically natural small parameters. Furthermore, Eq. (19) exhibits the link between the light neutrino masses and the dark sector, and can be used to relate the relevant scales for both sectors. In order to do that, let us recall again that in the ISS mechanism the light active neutrino mass scale m ν is of the order of with U νN the active sterile neutrino mixing. This mixing modifies the neutral and charged lepton currents, leading thus to numerous constraints. Cosmological observations [60][61][62][63] put severe constraints on neutral leptons with a mass below ∼ 200 MeV. Additional strong bounds for sterile neutrinos lighter that the mass of the W gauge boson are imposed from several laboratory searches, see for instance Refs. [64][65][66]. In the following, we focus on the case where the sterile neutrinos realizing the ISS mechanism are heavier than the electroweak scale, which are mainly constrained by electroweak observables and charged lepton flavor transitions [67]. A global fit analysis to these observables [68] imposes bounds on their mixings of the order |U νN | 2 10 −3 , which implies that µ O(1) keV in order to reproduce ∼ 1 eV neutrino masses.
On the other hand, Eq. (19) can be studied in two interesting limits, both related to the mass hierarchy between the lightest N R and ζ defining the DM candidate in our model. In the case where m N m ζ , the model contains a fermionic DM candidate with mass m N , and the µ scale of the ISS can be expressed, using Eq. (18), as In the opposite limit of m ζ m N , we have a scalar DM with mass m ζ , and where we have neglected the logarithmic dependence in the last step. Notice that the case m N ∼ m ζ leads to the same µ term than in the case m N m ζ , with just an additional factor of 1/2. From the above equations, we see that our model can easily generate the correct scale for the ISS µ term, with a DM candidate at the TeV scale, and relatively small y ζ and ∆m 2 ζ /m 2 ζ (and thus η µ ), which are protected by the accidental U (1) ζ symmetry.
In conclusion, our model connects the ISS µ term with the dark sector via the loop diagram in Fig. 1, in a similar fashion as the scotogenic model [57] does for neutrino mass generation. Notice however that in the original scotogenic model the dark sector generates directly small radiative masses for the active neutrinos, while in our model the dark sector generates only the Majorana masses for the sterile neutrinos ν S , which then induce masses for active neutrinos via Figure 2: Main channels contributing to the production of scalar DM ζ i , with i identifying the lightest state between the real and imaginary part of ζ, and j = R, I. Crossed diagrams are not shown, although taken into account in the analysis. f labels all the fermion fields coupled to the Higgs.
the ISS mechanism. In this sense, our model provides a scotogenic origin for the ISS µ term. A similar idea has been proposed, for instance, in Refs. [36,39,40,46], with different particle content and symmetries than ours, and without linking it to the gauging of U (1) B−L , as done in this work.

Dark matter
The residual Z 2 symmetry protects the lightest odd state of the dark sector, rendering it a viable candidate for DM. Depending on the mass hierarchy, we could have scalar or fermionic DM candidates if ζ or N R are the lightest states, respectively. The phenomenology of the WIMP paradigm of these two cases will be studied in the following.

Scalar dark matter
In this section, we consider the case where the DM candidate is the lightest component of ζ, either its real or its imaginary part. In the early Universe, it can be produced out of the SM thermal bath by the WIMP mechanism, via the different 2-to-2 scatterings presented in Fig. 2. DM can annihilate into SM Higgs bosons via the contact interaction, the s-channel exchange of a Higgs, or the t-and u-channel mediation of a ζ. It is also possible to have a scattering into SM charged fermions and vector bosons, mediated by the Higgs portal. Additionally, the new scalarsσ andχ could be in the final state or mediate the annihilation, as shown in the second row of Fig. 2. Finally, DM can annihilate into neutrinos via the t-channel exchange of N R (last  row). This latter channel exists not only for the DM annihilation but also for the coannihilation ζ R -ζ I , which is typically relevant if the mass splitting, see Eq. (13), is smaller that ∼ 20% [69]. However, it is interesting to note that both the annihilation and the coannihilation into active neutrinos depends on the active-sterile mixing angles (to the fourth power), which are suppressed at high temperature (i.e. at T 200 MeV) due to over dominating thermal masses of the active neutrinos [70][71][72], and therefore these channels are subdominant.
In the decoupling limit (λ φσ = λ φχ = 0), where h behaves as the SM Higgs, and if ζ is lighter than the other dark sector particles, the main annihilation channels correspond to the ones on the first row of Fig. 2, and the present scenario reduces to the scalar singlet DM model [73][74][75]. Figure 3 shows with a black line the values for the Higgs portal coupling λ φζ as a function of the DM mass m ζ required to reproduce the whole observed DM abundance. A full numerical computation has been performed using MicrOMEGAs [76,77]. 3 This scenario is constrained by collider measurements [78][79][80][81][82][83][84]. In particular, strong bounds on the invisible decay of the Higgs apply for m ζ 60 GeV, as one can have h → ζ R ζ R and h → ζ I ζ I . The corresponding decay width is given by A recent combination of searches with the ATLAS experiment found an upper bound for this branching ratio Br(h → inv) < 0.11 at 95% CL [85]. The corresponding excluded region is overlaid in brown in Fig. 3. The DM direct detection can further constrain the parameter space. The corresponding spin-independent DM-nucleon scattering cross section is given by [86][87][88][89][90] σ λ 2 φζ f n 4π where m n is the nucleon mass, µ ζ ≡ m ζ m n /(m ζ + m n ) is the DM-nucleon reduced mass, and f n ∼ 0.3 is the hadron matrix element. In Fig. 3 the recent limits from the PandaX-4T experiment [91] are shown in red. Additionally, indirect detection with γ-rays and antimatter could further constrain this scenario [92][93][94][95][96][97][98][99], however the limits are comparable to the ones coming from direct detection and thus we do not show them in Fig. 3. The combination of the direct detection bound together with the invisible decay width excludes a large fraction of the parameter space, letting unconstrained the narrow mass window close to the Higgs funnel (m ζ ∼ m h /2), and the region above the TeV, i.e., 2 TeV m ζ 30 TeV. Notice that larger masses would require non-perturbative λ φζ in order to not overclose the Universe.
Before closing this section, it is interesting to note that in the present scenario alternatives to the WIMP mechanism exists. For example, DM could also have been produced non-thermally in the early Universe via the FIMP mechanism [100][101][102][103][104]. In that case, much smaller Higgs portal couplings are required O(λ φζ ) ∼ 10 −11 with a broader mass range [105][106][107]. Additionally, if the Higgs portal coupling is suppressed and DM has sizable self-interactions due to a quartic coupling of order O(λ ζ ) ∼ 1, thermalization in the dark sector could have a strong impact on the DM abundance [108][109][110]. Moreover, DM could have been Hawking radiated by primordial black holes [111]. Finally, this model has also been studied in the framework of non-standard cosmologies [112][113][114][115].

Fermionic dark matter
In this case, the DM candidate is the lightest N R i state. In the early Universe, DM can annihilate into a) SM fermions and sterile neutrinos ν R and ν S mediated via the s-channel exchange of a Z , together with the annihilation into a pair of Z (first row in Fig. 4), or b) Z σ, Z Z , a couple of scalars or neutral fermions (second row). Finally, c) DM can coannihilate withσ or an active neutrino via the active-sterile mixing (last row). As in the case of scalar DM, here some processes are typically subdominant because the coupling λ φσ has to be very small to avoid tension with Higgs physics, and because the active-sterile mixing angles are suppressed at high temperatures.
If the DM N R is the lightest state of the dark sector, most of the channels described in Fig. 4 are kinematically closed, and therefore the annihilation into SM states via the exchange of a Z boson tends to be dominant. In that case the DM phenomenology is determined by three parameters: m N , M Z and the gauge coupling g B−L . Figure 5 shows contours for the values of the coupling g B−L required in order to successfully reproduce the whole observed DM abundance in the plane [m N , M Z ] (blue thick lines). The red dotted line corresponds to the resonant annihilation case (M Z = 2 m N ), and separates the regimes M Z 2 m N , where DM annihilates into SM pairs via the s-channel exchange of a Z , from the regime M Z 2 m N , where DM annihilates into a couple of Z dominantly via the t-channel mediation of N R . Notice that in the latter limit, also the first 3 diagrams of the second row are relevant, since we could Figure 4: Main channels contributing to the production of fermionic DM. f SM indicates all the SM fermions coupled to the Z and S = h,σ,χ, and ζ R,I . expect mσ to be of the same order of M Z . Nevertheless, achieving the hierarchy M Z m N for g B−L 1 requires large couplings y σ , leading to the non-perturbative regime shown in gray in Fig. 5, and thus this hierarchy is disfavored. Again, we have performed a full numerical computation with MicrOMEGAs.
In the limit of small momentum transfer, the Z exchange leads to a spin-independent scattering with nucleons. The corresponding cross section per nucleon is given by [116] σ g 4  exclude all the viable parameter space for Z lighter than 6 TeV. Furthermore, LEP measurements of e + e − → + − [118] were sensitive to heavy Z bosons, constraining the mass to be M Z > 6 g B−L TeV [119]. Finally, we note that the required gauge coupling g B−L should be higher than ∼ 0.7 in order to successfully reproduce the whole observed DM abundance. Nevertheless, too large g B−L couplings would violate perturbativity, leading to a well-defined allowed area in the parameter space that will be further probed in future experiments.

Conclusions
We have constructed a low-scale seesaw model with a local U (1) B−L symmetry, where some of the sterile neutrinos acquire a small Majorana mass radiatively, triggering then an inverse seesaw mechanism to generate tiny masses for the active neutrinos. In our proposed model, the SM gauge symmetry is supplemented by the spontaneously broken local U (1) B−L and discrete Z 4 symmetries, an extended scalar sector and additional fermionic singlets. The latter allows us to cancel the gauge anomalies of B − L and to realize an ISS mechanism for neutrino mass generation. The discrete Z 4 , spontaneously broken to a residual Z 2 , forbids the sterile neutrinos to acquire tree-level Majorana masses, which is instead generated from a one-loop diagram involving the dark sector particles. This loop suppression, together with the smallness of the couplings involved in the loop, provides a natural explanation for the smallness of the ISS µ term. Additionally, our model contains both scalar and fermionic DM candidates, whose stability is ensured from the residual Z 2 symmetry. We have shown that this model, in addition to providing an origin of the ISS mechanism, reproduces successfully the measured value of the DM abundance and is consistent with the constraints arising from the DM direct detection, invisible Higgs decays and LHC Z searches. The consistency with all these constraints restricts the mass of the scalar DM candidate to be either m h /2 or larger than about 2 TeV for values of the Higgs portal coupling of order unity. On the other hand, in the scenario of fermionic DM candidate, the aforementioned constraints set the masses of the fermionic DM component and the Z gauge boson to heavier than few TeV, and the U (1) B−L gauge coupling such that g B−L 0.7. Smaller couplings induce a DM overproduction that overclose the Universe. Alternatively, much higher couplings generate a DM underabundance, compatible with a multicomponent DM scenario. Future DM searches covering heavy DM candidates, further input from LHC and new detailed measurements of e + e − → + − at future lepton colliders will ultimately add new constraints on the Z mass and g B−L , as well as on the sterile fermion masses, rendering our scenario testable.

Acknowledgments
We would like to thank Miguel Escudero for bringing to our attention the LEP constraints on heavy Z bosons. NB received funding from the Spanish FEDER/MCIU-AEI under grant FPA2017-84543-P, and the Patrimonio Autónomo -Fondo Nacional de Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José de Caldas (MinCiencias -Colombia) grant 80740-465-2020. AECH is supported by ANID-Chile FONDECYT 1210378, ANID PIA/APOYO AFB180002 and Milenio-ANID-ICN2019 044. AECH is very grateful to the Laboratoire de Physique Théorique, Université Paris-Sud for hospitality and for partially financing his visit where this work was started. XM is supported by the Alexander von Humboldt Foundation. This project has received funding/support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 860881-HIDDeN, and from the DFG Collaborative Research Center SFB1258.

A Anomaly cancellation
Through the following conditions, we have checked the gauge and gravitational anomaly cancellation with A the anomaly coefficients appearing in the triangular amplitude involving gauge bosons (corresponding to the groups indicated in the subindices) in the external lines and SM fermions in the internal lines. Here B − L and Y refer respectively to the charges under U (1) B−L and U (1) Y , as given in Table 1.

B Stability and unitarity conditions
In order to determine the stability conditions of the scalar potential, one has to analyze its quartic terms because they will dominate the behavior of the scalar potential in the region of very large values of the field components. To this end, we introduce the following hermitian bilinear combination of the scalar fields and rewrite the quartic terms of the scalar potential as follows: − 2 λa 2 + λ σ b 2 + λ ζ c 2 + λ φσ + 2 λλ σ ab + λ φζ + 2 λλ ζ ac + λ φχ + λλ χ ad + λ ζσ + 2 λ σ λ ζ bc + λ ζχ + λ ζ λ χ cd + λ σχ + λ σ λ χ bd + λ ζ e 2 + h.c .
Regarding unitarity conditions, they are imposed to avoid that the scattering amplitudes involving physical scalars and longitudinal gauge bosons grow too much with energy. In practice, we consider the partial wave expansion, given for an amplitude M by M(θ) = 16π ∞ =0 (2 + 1) a P (cos θ) , (35) with P the Legendre polynomials and a the partial wave of order , and require that the partial waves are bounded from above for any of the possible 2 → 2 scatterings. This process can be simplified thanks to the equivalence theorem [122][123][124][125], as explained for instance, in Ref. [121]. The main idea is that the high-energy behavior of these scatterings can be described using the unphysical scalar particles. Moreover, these scatterings will be dominated by the dimensionless quartic couplings in Eq. (1), which contribute only to the = 0 partial wave, bounded to be |a 0 | < 1. Then, we can easily compute the high-energy scattering amplitudes in the unphysical basis, where all the interactions are just given by the quartic couplings, and then obtain the physical unitarity conditions by requiring that the eigenvalues of this S-matrix are lower than 16π.
The unitarity conditions imply that the modulus of the eigenvalues of S N and S C lie below 16π.