Signals of Leptophilic Dark Matter at the ILC

Adopting a model independent approach, we constrain the various effective interactions of leptophilic DM particle with the visible world from the WMAP and Planck data. The thermally averaged indirect DM annihilation cross-section and the DM-electron direct-detection cross-section for such a DM candidate are observed to be consistent with the respective experimental data. We study the production of cosmologically allowed leptophilic DM in association with $Z\, (Z\to f\bar f)$, $f\equiv q,\,e^-,\, \mu^-$ at the ILC. We perform the $\chi^2$ analysis and compute the 99\% C.L. acceptance contours in the $m_\chi$ and $\Lambda$ plane from the two dimensional differential distributions of various kinematic observables obtained after employing parton showering and hadronization to the simulated data. We observe that the dominant hadronic channel provides the best kinematic reach of 2.62 TeV ($m_\chi$ = 25 GeV), which further improves to 3.13 TeV for polarized beams at $\sqrt{s} = 1$ TeV and an integrated luminosity of 1 ab$^{-1}$.


Introduction
Dark Matter provides the most compelling explanation for many cosmological and astrophysical observations which defy an understanding in terms of luminous matter alone. However, in the absence of any direct evidence, the existence of such matter has, rightly, been a e-mail: Sukanta.Dutta@gmail.com b Corresponding Author, e-mail: bhartirawat87@gmail.com c e-mail: divyasachdeva951@gmail.com questioned and several alternatives proposed, including the modification of gravity at large distance scales. It should be appreciated, though, that starting with orbital velocities of stars in a galaxy or those of the galaxies themselves in a cluster [2], gravitational lensing [3,4], the dynamics of galactic (cluster) collisions [5] etc., the observations span a wide range of distance scales, and no single simple modification of gravitation can explain all, whereas dark matter (DM) can and does play an important role in understanding the data. Similarly, the fitting of the cosmological observables [6][7][8], requires that DM contributes about 23% to the energy budget of the universe, in contrast to only 4% contained in baryonic matter. Finally, postinflation, perturbations in the DM distribution, along with the gravitational perturbations, are supposed to have provided the seed for large-scale structure formation in the universe [7]. The constraints from the latter are rather strong. Indeed, while neutrinos would have been logical candidates for DM within the Standard Model (SM), a large energy component in neutrinos would have disrupted structure formation. For example, a recent study claims that if the equation of state for the DM be parametrized as p = w ρ, then the combined fitting of the cosmic microwave background, the baryon acoustic oscillation data and the Hubble telescope data restricts −9.0 × 10 −3 < w < 2.4 × 10 −3 , thereby clearly preferring a cold and dusty DM [9].
With the DM particle, by definition, not being allowed to have either electromagnetic 1 or strong interaction, the most popular 2 candidate is the weakly interacting massive particle (WIMP). Recently, the au-thors of reference [14] have proposed the cosmologically allowed super-light fuzzy dark matter of the order of 10 −22 eV which can possibly be testable in terrestrial neutrino oscillation experiments. And while its exact nature is unknown, several theoretical scenarios such as multi-Higgs models, supersymmetry, extra-dimensional theories, little Higgs models, left-right symmetric models all naturally admit viable candidates. Consequently, one of the most challenging tasks today is to identify the nature of the DM particle [15]. This, in principle, could be done in three kind of experiments. Direct detection can be achieved by setting up apparatii (very often ultra-cold bolometric devices) that would register the scattering of a DM particle off the detector material. While the DAMA experiment [16] did indeed claim the existence of a DM particle of mass ∼ 60 GeV from observed seasonal variation in the detector signal (originating, possibly, from a varying DM wind as the Earth traverses its path), subsequent experiments like CoGeNT [17], CRESST-II [18], XENON100 [19], PandaX-II [20,21] and LUX [22] have not validated this; rather, they have only served to impose bounds in the plane described by the DM particle mass and its coupling to nucleons. Indirect detection experiments, largely satellite-based, depend on the annihilation of a pair of DM particles in to SM particles which are, subsequently, detected. Although there, occasionally, have been claims of anomalies in the data, unfortunately the experiments have failed to validate each other's positive sightings, resulting, once again, in further constraints [23][24][25]. Finally, we have the collider experiments, wherein excesses (over the SM expectations) in final states with large missing momentum are looked for. It should be realized, though, that even if such an excess is established, a DM explanation would still only be an hypothesis, for the only statement that can be made with certainty is that the produced neutral and colour-singlet particle is stable over the detector dimensions.
It is the last mentioned approach that we assume in this paper. While an investigation of the nature of the DM particles needs an understanding of the underlying physics, we adopt, instead, a model independent approach. Eschewing the details of the underlying dynamics, we begin by postulating a fermionic DM particle, and consider four-fermi operators involving these and the SM fermions. The relative strengths of these operators would, of course, be determined by the underlying dynamics. We assume that the operators involving quark currents are subdominant, as could happen, for example, if the dynamics, at a more fundauntestable in the laboratory in the foreseeable future. Similar is the case for an ultralight gravitino DM candidate [12,13] mental level, involved a leptophilic boson. This immediately negates the constraints from the direct detection experiments [19][20][21][22] (as the dominant interactions therein are with the nucleons) as also bounds from the LHC [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].
Such an assumption also alters the conclusions (on the interrelationship between the DM-mass and the coupling strength) drawn from the deduced relic density. Starting with the Boltzmann equations describing the evolution of the particle densities, we derive the constraints on the same. As for the search strategies at a linear collider, the attention of the community, so far, has been largely commandeered by the final state comprising a single photon accompanied by missing energy [41][42][43][44]. While this continues to be an important channel, we augment this by other channels that are nearly as sensitive. Furthermore, we consider some novel effects of beam polarization.

Fermionic Dark Matter : A mini review
As we have already mentioned, rather than considering the intricacies of particular models, we adopt the more conservative yet powerful concept of an effective field theory. Assuming that we have a single-component DM in the shape of a Dirac fermion χ, we consider the least-suppressed (namely, dimension-six) operators involving a χ-current and a SM fermion-current 3 . This simplifying assumption rules out the possibility of resonances and co-annihilations affecting, significantly, the relic abundance of dark matter particle [45]. Furthermore, to reduce the number of possible operators, as also not induce flavour changing neutral current processes, we restrict ourselves to only flavour-diagonal currents. A convenient parametrization of such operators, for a single SM fermion ψ, is given by with the subscripts on the operators reflecting the Lorentz structure. The full interaction Lagrangian involving χ can, then, be parametrized as is a typical operator amongst those listed in eq.(1) and g f M N the corresponding strength. Λ refers to the cut-off scale of the effective theory.
Before we analyze any further, it is useful to consider the DM pair-annihilation cross section engendered by these operators [46]. Restricting ourselves to a single species 4 f , and denoting x i ≡ 2m 2 i /s, we have where N c = 1 (3) for leptons (quarks) and β χ,f are the center-of-mass frame velocities of χ (f ) with β i = √ 1 − 2 x i . Clearly, for non-relativistic DM particles, σ ann SS and σ ann SP are the smallest. Similarly, for m χ m f (as is the case except when f is the top-quark), we have These approximations would prove to be of considerable use in further analysis. Dark matter interactions, as exemplified by equation (1), have been well studied in the context of the LHC [26][27][28][29][30][31][32][33][34][35][37][38][39][40]. Clearly, the LHC would be more sensitive to operators involving the quarks. In this paper, we would be primarily interested in the orthogonal scenario, namely one wherein the DM current couples primarily to the SM leptons. To this end, we will assume that g q M N = 0 for all of the operators. As we shall see in the next section, this would render insensitive even the dedicated direct search experiments.

Relic Density of Leptophilic Dark Matter
Before we proceed further with our analysis, it behooves us to consider the existing constraints on such a DM candidate. This is best understood if we consider a single operator at a time, and, henceforth, we shall assume this to be so. In other words, all operators bar the one under discussion shall be switched off. As to the couplings g f M N , we consider two cases 5 , namely (a) Leptophilic: g l M N = 1 for all leptons l, and (b) Electrophilic: g e M N = 1 with all others vanishing. Concentrating first on the epoch where the DM (χ) has frozen out, let us consider their pair-annihilation to produce two light particles ( ), which, we assume, is in complete equilibrium, viz. n = n eq . The Boltzmann equation gives us: where a ≡ a(t) is the scale factor and · · · represents the thermal average. Here, v is the relative velocity of a pair of χ's, and σ is the annihilation cross section.
In the radiation era, which we are interested in, the temperature scales as a −1 . We can, therefore, rewrite the Boltzmann equation, in terms of new dimensionless variables as where H(m χ ) is the Hubble rate when the temperature equals m χ and is given by Here g is the effective number of the degrees of freedom at that epoch. To find the present density of DM particles, we need to find the solution of eq. (6) in terms of the final freeze out abundance, Y ∞ (at x = ∞). While this, can be done only numerically, it is instructive to consider an approximate analytic solution. At early times (T m χ ) reactions proceeded rapidly, and n χ was close to its equilibrium value, n eq ∝ T 3 , or Y ≈ Y eq . As the temperature falls below m χ , the equilibrium abundance Y eq = g/(2 π) 3 x 3/2 e −x is exponentially suppressed. Consequently, the DM particles become so rare that they do not find each other fast enough to maintain the equilibrium abundance. This is the onset of freeze-out. In other words, well after freeze out, Y Y eq and Integrating equation (8) from the freeze out temperature x f until very late times x = ∞, we get The energy density, for the now non-relativistic DM, ρ χ = m χ n χ , falls as a −3 . The number density n(a(T i ), T i ) in this post freeze out epoch at any given temperature i . The number of the degrees of freedom also changes with the temperature from g f (x f ) at the freezing epoch to that of today g 0 (x 0 ) ∼ 3.36 at the present day temperature T 0 . The present energy density of DM can, thus, be re-written as It is customary to parametrize ρ χ ≡ Ω χ h 2 ρ c , where ρ c = 1.05375 × 10 −5 h 2 GeV /c 2 cm −3 is the critical density of the universe and the Hubble constant today being expressed as H 0 = h × 100 km s −1 Mpc −1 .
We have, then, Using T 0 = 2.72548 ± 0.00057 K [47], the contribution of the Dirac fermion χ, with its two degrees of freedom, is While the WMAP 9-year data gives Ω DM h 2 = 0.1138 ± 0.0045, the Planck 2015 results [8] suggest a marginally different value, namely, Ω DM h 2 = 0.1199 ± 0.0022. This translates to a strict relation in the mass-coupling plane for χ. However, a conservative assumption would be that the relic density of the χ does not saturate the observed DM energy density, or Ω χ < Ω DM . This would impose a lower bound on the self-annihilation cross-section for the χ, or, in other words, an upper bound on Λ.

The annihilation cross-sections
In the calculation of the probability of a DM particle being annihilated by another one, it is useful to consider the rest frame of the first and re-express the aforementioned annihilation cross sections in terms of the velocity v of the second particle in his frame. Clearly, v = (2 β χ ) / 1 + 2 β 2 χ and, working with a small v expansion (relevant for non-relativistic DM), we have and s m 2 Clearly, the v-independent terms emanate from s-wave scattering alone, while the O(v 2 ) piece receives contributions from both s-wave and pwave annihilation processes.
The quantity of interest is not just σ ann , but the thermal average σ ann M N v since this provides a measure of the rate at which a DM particle will meet another with the appropriate velocity and get annihilated. Assuming χχ → e + e to be the dominant channel into which the χ's annihilate, we can estimate the thermallyaveraged value of σ(χχ → e + e − ) v.
The inferred value of Ω χ h 2 when ascribed to WIMPs with a mass of a few hundred GeVs, requires σ(χχ → e + e − ) v ≈ 3 × 10 −26 cm 3 s −1 . Taking cue from this, we analyze upper bounds on σ(χχ → e + e − ) v in section 5 from the sensitivity study of the DM signatures from ILC for the kinematically accessible DM mass range.
In Figure 1 we present the relic density plots for the interactions parametrized by the operators mentioned in section 2 corresponding to the Leptophilic and Electrophilic scenarios. We implement the model containing effective interactions of DM particle with the SM sector using FeynRules [49] and generate the model files for CalcHEP which is required for the relic density calculation in micrOmegas [50]. These results have been verified with MadDM [51,52].
The curves in Figure 1 imply that the correct amount of energy density requires the mass of dark matter particles to increase with Λ. We can infer this behaviour from (12) as Ω χ h 2 mainly depends on σv . From equations (13a)-(13d) we see that the thermal average in equation (12) i.e. σv is proportional to m 2 χ /Λ 4 for m χ values greater than ∼ 10 GeV. This implies that for a fixed value of Ω χ h 2 , Λ 4 /m 2 χ is a constant. This, in turn, translates to Λ ∝ √ m χ . WMAP data [48]. These curves are Obtained using micrOmegas-v-3.5.5 [50], the regions below the curves are allowed by the relic density bounds.

DM-electron scattering:
We now turn to direct detection search, where both electrophilic and leptophilic DM are scattered by (a) the bound electron of an atom and hence the recoiled electron is ejected out of the atom, (b) the bound electron and the recoil is taken by the atom as a whole and (c) the quarks, where the DM is attached to the loop of charged virtual leptons which in turn interacts with the quarks via photon exchange and hence the nucleus recoil is measured in the detector.
In this article we restrict our study for the DMelectron scattering only. We analyze the free electron scattering cross-sections with the DM for the scalar, pseudoscalar, vector and axial couplings respectively and are given as where v ∼ 10 −3 c, is the DM velocity in our halo. We observe that the cross section with the axial-vector couplings of DM dominates over all others. The pesudoscalar type coupling is highly suppressed due to its dependence on the fourth power of the electron mass and velocity. Respective cross sections for the case of bound electrons are likely to be enhanced than that of free electrons nearly by a factor of ∼ 10 5 [53]. Using the constraints from the relic density on the lower bound of the couplings for a given DM mass m χ , in Figure 2, we present the scattering cross section of the DM with free electrons via the most dominant channel with axialvector couplings. We compare our result with the recent 3σ limits from the DAMA/ LIBRA experiment [53] and 90% C.L. data from XENON100 [54]. However, even with the inclusion of the correction factor arising from the the bounded electrons, the contribution of the axialvector couplings of DM to the scattering cross-section are much below these experimental upper bounds. Indeed, the Xenon100 collaboration has constrained the σ(χe → χe) < ∼ 10 pb at 90% C.L. for 6 0.6 GeV < m χ < 1 TeV [54] for the case of axial-vector coupling. Expectedly, these constraints are much stronger (by nearly an order of magnitude) than the cross sections required to explain the DAMA/Libra results [16]. Interestingly, the strongest bounds emanate from solar physics. In general, the neutrino flux coming from DM annihilations inside the sun is proportional to the DM scattering cross section. But, working in the region where DM capture and its annihilations inside the star are in equilibrium, makes the neutrino flux independent of the DM annihilation. As the DM particles can be trapped in the solar gravitational field, their annihilation into neutrinos would modify the solar neutrino spectrum and, consequently, the SuperKamiokande measurements [53]. This translates to σ(χχ → τ + τ − ) < ∼ 0.1 pb and σ(χχ →  with the axial-vector interactions σ χ e − → χ e − is depicted for the DM mass varying between 25-400 GeV. Exclusion plots from DAMA at 90% and 3σ C.L. for the case of DMelectron scattering are also shown [53]. Bounds at 90% C.L. are shown for XENON100 from inelastic DM-atom scattering [54]. The dashed curves show the 90% CL constraint from the Super-Kamiokande limit on neutrinos from the Sun, by assuming annihilation into τ + τ − or νν [53]. versus the DM mass m χ using the relic density allowed lower bound on S-S, P-P, V-V and A-A couplings for DM. These cross-sections are compared with the median of the DM annihilation cross section derived from a combined analysis of the nominal target sample for the τ + τ − channels [23].

νν) <
∼ 0.05 pb which are relevant for the general leptophilic case with axial-vector coupling. Understandably, no such bound exists for the purely electrophilic case.

DM Annihilation:
Next we discuss the bounds from the indirect detection experiments, wherein we can directly use the cross sections given in equations (3a)-(3h) and fold them with the local velocity of the DM particle. The latter, of course, would depend crucially on the particular profile of the DM distribution that one adopts. In [23], the expected sensitivity is given as a limit on the thermally averaged DM annihilation cross section for the τ + τ − channel.
From equation (13a)-(13d), the thermally averaged annihilation cross-section of a leptophilic DM of mass m χ can be estimated using the relic-density allowed lower bound on the scalar, pesudoscalar, vector and axial-vector couplings of DM. The annihilation crosssection for the leptophilic DM is marginally higher than electrophilic DM and roughly are of the same order of magnitude for m χ ∼ 10 − 400 GeV. We compute the annihilation cross-section for the leptophilic DM with these lower bounds on DM couplings and depict in Figure 3. We compare our thermally averaged crosssection of leptophilic DM with the bounds from the current indirect detection searches [23].

Interesting Signatures at ILC
While DM particles can be produced in many different processes at a given collider, only a few of them are, potentially, of interest. Remembering that the DM particle has to be produced only in pairs, and that there must be at least one visible particle in the final state so as to register the event in a detector, the simplest process at a linear collider is, of course, e + e − →χχγ, and it has, rightly, been well studied in many different contexts [41][42][43]. Although the signature, namely monophoton with missing energy-momentum, can be masqueraded by both detector effects as well as a large and irreducible background emanating from e + e − → ν i ν i γ, the kinematic profiles are sufficiently different enough to merit the possibility of a discovery.
Given the exceeding simplicity of the aforementioned channel, attention has focussed on it almost exclusively. However, it is worthwhile to explore complementary channels, like, DM pair-production in association with on-shell Z or off-shell Z , e + e − → χχ + Z Z → jj, µ + µ − , e + e − (15) where j refers to a jet arising from any of u, d, s, c, b. In this paper, we shall not attempt to impose any btagging (or b-veto) requirements, and, hence, the b-jet is equivalent to any other light quark jet 7 . To effect Table 1: Accelerator parameters for each of the run scenarios considered in this paper. P − and P + represent the electron and positron polarizations respectively.

Unpolarized
Polarized respectively. The largest contribution to the first mode in eqn. (16a) accrues from e + e − → Z + γ * /Z * with Z → ν iνi and γ * /Z * splitting to two jets followed by the W + W − fusion channel e + e − → ν eνe W + W − → ν eνe jj and t-channel W exchange diagrams. Feynman diagrams analogous to both these sets also contribute to the two other processes given in (16b) and (16c) respectively. The relative strengths of these channels depend on a multitude of factors, such as the center-of mass energy, beam polarization (if any) and the kinematical cuts imposed. We use MadGraph5 [57] to perform the simulations for both SM background and the signal. We employ PYTHIA6 [58] to carry out the parton shower and we adopt an approach identical to theirs. Furthermore, we present, here, a much more refined analysis. 8 In all our analysis, the dijet (plus missing energy) signal is actually an inclusive one, in that we demand a minimum of two jets. hadronization. PGS [59] is used for Fast detector simulation with anti-k T algorithm for jet reconstruction. In accordance with technical design report for ILC detectors [60], we set the energy smearing parameters of the electromagnetic calorimeter and the hadronic calorimeter in the PGS card as The other accelerator parameters used in our analysis are also set according to the technical design report for ILC [61] and we tabulate them in Table 1 for convenience. Before we commence an analysis of the cross sections and the kinematic distributions we must impose basic (acceptance) cuts on transverse momenta (p T j/ ), rapidities (η j/ ) and separation (∆R jj/ ) of the visible entities (jets and µ ± , e ± ) along with the the missing transverse energy ( E T ), which commensurate with the typical detector capabilities on the one hand, and theoretical considerations (such as jet definitions) on the other.
p Ti ≥10 GeV where i = , j , -|η | ≤ 2.5 and |η j | ≤ 5,   Fig. 4: Normalized 1-D differential distribution of the cross-section w.r.t. kinematic observable P T j 1 , E T , sum of energy of jet pairs and ∆φ jj for the background process (shaded histograms) and the Z → j j associated DM pair production at the three representative values of M χ = 75, 225 and 325 GeV with √ s = 1 TeV, Λ =1 TeV and g l SS (g l P P , g l SP , g l P S ) =1. The bin width is 0.1 for the ∆φ ii distribution while it is 15 GeV for the E T and E j 1 + E j 2 .
The corresponding cross-sections with these basic cuts for the polarized and unpolarized beams of electron and positron are summarized in Table 2. Feynman diagrams consisting of different topologies and contributing to the background (16a) can be further suppressed by imposing the most-dominant invariant-mass cut of the visible particles 9 , and thus allowing Z + Z ( * ) final state with the on-shell Z decaying into the two visible entities, and the off-shell Z ( * ) going to a neutrino-pair. Within this approximation, the diagrams contributing to the SM background processes (16a)-(16c) have identical topologies and differ only in the coupling of the Z to the final state fermions. Consequently, this difference would manifest itself essentially only in the total rates with differences in the normalized distributions being of a subleading nature. This is true even for the signal events corresponding to these final states. Hence we have displayed the one-dimensional normalized differential cross-section for process with only jets as final state in Figures 4 and 5.

Analysis with one and two dimensional distributions:-
In order to augment the signal-to-noise ratio by imposing additional selection cuts on the kinematic observables, we recourse ourselves to the detailed study of the normalized one dimensional differential crosssection distributions w.r.t. for all the sensitive kinematic observables e.g. p T j/ l , E T , ∆φ ii and sum of energies from a jet or lepton pair (E i1 + E i2 ) for the signal with varying DM mass as well as varying effective DM pair -SM fermion pair couplings g f M N /Λ 2 at √ s = 500 GeV and 1 TeV.  The very structure of the DM interaction Lagrangian implies that, for a given final state, the normalized differential distributions for the SS, SP , P S and P P cases are very similar to each other, with the differences being proportional to the difference in mass of the SM fermions in the final state (i.e., e, µ or the light quarks). Similarly, the normalized distributions for the V V , V A, AV and the AA cases are, again, very similar to each other. Therefore, we have chosen to display the normalized histograms for the observables p T j/ l , E T , ∆φ ii and (E j1 + E j2 ) from signal processes at the three representative values of the DM masses m χ namely 75, 225 and 325 GeV in Figures 4 and 5 corresponding to the scalar and vector interactions of the DM bilinears with SM fermion pair at √ s = 1 TeV and Λ = 1 TeV. We list our observations as follows:-• In Figure 4, we observe that, as we increase the DM mass, the respective shape profile of the all kinematic observables become more and more distinct w.r.t. the SM contribution. This difference in shape profile for the scalar coupling with that of SM can be attributed to the chirality-flip interactions of the DM which becomes much more pronounced for higher values of m χ . On contrary, the shape profiles of all the kinematic distributions with vector coupling of DM in Figure  5, behave the same as those of SM Z → νν coupling. • We find that in the lower m χ region, the jet pair remains highly boosted as in the case of SM and hence, the peak of the differential cross-section distribution w.r.t. ∆φ jj coincides with that of SM background around 35 • . The peak however shifts towards higher values of ∆φ jj as we increase m χ . • We find that the energy profile of the visible pair E j1 + E j2 are sensitive to the choice of the DM mass and thus advocating that the upper limit on this kinematic observable can be used as a dynamic mass dependent selection cut to reduce the continuum background. This is realized by invoking a cut on the invariant mass of visible particles to be around This restriction translates to an upper limit on the variable E j1 + E j2 as a function of m χ : Since we are primarily interested inχχ production associated with on-shell Z → jj /ll, we shall use the condition on the total visible energy (18) as the Selection cut for the rest of our analysis. The signatures where a pair of leptons appear in the final state due to the alternative decay modes of Z boson, the differential cross-section distributions w.r.t. the observables ∆φ l + l − and E l + + E l − are similar to the observables ∆φ jj and E j1 + E j2 which are obtained from the respective E T + 2 jets. Therefore, we adopt the similar selection cut criteria on the E T and E l + + E l − .
Based on the study of 1-D differential cross-section distributions at the given luminosity we select the two most sensitive kinematic observables say E T and either ∆φ jj or ∆φ l + l − for the further analysis of two dimensional differential cross-section distributions. We analyse the efficiency of the Selection cut (18) through the double distributions of the simulated Z associated DM pair production events (for m χ values at 75, 225 and 325 GeV), which are then compared with the double differential cross-section distributions of the SM at the given integrated luminosity.
In Figure 6, we exhibit the normalized two dimensional differential cross-section with the Selection cuts in the E T and ∆φ jj plane for the background and the DM with m χ = 325 GeV. These 3D lego plots are generated from MadAnalysis 5 [62] at √ s = 1 TeV, Λ =1 TeV and g q SS = 1. The bin width is 15 GeV for the E T while it is 0.1 for the ∆φ ii distribution.
We observe that the suppression of background processes due to the implementation of selection cuts enhances the sensitivity of the signal. We study the signal efficiency (S/ √ S + B) at √ s = 1 TeV with the integrated luminosity 1 ab −1 and in the Table 3 we compare the 3σ reach of the cut-off Λ at a given mass m χ for the dominant signature e + e − → E T + 2 jets corresponding to the three representative masses of DM.

χ 2 Analysis
In this subsection we study the sensitivity of the cutoff scale Λ with the DM mass keeping the couplings of the DM bilinears with the SM fermionic pair to be unity. The χ 2 analysis is performed with the sum of the variance of the differential bin events at a given integrated luminosity due to the presence of the new physics (NP) contribution over the SM events.
We define the χ 2 for double distribution as 19) where N N P ij are number of differential events given by New Physics and N SM +N P ij are total number of dif-  Table 1. The contours at the left correspond to the unpolarized beam of e − and e + with an integrated luminosity of 500 fb −1 . The contours at the right correspond to 80% and 30% polarized beam for e − and e + respectively with an integrated luminosity of 250 fb −1 .  Table 1. The contours at the left correspond to the unpolarized beam of e − and e + with an integrated luminosity of 1 ab −1 . The contours at the right correspond to 80% and 30% polarized beam for e − and e + respectively with an integrated luminosity of 500 fb −1 .  Table 1. The contours at the left correspond to the unpolarized beam of e − and e + with an integrated luminosity of 500 fb −1 . The contours at the right correspond to 80% and 30% polarized beam for e − and e + respectively with an integrated luminosity of 250 fb −1 .  Table 1. The contours at the left correspond to the unpolarized beam of e − and e + with an integrated luminosity of 1 ab −1 . The contours at the right correspond to 80% and 30% polarized beam for e − and e + respectively with an integrated luminosity of 500 fb −1 .
ferential events for (ij) th grid in double distribution. δ sys is the total systematic error in the measurement. Although the systematic uncertainty inclusive of the luminosity uncertainty is considered to be of the order of 0.3% in the literature [61], We have considered a conservative total systematic error to be of order of 1%. Further, we find that the sensitivity of the cut-off Λ can be increased with the increase in the luminosity as the cut-off scales ∼ L 1/8 for a given χ 2 . Thus, for √ s = 1 TeV and Luminosity L = 1 ab −1 the sensitivity of the signal will improve by ≈ 10 %.
From the χ 2 analysis, we obtained the 3 σ contours in the m χ − Λ plane corresponding to all the three pro-cesses in Figures 7 and 8. We present the estimation of the 3σ sensitivity reach of the maximum value of the cut-off Λ max. that can be probed in ILC at √ s = 500 GeV and 1 TeV with an integrated luminosity L = 500 fb −1 and 1 ab −1 respectively in Table 3.
We compare our results for the unpolarized beams at √ s = 500 GeV and L = 500 fb −1 which are displayed in Figure 7 with Figures 4d and 4e of reference [56] corresponding to pseudo-scalar and axial-vector couplings respectively. We find that our analysis can provide a better kinematic reach ( m χ = 20 GeV ) for the pseudo-scalar coupling. We can probe Λ max. upto 0.98 TeV, which is further improved to 1.07 TeV with the  Table 3: Estimation of 3σ sensitivity reach of the maximum value of the cut-off Λ max. that can be probed in ILC at √ s = 500 GeV and 1 TeV with an integrated luminosity L = 500 fb −1 and L = 1 ab −1 respectively, through all visible channels of e + e − → Z + χχ. enhanced integrated luminosity of L = 1 ab −1 . However, for the case of axial-vector couplings we can only probe Λ max. upto 0.88 TeV as compared to 0.92 TeV given in Figure 4e of [56].
On the same note, we would like to compare our results for √ s = 1 TeV and L = 1 ab −1 which are displayed in Figure 8 with Figures 4d and 4e of reference [56] corresponding to the pseudo-scalar and axial vector respectively. We find that the sensitivity for the pseudo-scalar can be improved Λ max. = 2.62 TeV in comparison to 2.3 TeV given in Figure 4d of [56]. Similarly, for axial-vector couplings, we agree with Λ max. given in Figure 4e of [56].

Effect of Polarized Beam
We further study the sensitivity of Λ dependence on m χ with the polarized initial beams. The rate of pair production of the fermionic DM through Scalar (SS) and PseudoScalar (P P ) interactions can be enhanced by increasing the flux of right (left) handed electrons and right (left) handed positrons. Similarly V V and AA interaction of the fermionic DM can be enhanced by choosing the right(left) handed electron beam and left(right) handed positron beam. The background contribution from t-channel W exchange diagram can be suppressed significantly by choosing right handed polarized electron beam which then leaves us with the following choices of polarization combination w.r.t. helicity conserving(V V ,AA) and helicity flipping (SS,P P ) interactions respectively 1. + 80% e− and − 30% e+ 2. + 80% e− and + 30% e+ We exhibit the 99% C.L. contours corresponding to the polarized initial e − and e + beams in Figures 7 and  8 in the m χ -Λ plane for all possible visible signatures associated with DM pair production at √ s = 500 GeV and 1 TeV with an integrated luminosity of 250 fb −1 and 500 fb −1 respectively. We find that, for scalar and pseudoscalar (vector and axial vector) DM couplings, we can probe Λ max. values upto 1.1 TeV (0.89 TeV) with polarized initial beams at an integrated Luminosity of 250 fb −1 . Similarly, we observe that the analysis with √ s = 1 TeV and an integrated luminosity of 1 ab −1 the sensitivity of Λ max. can be improved to 3.13 TeV (2.05 TeV) for scalar and pseudoscalar (vector and axial-vector) DM couplings.
We perform the combined χ 2 analyses from all the three processes involving the DM pair production for both unpolarized and polarized initial beams. We compute the 99% confidence limit contours and plot them in Figures 9 and 10 for √ s = 500 GeV and 1 TeV respectively.

Comparison with the existing analysis
We present a comparison of our results with those that exist in the recent literature.

Sensitivity from mono-photon channel
First, we make an overall comparison with the other complementary dominant DM production mono-photon channel studied in the reference [41]. We have calculated and verified the cross sections and significance of various interactions in the mono-photon channel with the cuts mentioned in [41]. Our analysis shows that the sensitivity of Λ is higher than that of the mono-photon channel, when mediated by the scalar coupling of the DM, specially in the lower m χ region and behave as competitive DM production channel upto DM mass ≈ 300 GeV. It is important to mention that the enhancement in the sensitivity are obtained inspite of our conservative input for the systematic uncertainty ≈ 1% in contrast to that of ≈ 0.3% in the reference [41].

Unpolarized
Polarized  Table 4: Efficiency S/ √ B of the process e + e − → j j+ E T for m χ =10 GeV and Λ = 1 TeV at √ s = 500 GeV and an integrated luminosity L = 100 fb −1 , with unpolarized and polarized initial beams.

Results from mono-Z channels
Next, we compare our results with the existing analysis on the DM production channels in association with Z boson which decays to pair of light quarks, electrons and muons.
DM interactions mediated through the vector and axial vector couplings are compared with the results of the reference [55] in Table 4. We marginally improve the upper limit Λ max ( at m χ = 10 GeV) corresponding to the unpolarized initial beams. However, it should be noted that the authors of reference [55] have assumed α −1 em to be 137 at √ s = 500 GeV or 1 TeV instead of 127.9, thus leading to the gross under-estimation of the background events. Rectifying the error will further bring down the efficiency of the sensitivity of their analysis.
We also compare the significance of the mono-Z channel processes in the hadronic decay mode with those given in reference [56] at √ s = 500 GeV. We summarize this in Table 5 for some representative values of m χ and Λ. We find that the our cuts strategy and analysis shows enhancement in the signal significance w.r.t. the pseudo-scalar and axial-vector interactions of DM with charged leptons for both unpoalrized and polarized initial beams.

Summary and Outlook
The fact that a massive O(100 GeV) stable particle χ with near-weak scale interactions can constitute an at-  Table 5: Efficiency S/ √ S + B of the process e + e − → j j+ E T at √ s = 500 GeV and an integrated luminosity L = 100 fb −1 corresponding to the benchmark points mentioned in the reference [56] for both unpolarized and polarized initial beams.
tractive candidate for the Dark Matter in the universe (while satisfying all of competing constraints from the relic density, large-scale structure formation as well as myriad other cosmological and astrophysical observations) is well-known. While the correct relic abundance can be reproduced by ensuring the right annihilation cross-sections into any of the SM particles, direct detection of the DM candidate largely hinges on its having unsuppressed interactions with the first-generation quarks.
It, then, is of interest to consider the possibility that χ has very suppressed (if any) interactions with quarks and gluons. Not only would the direct-detection experiments have very little sensitivity to such particles, their pair-production cross-section at the LHC would be suppressed too. In such circumstances though, the DM must have significant couplings to at least some of the leptons and/or the electroweak gauge bosons.
In this article, we have investigated the first of the two possibilities. We obtained the upper bound on the cut-off Λ for the given range of electrophilic and/ or leptophilic DM mass m χ from the relic density constraints which is displayed in the Figure 1. Using these upper bounds on Λ for a given m χ , we estimated the thermally averaged annihilation indirect detection crosssection and leptophilic DM direct detection scattering cross-section and are displayed in Figures 3 and 2 respectively. We find that the present experimental limits in the respective searches not only favours the allowed parameter space from the relic density, but also constraints the DM model by providing the lower bound on the Λ for a given DM mass m χ .
If the couplings to the charged leptons be unsuppressed, the linear collider could play an interesting role in establishing the existence of such a χ. Using the allowed coupling to the electron, this has, traditionally, been undertaken using the mono-photon (accompanied by missing energy-momentum) channel. We have investigated, here, the complementary channel, namely e + e − → χχ + ff , where f is any of the light charged fermions (jets, electrons and muons) and exhibited the significance of these processes with basic kinematic cuts in Table 2.
We analysed the sensitivity of DM scalar and vector couplings through the one and two dimensional normalised differential kinematic distributions corresponding to the most dominant signature e + e − → j j+ E T in Figures 4, 5 and 6 respectively. 99 % C.L. contours for various sets of run parameters based on the χ 2 analysis of these differential distributions are shown in Figures  7, 8. Combining the χ 2 analysis, from all the three processes corresponding to √ s = 500 GeV and 1 TeV further improves the sensitivity of the cut-off for given DM mass and are shown in Figures 9 and 10 respectively. We find that specific choice of the initial beam polarisation enhances the sensitivity of the cut-off for the scalar and vector couplings of DM. Our analysis shows that the proposed ILC can probe the cosmologically allowed m χ − Λ region for the leptophilic and/ or electrophilic DM with higher sensitivity and can provide a direction towards the understanding of the nature of such DM.
We hope this study will be useful in studying the physics potential of the ILC in context to dark matter searches.