Radiative neutrino mass with ℤ3 dark matter: from relic density to LHC signatures

In this work we give a comprehensive analysis on the phenomenology of a specific ℤ3 dark matter (DM) model in which neutrino mass is induced at two loops by interactions with a DM particle that can be a complex scalar or a Dirac fermion. Both the DM properties in relic density and direct detection and the LHC signatures are examined in great detail, and indirect detection for gamma-ray excess from the Galactic Center is also discussed briefly. On the DM side, both semi-annihilation and co-annihilation processes play a crucial role in alleviating the tension of parameter space between relic density and direct detection. On the collider side, new decay channels resulting from ℤ3 particles lead to distinct signals at LHC. Currently the trilepton signal is expected to give the most stringent bound for both scalar and fermion DM candidates, and the signatures of fermion DM are very similar to those of electroweakinos in simplified supersymmetric models.


Introduction
Neutrino mass and nonbaryonic dark matter (DM) offer two pieces of robust evidence for the existence of physics beyond standard model (SM), but their origins remain mysterious. It would be appealing if they could be understood in the same framework. At low energies neutrino mass can be accommodated by a dimension-five operator in terms of the SM Higgs and lepton fields [1]. The operator can be realized at tree level in three different manners [2] which correspond exactly to the three types of conventional seesaws. Though simple enough, these models are difficult to test experimentally since they invoke very high energy scales or very weak couplings to SM particles in order to induce tiny neutrino masses. One way to alleviate this problem is to push the neutrino mass to a radiative effect of new physics which provides additional suppression. For this purpose, an (almost) exact discrete symmetry is usually required to forbid the generation of neutrino mass at a lower order. Such a symmetry can stabilize the lightest neutral member of all particles that transform nontrivially under the symmetry, and makes it a natural DM candidate.

JHEP05(2016)030
• New DM annihilation processes such as semi-annihilation (SE-A) [18] become available that allow for different numbers of DM particles to appear in the initial and final states. The processes can change significantly the evaluation of the DM relic density.
• Multi-component DM is possible. In this case, annihilation processes between different components lead to the so-called assisted freeze-out mechanism [28], which also influences the DM relic density.
In this paper, we focus on the first two features. We consider a specific Z 3 DM model that induces neutrino mass at two loops. The model was originally proposed in ref. [4], in which DM can be either a Dirac fermion or a complex scalar. Some phenomenological aspects of the model have been previously studied in ref. [29] with emphasis on the effect of SE-A processes on relic density and direct detection. Here we aim to implement a comprehensive analysis on DM properties and collider signatures. We will show that both SE-A and co-annihilation (CO-A) processes have significant effects on the evaluation of relic density while evading stringent constraints from direct detection. Moreover, the presence of many new decay channels of Z 3 particles induces a plenty of distinct signals at LHC for both scalar and fermion DM candidates which would be absent for Z 2 DM.
The rest of this paper is organized as follows. In section 2, we recall the model and discuss current experimental constraints on its parameters. Sections 3 and 4 contain the core content of this work, in which we systematically study DM properties and LHC signatures for both scalar and fermion candidates. In section 3, we explore the vast parameter space that survives the constraints from relic density and direct detection; all important annihilation channels will be presented and discussed in detail. In the following section 4, we first exhaust all decay patterns according to the mass spectra of new particles, and then analyze various LHC signatures and compare with the relevant LHC limits. Finally, section 5 is devoted to conclusions.

Model and constraints
In the model under consideration [4] a global and exact Z 3 symmetry is imposed to induce neutrino masses at the two-loop level through interactions with new particles charged under Z 3 . In one minimal version of the model, one introduces two scalars χ a (a = 1, 2) and one Dirac fermion S, both of which are neutral singlets of the SM gauge group, and one Dirac fermion doublet Ψ = (N, E) of hypercharge Y = −1. The fermions are assumed to be vector-like to avoid chiral anomalies. The new particles transform under Z 3 in the same way as χ a → χ a ω with ω = exp (i2π/3), while SM particles are neutral. The Yukawa and fermion mass terms involving new fields and the SM leptons F iL = (ν iL , iL ), iR and Higgs boson Φ = (G + , φ 0 ) are: whereΦ = iσ 2 Φ * . And the scalar potential is where M 2 and λ 3 are Hermitian, µ abc is complex and symmetric in the three indices, and λ ab;cd 2 = λ cd;ab 2 = (λ ba;dc 2 ) * is complex as well. Some of the phases in the above couplings can be removed by field redefinitions but there are still many physical ones. To make the number of independent parameters under control for our later numerical analysis, we will simply assume that λ ab;cd 2 = λ 2 , λ ab 3 = λ 3 , µ abc = µ, and that they are all real. There are some theoretical considerations that can be used to set a bound on the parameters in the scalar potential, such as perturbativity, unitarity [30], and Z 3 not to be spontaneously broken [22]. These constraints are easily respected in our numerical analysis. Since Z 3 is exact, new particles do not mix with SM particles but can mix among themselves. We assume that χ 1,2 are diagonalized by an angle α to χ L,H of masses M χ L ≤ M χ H . The electrically neutral fermions S and N also mix due to the Yukawa couplings z L,R by an angle β to N 1,2 of masses M N 1,2 . Our convention is that N 1 (N 2 ) is dominantly a singlet S (doublet N ) for small β but either mass order is possible. In terms of the mass-eigenstate fields the couplings will involve the mixing angles. For the Yukawa couplings, we simply replace the primed couplings by unprimed ones, e.g., x a L,R and h ai . For the scalar self-couplings, the angle α enters explicitly; e.g., the χ a χ b χ c coupling (now a, b, c = L, H) is proportional to µg abc , in which an index L (H) is associated with a factor of α − = cos α − sin α (α + = cos α + sin α), for instance, Therefore we can take M χ L,H , m N 1,2 , α, β, h ai , x a L,R , λ 2 , λ 3 , and µ as our input parameters.
The above interactions induce neutrino masses at the two-loop level [4] as shown in figure 1: where a, b, c = L, H refer to scalars χ L,H and k, l = 1, 2 to fermions N 1,2 . The loop functions are , , Our loop functions agree with ref. [16] which shares the same topology of Feynman diagrams in a different model, but the relative sign of the two terms differs from that in ref. [29] which computes neutrino mass from a Feynman diagram of same topology in another scenario of the Z 3 model. The induced 3 × 3 neutrino mass matrix has a degenerate structure implying a massless neutrino in either normal or inverted hierarchy. With a two-loop suppression factor of (4π) −4 ∼ 10 −5 , it is easy to accommodate a mass of order 0.1 eV for the other two neutrinos by assuming reasonable values of new parameters. For heavy masses of same order, the loop functions are of order 0.1. As will be shown below, the constraints from lepton flavor violating (LFV) transitions can be trivially fulfilled with h ai ∼ 0.01. This then requires µx a L,R sin 2 (2β) ∼ 0.1 GeV. As is well known, precise measurements of LFV transitions set strong constraints on relevant interactions. The diagram in figure 2 for the lepton radiative decay j → i γ yields the branching ratio [31,32]: where the loop function F (x) is given by (2.8) Currently, the most stringent limit comes from BR(µ → eγ) < 5.7 × 10 −13 (90% C.L.) [33], while the limits on τ decays are less stringent, BR(τ → eγ) < 3.3×10 −8 (90% C.L.) [ (2.9) For instance, when M E ∼ 200 GeV, the above become |h ae * h aµ | 0.0002, |h ae * h aτ | 0.12, and |h aµ * h aτ | 0.14. Thus, without requiring a special flavor structure we can assume safely a universal bound |h ai | 0.01. 1 The mass splitting between the charged fermion E and neutral ones N 1,2 will contribute to the custodial symmetry breaking measured by the parameter ∆T . Using the formulas in, e.g., refs. [36,37] and the fitting result |∆T | < 0.2 [38], we can set a bound on the mass splitting. In the most stringent case for a large mixing sin β ∼ 0.66, one gets |M N 2 −M N 1 | < 250 GeV. In the opposite case of a small mixing sin β < 0.1, |M N 2 − M E | is restricted to be less than a few GeV [39]. In our numerical analysis we will always work with a small β and assume N 2 , E are degenerate. Furthermore, for a light DM particle χ L or N 1 there are collider constraints on invisible Higgs decays, which will be examined in the next section.

Dark matter phenomenology
In this section, we investigate DM phenomenology of the Z 3 model. For this purpose, we generate the CalcHEP [40] file by using the FeynRules [41,42] package, which is used by micrOMEGAs4.1 [43] to calculate the DM relic density and DM-nucleon scattering cross section. We implement random scans over a vast parameter space (with a total of 3 × 10 5 samples), for the ranges or values of the input parameters shown in table 1. The constraints from relic density and direct detection are then imposed on each sample. For relic density, we use the combined Planck+WP+highL+BAO 2σ range, 0.1153 < Ω DM h 2 < 0.1221 [44]. As for direct detection, we adopt the currently most restrictive spin-independent limit provided by the LUX experiment [45]. 2 We make some comments before presenting our numerical results: • There are four new neutral particles in this model, χ L,H and N 1,2 , resulting in a rich annihilation pattern. In addition to the standard annihilation (ST-A) processes, the SE-A and CO-A processes play a crucial role in the DM relic density. The latter two processes are sensitive to the mass relations of new particles and should be thoroughly examined. We will discuss this issue in great detail in section 3.1. • For DM-nucleon scattering, it is sufficient to include tree-level contributions since oneloop terms are subleading. In the case of scalar DM, χ L interacts with quarks through the Higgs-portal λ 3 term, so that direct detection can set a stringent constraint on χ L χ † L → bb annihilation channel. This is a common feature of Higgs-portal models as we will discuss in section 3.2. For fermion DM, N 1 can scatter with quarks via the t-channel Z exchange due to the singlet-doublet mixing. Direct detection then imposes a very stringent bound on the mixing angle β. The additional t-channel Higgs exchange also contributes via the z L,R Yukawa terms but is subleading to the Z exchange. The angle β is therefore varied in a much narrower interval than α. The constraints from the Higgs and Z invisible decays and electroweak precision measurements also prefer a small β.

JHEP05(2016)030
• Since the quartic coupling λ 2 is only related to the DM self-interaction and has no further phenomenological effect, we assume a fixed value for it in the scan.
• According to our discussion in section 2, we assume a universal |h ai | 0.01 to avoid dangerous LFV processes, though a relatively large |h aτ | can be accommodated by a specific flavor structure.

Analysis of parameter space
We present our numerical scans in this subsection. Figure  • It is clear from the left panel that both scalar χ L and fermion N 1 DM samples can survive in the whole mass regions that we explored, but the scalar one has a much more number. This can be explained as follows. Since both CO-A and SE-A processes depend on specific mass relations, the survived samples resulting from them do not distinguish much between χ L and N 1 . The difference originates instead from ST-A processes. While χ L annihilates into gauge boson pairs and produces the correct relic density in a vast mass region, N 1 mostly annihilates into light and b quark pairs through the Higgs exchange, which can only give the correct relic density in a relatively low mass region.
• The right panel shows that survived samples tend to cluster in regions of small mass splitting for RD+LUX points, where CO-A and SE-A processes are generally dominant. Since ST-A channels are not sensitive to mass splitting and tend to cause a more scattered distribution in the plane, the clustering indicates that direct detection imposes a more severe constraint on ST-A processes. We will illustrate this feature explicitly in section 3.2.
In order to investigate the parameter space more comprehensively, the distributions of The direct searches for invisible Higgs decays have been carried out by ATLAS and CMS in the weak boson fusion (WBF) [46,47] and Zh associated production channels [47,48], with the 95% CL upper bounds on BR(h → invisible) of 28%(ATLAS), 65%(CMS) in the WBF channel and 75%(ATLAS), 83%(CMS) in the Zh channel, respectively. Alternatively, invisible Higgs decays also get constrained by fitting to visible Higgs decays [49][50][51].  The upper bound thus found is stronger, BR(h → invisible) < 25% (95% CL) [51], which will be employed in our following discussion. The decay width to scalar or fermion DM reads  with v ≈ 246 GeV. We take Γ vis = 4.07 MeV for the visible decay width at M h = 125 GeV and eq. (3.1) for the invisible one Γ inv . The scatter plot of the invisible decay branching fraction BR inv = Γ inv /(Γ vis +Γ inv ) is presented in figure 6 as a function of M DM for RD and RD+LUX survived samples, where the shaded area indicates the region excluded by the upper bound from ref. [51]. We found that for χ L DM, samples with M χ L < 55 GeV are totally excluded while for N 1 DM the corresponding bound can be relaxed to about 28 GeV.
One can also convert the upper bound on invisible Higgs decays into excluded regions in the [M χ L , λ hχ L ] ([M N 1 , β]) plane for χ L (N 1 ) DM. As shown in figure 7, the correlations among parameters manifest themselves more explicitly. In this manner, we obtain the most stringent bound λ hχ L 0.01 with M χ L < 55 GeV for χ L DM, or β 4 • with M N 2 − M N 1 = 500 GeV for N 1 DM in the most stringent case. We notice that these constraints are less stringent than from direct detection in the same mass regions, such that all of RD+LUX samples easily survive for either χ L or N 1 DM. Finally, N 1 DM also contributes to the invisible Z decay if kinematically allowed, The LEP 95% CL upper bound of 3 MeV [52] translates to β 17 • , which is weaker than from invisible Higgs decays. For light N 2 , the decays h → N 1N2 /N 2N1 , N 2N2 and Z → N 1N2 /N 2N1 , N 2N2 may also be allowed. These decay modes could provide more severe constraints [39], but are still weaker than the LUX constraint. Therefore we will not consider invisible decays as an individual constraint in our later analysis.
GeV is fixed, yielding the most stringent limit.

Analysis of annihilation channels
The aim of this subsection is to demonstrate the effects of various annihilation channels on relic density and direct detection, especially the crucial roles played by CO-A and SE-A processes. For this purpose, we first list all SE-A processes for both χ L and N 1 DM. As seen in figure  For χ L DM, we have the following observations: • As seen in figure 9, χ L DM has three ST-A channels. For light DM (M χ L < M h /2), it dominantly annihilates into bb, while for heavy DM (M χ L > M W ), the dominant annihilation processes are into gauge boson and Higgs pairs. Since DM annihilating through the Higgs portal type always tends to produce more gauge boson than Higgs pairs, the majority of samples is from the W + W − channel with rare samples coming from the hh channel. Furthermore, SE-A (CO-A) processes occur only when but all of them take a small fraction.
• For light DM, since annihilation cross section for the bb channel is suppressed by the Yukawa coupling of b, one first requires a relatively large λ hχ L to saturate relic density. As M χ L approaches M h /2, resonance enhancement and phase space suppression compete. Since the former dominates, the overall effect is to require a decline in λ hχ L . After M χ L climbs over the h resonance, the opposite takes place, resulting in the valley structure in the left panel of figure 10. This is indeed a common feature of Higgs-portal models. On the other hand, for heavy DM, the annihilation cross sections for the W + W − and hh channels are respectively proportional to the gauge coupling and Higgs self-coupling, so that relic density selects a narrow band in the [M χ L , λ hχ L ] plane.
• Upon imposing the direct detection constraint, most samples with the bb channel are excluded since λ hχ L as required by relic density is too large to evade the LUX bound for such light DM. The only exception is a DM mass near the resonance area, where a few samples survive due to a much smaller λ hχ L . In contrast, most of samples with SE-A and CO-A processes are safe in this case. This feature is mainly because, when relic density is determined by these two processes, a smaller λ hχ L is still allowed for the same order of DM mass, therefore alleviating the tension from direct detection. For instance, the mass interval M χ L ∈ [80,350] GeV is excluded by the LUX limit for the ST-A channel χ L χ † L → W + W − alone, but is allowed when the SE-A and CO-A channels χ L χ L → χ † L h and χ L χ † H → W + W − are taken into account. When M χ L > 350 GeV, all above channels satisfy the LUX bound, but SE-A and CO-A channels still keep smaller couplings.   For N 1 DM, we observe the following features: • Compared with χ L DM, N 1 DM has a more complicated annihilation pattern since more particles are involved in annihilation processes. As shown in figure 11, there are two ST-A channels in the RD survived samples, N 1N1 → dd, bb, both dominating for • Including the LUX limit, there are only five SE-A/CO-A annihilation channels that survive the combined RD+LUX constraints: figure 12. This is due to the similar reason as for χ L DM, i.e., they benefit from smaller couplings compared with ST-A channels, which breaks the tight correlation between relic density and direct detection.

Comment on gamma-ray excess from the Galactic Center
While a complete analysis on indirect detection constraints is beyond the scope of this paper, we discuss briefly one of the most interesting anomalies in DM searches, namely a possible gamma-ray excess from the Galactic Center (GCE). The excess has been reported by a series of theoretical analyses using public data from the Fermi-LAT since 2009 [53][54][55][56][57][58][59].  Very recently, the Fermi collaboration has also released their own analysis [60]. This has attracted great attention in both astrophysics and particle physics communities. When using a ST-A process to interpret the excess, the spectrum is best fit by bb final states for a DM mass of 30 − 50 GeV with σv bb ∈ [1.4, 2] × 10 −26 cm 3 s −1 [53], and the morphology of DM density distribution matches the canonical Navarro-Frenk-White (NFW) halo profile. The τ + τ − , qq and cc (gg, W + W − , ZZ, hh and tt ) final states with a lighter (heavier) DM mass and a slightly different annihilation cross section are also acceptable [61][62][63][64]. Furthermore, it does not conflict with current limits from dwarf spheroidal, antiproton and CMB observations when taking into account uncertainties in the DM halo profile and propagation model [64][65][66][67]. As usual, the excess can also be incorporated by astrophysical phenomena, including millisecond pulsars or unresolved gamma-ray point sources [68][69][70][71][72]. However, astrophysical interpretations encounter some challenges on matching the spectrum and morphology of the excess. In any case, GCE has triggered extensive model building studies in both general and specific frameworks . These models can be divided into two scenarios from a model-independent viewpoint: one-step direct annihilation and multi-step    cascade annihilation [63,76,77]. In the first scenario, DM annihilates directly into SM final states, so that its mass and cross sections are tightly bounded with the resulting photon spectrum. More critically, this scenario usually suffers from stringent constraints from direct detection and collider searches on DM or exchanged particles. On the other hand, in the second scenario, DM annihilates into lighter mediators which subsequently decay to SM particles. Since cascade decays modify observed signals of DM annihilation, shift SM final states (and thus the resulting photons) to lower energies and broaden their spectra, the corresponding parameter space will be considerably extended and could evade bounds from direct detection.
In the Z 3 model under study, ST-A channels also face the same difficulty mentioned above. In figure 13 we plot the distribution of σv bb for all survived samples, except for RD+LUX samples in the case of N 1 DM, which are entirely excluded by the LUX constraint. We see that the parameter region consistent with GCE is excluded by direct detection. In order to avoid this conflict, some recent papers proposed a class of DM models with a local Z 3 symmetry [83,110,111], which often arises as a remnant of a spontaneously broken hidden gauge symmetry. The GCE may then be explained using SE-   A channels associated with new Higgs/gauge bosons. More interestingly, as pointed out in refs. [108,109], singlet models with a global Z 3 symmetry can also fit the GCE signal when taking into account SE-A contributions properly. In such models, DM candidates can be either a scalar or a two-component scalar and fermion. It has been shown that the GCE signal can be accommodated in either case when the DM mass is close to the SM Higgs so that the produced single Higgs through SE-A processes is nearly at rest. This mechanism also works for the Z 3 model under consideration, and the relevant SE-A channels correspond to χ L χ L → χ † L h and N 1 χ L →N 1 h. However, since the model content here is richer, the parameter space required by GCE could be very different. A comprehensive and highly efficient analysis of this issue would employ the MCMC method, which we hope to address in the future.

LHC phenomenology
In carrying out the LHC study of the Z 3 model, we use MadGraph5 aMC@NLO [112,113] to calculate the production cross sections of Z 3 particles with CTEQ6L1 [114] parton distribution functions (PDFs). The leading contributions under consideration are the pair and associated production of the doublet fermions via the s-channel Drell-Yan process: The total cross sections of these processes are plotted in figure 14 as a function of the mass M Ψ , where an overall K-factor of 1.2 is applied to both 13 TeV and 14 TeV cases [115]. the typical production rate of electroweakinos (charginos and neutralinos) in the minimal supersymmetric standard model (MSSM) [116].

Decay properties
To prepare for the study of LHC signatures, we discuss in this subsection the decay properties of Z 3 particles. Figure 15 shows the decay patterns for all nine cases allowed by DM considerations, and the decay branching ratios for all Z 3 particles are shown in figures 16-20. We assume E and N 2 to be degenerate to reduce the number of parameters, which corresponds to a degenerate fermion doublet in the limit of no mixing. These decay patterns not only affect the DM properties discussed in the previous section, but also have a great impact on the LHC signatures. Cases AI-III correspond to fermion DM, for which we only consider the case M N 1 < M N 2 due to severe constraints on the opposite mass order from direct detection [117]. Cases BI-CIII correspond to scalar DM with M N 1 < M N 2 in cases BI-BIII or M N 1 > M N 2 in cases CI-CIII. From these decay patterns, we see clearly that the decays of the fermion doublet are heavily dependent on the mass spectrum of the Z 3 particles. Thus in the following studies for each Z 3 particle, we will choose several mass spectra to illustrate such an impact. The decay channels can be classified into three categories according to interactions via (1) the gauge coupling (e.g., E − → W − N 1 ), (2) the Yukawa coupling (e.g., E − → χ L − ), or (3) the scalar self-coupling (e.g., χ H → χ L h). Decays like E − → W − N 1 in category (1) are possible due to the mixing between the singlet and doublet neutral fermions determined by the angle β. As mentioned previously, in the case of fermion DM, β is tightly constrained by direct detection, while in the case of scalar DM a large β is still allowed. For simplicity, we will choose a small β in both cases in our following discussion and other relevant parameters, as follows: We will take several sets of Z 3 particle masses to illustrate different decay patterns. We first discuss the decays of the heavy charged fermion E − . There are three decay channels: The branching ratios of E − as a function of M E are presented in figure 16 for three cases of Z 3 particle spectra. Case (a) corresponds to fermion DM, while cases (b) and (c) correspond to scalar DM. In case (a), the decay channel E − → W − N 1 is dominant in the whole mass region. BR(E − → χ L − ) reaches maximum 0.1 around M E = 400 GeV, while BR(E − → χ H − ) is a little bit smaller due to phase space suppression. In cases Because of the mixing between the neutral fermions, N 2 has more decay channels than E: In addition to N 2 → ZN 1 , N 2 can decay into N 1 through emission of h, χ L , and χ H . More interestingly, the decay N 2 → χ L ν is totally invisible at colliders in the case of scalar DM, which will intensively contribute to the signature of mono-jet, -γ, and -Z [118]. For the same cases of Z 3 spectra as in the discussion of E − , the branching ratios of N 2 as a function  of M N 2 are plotted in figure 17. In case (a) for fermion DM, N 2 → ZN 1 is dominant before N 2 → hN 1 is opened, and BR(N 2 → ZN 1 ) ≈ BR(N 2 → hN 1 ) ≈ 0.5 soon after the latter is opened. The branching rations of other decay channels are always smaller than 0.1. With the choice of h ai = sin β = 0.01, we have approximately BR(N 2 → χ L,H ν) ≈ BR(N 2 → χ † L,HN 1 ), for all three sets of masses. In case (b) for scalar DM, N 2 decays dominantly into χ L ν in the low mass region below 400 GeV, and into ZN 1 /hN 1 in the high mass region above 600 GeV. In the intermediate mass region around 500 GeV, the four decay channels N 2 → χ L ν, χ † LN 1 , ZN 1 , and hN 1 are comparable with each other. In case (c), with M χ H lighter than M N 1 , BR(N 2 → χ H ν) could reach over 0.3 before ZN 1 is open.
Although the direct production rates of N 1 , χ L , and χ H are small at colliders, they can be produced via the cascade decays of E, N 2 and subsequent decays into lighter particles. Possible promising signatures might occur in certain cascade decay chains, thus we also present the decay channels of these singlet particles for completeness. We first discuss the decays of N 1 , which happen only in the case of scalar DM: In figure 18, the branching ratios of N 1 is displayed as a function of M N 1 for three cases of Z 3 particle spectra. The decay N 1 → χ L ν is totally dominant in the low mass region In the high mass region where all channels are opened, the three channels N 1 → W + E − , ZN 2 , hN 2 dominate, and have the approximate relations, (4.6) due to the Goldstone nature of W, Z [119]. In contrast to N 1 , χ L can only decay in the case of fermion DM. Being only mediated by Yukawa couplings, the decay channels of χ L are: (4.7) In figure 19, we show the branching ratios of χ L as a function of M χ L for two sets of Z 3 particle masses. In the low mass region, the only allowed decay is χ L → N 1 ν. In the high mass region above 2M N 1 , the decay χ L →N 1N1 becomes dominant. Since the mixing angle β must be tiny in the case of fermion DM, the channels χ L →N 1N2 ,N 2N2 are always negligible. The decay channels χ L → N 2 ν, E − + depend heavily on the mass relations between M N 2 and M N 1 . For instance, if M N 1 < M N 2 < 2M N 1 as in case (a), both can be the main decay channels with an approximately equal branching ratio of ∼ 0.5 in the range between M N 2 and 2M N 1 . On the other hand, if M N 2 > 2M N 1 as in case (b), neither of them dominates. The scalar self-interactions result in a richer decay pattern for the heavier scalar χ H than the lighter χ L : are illustrated in figure 20 for six cases of Z 3 particle spectra. Cases (a)-(c) correspond to scalar DM, and cases (d)-(f) to fermion DM. Similar to χ L , χ H →N 1N1 is the only dominant decay in the high mass region above 2M N 1 . But in the mass region below 2M N 1 , the decays of χ H can be quite different from χ L . For scalar DM, χ H decays into N 1 ν (cases (a) and (b)) or N 2 ν/E − + (case (c)) in the low mass region, depending on which of N 1 and N 2 is lighter. In the intermediate region between M χ L + M h and 2M N 1 , the cascade decay χ H → χ L h dominates. Such decay channels play a very important role in the detection of scalar interactions at collider. And once allowed, the branching ratio of χ H → χ † L χ † L could reach about 0.2, which is the dominant invisible decay of χ H . Furthermore, for a large µ, e.g., µ = 100 GeV, the invisible decay χ H → χ † L χ † L is expected to be even larger than χ H → χ L h. For fermion DM, χ H can only decay as χ L into N 1 ν in the low mass region. Case (d) is most interesting among all three, where the four main decay channels This special case requires the mass relation M N 2 < M χ L + M h < 2M N 1 to be satisfied. If not, χ L h will be the main decay channel for a heavy N 2 as shown in case (e), or N 2 ν/E − + take over for a heavy χ L as shown in (f). If both N 2 and χ L are relatively heavy, χ H will decay the same way as χ L as shown in case (b) of figure 19.

LHC signatures and constraints
After the systematic study on the decay properties of the Z 3 particles in section 4.1, we now address their possible signatures at LHC. To a large extent, the LHC phenomenology is -21 -

JHEP05(2016)030
governed by the fermion doublet decays, since they can be pair or associated produced. The various decay channels of N 2 and E ± as well as the cascade decays of other Z 3 particles will lead to characteristic collider signatures. For instance, the final states of Z 3 particles will always have missing transverse energy E T due to the existence of DM. At the same time the most interesting and easiest way to detect signatures of neutrino mass models usually involves multi-lepton final states [120][121][122][123][124][125], and so is expected for the two-loop radiative neutrino mass model under consideration. Furthermore, with a Higgs boson h [126,127] in the decays of N 2 and χ H , it is also interesting to probe signatures with this h. Therefore we will explore the LHC signatures involving multi-( = e, µ) plus E T with or without a Higgs boson h. These signatures are naturally classified in terms of the number of leptons (up to four) in the final states.
Before delving into details, we make some comments on the choice of parameters in our collider analysis. We illustrate our results in figures 21-24 with the values of parameters that respect the constraints established in previous sections, and provide the DM results at a few benchmark points as shown in tables 3-5. In both cases of fermion and scalar DM, we choose the DM particle mass to be 150 GeV.
For fermion N 1 DM, we have to discriminate between two mass regions of N 2 (or simply Ψ since the singlet-doublet mixing is small) according to which SE-A channel of N 1 dominates. When M N 2 < 300 GeV, the channel N 1 N 1 → E + − dominates and is sensitive to the universal Yukawa couplings h ai . The other parameters are less relevant and can be given a fixed and save value. To be specific, we assume where we choose three typical values of h ai to illustrate their potential impact on the multisignatures. We show in table 3 the benchmark points corresponding to those values of h ai . One sees that a larger h ai requires a heavier N 2 to fulfil the relic density constraint since the cross section for N 1 N 1 → E + − increases monotonically with h ai . Note also that a larger universal h ai tends to touch the edge of the LFV bounds but this can be avoided simply by arranging a slight flavor structure in the eµ sector without modifying the multisignatures as a sum over = e, µ events.
When M N 2 > 300 GeV, the above channel is kinematically closed and another one, N 1 N 1 → χ † L h takes over. As we discussed in section 3.2, this channel is mediated by the exchange of an s-channel χ L,H or a t-channel N 1,2 . Fixing M χ H = 2M N 1 , the χ H exchange dominates and it is largely dependent on the couplings x a L,R while not sensitive to other parameters. Thus for M N 2 > 300 GeV, we assume instead  Table 3. Three benchmark points for N 1 DM with parameters given in eq. (4.10) when 0.34 0.1182 3.9 × 10 −11 Table 4. A benchmark point for N 1 DM with parameters given in eq. (4.11) when M N2 > 300 GeV.
For scalar χ L DM with M χ L = 150 GeV, the most effective annihilation channel to fulfil the DM constraints is χ L χ L → χ † L h. Since this channel is mediated by an s or a t-channel χ L,H , a relatively large µ and an appropriate λ 3 help to enhance its rate while other parameters can vary in a wide range. We choose A benchmark point is shown in table 5.

Signatures for N 1 DM
We first highlight the signatures appearing in the case of fermion DM. Since we concentrate on the multi-lepton signatures, we will consider the leptonic decays of the gauge bosons W and Z. The possible signatures are listed as follows.
(F1) 0 2h. This signature of no leptons and a pair of Higgs bosons [128,129] has a large E T , which can be used to suppress the SM background. The production mechanism is with h → bb/γγ. The same signature is also expected in supersymmetric (SUSY) and canonical seesaw models [130]. With BR(N 2 → hN 1 ) ≈ 0.5 in our benchmark scenario, the production rate of this signature is a quarter of σ(N 2N2 ). A search for this signature in the SUSY scenario has been performed by CMS [131] in gauge-mediated SUSY-breaking model -23 -

JHEP05(2016)030
where the lightest superparticle (LSP) is gravitino and the next-to-lightest superparticle is higgsino. For nearly massless LSP, there is no exclusion limit for N 2 up to 500 GeV if one matches N 2 − N 1 to the higgsino-gravitino system; and the sensitive mass region is M N 2 > 200 GeV for BR(N 2 → hN 1 ) > 0.5. However, for 14 TeV LHC with an integrated luminosity L = 3000 fb −1 , we may have a chance to probe this signature for a small production rate down to ∼ 0.1 fb or m N 2 up to 800 GeV [130].
(F2) 1 1h. This signature follows from the associated production of the doublet fermions: with h → bb/γγ, exactly as in the chargino-neutrolino system in SUSY models. If we further consider h → W W * → ± νqq , this channel can also produce the like-sign dilepton signature ± ± . Searches for this signature have been carried out by CMS [131,132] and ATLAS [133,134]. Again, matching E ± N 2 (N 1 ) toχ ± 1χ 0 2 (χ 0 1 ) and assuming BR(E ± → W ± N 1 ) ≈ BR(N 2 → hN 1 ) ≈ 100%, the limits on M N 2 have been set to 200 GeV by CMS [132] and 250 GeV by ATLAS [133] for massless N 1 . But as discussed in section 4.1, BR(E ± → W ± N 1 ) ≈ 100% and BR(N 2 → hN 1 ) ≈ 50% for the model considered here, one expects the limits on M N 2 to be relaxed.
(F3) 2 (\ Z). In this signature the two opposite-sign leptons are required not to reconstruct a Z boson. Such events are produced as with W ± → ± ν. In figure 21, the cross section of this signature at 13 (14) TeV LHC is presented for three typical values of h ai . For a fixed DM mass M N 1 = 150 GeV, the DM constraints are accommodated by different annihilation channels as discussed above according to whether M N 2 < 2M N 1 or opposite. We therefore have plotted the cross section using the parameters in eq. (4.10) for M N 2 < 300 GeV (solid curves) and in eq. (4.11) for M N 2 > 300 GeV (dashed ones). One sees that the cross section for a relatively light N 2 almost does not vary with h ai , but drops with the increase of h ai for M N 2 > 300 GeV. For instance, it is approximately two orders of magnitude smaller at h ai = 0.05 than at h ai = 0.01 in the latter case. Due to the large SM background from dibosons (W W ) and top quarks (mainly from tt and W t), constraints on this signature are relatively loose. The current LHC limits are sensitive to this signature in the mass region 100 GeV < M E < 180 GeV and M N 1 < 30 GeV which are based on 2 + E T searches of direct production of electroweakinos and sleptons [135]. But as discussed in section 3, such a low mass can hardly pass the constraints from DM. A brief discussion on such a signature with a much heavier N 1 at LHC has been performed in ref. [39].
(F4) 2 (Z)2j(Z). In this signature a pair of opposite-sign leptons is required to reconstruct a Z boson while a pair of jets is required to reconstruct a second Z boson. This ZZ signature comes from the decays of a neutral pair: with one Z → + − and the other Z → qq. There can also be fake contributions coming from W Z and hZ decays. Current LHC limits for this signature also come from direct searches of electroweakinos and sleptons [131,132]. Assuming M N 1 = 1 GeV, the most sensitive mass region is M N 2 > 250 GeV for BR(N 2 → ZN 1 ) > 0.5; and M N 2 < 380 GeV is likely excluded by CMS for BR(N 2 → ZN 1 ) ≈ 1 [131]. With a much heavier N 1 and BR(N 2 → ZN 1 ) = 0.5 in this model, the exclusion limits are considerably weakened.
(F5) 2 (Z)1h. This Zh signature also comes from the decays of a neutral pair N 2N2 : with h → bb/γγ. The production rate of this Zh signature is twice as large as ZZ above in our benchmark scenario for case (a) of figure 17. Analogously to previous signals, we found that most relevant LHC limits come from ref. [131]. The most sensitive mass region is 160 GeV < M N 2 < 430 GeV with BR(N 2 → hN 1 ) ∈ [0. 45, 0.85]. For the signature dominated by the bb channel, no exclusion limits are set due to large tt backgrounds.
(F6) 3 (Z). The production mechanism for this trilepton signature is The cross section for the 3 + E T signature at 13 (14) TeV LHC is shown in figure 22.
It is comparable with that the di-lepton signature in eq. (4.15) due to a relatively large production rate of E ± N 2 . But with a much cleaner background, this signature is expected to be the most promising one and to set the most stringent constraints in the mass region M N 2 250 GeV. Once again, current limits for this signature have been set by ATLAS [136] and CMS [132,137] from direct searches for electroweakinos. A recasting work [138] based on ATLAS limits has been performed in the gaugino-higgsino sector in MSSM with binolike DM and decoupled sfermions. We can transfer their recasting limits to our signal. Instead of M N 2 > 370 GeV set by ref. [136], recasting shows that the ATLAS limits are sensitive in the mass region M N 2 270 GeV and M N 1 75 GeV [138]. In addition, a combined analysis on the 2 and 3 signals by ATLAS [135] shows that M N 2 > 425 GeV. However, most of current limits are based on simplified models and can be significantly relaxed with different spectra, decay chains and branching ratios.
(F7) 4 (ZZ). This signature requires two pairs of opposite-sign dilepton to reconstruct the Z pair. It results from the process with both Z → + − . Although this signature is very clean, its production rate is suppressed due to the small leptonic branching ratio of Z. For this signal, the constraint from CMS searches [131] is less stringent than from the 2 2j signature discussed above. In summary, for fermion DM, the current most stringent LHC limit comes from the 3 signal resulting from W Z bosons. At upcoming LHC run II, other signatures such as 4b from hh, 2 2b from hZ, and 2 2j from ZZ are expected to have better sensitivity than this one in the high mass region. More importantly, noting the similarity of signals between fermion DM in the Z 3 model and electroweakinos/sleptons in SUSY models, it is very interesting to recast search limits on the latter to this scenario and examine their interplay with DM constraints. For this purpose, a detailed simulation and recasting is necessary using the tools already available [139][140][141]. We hope to come back to this in another work.

Signatures for χ L DM
Now we turn to the signatures related to scalar DM. A distinct decay mode of N 2 in this scenario is N 2 → νχ L , where both ν and χ L are invisible at colliders. This results in various mono-X (X = j, γ, W, Z, h, ) signatures at LHC. In what follows, we first discuss these mono-X signatures, and then the signatures of multi-leptons plus E T with or without h. (S1) 1j. This mono-jet signature is extensively studied in DM searches at LHC. It proceeds as pp → N 2N2 + j → νχ L νχ † L + j, (4.20) and in the low mass region M χ L < M h /2, the following signal channel should also be considered: The second process depends on the coupling λ hχ L , and according to ref. [142], the 14 TeV LHC with 300 fb −1 luminosity has the ability to probe |λ hχ L | > 8 × 10 −3 . The mono-jet searches by both CMS [143] and ATLAS [144] are based on the effective field theory approach to weakly interacting massive particles of DM, where only the DM pair contributes to the signature E T . Differently from this, the signature E T in this Z 3 model is also contributed by the neutrino pair as shown in eq. (4.20). Since N 2N2 can be copiously produced through the Drell-Yan process, we expect that there could be tight constraints from the mono-jet signature. Moreover, the mono-γ [145,146] and mono-W/Z [147,148] signatures are also possible at LHC. Although such signatures are less promising than mono-jet, they can be used as a diagnostic tool of the underlying models [118].
(S2) 1h. This is the so-called mono-h signature at LHC [149,150], which has attracted attention since the Higgs discovery [126,127]. The signature arises from L , when χ L,H are both lighter than N 2 . Searches for the signature have been recently published by ATLAS in the h → γγ [151] and h → bb [152] channel. The upper limit on the cross section is 0.7 fb for γγ and 3.6 fb for bb with E T > 90 GeV and E T > 150 GeV respectively. Similar to the mono-jet signature, this mono-h also has a pair of neutrinos contributing to E T . Since χ H must be 125 GeV heavier than χ L , BR(N 2 → νχ H ) should be always smaller than 0.5, but on the other hand BR(χ H → hχ L ) is totally dominant. Therefore, this signature is also promising.
(S3) 2h. This double Higgs plus E T signature is also produced in the case of scalar DM with both χ L . The searches by CMS [131] are also applicable here. Differently from the case of fermion DM, the h-pair now comes from the cascade decay of χ H and thus their sequential decay products bb/γγ are expected to be less energetic.
(S4) 1 . This signature can be regarded as a mono-with a large E T at LHC, and it arises from As shown in figure 14, the production rate of E ± N 2 is the largest at LHC. For both χ H and N 1 heavier than N 2 , E ± → ± χ ( †) L and N 2 → νχ L are totally dominant. The monosearch has been performed by CMS [153]. With both ν and χ L contributing to E T , we expect severe constraints on an electroweak-scale N 2 .

JHEP05(2016)030
(S5) 1 1h. This signature is quite similar to the W h signature in the fermion DM case. The production mechanism is L . The searches for the W h signature [131][132][133][134] can be applied to set a constraint on this signature as well. But differently from the fermion DM case, the branching ratios of E ± → ± χ ( †) L,H and N 2 → νχ L,H can be varied by tuning M χ H and the corresponding Yukawa couplings h ai .
(S6) 1 2h. This signature can only be produced in the case of scalar DM, and thus can be used to distinguish the character of DM at LHC. It follows from the process A similar signature has been studied in the context of type-II seesaw [154], where the lepton comes from an off-shell W . The additional and E T provides more efficient cuts than the pure Higgs pair to suppress the background, hence this signature is within the reach of LHC for a light N 2 [154].
(S7) 2 (\ Z). Differently from the fermion DM case, the lepton pair is produced from direct decays of E ± , 27) and is expected to be much more energetic for a large mass splitting between N 2 and χ L than from the W pair in the fermion DM case. This will lead to a more stringent constraint at colliders. The cross section at 13 (14) TeV is depicted in figure 23 for a universal Yukawa coupling h ai , so that BR(E ± → e ± χ L,H ) = BR(E ± → µ ± χ L,H ) = BR(E ± → τ ± χ L,H ). Contrary to the fermion DM case, the cross section now increases with h ai . The search for this signature by ATLAS [135] has excluded the mass of E ± between 160 GeV and 310 GeV with M χ L = 100 GeV for a simplified model.
(S8) 2 (\ Z)1h. Though sharing the same final state as the hZ signature in the case of fermion DM, the lepton pair is from the direct decays of E ± , L . As shown in case (c) of figure 16, BR(E ± → ± χ † H ) can reach over 0.3, so the production rate for this signature could be promising. Since the h → bb channel suffers from quite large background, we expect the h → W W/ZZ/γγ/τ + τ − channels to enhance the observability. (S9) 2 (\ Z)2h. As far as we know, the + − hh + E T signature has been seldom studied in previous papers. To have a pair of h in the final state, we require two χ H s in the cascade which further cascade decay as χ L . Since BR(E → χ H ) ≈ 0.3 and BR(χ H → hχ L ) ≈ 0.8 − 1, the cross section of this signature is roughly one-tenth of σ(E + E − ). On the other hand, the backgrounds such as ZZhh, W W hh, ttjj, etc., are relatively small. So this signature may also be promising at LHC.
(S10) 3 (\ Z). The trilepton signature is also possible in the case of scalar DM following the production pp → E ± N 2 → ± χ L + νχ H , ± χ H + νχ L , (4.30) and decays χ ( †) L mediated by an off-shell E ± . To have a relatively large branching ratio for the decays, the mass splitting between χ H and χ L must be less than M h . The theoretical cross section for the signature is plotted in figure 24, and it can be about ten times larger than that from W Z in eq. (4.18) for the fermion DM case. The same final state has been searched for by CMS [132] and ATLAS [136] for sleptons lighter than charginos and neutrolinos, with an exclusion limit on M N 2 up to about 700 GeV. But these constraints cannot be applied directly to the signature here, mainly because of the softness of the dilepton from χ (S11) 4 (\ Z). There are two processes contributing to this signature L as well. The first process has one pair of energetic leptons from the direct decays of E ± , while all leptons in the second are expected to be soft. The search for this signature has been carried out by ATLAS based on the simplified versions of R-parityconserving, R-parity-violating, and general gauge-mediated SUSY breaking models [155]. With appropriate matching of particles and decay chains, we obtain that M N 2 < 600 GeV with M χ L < 100 GeV has been excluded by the direct search [155]. For M N 2 < 500 GeV, the exclusion limit on M χ L of this 4-lepton signature is comparable with that of the trilepton signature. But for the same reason as discussed for the trilepton signature, the constraint cannot be taken for granted before a detailed simulation is performed.
To summarize the case of scalar DM, the most stringent constraint is also expected to come from the 3 signature. More interestingly, we find that various mono-X (X = j, γ, W, Z, h, ) signatures appear in this case, and differently from the current searches [143-148, 151, 152], missing transverse energy involves both scalar DM χ L and neutrinos. A detailed simulation and recasting of these mono-X and multi-signatures with or without h is necessary to clarify the feasibility of testing the scalar DM scenario at LHC.
Before ending this section, we briefly discuss how to distinguish between the collider signatures of Z 2 and Z 3 DM models. Based on the method developed in refs. [25][26][27], the two symmetries can be potentially discriminated by using multiple kinematical edges and cusps. The basic idea is that the cascade decay of a Z 3 particle can result in two visible particles that are separated by a DM particle. Such a decay chain involves a triple coupling of Z 3 particles which is absent in the Z 2 case. But in the minimal case with only two Z 3 scalars (χ H,L ), the desired decay chain is hard to realize. For that purpose, we may introduce a third scalar χ 3 . Then a concrete example would be the decay chain, E − → − χ 3 → − χ † L χ † H → − χ † L hχ † L , assuming the mass hierarchy M E > M χ 3 > M χ H > M χ L and suitable mass splitting. The charged lepton and Higgs boson h are then separated by the DM particle χ † L , which then results in a cusp in the distributions of the kinematical variables M h (invariant mass of the h system) and M 2 h (invariant mass squared) [25].

Conclusion
We have made a comprehensive analysis on the phenomenology of a Z 3 DM model that generates neutrino mass at two loops. We have examined in great detail its properties in relic density, direct detection and LHC signatures. For indirect detection, we also briefly discussed the GCE issue. To conclude, we summarize the key features separately for the scalar and fermion DM as follows.
For the scalar χ L DM, there are three ST-A channels χ L χ † L → bb, W + W − , hh, and three SE-A/CO-A channels χ L χ L → χ † L h and χ L χ † H , χ H χ † H → W + W − . The χ L χ † L → W + W − channel can satisfy both relic density and direct detection constraints in a vast mass region and thus gives the dominant contribution in the parameter space. Upon imposing the direct detection constraint, the bb channel is almost excluded while most of SE-A and CO-A processes survive. This is due to the fact that λ hχ L required by relic density is considerably relaxed for the SE-A/CO-A channels, thus alleviating the tension from direct detection. Concerning the LHC constraints, the 3 signal is expected to give the most stringent bound. Moreover, various mono-X (X = j, γ, W, Z, h, ) signatures are different from those in current LHC searches since missing transverse energy now involves both scalar DM and neutrinos. A detailed simulation and recasting of these mono-X and multi-signatures with or without h will be helpful to test the scalar DM scenario at LHC.
If the lighter of neutral fermions (N 1 ) plays the role of DM, the direct detection requires it to be an almost singlet with a mixing angle β < 2 • . Compared with χ L DM, it has more annihilation channels, including two ST-A channels N 1N1 → dd, bb and eight SE-A/CO-A channels N 1 χ L →N 1 h, E + W − ,N 1 Z; N 1 N 1 → χ † L h, E + − ,N 1 ν; N 2 E + → tb, ud. However, only five SE-A/CO-A channels (N 1 χ L →N 1 h; N 1 N 1 → E + − , χ † L h; N 2 E + → ud, tb) survive the LUX constraint, due to the same reason as for scalar DM. Interestingly, the LHC signatures of fermion DM are very similar to those of electroweakinos in simplified SUSY models. Currently, the 3 signal resulting from W Z bosons provides the most severe bound. At upcoming LHC run II, other signatures such as 4b from hh, 2 2b from hZ, and 2 2j from ZZ may be more promising in the high mass region. To make accurate estimation, it is necessary to recast current search limits on electroweakinos/sleptons and combine them with DM constraints. Finally, this model can also schematically explain the GCE observed by Fermi-LAT when taking into account contributions from SE-A processes for appropriate DM mass. The corresponding annihilation channels are χ L χ L → χ † L h for χ L DM and N 1 χ L →N 1 h for N 1 DM. A comprehensive analysis of this issue based on the MCMC method deserves a separate work.