Characterizing Higgs portal dark matter models at the ILC

We study the Dark Matter (DM) discovery prospect and its spin discrimination in the theoretical framework of gauge invariant and renormalizable Higgs portal DM models at the ILC with $\sqrt{s} = 500$ GeV. In such models, the DM pair is produced in association with a $Z$ boson. In case the singlet scalar DM, the mediator is just the SM Higgs boson, whereas for the fermion or vector DM there is an additional singlet scalar mediator that mixes with the SM Higgs boson, which produces significant observable differences. After careful investigation of the signal and backgrounds both at parton level and at detector level, we find the signal with hadronically decaying $Z$ boson provides a better search sensitivity than the signal with leptonically decaying $Z$ boson. Taking the fermion DM model as a benchmark scenario, when the DM-mediator coupling $g_\chi$ is relatively small, the DM signals are discoverable only for benchmark points with relatively light scalar mediator $H_2$. And the spin discriminating from scalar DM is always promising while it is difficult to discriminate from vector DM. As for $g_\chi$ approaching the perturbative limit, benchmark points with the mediator $H_2$ in the full mass region of interest are discoverable. And the spin discriminating from both the scalar and fermion DM are quite promising.


I. INTRODUCTION
Since the existence of Dark Matter (DM) is confirmed by many astrophysical observations [1], identifying the properties such as their masses and spins and couplings of the DM is one of the most important tasks of particle physics. The most often considered thermal DM candidate is the weakly interacting massive particle (WIMP), which has the mass around O(100) GeV and interacts with Standard Model (SM) particles via the electroweak force. Thus it can be produced directly in collider experiments in principle. The DM signal at colliders can be probed as the momentum imbalance at the detector if it is produced with recoiling against visible objects. Probing the DM signals at colliders could elucidate the particle physics properties of DM without suffering from astrophysical uncertainties thus becomes one of the main object of the current and future colliders.
There are three theoretical frameworks [2] that are used for describing the DM phenomena at the colliders, each has its own advantages and limitations: • The DM effective field theory (EFT) [3][4][5] is the low energy approximation of a renormalizable theory after integrating out the heavy particle that mediates the DM-SM particles interactions. The number of free parameters in the EFT is minimal, i.e. only two parameters are relevant for each operator, the coefficient of the effective operator and the DM mass. However, the EFT descriptions of DM interactions are valid only when momentum transfer is much smaller than the mass of the mediator such as in DM direct detection experiments. While in the collider experiments, where the momentum transfers can be quite high, the kinematic distributions that are predicted by the underlying UV completion are not correctly captured by the EFT [6][7][8][9][10], especially in the region with light mediator or heavy DM.
• In the DM simplified model [11][12][13], the DM is neutral under the SM gauge group and interacts with the SM particles via the portal of a single particle. Models in this class usually contain 5 free parameters: DM mass m χ , mediator mass m med , DM-mediator coupling g χ , SM particle-mediator coupling g SM and the mediator decay width Γ med . Considering mediators of different masses makes it possible to consider different kinematic distributions that cannot be mapped onto effective operators, thus providing a more general framework for describing the DM phenomena. However, simplified DM models with a single scalar mediator often violate the SM gauge symmetry [14,15] and perturbative unitarity, thus may become invalid for describing some sort of UVcomplete models.
• There are growing interests in second generation simplified DM model that respects the SM gauge symmetry and renormalizability [14][15][16][17][18][19][20][21][22][23][24][25]. Among them, the simplest ones are singlet DM extension of SM with Higgs portal. In these models, depending on the DM spin, the gauge invariant DM-SM interactions may require at least two mediators. Even though in the parameter region where only the contribution of one mediator dominates, the prediction of this model coincide with that of the simplified model with a single mediator. We have shown in Ref. [15] that the interference effect between the two mediators can affect the exclusion bounds considerably in some parameter space 1 . While the models in this class give more realistic predictions regarding a UV completion, there will be ad hoc constraints from many experiments, which may be quite specific and only applicable to certain UV completed models. For example, for singlet fermion DM extension of SM, the constraints and the prediction of the model with h SM + singlet scalar portal is quite different from those of the model with two-Higgs-doublet portal.
All the above frameworks have been widely used in studying DM phenomenology at colliders. Each case contains quite a lot of possible operators/models that describe the nature of DM and its couplings. If any excess is observed within a given theoretical framework, it will be important to ask which operator/model provides the best description, i.e. characterize the DM properties. There are several studies [26][27][28] devoted to distinguishing the DM EFT operators and its spin at the LHC. In the framework of DM simplified model with single mediator, many current works [29][30][31][32][33][34][35][36][37] are mainly focused on distinguishing the spin of the mediator and identifying the coupling forms between the mediator and SM particles. Because here the DM are dominantly produced by the decay of the on-shell mediator, those visible final states does not carry any information of the DM nature.
In this work, based on the gauge invariant and renormalizable DM models with Higgs portal, we will study the fermion DM (FDM) discovery prospects and its spin discriminations against scalar DM (SDM) and vector DM (VDM) at the ILC. A very preliminary study along the similar direction has been given by one of us in Ref. [38], where the detector effects were completely ignored in discussions of the DM discovery and only qualitative arguments were given regarding the spin discrimination. By curing these two problems, we find the hadronic channel of DM production provides a better sensitivity for DM discovery than the leptonic channel. Taking the FDM model as a reference model, the FDM with the coupling in a wide range can be discovered in the hadronic channel when the second mediator is relatively light. In this region, the spin discriminating from SDM is always quite promising, because the SDM model is intrinsically different from the FDM model with only one mediator being introduced. However, the spin discriminating from VDM is much more difficult, which becomes possible only in the region where the coupling between the DM and the second mediator is approaching the perturbative limit. This paper is organized as follows. In Sec. II, we briefly describe the renormalizable and gauge invariant Higgs portal DM models for scalar, fermion and vector DM. Their complete Lagrangians as well as the interaction Lagrangians that are relevant to the DM search at collider are provided. Sec. III details the analysis for the DM discovery and the strategy for the DM spin discrimination based on a benchmark scenario. Similar methods are then applied to the leptonic channel of the benchmark scenario and the hadronic channel with different couplings in Sec. IV and Sec. V, respectively. Then we conclude the work in Sec. VI.

II. HIGGS PORTAL DM MODELS
In this section, we define the Higgs portal DM models with SM gauge invariance and renormalizability, where the DMs are scalar, fermion and vector particle, respectively.
The SDM model can be constructed by simply introducing a new scalar S in addition to the SM [39][40][41] where H is the SM Higgs doublet and S is assumed to be odd under a Z 2 symmetry and thus becomes a DM candidate. After the electroweak (EW) symmetry breaking H → (0, (v h + h)/ √ 2) T and assuming S = 0, we can write down the interaction Lagrangian for DM production at the ILC as In this model, the DM can only be pair produced through the SM Higgs (h) mediation. The simplest Higgs portal singlet FDM model with SM gauge invariance and renormalizability contains a SM singlet Dirac fermion DM χ and a real singlet scalar mediator S 2 in addition to the SM particles [16,17]: where the singlet scalar S can not have direct renormalizable couplings to the SM particles due to the SM gauge symmetry and the singlet Dirac fermion χ is assumed to be odd under a Z 2 dark parity χ → −χ. When both scalar fields H and S develop nonzero vacuum expectation values (VEV), v h and v s , so that The interaction Lagrangian of interest can be written in the mass eigenstates as In contrast to the SDM model, there are two scalar bosons that mediate the DM production in the fermion DM model. The interference effects between two mediators can lead to interesting applications to DM searches at colliders [14,15]. If the H 1 is assumed to be the 2 Here the singlet scalar S is different from the singlet scalar DM defined in Eq. (II.1), although we use the same notation. In the FDM case, there is no Z 2 symmetry (S → −S) so that S cannot be a DM candidate, and S is a messenger between the dark sector and the SM sector through the Yukawa coupling (y χ -term) in Eq. (II.3).
As for constructing a renormalizable and gauge invariant model for vector (VDM), we need to introduce an abelian dark gauge group U (1) X and a dark Higgs field Φ [23,49]: where the U (1) X charge of Φ will be taken as Q Φ = 1 throughout the paper. In this model, a Z 2 symmetry (V µ → −V µ ) and charge conjugation symmetry have been imposed by hand, thereby forbidding the kinetic mixing between V µ and the SM U (1) Y gauge boson and making the vector boson V µ stable. It can also be implemented by some unbroken local dark gauge symmetry as proposed in Ref. [50].
Similarly to the FDM model with Higgs portal, there are two scalar mass eigenstates (H 1/2 ) that are originated from the mixing of SM Higgs h and dark Higgs φ, with the mixing angle given by Then, the interaction Lagrangian that is relevant to the collider study can be written as So far we have derived the relevant interaction Lagrangians for scalar, fermion and vector DMs with Higgs portal in Eqs. (II.2), (II.7), (II.10) respectively. Note that there is only one scalar mediator (h) in the scalar DM model, while there are two scalar mediators (H 1/2 ) in fermion and vector DM models. The difference in the number of mediators can lead to quite different kinematic distributions, which can be used to discriminate scalar DM model against fermion/vector DM models. On the other hand, distinguishing fermion DM models from vector DM models is more involved. First of all, if the DM production is dominated by on-shell H 1/2 production with subsequent invisible decay, it will be impossible to observe any differences in the final state distribution. The spin discrimination between fermion and vector DM is possible only if the off-shell contributions become important. Then, given the same decay width of H 1/2 , the fermion and vector DM model will predict different DM production rate as well as final state kinematics.

III. A BENCHMARK STUDY
At the ILC, the Higgs portal DM is dominantly produced through the Higgs-strahlung process where D = S, χ, V µ for scalar, fermion and vector DM, respectively. The Z boson can decay either leptonically or hadronically. We will show later that the leptonic mode which is suffering from branching ratio suppression is less sensitive than the hadronic mode. In this section, we will focus on the discovery prospect of the hadronic mode of fermion DM and discuss its spin discrimination against vector/scalar DM. Note that for scalar DM, only one mediator H 1 = h SM is introduced. In order to guarantee sufficient DM production rate at colliders while consistent with current measurements, the relevant parameters for the fermion DM production are chosen as Four benchmark points with different m H 2 = (200, 300, 400, 500) GeV will be studied, which are denoted as FDM200, FDM300, FDM400 and FDM500, respectively. For each benchmark point, we assume that the decay width for heavier scalar H 2 into the H 1 pair is negligible 3 . Then we can express the minimal decay width of H 2 as where f is the SM fermion, V = Z, W and δ V = 1(2) for Z(W ± ).
To study the spin discrimination, the parameters for the vector DM production are chosen accordingly: and the g V is chosen such that the total decay width of H 2 is the same with that in the fermion DM case, since one can rely on other method to measure the total decay width of H 2 . We give the total decay widths of H 2 for four benchmark points in FDM model as well as the corresponding g V of VDM model in Tab  The scalar DM model is much simpler, since there are only two parameters: m S and λ HS . In studying the spin discrimination, we will fix m S = 80 GeV and take appropriate λ HS such that the number of signal events after all selections are kept the same as that of each benchmark point of the FDM model. However, changing the λ HS can only lead to total rescaling of the cross section and will not affect the differential distribution of kinematic variables. In the following, we will fix λ HS when discussing the kinematic shapes without loss of generality. The SM processes with any species of neutrino in the final state could mimic the DM signal. The dominant SM background processes to Eq. (III.1) are shown in Fig. 1. Among them, the first and the second diagram (including three species of neutrino) give similar amount of contributions, while the third one is slightly smaller. At the ILC with √ s = 500 GeV and unpolarized beams, the total production cross section including the interference effects between different diagrams is 219 fb. Since the left and right handed fermions have different electroweak charges, the background cross section, especially the contribution from vector boson fusion (VBF) process, strongly depends on the beam polarization. The ILC will be able to provide highly polarized electron beam (80%) and moderately polarized positron beam (30%) [51]. The background cross sections with respect to the varying beams polarizations are plotted in the left panel of Fig. 2, where we have used the positive sign for right handed polarization and negative sign for left handed polarization. We can see that the background cross section is largest for electron polarization P e − = −80% and positron polarization P e + = 30%, while it is smallest for (P e − , P e + ) = (80%, −30%). Meanwhile, the cross sections of signal processes also mildly depend on the beam polarization. Taking the benchmark point of FDM200 for illustration, the signal to background ratio with respect to the varying beams polarizations are given in the right panel of Fig. 2, where the values have been normalized to unit at (P e − , P e + ) = (0, 0). It can be seen that the signal to background ratio can be either reduce by a factor of 0.7 or enhanced by a factor of ∼ 3 comparing to the value at (P e − , P e + ) = (0, 0). Although polarized beams improve the sensitivity, we report the results with the unpolarized beam in this work. In this work, the cross sections and events for signal and background are generated by MadGraph5_aMC@NLO_v2.4.3 [52]. The Pythia6 [53] is used for parton showering and hadronization. The final state jets are clustered using the Fastjet [54]. We also include the detector effects by using Delphes_v3.4.1 [55] with input of ILD card [56]. The track momentum and calorimeter energy resolutions of the card are listed in Tab. II. It should be noted that a more realistic detector simulation should also consider the energy spectra of income beams, the effect of which is neglected in our simulation.

A. Features of DM spin
For our signal processes at the ILC, the 4-momentum of the DM pair system can be solved as where the √ s is the collision energy and E Z ( p Z ) is the energy (momentum) of the Z boson. Therefore the invariant mass of the DM system is an observable at the ILC: The differential cross section with respect to m 2 DD for scalar, fermion and vector DM production have been calculated in Ref. [38]. It can be factorized as an off-shell mediator production and decay: where t ≡ m 2 DD and D = S, χ, V for scalar, fermion, vector DM respectively. The off-shell mediator production cross section is universal for all DM spins. In above equation, being the weak mixing angle is the averaged spin factor for initial electron and positron; which is different from spin to spin shows the spin dependent behaviour of the differential cross section: We can see from above that different DM spins can lead to different collision energy √ s dependence of the production cross sections and different distributions of the DM pair invariant mass m DD . Especially the threshold behaviors (t 4m 2 DD ) or the large-t bahaviors clearly depend on the DM spin. In Fig. 3, we show the DM total production cross section in SDM, FDM and VDM models by integrating over t in Eq. (III.7). The cross sections of benchmark points in FDM and VDM increase faster than that in SDM, due to the contributions from the second mediator. Comparing FDM and VDM, we can find that the VDM has slightly larger cross section than FDM when the m H 2 200 GeV, while it can have much smaller cross section for heavy H 2 . The differences are largest when the collision energy is relatively small √ s ∼ [400, 500] GeV. In the following discussion, we will study the collider phenomenology with fixed √ s = 500 GeV, so that FDM and VDM may possibly be distinguished by their production rate directly. For fixed √ s, a powerful spin discriminator at the ILC is the invariant mass of the DM pair m DD ≡ √ t. We plot the m DD distributions for signals with different DM spins as well the background both at parton level (left panel) and at detector level (right panel) in Fig. 4.
At parton level, the m DD for SM background corresponds to the invariant mass of the neutrino pair in the final state, since they will mimic the missing energy from the DM pair at detector level. As we have discussed before, there is a large fraction of background events in which the neutrino pair is produced from Z boson decays. Thus the m DD will show a sharp peak at m Z which is a SM background. The m DD is usually quite large for the VBF background process (first panel in Fig. 1), which gives another broad peak at m DD ∼ 400 GeV. In the SDM model, the DM with m S = 80 GeV is pair produced through the off-shell SM Higgs mediation. The m SS will peak at 2m S and decrease as 1/m 4 SS with increasing m SS . In FDM and VDM models, there is another resonant enhancement at m DD ∼ m H 2 because of the existence of the additional scalar mediator, especially when the mass of H 2 is relatively  light and decay width of the H 2 is small. This explains the clear peaks for FDM200 and VDM200. The peaks become much broader for m H 2 = 300 GeV since the decay width of H 2 is large. As the on-shell H 2 production is (almost) kinematically closed for m H 2 = 400/500 GeV, the peaks no longer exist. The FDM and VDM also show distinguishable structures in the m DD distributions. When the second scalar mediator is light the VDM has more events in the small m DD region than the FDM while this becomes opposite when the second mediator is heavy. The features at parton level can be smeared out to some extent by the detector effects. First of all, the momenta of DMs/neutrinos are not observables. One can only calculate the m DD from Eq. (III.6) by using the momentum of the Z boson, which is identified as the vector sum of the momenta of two leading jets. In some cases, only one of two jets from the Z boson decay is reconstructable at the detector (p T (j) > 20 GeV and |η(j)| < 3.0). These events will be dropped. The detector level distribution of the m DD is given in the right panel of Fig. 4. We can see that the peaks are broadened and the edges get ambiguous. In particular, for the background process, the peak at Z boson mass is almost disappearing and the distribution of m DD is quite flat, rending the discovery of signal processes difficult. The edges for signal distributions at 2m D and √ s − m Z are less steep. Nevertheless, we are still able to observe distinguishable distributions between signal and background as well as between signals with different DM spins. These features can be used to search and characterize the signal as will be discussed in the following.

B. Discovery prospect of FDM and spin discriminating power
A signal has to be discovered with high significance before being characterized. In this section, based on the benchmark scenarios that we have set at the beginning of this section, we will study the discovery prospects of the FDM and discuss its spin discriminating power against SDM and VDM at the ILC with L = 1000 fb −1 and √ s = 500 GeV.
In the event reconstruction, leptons are required to have p T ( ) > 10 GeV, |η( )| < 2.4 4 and be isolated which means the scalar sum of transverse momenta of all particles with p T > 0.5 GeV that lie within a cone of radius R = 0.5 around the e(µ) is less than 12%(25%) of the transverse momentum of the e(µ). Jets are reconstructed from particle flow objects from Delphes using the anti-kt jet clustering algorithm [57] with a radius parameter R = 0.5.  Table III, where we have taken into account the Z boson hadronic decay branching ratio. It can be seen that the total cross section decreases quickly with increasing the mediator mass. The preselection efficiency is relatively flat (∼ 0.7 − 0.8) and is smallest when m H 2 = 400 GeV. This is because for each event, the DM pair recoil energy (E miss T ) is in inverse proportion to the invariant mass of dark matter pair (m DD ). As can be seen clearly in the right panel of Fig. 4, the distribution of m DD is hardest for FDM400, while it is decreasing with either larger or smaller m H 2 .
On the other hand, the production cross sections of the SM background before and after the preselection are found to be 219 fb and 109.1 fb, respectively, which are typically more than two order of magnitude larger than that of our signals. Such small signals can be easily 4 It would be more conventional to use variables of momentum p and polar angle θ at electron positron collider, which is, however, not supported in Delphes yet. We will follow the notation as in Delphes ILD card with selections applied to p T and η throughout the paper. It has to be noted that such a choice will not bring much differences into our final results because of the following reasons: (1) θ is simply given by θ = 2 arctan(e −η ); (2) the p T and p are highly correlated, they have similar sensitivities in signal and background discrimination. 5 E miss hidden in the background with relatively large uncertainty. One would rely on more refined cuts to improve the signal-to-background ratio as well as the signal significance.
From the left panel of Fig. 4, we know the m DD can play an important role in signal and background discrimination. Moreover, in signal processes, the DM pair is produced with recoiling against a Z boson which decays into two detectable jets. The two DM particles are flying along the similar direction. While in the background process, in particular the first and third diagrams in Fig. 1, the momenta of two neutrinos are unlike to align with each other leading to a cancellation in missing transverse momentum. As a result, both the missing transverse energy (E miss T ) and the transverse momentum of the Z boson (p T (Z)) and the leading jet (p T (j 1 )) get softened for the background, as being demonstrated in the Fig. 5. We note that the distributions of E miss T , p T (Z) and p T (j 1 ) are highly correlated: hardest for VDM200 and SDM; softest for VDM400 and SM background. Another useful and less correlated discriminator is the azimuthal angle separation between the p miss T and the momentum of the closer jets: In the signal process, the DM pair is flying around the opposite direction of an energetic Z boson, which decays to two collinear jets. The ∆φ min is distributed toward ∼ π. As for background processes, where the Z boson energy is much smaller, the ∆φ min distribution is flatter. We will adopt the BDT method [58,59] that takes into account all the above variables as well as the transverse momentum of the second leading jet (p T (j 2 )) and the invariant mass of jet pair (m jj ) in order to discriminate each signal benchmark point against the SM background. The BDT method uses a 100 tree ensemble that requires a minimum training events in each leaf node of 2.5% and a maximum tree depth of three. For each benchmark point, it is trained on the half of the preselected signal and backgrounds events and is tested over the rest of the events. To avoid overtraining, the Kolmogorov-Smirnov test [60] in the BDT training is required to be greater than 0.01. After the BDT training, one can assign a BDT response to each event, which is usually larger for signal than for background. Distinguishable distributions of BDT response for signal and background can be obtained by taking into account a large number of events. Then, a cut on the BDT distribution can help to improve the signal purity. We plot the signal-tobackground ratios (N S /N B ) and the signal significances (N S / √ N S + N B ) with varying cuts on the BDT distributions for FDM benchmark points in Fig. 6. We can see that the cuts on BDT can improve the N S /N B by at least one order of magnitude, while improvements on the signal significance are only significant for benchmark points with relatively light mediator mass.
The corresponding cut on BDT for each benchmark point in FDM model that maximizes the signal significance is given in the Table III, where we also provide the numbers of signal and background events and the signal significance after the BDT cut. We find that detections on the benchmark points of FDM200 and FDM300 can be made at 3-σ level at the ILC with collision energy of √ s = 500 GeV and integrated luminosity of 1000 fb −1 . This would allow us to perform the spin discrimination for those two benchmark points.
The procedure of the spin discrimination can be described as the following. Firstly, events are simulated and production cross sections are calculated for benchmark points in SDM model (SDM200, SDM300) and in VDM model (VDM200, VDM300). The SDM200 (SDM300) denotes benchmark point in SDM model that has the same signal yields after the event selection as the FDM200 (FDM300) and the VDM200 (VDM300) denotes the benchmark point in VDM model that has the second mediator mass of 200 (300) GeV. Next, after the event reconstruction, the same preselection cuts as for FDM are applied. The cross sections as well as the preselection efficiencies for those benchmark points are provided in the Table IV. Note that the preselection efficiencies for SDM200 and SDM300 are the same, since the only free parameter λ HS in SDM model can not change the kinematic features of the final state. Then, we apply the BDT that has been trained on the benchmark point FDM200 (FDM300) to the corresponding benchmark point SDM200 (SDM300) and VDM200 (VDM300). Finally, we apply the BDT cuts as given in the fourth row of Table III to the corresponding benchmark points in SDM and VDM model. The event numbers at L = 1000 fb −1 for those benchmark points are given in the fourth row of   The survived events are used to plot the distributions of m DD for different models. In Fig. 7, we give the 5-bin distributions of m DD after applying the BDT cut for signals of different DM spin adding to the SM background. We can observe that the m DD distributions of benchmark points in FDM and SDM model have visible difference, while that of benchmark points in FDM and VDM are almost the same. To assess the degree of difference between the benchmark points in FDM and SDM, we construct the χ 2 statistic ) is the number of FDM (SDM) signal plus background events in the i-th bin and the i runs over five bins of the histograms in Fig. 7. The δχ 2 value is compared to the χ 2 distribution with 4 degrees of freedom to calculate the p-value, which can be further transformed to the significance level (S) from a Gaussian distribution. The S for each benchmark point in SDM model is given in the fifth row of Table IV. Both benchmark points in SDM model can be distinguished from the benchmark points in FDM at significance level of more than 2-σ. We note that the number of events after the BDT cut contains not only the information of normalization of the m DD distribution but also the information of its shape, since the BDT selection used the m DD distribution. Therefore, for discriminating FDM and VDM, the significance level will be simply estimated by S = N FDM ) is the number of FDM (VDM) signal events as given in Table III (Table IV), N B is the number of background events after applying BDT cut. We find both benchmark points in VDM model can only be distinguished from the benchmark points in FDM with significance level below 1-σ.

C. DM properties of benchmark points
In this subsection, we will briefly discuss the DM relic density [1] and DM direct detection bound [61] for our benchmark points 6 . These values are calculated numerically by micrOMEGAs [63] with the CalcHep/CompHEP [64] model files that are written by FeynRules [65,66]. For all benchmark points, the DMs are dominantly annihilatting into W W * through scalar mediator(s) where W * is the off-shell W boson. Due to the relatively large couplings between the mediator and DMs being chosen, the relic abundances of our DM particles are always below the measurement (Ωh 2 0 = 0.1198) as can be seen in Tab. V, rendering our DM particle only as a component of a full DM sector. Among DM spins, the fermion DM has suppressed s-wave annihilation, thus largest relic density.
In comparison between the DM-proton scattering cross section in our model and the LUX constraint, the cross section (σ SI p ) calculated in micrOMEGAs should be rescaled by a factor of Ωh 2 0.1198 with Ωh 2 being the calculated relic density of each benchmark point. According to Ref. [61], the current LUX measurement has excluded σ SI p · Ωh 2 0.1198 > 1.4 × 10 −10 pb for m DM = 80 GeV, which means all of our benchmark points should have been excluded already. However, the direct detection limits rely on assumptions about the local dark matter density and velocity distributions, which are expected to vary from the standard assumptions used in the experimental results [67][68][69][70]. Moreover, if there is indeed a DM sector, our DM particle can either decay or annihilated into other dark particles, so that the direct detection can be evaded. It should be noted that those modifications will not lead to any effects in the collider phenomenology of DM searches.

IV. THE LEPTONIC CHANNEL
As we have seen in Fig. 4, the hadronic channel is suffering from the large uncertainty in jet momentum measurement, leading to smearing effects in the m DD distributions. On the other hand, much better lepton (e/µ) momentum resolution of the leptonic channel may help to improve the discovery sensitivity as well as the spin discriminating power.
However, the main drawbacks of the leptonic channel are its small production cross section and relatively large SM background. The Z boson in the leptonic channel is required to decay into electron or muon pair, the decay branching ratio of which is around one order of magnitude below that of hadronic mode: Br(Z → + − ) = 6.7% with = e, µ, Br(Z → qq) = 69.9% with q = u, d, c, s, b. Moreover, aside from the background processes listed in Fig. 1 with j being replaced by , there are new SM backgrounds such as the single W and W boson pair productions where the W bosons are decaying leptonically. The total production cross section of the SM process e + e − → νν is 505 fb at the √ s = 500 GeV ILC, which we find is dominated by the contributions from processes with W boson in the final state. In Fig. 8, we plot the m DD distribution for the leptonic channels of signals and background at parton level (left panel) and detector level (right panel). We can find that the shapes of m DD distributions are largely unaltered after taking into account the detector effects, i.e. peaks are sharp and edges are steep even at the detector level. Comparing to the Fig. 4, the main features of signal distributions are kept the same as that in the hadronic channel, since the two channels only differ in the Z boson decay final state. As for background, the Z peak in the leptonic channel is less notable because the processes with W in the final state are dominating. We note that in some events, only one of the two leptons in the final state is reconstructable at the detector (p T ( ) > 10 GeV and |η( )| < 2.5). Those events are corresponding to those with m DD = 0 GeV in the right panel of Fig. 8.
Events for the leptonic channel are reconstructed with the same method as adopted for the hadronic channel. The candidate signal events are selected with the following preselection cuts: (1) exactly two opposite sign same flavor leptons in the final state; (2)  GeV. The total cross sections of the leptonic channels of benchmark points in the FDM model and their preselection efficiencies are given in Table VI. We also find the corresponding preselection efficiency of background is ∼ 0.029 which is much smaller than that of signal. Nevertheless, after the preselection, the production rates of our signals are still around 2-3 order of magnitude smaller below that of the background. To increase the signal significance, we follow the similar strategy as in the hadronic channel, i.e. adopting the BDT method. The discriminating variables that are used in the leptonic channels are where the ∆r ( , ) ≡ (∆η ( , )) 2 + (∆φ ( , )) 2 is the angular distance between two leptons and ∆φ min ≡ min i=1,2 ∆φ p miss T , p ( i ) is the azimuthal angular separation between the missing transverse momentum and the closer lepton.
After training the BDT on each benchmark point in the FDM model, we can obtain the distributions of BDT response for signal and background. The cut on the BDT distributions is chosen such that the signal significance (N S / √ N S + N B ) of each benchmark points is maximized. The corresponding BDT cuts, the number of signal and background events as well as the signal significance after BDT cuts are given in Table VI. Only the FDM200 is discoverable at the ILC with √ s = 500 GeV and L = 1000 fb −1 . For all benchmark points, the signal significances of the leptonic channel are 2-3 times smaller than those of the hadronic channel.
We can also discuss the spin discriminating of the FDM against SDM and VDM for the benchmark point FDM200. The production cross sections and preselection efficiencies of benchmark points SDM200 and VDM200 are given in the second and third row of Table VII. As in the hadronic channel, the significance levels (S) of spin discriminations between FDM and SDM and between FDM and VDM are calculated with two different methods. The results are given in the fifth row of Table VII. The FDM200 can be distinguished from SDM200 at significance level of around 2-σ, while it is impossible to be distinguished from VDM200. We can conclude that the hadronic channel provides better sensitivities in both signal discovery and spin discrimination than the leptonic channel.  Table IV, with the leptonic channel instead.

V. VARYING THE COUPLING IN THE HADRONIC CHANNEL
So far, we have studied the benchmark points with g χ = 3 in the FDM model. In this section, we will survey the discovery and spin characterizing prospects of benchmark points with g χ = 1 and g χ = 10 in the FDM model, while keep sin α and m χ unchanged. For each g χ , four different choices of m H 2 = (200, 300, 400, 500) GeV will be considered. As have been done for the g χ = 3 case, the corresponding benchmark points in VDM model are chosen such that the decay widths of H 2 are kept the same as the ones in the FDM model. We note that the branching ratio of H 2 → H 1 H 1 is assumed to be negligible in calculating the decay width of a H 2 . Benchmark points in the SDM model are chosen with the criterion that the signal yields after the event selection for signal process is the same with that of benchmark points in the FDM model by tuning the free parameter λ HS .
The most important effect of changing the g χ is that the total decay widths of the H 2 become different in the FDM and VDM models. As shown in Fig. 9, for FDM and VDM, the peaks in the m DD distribution are quite sharp when the g χ = 1. Especially, when m H 2 = 400 GeV, the contribution from the on-shell H 2 is still dominating even with the small kinematic phase space. This is in contrast to the Fig. 4 where the decay width of H 2 is much wider rendering the disappearance of the H 2 peak. We note that differences in the distributions of m DD between the FDM and VDM only occur in the off-shell H 2 processes. Otherwise, it is simply the on-shell H 2 production with subsequent invisible decay, which leaves no information of DM spin in the visible products. This explains why the m DD distributions for FDM and VDM almost overlap when H 2 is light, while the difference becomes visible in the region m H 2 300 GeV where the off-shell contribution is sizable. For g χ = 10 which is close to the perturbative limit, the decay width of H 2 is so wide that the off-shell H 2 contribution is important when m H 2 200 GeV and is dominant when m H 2 300 GeV. Then, it is possible to distinguish the FDM against VDM in the full range of m H 2 . From the right panel of Fig. 9, we can also see that m DD distributions for FDM (VDM) with m H 2 300 GeV are almost identical, because the signal events are occupying the lower side of the off-shell H 2 propagator irrespective of the H 2 mass and decay width. The signals are searched with the same method as used for benchmark points with g χ = 3. We will only discuss the hadronic channel, since we have shown that it has better sensitivity than the leptonic channel. We first list the production cross sections of the benchmark points in the hadronic channel and the corresponding preselection efficiencies in Table VIII.  Compared to Table III, we can find that all benchmark points in the FDM have similar total production rate when the H 2 is relatively light. While for m H 2 300 GeV, the production cross section increases with the coupling g χ . The increase is more dramatic for heavier H 2 . Eventually, the signal production cross sections are approaching to the same value when g χ is close to the perturbative limit due to the dominance of the off-shell H 2 contribution. The preselection efficiency for most of the benchmark points are similar, i.e. between 0.7-0.8, except for the FDM400 with g χ = 1. For this benchmark point with g χ = 1, the final state particles are a Z boson (m Z = 91.2 GeV) and an almost on-shell H 2 (m H 2 = 400 GeV), rendering the kinetic energy of final states quite small, E kin ∼ O(10) GeV. The preselection condition E miss T > 50 GeV can cut out a large number of events. The same BDT method that has been used in subsection III B for benchmark points with g χ = 3 is also adopted here. The BDT is trained on the preselected events of each benchmark point with given g χ and m H 2 in the FDM model and the SM background. A cut on the BDT responses of signal and background can be applied later to improve the signal significance. The BDT cut for each benchmark point that maximizes the signal significance (N S / √ N S + N B ) is given in Table VIII. We can find that at the ILC with √ s = 500 GeV and L = 1000 fb −1 , for g χ = 1, only the benchmark points FDM200 and FDM300 can be discovered at more than 3-σ level while for g χ = 10, all of the benchmark points can be discovered with signal significance great than ∼ 8-σ.
The production cross sections and the preselection efficiencies of benchmark points in SDM and VDM models corresponding to those in FDM model with g χ = 1 and g χ = 10 are listed in Tables IX and X, respectively. For the case of g χ = 10, the benchmark points FDM200 FDM300 FDM400 FDM500 in VDM model has much larger (smaller) production cross section than those in FDM model when the H 2 is light (heavy). So that it is possible to distinguish FDM and VDM even by using the production rates of signal alone. The number of signal events for each benchmark point in the SDM and VDM model after applying the BDT cut as well their significance level S of spin discrimination are calculated with the same strategy as introduced in subsection III B. In the case of g χ = 1, we can see in Table IX that only benchmark points SDM200 model can be distinguished from FDM model with S > 3, while it is impossible to discriminate the FDM benchmark points against the VDM benchmark points. When the g χ is close to the perturbative limit, the spin discrimination is quite promising as given in Table X. The DM spin of our benchmark points with H 2 in the full mass region of interests can be identified with high significance level. Owning to the considerable difference in the production rate between the FDM and VDM, the VDM has better discriminating power against FDM than the SDM.  Table IV, but the benchmark points have been changed to those which are corresponding to benchmark points with g χ = 10 in FDM model.

VI. CONCLUSION
In this paper, we have considered DM discovery prospect and its spin discrimination at the ILC in the theoretical framework of gauge invariant and renormalizable Higgs portal DM models for the first time. The gauge invariances of the FDM model and the VDM model require another new scalar field (in addition to the SM Higgs boson) that mediates the DM and SM particles interaction, while the gauge invariant SDM model only needs one medatior, the SM Higgs boson.
Taking the FDM model with g χ = 3 as a benchmark scenario, we study the discovery prospects and spin discriminating powers of both its hadronic channel and leptonic channel at the ILC with √ s = 500 GeV and L = 1000 fb −1 . In the hadronic channel, we first employ the BDT method with input of a few discriminative kinematic variables such as the DM pair invariant mass m DD and the azimuthal angular separation between the missing transverse momentum and the closer jet ∆φ min to improve the signal sensitivity. We find the benchmark points with m H 2 300 GeV can be probed at more than 3-σ level. For those discoverable benchmark points in the FDM model, the spin discriminating against SDM can be made with 3-σ level, due to the intrinsic difference between the FDM model and the SDM model, i.e. the FDM model contains two mediators while the SDM model only gets one. However, the spin discriminating against VDM is almost impossible, with the significance level below one for all discoverable benchmark points. The leptonic channel is also considered with the similar strategy. We find that the leptonic channel has worse discovery potential than the hadronic channel. Only benchmark points of FDM model with the mediator mass m H 2 200 GeV is discoverable. As with the hadronic channel, the spin discrimination between FDM and SDM can be made while it is quite difficult to distinguish FDM and VDM.
We also survey the discovery and the spin characterizing prospects of the benchmark points in the FDM model with varying g χ . Choosing smaller g χ does not reduce the DM production cross section in benchmark points with small m H 2 much as long as the H 2 → χχ branching is dominating. Furthermore, the smaller g χ which gives narrower decay width of H 2 will increase the difference between the m χχ distributions of the FDM and the SDM models. Thus benchmark points with g χ = 1 even have better signal significances and spin discriminating powers than those with g χ = 3. As for benchmark points with g χ approaching the perturbative limit, the off-shell H 2 contribution becomes quite important, leading to the increased production rate especially for those with heavy H 2 . We find that the benchmark points with H 2 in the full mass region of interest are discoverable. The spin discriminating against both the SDM and VDM are quite promising.
It should be noted that for FDM/VDM comparison throughout the work, the bench-mark points of VDM are chosen such that the decay widths of H 2 are kept the same as the ones in the FDM model. This can be possible provided that the decay width of H 2 can be measured elsewhere. Then, the normalization of m DD distribution become an important handle for FDM and VDM discrimination. We also considered the FDM/VDM comparisons without the information of normalization and find the discrimiantions are impossible except for the cases of g χ = 10. The S calculated from Eq. III.13 are 1.07, 1.24, 1.56 and 1.48 for FDM200/VDM200, FDM300/VDM300, FDM400/VDM400 and FDM500/VDM500, respectively.