Singlet Majorana fermion dark matter: a comprehensive analysis in effective field theory

We explore a singlet Majorana fermion dark matter candidate using an effective field theory (EFT) framework, respecting the relations imposed by the standard model $SU(3)_C \times SU(2)_L \times U(1)_Y$ gauge invariance among different couplings. All operators of dimension-5 and dimension-6, forming a complete basis, are taken into account at the same time, keeping in view ultraviolet completions which can give rise to more than one operator at a time. If in addition CP-conservation is assumed, the remaining parameter space, where an EFT description is valid, is found to be rather restricted after imposing constraints from relic abundance, direct, indirect and collider searches. On including the CP-violating dimension-5 operator, $(\overline{\chi}i \gamma_5 \chi) (H^\dagger H)$, a significantly larger parameter space opens up. We use the profile likelihood method to map out the remaining landscape of such a DM scenario. The reach of future searches using ton-scale direct detection experiments, an $e^+ e^-$ collider like the proposed ILC and limits from future gamma-ray observations are also estimated.


Introduction
The evidence for non-baryonic dark matter (DM), which makes up more than 80% of the matter content in the Universe, is still entirely from its gravitational interactions and there is no convincing evidence so far of any couplings of the DM particle with the standard model (SM) sector [1,2]. Dedicated efforts for more than a decade in underground DM direct detection experiments as well as the search for DM pair-annihilation products in gamma rays, antiparticles and neutrinos have now considerably restricted the possible strength of such couplings. In fact, the dual requirements of obtaining a thermal relic abundance as required by the WMAP and Planck measurements, and a low enough spin-independent (SI) scattering rate with nuclei to be consistent with the impressive bounds from experiments like XENON100 and LUX, have cornered several possible models for particle DM.
From the theoretical viewpoint there are two possible roadmaps to explore weakly interacting massive particles (WIMP). One is to adopt a beyond SM (BSM) scenario which is motivated by some other particle physics considerations (for e.g., the naturalness problem of the electroweak scale) and can furnish a viable DM candidate.
The most well-known example in this class is the lightest neutralino in the R-parity conserving minimal supersymmetric extension of the SM. The alternative possibility is to study DM interactions with the SM sector in an effective field theory (EFT) set-up 1 . Needless to say, these two approaches are not completely independent. If only the DM particle is much lighter than the other new states in the model, at energies below the mass scale of the heavier states, any such BSM scenario is amenable to an EFT description. In this paper, we shall focus on the latter approach. To begin our discussion of an EFT for the DM particle, we first need to fix its spin and its quantum numbers under the SM gauge group of SU(3) C × SU(2) L × U(1) Y . We take the DM to be a spin-1/2 Majorana fermion, and a singlet under the SM gauge interactions.
The low-energy EFT of Majorana fermion DM has been studied on several occasions in the context of direct, indirect and collider experiments . The purpose of this paper is to improve upon the previous studies, update them in the light of recent data in all frontiers, and to perform a complete profile likelihood analysis of 1 In this study, an EFT is described by all the renormalizable and a complete set of higherdimensional operators (upto a given dimension) respecting the symmetries of the SM, with arbitrary low-energy co-efficients. This does not include matching to ultraviolet complete theories or running of the operator coefficients as in the original considerations of EFT's [3].
the EFT to determine its currently allowed parameter space. The directions in which our study goes beyond the previous approaches to the problem are the following: 1. We consider all operators of dimension-5 and dimension-6 allowed by the symmetries of the theory at the same time in our analysis. Not only is this justified in the spirit of an EFT, generically UV completions can lead to the presence of more than one operator at a time in the low energy effective theory. Therefore, we work with a complete basis of all dimension-5 and dimension-6 operators.
As is well-known, often there are redundancies among different higher dimensional operators which are related by the equations of motion (EOM). We have enumerated such operators as well, and have demonstrated how they can be eliminated with the help of the EOM.
2. A consequence of the above consideration is that given the DM mass, the relic density (Ωh 2 ) requirement does not fix the co-efficient of an operator, neither does the Ωh 2 < 0.1 condition put a lower bound on the same. All operators contribute additively to the annihilation cross-section at freeze-out, modulo interference effects which can be negative. These interference effects are entirely missed if we only consider one operator at a time. 3. We write down the effective Lagrangian by respecting the SM gauge invariance under SU(3) C × SU(2) L × U(1) Y . As a result, the relations between different DM couplings, as implied by the SM gauge invariance are taken into account. 4. A complete likelihood analysis is performed by including the requirements of relic density (from Planck [25]), direct detection (the recent LUX [26] and XENON100 [27] results), indirect detection (Fermi-LAT dwarf spheroidal observations [28] and IceCube [29] limits) and collider (Z-boson invisible width [30] and monophoton search limits from LEP [31] and Higgs invisible width and monojet+ missing energy search limits from LHC [32]). We also properly take into account the uncertainties coming from each of these measurements in our likelihood.
The remaining sections of the paper are organized as follows. In section 2, we describe the basis of dimension-5 and 6 operators considered by us, electroweak symmetry breaking (EWSB) effects and the range of validity of the EFT calculations. The experimental constraints from direct, indirect and collider searches are detailed in section 3, which also contains specifics about the construction of the likelihood function including the treatment of uncertainties. Section 4 is devoted to the role of individual operators in determining the thermal component of the DM relic abundance. Our results in the CP-conserving scenario are discussed in section 5, followed by the prospects of future experiments. The case including the CP-violating dimension-5 operator is discussed in section 6. Section 7 summarizes our findings. Finally, the construction of the likelihood function for direct detection experiments, and the validation of the LHC monojet simulation framework can be found in Appendices A and B respectively.

The effective field theory framework
In this section, we describe our EFT framework, in which the low-energy degrees of freedom consist of the SM particles and the Majorana fermion DM field χ. All interactions of the DM field with the SM sector are encoded by higher-dimensional operators, with a suppression scale Λ, which is the mass scale of the heavy fields integrated out to obtain the following low-energy Lagrangian: where, L SM is the renormalizable SM Lagrangian and L n represents higher-dimensional operators of dimension-n. In this study, we focus on the dimension-5 and dimension-6 operators only. In order to ensure the stability of the DM particle, we need to impose a Z 2 symmetry, under which the DM field is odd, and the SM fields are even. Therefore, any interaction term involving χ has to be at least bilinear in this field. Since this bilinear itself has a mass dimension of 3, the lowest dimension gauge-invariant operators that can be written down involving the interaction of a Majorana fermion DM and the SM sector are of dimension-5 and involve the Higgs doublet bilinear H † H: where, the CP-conserving operator is and the CP-violating one is given by As mentioned in the introduction, we shall first consider the CP-conserving case, and hence drop L CPV 5 from the dimension-5 Lagrangian. Subsequently, we shall discuss the consequences of adding the CP-violating operator.
After electroweak symmetry breaking, both the operators in Eqs. 3 and 4 lead to additional mass terms for the DM field. Moreover, the mass term coming from L CPV 5 is purely imaginary. As a result, we obtain a complex mass for the DM field: However, one can perform the following chiral rotation of the DM field to make the mass term real and positive 2 : where, Such a field redefinition keeps the kinetic as well as the dimension-6 current-current interactions of the same form, but can of course mix the scalar and pseudo-scalar interactions once CP-violation is allowed. However, since we keep both the co-efficients g S and g P S arbitrary, such an operator mixing is taken into account. Similarly, in our subsequent analysis, the DM mass is taken as a free parameter, m χ , as the bare mass term M χ is also a priori arbitrary.
The complete set of SU(3) C ×SU(2) L ×U(1) Y invariant dimension-6 four-fermion operators are the following: where, the sum over i runs over the three generations, and we have assumed flavourblind couplings to the SM fermions with the same gauge and global quantum numbers. Q i L , u i R , d i R , ℓ i L and e i R represent the left-handed (LH) quark doublet, righthanded (RH) up-type quark, RH down-type quark, LH lepton doublet and RH lepton fields respectively (the RH fields are singlets under SU(2) L ). At dimension-6, there is an additional operator involving the derivative of the Higgs field [18,20] L Higgs 2 While this manuscript was in preparation, similar observations were made in Ref. [19] After EWSB, apart from the χχhhZ and χχhZ couplings, L Higgs 6 leads to the following 3-point interaction of χ with the Z-boson: where, g is the SU(2) L gauge coupling, θ W is the Weinberg angle, and v = 246 GeV is the vacuum expectation value of the Higgs field. This term plays a crucial role in determining the low DM mass region in the m χ − Λ plane that satisfies the relic abundance requirement.
For a Majorana fermion DM, the vector current (χγ µ χ) and dipole moments (χσ µν χ, χσ µν γ 5 χ) vanish identically. While the so-called anapole moment term can exist [21], this operator can be written as a linear combination of the four-fermion operators in Eq. 8, by using the EOM of the gauge field (in this case, the photon field F µν ) as follows: where, Q f is the electric charge of f . Similarly, any operator involving the derivative of the DM field can be eliminated using the EOM of the DM field, which deviates from a free-field equation only with terms proportional to Λ −1 or lower.
Since we include only terms upto dimension-6 in the EFT, there is an implicit assumption that all operators are suppressed by a similar scale Λ, and therefore the contribution of dimension-7 terms in DM phenomenology is sub-leading. For the dimension-7 terms involving the SM fermions, for e.g., (χχ)(Q L Hd R ) this then readily justifies dropping these operators, since the dimension-6 operators in Eq. 8 as well as the Higgs-exchange-induced Yukawa couplings via the scalar operator in Eq. 3 should lead to much stronger interactions with the SM fermions. Although we also drop the dimension-7 couplings with the gauge field strength tensors G a µν , namely, χχG a µν G aµν (here the gauge field may belong to any of the SM gauge groups), the effective DM coupling to gluons via heavy quark loops is taken into account while considering the spin-independent direct detection rates. Apart from this case, the effect of these couplings is expected to be much smaller than the dimension-5 and 6 terms, since not only are they suppressed by 1/Λ 3 , the coupling constant is also of the order of 1/16π 2 , since the DM particle being a gauge singlet, can couple to gauge boson pairs only via loop diagrams.
Generically, an EFT description is valid as long as Λ >> m χ , and therefore, it is justified to integrate out the heavy fields with mass of the order of Λ 3 . However, in order to get a concrete idea about the minimum value of Λ to be considered while computing DM observables, let us take the example of a heavy particle X (scalar or vector) which can mediate in the s-channel annihilation of a DM-pair to a pair of SM particles. An EFT description in this case boils down to replacing by a constant mass scale the product of couplings of X to the DM-pair (g DM X ) and to the SM-pair (g SM X ), and its propagator denominator: In an weakly coupled underlying theory, g DM X ∼ g SM X ∼ O(1), and henceforth we shall assume this to be the case. Therefore, when matching the UV theory to the EFT, we obtain the relation Λ = m X for g DM X × g SM X = 1. Now, consider the pair annihilation rate of χ in the early universe, which is relevant for determining its current abundance. In this case, if v is the relative velocity between the DM particles, then the centre of mass energy squared is given by Therefore, the 2nd term in the propagator expansion in Eq. 13 now reads Ignoring the v 2 piece (which is smaller than the leading term by a factor of ∼ 0.025), and comparing this term to the leading term of −1/m 2 X , we observe that if we assume the minimal requirement of not producing X on-shell in a process, i.e., m X > 2m χ , then the 2nd term in the expansion in Eq. 13 is also of the order of −1/m 2 X , which is the same as the leading term. Therefore, in such a case, m X > 2m χ is a very poor approximation of the full theory. If on the other hand, we assume m X > 10m χ , the 2nd term in Eq. 13 can at most be −1/(25m 2 X ), and hence it contributes only 4% of the leading term (at the amplitude level). Therefore, at the level of crosssections, the error will be less than 8%. With such considerations, we find it justified to consider Λ ∼ m X > 10m χ . There can be modifications to this argument if the underlying theory is not weakly coupled or we have a heavy particle mediating in the t-channel etc. Keeping such modifications in mind, while presenting our results, we have separately indicated the regions where Λ > 10m χ and where 2m χ < Λ < 10m χ .
Apart from the computation of relic density, a careful choice of the scale Λ is also necessary for the high-energy collider experiments like LEP and LHC, in order for the EFT to be a valid description. We shall discuss the choice of such scales in Sec. 3. For the direct detection experiments, the momentum transfer in the relevant processes is in the MeV scale. Therefore, they do not lead to any further constraint on the region of validity for the EFT. Finally, as far as the indirect detection experiments are concerned, the conditions should be similar to the ones obtained above for the annihilation processes in the early universe. The only difference between DM-pair annihilation rates (σv) in the early and the present universe is the DM relative velocity v, which was already a small effect in our estimates above for the early universe.

Experimental constraints and the likelihood function
In this section, we provide the details of the experimental constraints employed in our analysis: originating from cosmological, astrophysical, laboratory-and colliderbased searches. We also briefly describe the profile-likelihood method used in our study and the various uncertainties in the observables that enter into the likelihood.

Profile likelihood method
In this study, we employ the profile-likelihood (PL) approach [33] to explore high probability regions of the multi-dimensional EFT parameter space. PL is a statistical method motivated in parts by both the Bayesian and the frequentist approaches.
It treats the unwanted parameters as nuisance parameters as in the Bayesian theory.
However, unlike in a Bayesian approach, in which one marginalizes over all unwanted parameters, the PL method takes the frequentist's concept of maximizing the likelihood along the directions one is profiling over. In other words, if a model has an n-dimensional parameter space, and we are only interested in p of those dimensions, then the PL is obtained by maximizing the likelihood over the (n − p) dimensions we are not interested in. Therefore, unlike in the marginal posterior in Bayesian theory, the prior does not contribute to the PL. However, it is very difficult to cover the full multi-dimensional parameter space with finite samples in a numerical scan.
Thanks to the advantage of prior-independence, it is nevertheless possible to combine several fine-grained scans focused on particular regions of the parameter space.
For example, even though in a region like m W < m χ < m t , the DM particle χ can achieve the required relic density via its pair-annihilation to the W + W − final state, this solution spans only a small volume in the whole hyperspace. A focused scan of such regions, therefore, becomes a necessity.
Our results will be primarily described in the relevant set of two-dimensional parameter regions which are in best agreement with all current experimental data, an example being the (m χ , Λ) space in 68% (1σ) and 95% (2σ) confidence intervals. After profiling over the rest of the parameters, one can write down confidence intervals in the (m χ , Λ) plane as an integral of the likelihood function L(m χ , Λ) where, the normalization in the denominator is the total probability with R → ∞ and R is the smallest area bound with a fraction ̺ of the total probability. If the likelihood can be modelled as a pure Gaussian distribution, the 68% (95%) confidence intervals in a two dimensional parameter space corresponds to −2 ln(L/L max ) = 2.30 (5.99), where L max is the maximum value of the likelihood in the region R.
Hereafter, for convenience, we introduce the variable χ 2 = −2 ln(L). We note that χ 2 is not exactly the same as in the usual chi-squared analysis unless the likelihood is described by a pure Gaussian.
The experimental constraints employed in our analysis, along with their central values, 1σ experimental uncertainties, and functional form of the likelihood functions are shown in Table 1. The details of most of the constraints and their likelihoods are provided in the relevant sections referred to in the Table. Our numerical scan of the parameters in the EFT span the following ranges: 10 GeV ≤ m χ ≤ 5 TeV We use a flat prior for all the operator co-efficients in the range |g i | < 1. For m χ and Λ, we combine both flat and log priors to obtain a maximal coverage of the whole  parameter space. As mentioned earlier, the co-efficients of the effective operators are taken to be |g i | < 1 , since we assume an weakly coupled UV completion to the EFT.

Relic abundance
From the relative heights of the acoustic peaks in the cosmic microwave background, the Planck experiment has measured the cold dark matter density to an accuracy of 3% [25]: If in addition, the WMAP polarization (WP) data at low multipoles is included, the above number on the 1σ error changes slightly: Although the difference between the two is not very significant, for definiteness, we use the value in Eq. 19 in our analysis. The likelihood function is taken as a Gaussian, and apart from the experimental error bar quoted above, an additional theoretical uncertainty of 10% in the computation of Ω χ h 2 has been assumed. The above number on the relic abundance is taken only as an upper bound, i.e., we demand Ω χ h 2 ≤ Ω c h 2 , such that the DM candidate χ does not overclose the universe.
However, we do not assume the existence of some other DM candidate making up for rest of the required relic abundance. Therefore, if for a parameter point Ω χ h 2 < Ω c h 2 as computed within the EFT, then a non-thermal production of χ should give rise to the additional required DM density. Such a non-thermal mechanism is not described by the EFT, but can exist in the UV completion. For example, in a supersymmetric theory, the late-time decay of gravitinos or moduli fields can produce a neutralino DM. Since such a gravitino or moduli field interacts very feebly with the DM field, it does not affect the DM phenomenology otherwise. To summarize, even though we accept parameter points which satisfy Ω χ h 2 ≤ Ω c h 2 , the DM particle χ is assumed to have the relic density of Ω c h 2 in the present universe, produced with a combination of thermal and non-thermal mechanisms.
The relic density 4 , as well as all other observables in DM experiments have been computed using the code micrOMEGAs [34] with the input model files for CalcHEP [35] generated using FeynRules [36]. However, in many cases, we replace the default parameters used in micrOmegas for astrophysical and nuclear physics inputs as described in the following subsections.

Spin-independent and spin-dependent scattering with nuclei
In the non-relativistic limit, relevant for spin-independent (SI) DM scattering with nuclei, the only Majorana DM bilinear that plays a role is the scalar one, χχ. The pseudo-scalar bilinear χγ 5 χ in L CPV 5 vanishes in the zero DM velocity limit, while the axial vector current χγ µ γ 5 χ leads to spin-dependent (SD) scattering. Therefore, only the scalar operator L CPC 5 is constrained by the SI scattering limits and the fourfermion interactions with SM quarks as well as the DM interaction with the Z-boson are constrained by the SD limits. For the recent-most bounds on these scattering cross-sections, we have used data from the LUX experiment [26] (for SI) and the XENON100 experiment (for SD) [27].
The event rate in direct detection experiments suffers from uncertainties coming from astrophysical and nuclear physics inputs. The astrophysical uncertainties originate from our lack of precise knowledge of the DM local density (ρ ⊙ ) as well as its velocity distribution in the rest frame of the detector, notes the DM velocity in the galactic rest frame, and v E represents the motion of the earth with respect to the galactic frame. At present, the experimental collaborations present their limits on SI and SD cross-sections by fixing the astrophysical inputs to specific values, and by assuming f ( v + v E ) to be a truncated Maxwell-Boltzmann distribution, with two additional parameters, the velocity dispersion v 2 and the galactic escape velocity v esc . In order to take into account the uncertainties of all the astrophysical parameters, we adopt the method developed in Ref. [37]. Given a choice of the DM density profile in the halo ρ χ (r), and a mass model for the Milky Way, Eddington's inversion formula [38] is used in Ref. [37] to determine the phase- The two primary assumptions in this approach are that the DM particle χ makes up the entire DM component of the Universe, and that the DM distribution is spherically symmetric. We adopt the phase-space density factor and its associated 2σ error bars as computed in Ref. [37], to which we refer the reader for further details 5 . The Burkert profile is chosen for ρ χ (r), which tends to give a slightly larger velocity dispersion compared to the NFW and Einasto profiles [37].
The nuclear physics uncertainties in SI and SD scattering rates stem from the corresponding nuclear matrix elements N|qq|N and N|qγ µ γ 5 q|N respectively.
Recent progress in lattice QCD calculations predict a rather small value for the strange quark content of the nucleon, f T s , and the results from different lattice simulation groups seem to have converged on this fact. Similarly, the pion-nucleon sigma term Σ πN , entering the SI rates along with f T s , has also been determined by lattice calculations rather precisely. For SD scattering rates, the nuclear physics inputs are encoded by ∆q n (q = u, d, s), which gives the fraction of spin due to each quark in the neutron (the corresponding numbers for proton are related by isospin rotation). The values of these nuclear physics inputs used in our calculations, along with their uncertainties, are listed in Table 2.
For further details on the construction of the likelihood function for LUX (SI) and XENON100 (SD), we refer the reader to Appendix-A.

DM annihilation in the present universe
Pair annihilation of DM particles in the present universe can lead to several SM final states. In our analysis, we include constraints from gamma ray and neutrino searches.
Hadronic nuisance parameters f T s 0.043 ± 0.011 [40] ∆u n −0.319 ± 0.066 [41] ∆d n 0.787 ± 0.158 [41] ∆s n −0.020 ± 0.011 [41]  Since the constraints obtained from the observation of gamma rays originating at Milky Way satellite galaxies are more robust than those obtained from Galactic Centre observations, we consider only the former 6 .
In this connection, it should be mentioned that DM pair annihilations can also lead to positron and anti-proton signals. However, astrophysical positron backgrounds are not yet known precisely enough to use them as robust constraints. As for anti-protons, although constraints can be obtained, but they depend strongly on the used propagation model for anti-protons under the galactic magnetic fields.
Since we want to determine as robust a limit on the parameter space as possible, even though relevant in certain cases, we abstain from using the positron and anti-proton data in this study.

Gamma ray observations
In the EFT setup considered here, the gamma ray line signal can either appear from loop-level processes involving SM fermions, or from dimension-7 operators which are suppressed for reasons explained in Sec. 2. Therefore, as far as indirect detection signals are concerned, the continuum gamma ray observations can put constraints on our scenario, although the scalar and axial-vector DM currents lead to annihilation rates in the present Universe which are p-wave suppressed. Hence the gamma ray observations are mostly relevant for the pseudo-scalar operator in our EFT setup.
We use the constraints obtained by Fermi Large Area Telescope (Fermi-LAT) observation of diffuse gamma rays from the milky way satellite galaxies (dwarf spheroidal galaxies (dSphs)) [28] by including the eight classical dSphs in our analysis, since the dark matter distribution in the classical dSphs is measured with a higher accuracy from the velocity dispersion of the luminous matter [44]. We combine the 273 weeks' Fermi-LAT data ( from 2008-08-04 to 2013-10-27) using the Pass-7 photon selection criterion, as implemented in the Fermi Tools. The J-factors for the dSphs included are taken from However, we predetermine, in a model independent manner, the likelihood map in the residual flux -energy plane, by combining the data for the eight classical dSphs.
Here, the residual flux refers to the background subtracted gamma ray flux, scaled by the J-factor. For details on the statistical analysis, we refer the reader to Ref. [45] 7 .
We include data in the whole energy range as observed by Fermi, namely, 200 MeV to 500 GeV. In the dSphs likelihood function used by us, the J-factor uncertainty is modelled by a Gaussian function, while the statistical uncertainty in the number of observed photons is modelled by a Poission distribution, with parameters as given in Ref. [45].

Neutrino telescopes
The IceCube neutrino telescope has been looking for muon neutrinos originating from DM annihilations in the Sun. The 317 days' data collected with the 79-string IceCube detector (including the DeepCore subarray) during the period June 2010 to May 2011 is found to be consistent with the expected atmospheric backgrounds [29], and therefore leads to bounds on DM annihilation rates in the Sun to final states involving neutrinos. Assuming that the DM capture and annihilation rates in the Sun are in equilibrium, these bounds can then be translated to limits on SI and SD scattering cross-sections of the DM particles with proton (σ SD p−χ ). While for SI scattering rates the XENON100 and LUX constraints are much stronger, for σ SD p−χ scattering with m χ > 35 GeV, the IceCube limits are stringent and competitive with ground-based experiments. We should remark here that, due to the higher spin 7 We thank Qiang Yuan for providing us the likelihood map using the 273 weeks' data, and Xiaoyuan Huang for careful cross-checks. expectation value of the neutron group compared to the proton group in XENON nuclei, the bounds on σ SD n−χ from XENON100 is stronger than the bound on σ SD p−χ . This is the reason the IceCube limits on σ SD p−χ are important even after including the XENON100 limits on SD scattering rates, since in the EFT, the σ SD p−χ can be enhanced compared to σ SD n−χ in certain regions of the parameter space (for example, with |g LQ − g Ru | > |g LQ − g Rd |).
As is well-known, the energy spectrum of neutrinos produced from the annihilation channel χχ → W + W − (or, χχ → τ + τ − for m χ < M W ) is harder than that produced from the channel χχ → bb. Hence, the former leads to stronger bounds on σ SD p−χ , assuming a 100% annihilation to either of the channels. Since the limit from the bb channel is the weakest, any parameter point in the EFT, with any branching ratio to the bb final state, has to at least satisfy this limit. Therefore, in our analysis, we reject parameter points which lead to σ SD p−χ higher than the 95% C.L. IceCube bound as obtained assuming the bb mode of DM annihilation.

Z-boson invisible width
For m χ < M Z /2, the interaction term in Eq. 10 will lead to the invisible decay Z → χχ. The decay width of the Z-boson has been precisely measured at the LEP experiment, and apart from the width originating from Z → νν, the 95% C.L. upper bound on the invisible width of the Z-boson is given by [30] Γ Z inv < 2 MeV (95% C.L., LEP).
This, therefore, acts as a powerful constraint on the dimension-6 operator in Eq. 9 for a light DM. The likelihood function is taken as a Gaussian with parameters as shown in Table 1.

Mono-photon search
Single photon events were looked for at the LEP collider to search for signatures of graviton production, the null results of which leads to bounds on the radiative process e + e − → χχγ. In this paper, we use the limits from the DELPHI collaboration obtained using the 650 pb −1 of LEP2 data with centre of mass energy in the range 180 − 209 GeV [31]. We compute the relevant cross-sections using MadGraph5 We have used Γ h SM Total = 4.21 MeV for m h = 126.0 GeV in our numerical analysis. The Higgs invisible search constraints the CP-even and -odd scalar operators considerably, in the DM mass range in which the decay mode is kinematically allowed. upper bounds in the m χ − Λ plane for the axial vector current-current four fermi interaction with quarks. The agreement found is accurate to within 5%. We briefly discuss this validation procedure and our MC setup in Appendix-B, where the LHC likelihood function is also described.
The validity of an EFT approach at a high-energy hadron collider like the LHC has recently been discussed at length [22,23]. Since rather strong jet p T and E T / cuts are used by the collaborations in order to reduce the very large SM backgrounds, the subprocess centre of mass energy involved is also large. In such a case, if the suppression scale of the operators Λ (or equivalently the mediator mass in a weakly coupled UV completion) is comparable to the subprocess COM energy, the mediator particles can be produced on-shell, and an EFT description breaks down. As an order of magnitude estimate, one can consider the minimum value of the cut-off scale, in order for the EFT to be valid as [22] Λ > 4m 2 where, for the CMS analysis, E T / min = 400 GeV. For low DM masses, the above requirement is dominated by the E T / cut. Comparing this order of magnitude bound to the condition imposed by us (Λ > 10m χ ), we can see that, roughly for m χ > 41 GeV, we can safely use the EFT framework at the LHC when Λ > 10m χ . For 10 GeV < m χ < 40 GeV, we find that even if Λ > 10m χ , the expected cross-section after the CMS cuts is so large that even if the EFT limit is stronger than that given by an UV complete theory including the mediator particles, one can still exclude these points in the UV theory as well.

Relic density: the role of individual operators
The requirement that Ω χ h 2 (Thermal) ≤ Ω c h 2 , places important lower bounds on the operator coefficients g i /Λ and g i /Λ 2 . Although all the couplings enter together   Ω χ h 2 is observed in Fig. 1. This essentially implies that for those m χ the requirement on Ω χ h 2 can be satisfied even for very large values of Λ or for very small values of the coupling constants g i . We now briefly discuss the case for each operator separately.

Dimension-5 Higgs portal operators
With only the CP-conserving Higgs-portal operator L CPC

Dimension-6 four-fermion operators
In Fig 1 we also

Allowed region in m χ − Λ plane
In the EFT analysis including all operators upto dimension-6 simultaneously, the most important two-dimensional likelihood map is in the m χ − Λ plane. In Fig. 2 (left panel) we show the 68% (yellow) and 95% C.L. (blue) allowed regions in the m χ −Λ plane, whereby the condition Λ > 10m χ has been imposed. For completeness, 8 An additional large s-wave contribution to σv F.O. , and therefore a drop in Ω χ h 2 , is expected after the tt annihilation mode opens up. However, this is not visible in Fig. 1, since for this figure we have set g LQ = g Ru = 1, and therefore only the vector currenttγ µ t survives, which does not lead to an m 2 t enhanced term in σv F.O. . When g LQ = g Ru such a feature is observable due to the additional axial vector current, see, for e.g., Fig. 2 (left panel).
we also show the 95% C.L. (grey) allowed region with 2m χ < Λ < 10m χ , although as discussed before, in this region the validity of the EFT from the point of view of relic abundance and collider physics computations is questionable. Finally, the region with Λ < 2m χ , where an EFT analysis cannot be applied, is shown in pink.
As mentioned in Sec. 3.1, the two-dimensional allowed regions are obtained after profiling out rest of the parameters in the EFT, which, in this case, are the O(1) coefficients g i for the different operators. The shape of the allowed contours is mostly determined by the requirement that Ω χ h 2 (Thermal) < Ω c h 2 (Planck) as discussed in detail in the previous subsection.
We now discuss the allowed and ruled out regions of DM mass m χ by dividing them in a few windows:  to m 2 t , and therefore slightly larger values of Λ are also allowed.
It is interesting to note that, when the condition Λ > 10m χ is imposed, there is an upper limit on the allowed value of m χ , mainly coming from the relic density constraint. For Λ > 10m χ , this upper limit is m χ ≃ 300 GeV, while for Λ > 2m χ , it extends upto m χ ≃ 3 TeV, as can be seen in Fig. 2.

Future prospects
Having determined the currently allowed parameter space for the CP-conserving scenario in the m χ − Λ plane in the previous subsection, we now discuss the predictions for future experiments in this region. We focus on the future direct detection experiments (both SI and SD), as well as an e + e − collider experiment like the proposed International Linear Collider (ILC) [49]. In the CP-conserving case, the predictions for DM annihilation rates in the present Universe, σv 0 , are rather low in our EFT framework (below 10 −27 cm 3 sec −1 ), since most of the operators lead to annihilations which are p-wave suppressed. Therefore, the reach of future gamma-ray observation experiments is not very high in this case, except for m χ > m t , where one can obtain a rate as high as σv 0 ∼ O(10 −26 ) cm 3 s −1 due to the additional s-wave contribution.
In Fig. 3 (top-left) we show σ SI p as a function of m χ as found in the allowed EFT parameter space. The blue (yellow) shaded region is allowed at 95% (68%) C.L. after profiling over all the other parameters except m χ with the condition Λ > 10m χ . If the last condition is relaxed to Λ > 2m χ the additional grey shaded region also survives at 95% C.L. There is a small difference between the 90% LUX constraint (red solid line) and the 95% C.L. region obtained in our scan after including the LUX limits, primarily because the astrophysical parameters related to the DM local density and the DM velocity distribution are kept fixed in the LUX analysis, while they have been profiled out in our analysis including the errors as discussed in Sec. 3.   After a particular detector design is adopted, the overall efficiencies will be modified, and the upper bound will be somewhat shifted. Therefore, Fig. 3 gives us the order of magnitude expectation of σ(e + e − → χχγ) in the allowed region of the EFT, and we see that cross-sections in the range of 2 pb to 10 −3 fb or lower can be possible.
The upper bounds from future ILC measurements can be around O(fb) and thus can cover a significant range of the possible cross-sections. It is clear that the resonance regions (Higgs or Z) may not yield high signal rates, as small couplings suffice to satisfy the relic abundance requirement.
Finally, we show in Fig. 4 a combined view of the different possible future constraints, and the EFT parameter space in the m χ − Λ plane that survives at 95% C.L. after applying each of them individually. The grey points marked as δχ 2 < 5.99 represent the currently allowed parameter space. We first apply the projections of the XENON-1T experiment for SI and SD direct detection rates, and the yellow points are found to survive the cuts. The prospects of the ILC are considered next, whereby we include the expected constraints from a very precise measurement of the Higgs boson invisible branching ratio and from the mono-photon+E T / search. The sensitivity of the Higgs invisible branching ratio determination is expected to reach 0.7% level in the ZH production mode after including the Z → qq decay channel [54].
For the γ + E T / search in e + e − collisions, we have taken into account collisions with √ s = 250, 500 and 1000 GeV, in each case assuming an integrated luminosity of 1 ab −1 . The reason for taking into account all the centre of mass energies is the fact that for a given √ s, we can only reliably compute the production cross-sections in an EFT for Λ > √ s. Therefore, only by taking into account the three different energies, we can cover an wide range of values of Λ. The green points remain after imposing the ILC projections, a large fraction of which are in the Higgs resonance region. However, this region can also be covered by the projected sensitivity of the LZ experiment, which is applied next, and the blue points survive the LZ projections of SI and SD scattering rates.
For Λ > 10m χ , the finally surviving points after considering all the future ex-periments, are concentrated in two regions. One of them is the Z-resonance region, and the other one is the region where the DM is heavier than 100 GeV. It is interesting to note that if the ILC can be operated at a centre of mass energy close to the Z-pole (the so-called Giga-Z option), then a significant part of the remaining points in the Z-resonance region can also be covered. For 2m χ Λ 10m χ , a large number of points survive for m χ 200 GeV, which cannot be covered by the future experiments considered in our study.
The high-energy run of the LHC at 13 TeV, including its upgrade to a highluminosity phase, may be able to cover part of the last mentioned region above using the mono-jet+E T / channel. However, usually we require strong kinematical cuts at 13 TeV LHC in order to suppress the very large background coming from the Z(→ νν) + jets process, which essentially translates to a requirement of high subprocess centre of mass energies in the events (of the order of 1 TeV or higher). Since Λ in this region also lies in the ballpark of a TeV, it is difficult to estimate the future LHC reach in this part of the parameter space within the EFT. A proper estimate can only be made by using suitable simplified models representing the effective operators (with s-or t-channel mediators between the DM and the SM sector). Such a study is, however, beyond the scope of the present article.  experiments, in the χχ → W + W − channel, as estimated in Ref. [56]. For m χ < M W , the χχ → bb channel is relevant, which gives rise to similar constraints as well. This sensitivity line is only indicative of the future reach, and will vary depending upon the specific annihilation channel(s) relevant in different parameter regions of the EFT.
In Fig. 6 (right panel) we show the surviving parameter points in the CP-violating scenario, in the m χ − Λ plane. The various constraints imposed in this figure are discussed in the context of the CP-conserving scenario (Sec. 5.2). In addition to the Z-resonance region and a much larger bulk region with m χ > 100 GeV, an additional set of points in the Higgs-resonance region also survives all the future experiments considered, since the CP-violating operator is not constrained by SI scattering rates, which, for the CP-conserving case, could completely cover the Higgs-resonance part.
If indeed a CP-violating operator exists in the DM sector, it is then imperative to find out other possible experiments which can probe such a coupling.

Summary
To summarize, we studied a standard model singlet Majorana fermion DM candidate in the framework of an effective field theory, by taking into account the presence of a complete basis of gauge invariant operators of dimensions 5 and 6 at the same time.
The profile likelihood method was used to determine the currently allowed region of parameter space at 95% C.L., and the span of different DM related observables as a function of its mass. Considering an weakly coupled ultra-violet completion for the EFT, we find that imposing the condition Λ > 10m χ makes the EFT an excellent approximation to possible UV complete models. For completeness, we have also computed the confidence intervals in the region 2m χ < Λ < 10m χ . Including the Among the presently allowed parameter region, surprisingly enough, it seems viable that the Higgs resonance region will be completely covered by the combined efforts of future ton-scale direct detection experiments like XENON1T and LUX, as well as the measurement of the Higgs invisible branching ratio at the proposed International Linear Collider. The main difference between the Higgs and Z-resonance regions is that the former leads to SI scattering rates, while the latter to SD ones.
Therefore, the direct detection experiments will probably not be able to cover whole of the Z-resonance region with their currently projected sensitivities. Here also, the ILC can play a major role, when we use the mono-photon search channel for DM pair-production. For DM heavier than the top mass, where the four-fermion operators play a major role, there is a large set of allowed points if 2m χ < Λ < 10m χ . If the couplings to quarks are flavour universal, the high energy and high luminosity LHC runs should have an impact here. However, it is difficult to judge it within an EFT, especially if only the weak condition Λ > 2m χ holds. We leave a detailed study of such a region within realistic simplified models with s-and t-channel mediator exchange for a future work.
If the CP-violating operator is present, the allowed parameter region is larger, since DM pair annihilation to light quarks are no longer p-wave suppressed and contribute to Ω χ h 2 (Thermal). DM masses of upto about 1 TeV now become viable even with Λ > 10m χ . The same reason also gives us a hope of testing at least part of the parameter space using gamma-ray observations from DM annihilation to quarks and leptons. We presented some estimates of this indirect detection capability, and further detailed studies are necessary here as well. It also seems necessary to pursue new avenues to probe the CP-violating coupling itself, since it is not testable using low-energy scattering experiments with nuclei.
Appendix A: The direct detection likelihood As discussed in Sec. 3.3, the LUX collaboration obtained the bounds on σ SI p by fixing the astrophysical inputs to specific values and by taking the velocity distribution of DM to be a truncated Maxwellian [26]. Since we use the DM phases-space density computed in Ref. [37] along with its 2σ error bars, we first rescale the LUX bounds accordingly. For calculating the time-averaged number of events we adopted the Helm form factor [57] for σ SI p as used by LUX, together with the energy resolution function and efficiency factors for the LUX detector [26]. In the following, the rescaled values of the LUX 90% C.L. bounds as a function of the DM mass are denoted by σ SI p,90% . A similar approach is applied to rescale the XENON100 [27] bounds on σ SD n as well, with the corresponding form factors relevant for SD scattering [58].
In constructing the likelihood, the nuclear physics inputs discussed in Sec. 3.3, namely f T s and Σ πN for σ SI p , and ∆q n,p for σ SD n , are treated as nuisance parameters. Together with the rescaled 90% C.L. limit from LUX for σ SI p and from XENON100  The likelihood function for the LHC monojet search can be expressed as follows where, b is the expected number of background events, s is the expected number of signal events for a given parameter point and N obs represents the number of events observed by CMS after employing the kinematic selection criterion described above. The systematic uncertainty is taken into account by convoluting the Poission likelihood function with a Gaussian with mean b (which is treated as a nuisance parameter and profiled out) and variance σ b .