Radiative Neutrino Mass with $Z_3$ Dark matter: From Relic Density to LHC Signatures

In this work we give a comprehensive analysis on the phenomenology of a specific $\mathbb{Z}_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 $\mathbb{Z}_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.


I. 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.
The simplest discrete symmetry is a Z 2 parity. However, if it appears as a remnant of a broken gauge group (U (1) X ), other Z N s are also possible in general [16]. The possibility with N > 2 has been investigated in Refs [4,[17][18][19][20][21][22]. Compared with the Z 2 case, DM with Z N symmetry has the following distinct features: • New DM annihilation processes such as semi-annihilation (SE-A) [17] 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.
• DM particles have new interesting decay modes that result in richer phenomenology and distinguishable signatures at colliders [23,24].
• Multi-component DM is possible. In this case, annihilation processes between different components lead to the so-called assisted freeze-out mechanism [25], 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. [26] 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 Sec. II, we recall the model and discuss current experimental constraints on its parameters. Sections III and IV contain the core content of this work, in which we systematically study DM properties and LHC signatures for both scalar and fermion candidates. In Sec. III, 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 Sec. IV, 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, Sec. V is devoted to conclusions.

II. 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 [27], and Z 3 not to be spontaneously broken [21]. 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.
Feynman diagram for neutrino mass.
The above interactions induce neutrino masses at the two-loop level [4] as shown in Fig. 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 , with Our loop functions agree with Ref. [15] 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. [26] 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 Fig. 2 for the lepton radiative decay j → i γ yields the branching ratio [28,29]: 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. [33,34] and the fitting result |∆T | < 0.2 [35], 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 [36]. 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.

III. DARK MATTER PHENOMENOLOGY
In this section, we investigate DM phenomenology of the Z 3 model. For this purpose, we generate the CalcHEP [37] file by using the FeynRules [38] package, which is used by micrOMEGAs4.1 [39] 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 I. 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 [40].
As for direct detection, we adopt the currently most restrictive spin-independent limit provided by the LUX experiment [41] 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 Sec. III A.
• For DM-nucleon scattering, it is sufficient to include tree-level contributions since one-loop 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 Sec. III B. For fermion DM, 1 We recall that in the case of fermion DM with a Z2 symmetry it is hard to provide correct relic density with such small Yukawa couplings [32]. 2 Since the exclusion limit varies with the DM particle mass MDM, we interpolate the LUX data to obtain the corresponding exclusion limit for each randomly generated MDM. 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 β.
• 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 Sec. II, 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. The SM Higgs has the mass M h = 125 GeV.

A. Analysis of parameter space
We present our numerical scans in this subsection. Figure  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 III B.   The direct searches for invisible Higgs decays have been carried out by ATLAS and CMS in the weak boson fusion (WBF) [42,43] and Zh associated production channels [43,44], 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 [45][46][47]. The upper bound thus found is stronger, BR(h → invisible) < 25% (95% CL) [47], 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. 10 for the invisible one Γ inv . The scatter plot of the invisible decay branching fraction is presented in Fig. 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. [47]. 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 Fig. 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 [48] 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 [36], 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.

B. 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 Fig. 8, a χ L pair can annihilate into χ † L,H h, N 1,2 ν and E + − final states via the sor t-channel exchange of χ † L,H . Similarly, an N 1 pair annihilates into χ † L,H h,N 1,2 ν and E + − final states via the exchange of an s-channel χ † L,H or a t-channel N 1,2 . Therefore, the s-channel annihilation may dominate when As we will show, the t-channel annihilation can also dominate in some regions. Moreover, CO-A processes are important in this model, which occur and even dominate in the case of  Table II. We can gain some useful information from the figures and table.
For χ L DM, we have the following observations: • As seen in Fig. 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 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 op-posite takes place, resulting in the valley structure in the left panel of Fig. 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 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 Fig. 11, there are two ST-A channels in the RD survived samples, N 1N1 → dd, bb, both dominating for M N 1 < M h /2. For SE-A processes, Finally, for CO-A processes, • Including the LUX limit, there are only five SE-A/CO-A annihilation channels that survive the combined RD+LUX constraints: Fig. 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.

C. 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 [49][50][51][52][53][54][55]. Very recently, the Fermi collaboration has also released their own analysis [56]. 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 [49], and the morphology of DM density distribution matches the canonical Navarro-Frenk-White (NFW) halo profile. The τ + τ − , qq and cc (gg, including millisecond pulsars or unresolved gamma-ray point sources [64][65][66][67][68]. 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 [59,72,73]. 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 Fig. 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 [79,105,106], which often arises as a remnant of a spontaneously broken hidden gauge symmetry. The GCE may then be explained using semi-annihilation channels associated with new Higgs/gauge bosons. More interestingly, as pointed out in Ref. [103,104], 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.

IV. LHC PHENOMENOLOGY
In carrying out the LHC study of the Z 3 model, we use MadGraph5 aMC@NLO [107] to calculate the production cross sections of Z 3 particles with CTEQ6L1 [108] 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 Fig. 14 [110].

A. 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 Figs. 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 [111]. 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 selfcoupling (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: sin α = 0.1, sin β = 0.01, λ 2 = λ 3 = 0.1, h ai = 0.01, x a L,R = 1, µ = 10 GeV.
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 Fig. 16 for three cases of  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 [112]. 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 Fig. 17. In case (a) for fermion DM, 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  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 Fig. 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 for all the three cases. In case (b), the decay 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, due to the Goldstone nature of W, Z [113].
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: In Fig. 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 The scalar self-interactions result in a richer decay pattern for the heavier scalar χ H than the lighter χ L : Among these, χ H →N 1N2 ,N 2N2 are severely suppressed by the tiny mixing angle β. Note that these decay channels can become relatively important in the case of scalar DM, where the mixing angle β could be larger. The branching ratios of χ H as a function of M χ H are illustrated in Fig. 20 for six cases of χ H → N 1 ν, N 2 ν/E − + , χ L h,N 1N1 become dominant sequentially as M χ H increases. 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 Fig. 19.

B. LHC signatures and Constraints
After the systematic study on the decay properties of the Z 3 particles in Sec. IV A, we now address their possible signatures at LHC. To a large extent, the LHC phenomenology is 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 involve multi-lepton final states [114][115][116][117][118], and so is expected for the two-loop radiative neutrino mass model under consideration. Furthermore, with a Higgs boson h [119,120] 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.

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 [121,122] 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 [123]. 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 [124] in gauge-mediated SUSY-breaking model 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 [123].
(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 [124,125] and ATLAS [126,127]. 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 [125] and 250 GeV by ATLAS [126] for massless N 1 . But as discussed in Sec. IV A, 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 Fig. 21, the cross section of this signature at 13 (14)  are based on 2 + E T searches of direct production of electroweakinos and sleptons [128]. But as discussed in Sec. III, 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. [36]. (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 [124,125]. 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 [124].
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 Fig. 17. Analogously to previous signals, we found that most relevant LHC limits come from Ref. [124]. The most sensitive mass region is 160 GeV < M N 2 < 430 GeV with BR(N 2 → (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 Fig. 22. It is comparable with that the di-lepton signature in Eq. (24) 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 [129] and CMS [125,130] from direct searches for electroweakinos. A recasting work [131] based on ATLAS limits has been performed in the gaugino-higgsino sector in MSSM with bino-like DM and decoupled sfermions. We can transfer their recasting limits to our signal. Instead of M N 2 > 370 GeV set by Ref. [129], recasting shows that the ATLAS limits are sensitive in the mass region M N 2 270 GeV and M N 1 75 GeV [131]. In addition, a combined analysis on the 2 and 3 signals by ATLAS [128] 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 [124] 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 [132][133][134]. 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 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. [135], the 14 TeV LHC with 300 fb −1 luminosity has the ability to probe |λ hχ L | < 6 × 10 −3 . The mono-jet searches by both CMS [136] and ATLAS [137] 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. 29. 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-γ [138,139] and mono-W/Z [140,141] 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 [112].
(S2) 1h This is the so-called mono-h signature at LHC [142,143], which has attracted attention since the Higgs discovery [119,120]. The signature arises from with χ L , when χ L,H are both lighter than N 2 . Searches for the signature have been recently published by ATLAS in the h → γγ [144] and h → bb [145] 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 [124] 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.
This signature can be regarded as a mono-with a large E T at LHC, and it arises from As shown in Fig. 14, the production rate of E ± N 2 is the largest at LHC. For both χ H and N 1 heavier than L and N 2 → νχ L are totally dominant. The mono-search has been performed by CMS [146]. With both ν and χ L contributing to E T , we expect severe constraints on an electroweak-scale N 2 .
(S5) 1 1h This signature is quite similar to the W h signature in the fermion DM case. The production mechanism is with χ L . The searches for the W h signature [124][125][126][127] 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 with both χ A similar signature has been studied in the context of type-II seesaw [147], 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 [147].
(S7) 2 (\ Z) Differently from the fermion DM case, the lepton pair is produced from direct decays of 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 Fig. 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 [128] 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 ± , with χ L . As shown in case (c) of Fig. 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 decays of E + E − , which further cascade decay as χ  (S10) 3 (\ Z) The trilepton signature is also possible in the case of scalar DM following the production 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 Fig. 24, and it can be about ten times larger than that from W Z in Eq. 27 for the fermion DM case. The same final state has been searched for by CMS [125] and ATLAS [129] 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 χ L . A recasting of it on the LHC searches would reveal a more realistic constraint. (S11) 4 (\ Z) There are two processes contributing to this signature with χ ( †) 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-parity-conserving, R-parity-violating, and general gauge-mediated SUSY breaking models [148]. 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 [148]. 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 [136-141, 144, 145], missing transverse energy involves both scalar DM χ L and neutrinos. A detailed simulation and recasting of these mono-X and multisignatures with or without h is necessary to clarify the feasibility of testing the scalar DM scenario at LHC.
Before ending up 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. [23,24], 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 (energy of the h system) and M 2 h (invariant mass squared) [23].

V. 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.