Dark-matter-spin effects at future $e^{+} e^{-}$ colliders

We discuss possibility to detect spin 0, 1 and 1/2 dark matter (DM) at future $e^{+} e^{-}$ colliders. The models considered here are simple, consistent and renormalizable field theories, that provide correct DM abundance and satisfy direct detection, indirect detection and collider constraints. The intention of this paper was to verify to what extend it might be possible to disentangle models of different DM spins by measurement of the cross section for $e^{+} e^{-} \to Z + \ldots\,$ at future $e^{+} e^{-}$ colliders. We specialize to the case of the ILC operating at $\sqrt{s} = 250~\text{GeV}$, however our results apply as well for the FCC-ee and the CEPC colliders. For each model the cross section maximized with respect to parameters was calculated and compared to the expected 95% CL cross-section limits estimated for the ILC. It turned out that near $2 m_{\text{DM}}\simeq m_{1,2}$ resonances, where $m_{1}$ and $m_{2}$ are the SM Higgs boson and a non-standard Higgs boson masses, respectively, there exist substantial regions where the models are testable. A special attention has been payed to calculation of the cross section in the region where $m_{1}\simeq m_{2}$.


Introduction
In spite of the Higgs-boson discovery at CERN's Large Hadron Collider (LHC) by the ATLAS [1] and CMS [2] collaborations, the underlying theory of fundamental interactions is still missing since the Standard Model (SM) does not provide a candidate for dark matter (DM), while its existence has been confirmed by many independent experiments (see e.g. [3][4][5][6][7][8][9][10][11]). In this project we are going to discuss minimal extensions of the SM that describe dark matter of various spins (0, 1, 1 /2) in a framework of a consistent, renormalizable quantum field theory. Even if the ultimate theory of DM will prove to be non-minimal, it is reasonable to expect that the minimal models discussed here will capture its major low-energy properties. Our intention is to verify to what extend future e + e − colliders operating near √ s = 250 GeV: the Future Circular Collider (FCC-ee) [12][13][14], the Circular Electron Positron Collider (CEPC) [15] and the International Linear Collider (ILC) [16][17][18], could be useful for detecting DM in the process of mono-Z production, e + e − → Z + · · · . Our strategy is to impose existing constraints on simple models of pseudoscalar (pGDM), vector (VDM) and fermion (FDM) dark matter and determine regions of parameters in which DM-production cross section at e + e − colliders is maximal. Then we compare the maximized predictions with the expected 95% CL cross-section limits at the ILC, assuming that it will provide a satisfactory estimate for the other colliders as well. That way we are trying to verify whether the future electron-positron colliders operating in the vicinity of √ s = 250 GeV could be used to test theories of DM. DM production at future e + e − colliders has already been discussed in the literature, see [19,20]. However, our approach has another motivation, also the models adopted here are not the same. The goal of this project is different as well.
The paper is organized as follows: after the introduction in Sec. 1, in the subsequent Secs. 2, 3 and 4 we describe the pseudo-Goldstone, vector and fermion dark matter models, respectively. Sec. 5 is devoted to the constraints on dark matter scenarios that are adopted in the paper. In Sec. 6 we calculate the cross section for the e + e − → Z + · · · process and discuss subtleties of the mass degenerate case, m 1 m 2 . The next section, Sec. 7, is to review the expected sensitivity to this process at the ILC. Sec. 8 contains our numerical results with determination of regions in the parameter space that could be tested at the FCC-ee, CEPC and ILC. In the final section, Sec. 9, we summarize our findings. In appendices we collect results concerning the Higgs boson decay widths and 2-point 1-loop scalar Green's functions.

Pseudo-Goldstone dark matter
In spite of the fact that the minimal model of scalar (spin zero) DM [21,22] assumes merely an addition of a real scalar field odd under a Z 2 symmetry, here we are going to consider a model (pGDM) that requires an extension by a complex scalar filed S. The model is in some sense very similar to vector and fermion dark matter models that will be discussed here as well, so it is worth to compare all of them. In order to stabilize a component of S we require an invariance under DM charge conjugation C : S → S * , which guarantees stability of the imaginary part of S, A ≡ Im S/ √ 2. The real part, φ S ≡ Re S/ √ 2, is going to develop a real vacuum expectation value (vev) φ S = S = v S / √ 2. 1 Therefore, φ S will mix with the neutral component of the SM Higgs doublet H, in exactly the same manner as it happens for the VDM or the FDM. In order to simplify the potential we impose in addition a Z 2 symmetry S → −S, which eliminates odd powers of S. Eventually, the scalar potential reads: with µ 2 real, as implied by the C symmetry. Note that the µ 2 term breaks the U (1) explicitly, so the pseudo-Goldstone boson A is massive. In the limit of exact symmetry, A would be just a genuine, massless Goldstone boson. Since the symmetry-breaking operator µ 2 (S 2 + S * 2 ) is of dimension less that 4, its presence does not jeopardise renormalizability even if non-invariant higher dimension operators were not introduced, see for instance [23]. Note that dimension 3 terms are disallowed by the Z 2 and gauge symmetries. In other words, we can limit ourself to dimension 2 U (1)-breaking terms preserving the renormalizability of the model. The freedom to introduce solely the soft breaking operators offers a very efficient and economical way to generate mass for the pseudo-scalar A without the necessity to introduce dimension 4 terms like S 4 or |S| 2 S 2 , and keeping the renormalizability of the model. It is also worth noticing that the Z 2 symmetry S → −S is broken spontaneously by v S and, therefore, φ S , the real part of S, is not stable, making A the only DM candidate. The scalar fields can be expanded around the corresponding generic vev's, v for H and v S for S, as follows: The global minimum of the potential with corresponding value of the potential and the scalar mass-squared matrix read: in the basis (φ H , φ S , A). Note that the third spin-zero state A does not mix with the former ones.
Conditions necessary to guarantee the asymptotic positivity of the potential and the global minimum at (v H / √ 2, v S / √ 2) with non-zero vevs will be discussed in Sec. 5.5. It is worth to notice that in the vector DM model considered in the following section, A becomes a genuine Goldstone boson (µ 2 = 0) and disappears as a longitudinal component of the massive DM vector X.
There are two mass eigenstates, h 1 and h 2 , in this model. The mass matrix (2.5) can be diagonalized by the orthogonal rotation matrix R −1 acting on the space spanned by the two CP-even scalars φ H and φ S : We assume hereafter that h 1 is the 125.09 GeV boson observed at the LHC. Moreover, since sin α appears in calculations of Secs. 5 and 6 in the second power, we will assume without losing generality that the sign of κ is chosen in such a way that sin α > 0. We choose as independent parameters of the model the set: v S , sin α, m 2 and m DM = m A . Together with v = 246.22 GeV and m 1 = 125.09 GeV this set is sufficient to determine all the 6 parameters of the potential; relevant relations will be presented in Sec. 5.5. As it will be seen later, scalar potentials in other theories discussed in this work could be also parametrized in terms of the same parameters, allowing for meaningful comparison between the models. 2 Vertices relevant for the calculation of annihilation cross section in the pGDM model have been collected in Fig. 1. Similar models have been considered in a more general context including a possibility of fast first order phase transition in [24][25][26].

Vector dark matter
The next model that we want to compare with the pGDM is the popular vector DM (VDM) model [27][28][29][30][31][32] that is an extension of the SM by an additional U (1) X gauge symmetry and a complex scalar field S, whose vev generates a mass for the corresponding gauge field. The quantum numbers of the scalar field are None of the SM fields are charged under the extra gauge group. In order to ensure stability of the new vector boson a Z 2 symmetry is assumed to forbid U (1)-kinetic mixing between U (1) X and U (1) Y . The extra gauge boson X and the scalar field S transform under the All other fields are neutral under the Z 2 . The vector bosons' masses are given by: where g and g are the SU (2) and U (1) gauge couplings, while, as in the previous model, v and v S are the vev's of H and S, respectively: ( H , S ) = 1 √ 2 (v, v S ). 3 The scalar potential for this model is given by It is easy to find solutions of the potential minimization conditions for the scalar fields: Both scalar fields can be expanded around corresponding vev's as follows The mass-squared matrix M 2 for the fluctuations (φ H , φ S ) is identical as the 2 × 2 block of the mass matrix for the pGDM model (2.5), so that the diagonalization 2.6 and relation (2.7) remain applicable.
Conditions for existence of non-zero vevs, globality of the minimum and asymptotic positivity of the potential will be discussed in Sec. 5.5. The input parameters adopted here are: v S , sin α, m 2 and m DM = m X .
Vertices relevant for the calculation of annihilation cross section in the VDM model have been collected in Fig. 2. Figure 2. The vertices relevant for the VDM model.

Fermion dark matter
In the case of minimal fermion DM, the gauge group remains the standard one, i.e.
The DM candidate χ (left-handed Dirac fermion) is introduced together with a real scalar S that is necessary to mediate DM interaction with the SM. The extra states are charged under Z 4 : S → −S while χ → iχ. The resulting symmetric Lagrangian reads: where χ c ≡ −iγ 2 χ * . Note that the above potential is the same as in the VDM case (see 3.4), up to normalization of singlet mass and its couplings. The positivity conditions of the potential remain, of course, the same for this model as for the previous two since all the potentials have the same asymptotic behavior.
We parametrize fluctuations of scalar fields as follows: with v and v S being the vev's of the neutral component of the doublet H and the singlet S, respectively, determined by (3.5).
After SSB, relevant parts of the Lagrangian take the following form: Here, as in the models discussed earlier, there are two physical (mass eigenstates) scalar degrees of freedom, h 1 and h 2 , that are linear combinations of φ H and φ S . Note that because of appropriate normalization of terms involving S in the potential (4.2) the mass matrix and its diagonalization remain the same as in the other models. It is convenient to use the analogous input parameters to discuss this model: v S , sin α, m 2 , and m DM = m ψ .
Positivity and minimization conditions for this model will be discussed in Sec. 5.5.

Astrophysical and other constraints
Hereafter we are going to allow for resonant DM annihilation process, so we will adopt the Breit-Wigner propagators for mediating particles, i.e. the Higgs bosons h 1,2 . Γ 1,2 will denote the total width of h 1,2 , respectively.

Dark matter abundance
The thermally averaged cross section for DM annihilation into a SM fermion-anti-fermion pair, σ(DM DM →f f ), reads: 4 with n c = 1(3) for f being lepton (quark) and the variable X defined 5 as The DM abundance observed by the Planck Collaboration [10], constraints the annihilation cross section at the freeze-out temperature by σv freeze out = (n + 1) · 2.2 · 10 −26 cm 3 s −1 = (n + 1) · 1.9 · 10 −9 GeV −2 , what corresponds to the current value of annihilation cross section equal to σv now = (T 0 /T f ) n · σv freeze out = (T 0 /T f ) n · (n + 1) · 1.9 · 10 −9 GeV −2 , where T 0 is the present CMB temperature while T f ∼ m DM /25 is temperature at the moment of freeze out. Value of n is 0 for the bosonic models (pGDM, VDM) and 1 for the FDM.
Hence, keeping only the leading (bb) contribution in eq. (5.1), we obtain the following constraint .
(5.6) 4 Other final states are not accessible kinematically for mass ranges considered here. 5 Note that at the tree level X reduces to κ 2 .

Dark matter indirect detection
Since we fix the DM abundance to its observed value (5.3), the present annihilation cross section is also fixed by (5.5), so that it remains to be a function of m DM only. Therefore, the limit on the present annihilation cross section, for instance from Fermi-LAT [33], implies a lower limit on DM mass. For the bosonic models (pGDM, VDM), adopting data for the bb final state, one obtains m DM > ∼ 20 GeV. Hereafter, we will consider this region only. In the case of the FDM, the extra suppression by T 0 /T f implies that the cross section is by a factor of 10 −11 − 10 −13 smaller than for the bosonic models and, therefore, there is no constraint on m ψ .

Dark matter direct detection
The DM direct detection (DD) experiments impose severe constraints on the parameter space of DM models. In the models discussed here the spin-independent cross sections for the DM-nucleon scattering are given by where m N denotes nucleon mass and µ is the reduced mass for the DM-nucleon system while for the form factor we have adopted f N 0.3 GeV. Widths and momentum transfer in the denominator have been neglected as much smaller than masses. It turns out that in the pGDM model the cross section vanishes [34][35][36] in the limit of zero momentum transfer, so 1-loop calculations are needed. The 1-loop results are encoded above through the factor containing A, defined according to [35] 6 as where the functions C and D are defined in Appendix B. In eq. (5.9), the sign of sin α is relevant, what seems to contradict our statement that chosing sin α > 0 does not spoil generality of our considerations. However, the only place where we use results of eq. (5.7) for the pGDM is the comparison in Fig. 4. Regardless of the sign of sin α, the conclusion that X (DD) is orders of magnitude larger than X (Ω DM 0 ) remains true, and hence, we do not have to consider the sin α < 0 case separately.
For practical purposes, the XENON1T limit [37] for m DM 40 GeV can be parametrized as follows Hence, X is constrained from above by DD limit: .
It turns out that in considered range of parameters, in the case of the pGDM, the DD upper bound on the value of X is always higher than the value corresponding to the correct relic density, see Fig. 4. Therefore, in the case of the pGDM, the DD constraint does not limit the range of (m 2 , m DM ). . Comparison between the DD upper bound for the value of X (denoted by X (DD), see eq. (5.11)) and the value providing correct relic density (denoted by X (Ω DM 0 ), see eq. (5.6)), in the case of the pGDM. Since the upper bound is always higher than the required value, the DD constraint does not affect the range of (m 2 ,m DM ) in the case of this model.

Collider constraints
The mixing angle α is constrained from the measurement of the SM signal strength µ LHC . The latest LHC bound is µ LHC = 1.09 ± 0.11 which amounts to sin 2 α < 0.13 at the 2σ CL [38]. Hereafter we will adopt a bit stronger limit sin α < 0.30.
When the DM mass is smaller than half of the SM-like Higgs boson h 1 , m DM < m 1 /2, the Higgs invisible decay provides another constraint on DM scenarios. In the models discussed here, the width for invisible decays are as follows .
Current LHC measurements [39] provides the following limit on the invisible branchig ratio: at the 95% CL.

Theoretical constraints
In order to ensure that the leading order calculations adopted here are meaningful, we impose the following perturbativity conditions on the U (1) X gauge coupling in the VDM model and the Yukawa coupling in the FDM model: g X < 4π and y x < 4π. Both of them correspond to v S > m DM 4π . In the pGDM model, the AAh i coupling is proportional to m 2 i /v S (cf. Fig. 1), therefore we also require m i /v S < 4π (i = 1, 2). It is interesting to note that there exist regions (e.g. m 2 ∼ m 1 ) in the parameter space where the proper abundance of DM requires small v S . In these regions some quartic couplings might be too large (nonperturbative), since λ S ∝ m 2 i /v 2 S and κ ∝ (m 2 1 − m 2 2 )/(vv S ), see Fig. 11. Therefore we also impose the conditions: λ S , |κ| < 4π. Summing up, the conditions adopted here in order to ensure perturbativity within considered models are Let us now consider conditions for stability of the vacuum state. Scalar potentials of the models read (see eqs. (2.1), (3.4) and (4.2)): To ensure asymptotic positivity of all the potentials, the following conditions must be satisfied: Vacuum expectation values of the scalar fields H, S are denoted as follows In each case, v, v S = 0 must minimize the value of the potential. The corresponding point in the (H, S) space is a critical one if and only if in the case of the pGDM and in the case of the VDM and the FDM.
To ensure that the critical point is a strict minimum, we demand the second derivative of the potential to be positive definite, therefore .

(5.24)
Hence, assuming v 2 , v 2 S and λ H are positive, the following condition must hold Positivity of the vev's squared requires (cf. eqs. (5.21) and (5.22)) in the case of the pGDM and in the case of the VDM and the FDM. Let us check when the points given by eqs. (5.21) and (5.22) are global minima. In the case of the pGDM, due to the presence of the µ 2 (S 2 + S * 2 ) term, in principle the phase of the vacuum expectation value of S could be relevant. Hence, let us assume that S = (v S + iv A )/ √ 2. Now, we have to minimize the potential with respect to v, v S and v A . There are six critical points of the potential, namely Assuming the asymptotic positivity conditions (5.18), the strict-minimum condition (5.25) and positivity of m 2 DM = −4µ 2 , minimum (5.33) is always smaller than (5.29) and (5.30). To ensure that minimum (5.33) is smaller than (5.32), the following additional condition must hold: Value of (5.31) is obviously greater than (5.30) and, therefore, greater than (5.33) if Both of these conditions are checked for considered region of parameter space at the end of this subsection.
In the case of the VDM and the FDM, we can assume that S is purely real without losing generality. Therefore, V is minimized with respect to v and v S . The critical points are We can express parameters of the potential in terms of the input parameters: m 1 , m 2 , v, v S , m DM and sin α as follows: (5.43) It appears that the stability and positivity conditions (5.18), (5.25) and (5.26) expressed in terms of the input parameters are automatically satisfied: In fact, our choice of the input set implicitly assumes that coefficients of V are such that v 2 , v 2 S , m 2 1 , m 2 2 > 0. The global-minimum conditions (5.34) and (5.35) for the case of the pGDM are expressed in terms of the input parameters as follows: It can be numerically shown (see Fig. 5) that in the considered range of masses theses conditions are always satisfied. 6 Production of DM pairs at future e + e − colliders The DM models can be tested at e + e − collider experiments. In particular, these experiments allow for the copious production of DM states associated with a Z boson, what is referred to as so called Higgsstrahlung process or mono-Z emission [19,20,[40][41][42][43][44], see  Figure 6. Feynman diagram for considered channel of DM production. DM denotes the dark particle that is either A, X or ψ.
We assume that the energy of the Z boson can be reconstructed from data, therefore allowing for determination of the recoil mass ( Q 2 ), corresponding to invariant mass of the dark particles.
The differential cross section for DM pair production at e + e − colliders reads where the parameter X is defined in eq. (5.2) and is the cross-section for the e + e − → Zh SM process, with mass of the SM Higgs particle equal to Q 2 . Here, λ(a, b, c) denotes the Källén function, defined in Appendix A, and g V , g A stand for the vector and axial coupling, respectively. 7 The above result has been obtained by adopting the standard Breit-Wigner propagators for the virtual/real Higgs bosons h i .
Note that in the limit of m 2 → m 1 the cross section dσ/dQ 2 (6.1) seems to be amplified for h i being on-shell, i.e. Q 2 → m 2 1,2 . This is a surprising observation since, on the other hand, the second relation in (5.40) between masses and the portal coupling κ implies that in the limit m 2 → m 1 , whenever v S = 0, κ → 0 so that the dark sector decouples in each model discussed here. Therefore, all cross sections for DM production or annihilation from the SM must vanish in this limit. This notice definitely deserves a careful explanation.
Let's investigate the parameter X . First, it is easy to see that lim m 2 →m 1 X = sin 2α From (5.40) one finds that if v S = 0 then the limit m 2 → m 1 implies κ → 0 and λ H v 2 − λ S v 2 S → 0. Therefore, according to (2.7) tan 2α is undefined. For instance, for fixed λ H , v and v S it is easy to see that, approaching the limiting point (λ H (v/v S ) 2 , 0) in the (λ S , κ) plane, one can get α = 0, α = π/4 or α = 1/2 arctan(v/v S ), choosing the corresponding trajectories: κ = 0, λ S = λ H (v/v S ) 2 or κ = −λ S + λ H (v/v S ) 2 , respectively. Since in the limit m 2 → m 1 neither sin α → 0 nor Γ 1 → Γ 2 , so X does not vanish, in spite of justified arguments mentioned above. The solution to this puzzle lies in the fact that for m 2 → m 1 also off-diagonal (i = j) Higgs boson self-energies are relevant and should be resummed so the naive, diagonal, Breit-Wigner propagators are not appropriate. To illustrate this point let us consider the e + e − → ZXX process. The matrix element reads In the case of polarized beams, Figure 7. The Higgs-boson mediators with their self-energies.
By ∆ we denote the propagator contracted with the mixing matrix. From [45] (see also [46]), the contracted propagator can be expressed explicitly as In magnitude, all of them are comparable to m i Γ i . Results for Π ij are collected in Appendix A.

(6.6)
It is easy to see that the above propagator could be rewritten (dropping terms proportional to Γ 1,2 in the numerator) as which reduces to the standard Breit-Wigner propagator. This simplified result has to be replaced by the full formula only when |m 1 − m 2 | is comparable to the widths. In order to investigate the case m 1 ∼ m 2 one has to calculate Π ij . The explicit calculation (see Appendix A) confirms that Hence, the full propagator (6.5) vanishes in the limit m 1 = m 2 , exactly as it should. However, in practice, the region |m 1 − m 2 | < ∼ Γ 1,2 is so narrow that the naive Breit-Wigner approximated resummation (6.7) could be adopted, keeping in mind that exactly on the diagonal m 1 = m 2 cross sections do vanish. 7 Constraints expected from future e + e − colliders Standard Model Higgs boson production in the Higgsstrahlung process is considered as a "golden channel" for a model independent determination of the Higgs boson properties at future e + e − colliders. By reconstructing the produced Z bozon, Higgsstrahlung events can be selected with high efficiency independently on the Higgs boson decay.
Largest sample of events can be selected when both leptonic and hadronic decay channels of the Z boson are considered. Reconstructing just the Z boson is of particular interest when we look for rare processes involving the Higgs boson, for instance possible decays into DM states. Events with mono-Z production, and no other activity in the detector, can be considered as candidate events for the invisible Higgs boson decays, if the recoil mass, Q 2 , reconstructed from energy-momentum conservation, is consistent with the Higgs boson mass. Highest sensitivity to invisible decays of the 125 GeV boson is expected at √ s 250 GeV, corresponding to the maximum of the Higgsstrahlung cross section. The main background processes that limit the sensitivity at this energy range are the production of ZZ and W + W − pairs, as well as single Z production via the W W fusion, e + e − → ν eνe Z. For the Z-boson pair production with one boson decaying into neutrinos, the final state reconstructed in the detector is identical to the one expected for the invisible Higgs boson decays and the recoil mass can be significantly overestimated due to beams spectra, 8 or large initial state radiation. For hadronic Z-boson decays, also detector resolution effects, dominated by the jet energy resolution, are very important. Nevertheless, due to branching fraction much larger than in the leptonic case, the expected limits on the invisible decays of the 125 GeV Higgs boson are dominated by hadronic Z decay measurements. For 2000 fb −1 of data collected at 250 GeV ILC, the expected limit on the invisible branching fraction is 0.23%, when combining hadronic and leptonic channels [47]. Similar sensitivity is expected also for other future e + e − collider projects [48].
The Higgsstrahlung analysis can be extended to the search for production of a generic scalar of arbitrary mass, assuming it is produced in association with the Z boson, as described in the previous section. The analysis procedure is the same as for the 125 GeV Higgs boson, only the event selection criteria have to be tuned to the considered scalar mass. The cleanest sample of Higgsstrahlung events is obtained when selecting Z boson decaying into muons, as the invariant mass of the µ + µ − pair can be reconstructed with sub-GeV precision and the background levels are significantly smaller than for the hadronic channel. This channel gives the best sensitivity to the production of light scalars, below 125 GeV, as the hadronic background levels increase rapidly towards low recoil masses and superior recoil mass reconstruction in muon channel allows for much better suppression of nonresonant background. No assumptions are made on the scalar decay modes or branching ratios. The expected number of events due to SM background processes remaining after the optimized selection cuts and the corresponding signal selection efficiency can be used to extract the expected cross section limit on the new scalar production as a function of its mass. Shown in Fig. 8 are the 95% CL exclusion limits expected for the ILC running at √ s = 250 GeV, normalized to the cross section for the SM Higgs boson production of a given mass [18,49,50]. Presented results assume Z-boson identification by its µ + µ − decays only. In the frequentist approach, the limit value is defined as the signal production cross section which, with probability of 95%, would result in the observed number of events higher than the SM expectation. The sensitivity is weakest for the scalar mass close to the mass of the Z boson, due to the background from Z-boson pair production (with one Z decaying into muons). Scalar masses up to the kinematic limit of √ s − m Z ∼ 155 GeV can be probed at 250 GeV.
To extend the limits towards higher scalar masses, e + e − collider running at higher energies is needed. If the new scalar is heavier than 125 GeV and it is expected to decay predominantly in invisible channels, the cross-section limits can be improved by considering hadronic Z boson decays. This gives factor of 20 increase in the expected signal event statistics (compared to the Z → µ + µ − decay channel) with only moderate increase in background levels, as the mono-Z signature allows for efficient suppression of SM background processes [51].

Numerical results
Here we will apply the strategy described in earlier sections to investigate how large the total cross section for Z and DM production at an e + e − collider could be. In order to maximize the cross section we will focus on colliders running at the CoM energy close to √ s = 250 GeV, while drawing figures we specialize to the case of the ILC at exactly √ s = 250 GeV [18,49]. The cross section depends on four independent variables: m 2 , m DM , sin α and v S . Instead of v S one can use X defined by (5.2), which is fixed by the relic abundance. Then, for each point (m 2 , m DM ) in our plots, Figs. 9-13, we choose such a value of sin α ≤ 0.3 that maximizes the cross section. Due to the resonant enhancement the cross section is much greater in area where at least one of on-shell mediators (h 1,2 ) can decay into pair of DM particles. As seen from (6.1), the differential cross section is maximized when the two Higgs bosons are on-shell at the same missing invariant mass Q 2 m 1 m 2 . Therefore, the total cross section is largest when m 1 m 2 . Hence, the maximum appears in the lower-left quarter as close to the diagonal m 1 = m 2 as allowed by the invisiblebranching-ratio condition. 9 In the case of vector and fermion DM models, the direct detection limits on the DM-nucleon cross section are very strong, a consequence of that is that couplings between DM and the SM (parametrized by X ) must be severely suppressed. Therefore, in general, in order to satisfy the DD constraint and at the same time provide appropriate DM abundance, the early-Universe DM annihilation must occur in a vicinity of a resonance, i.e. either 2m DM − m 1 0 or 2m DM − m 2 0. For the pGDM, because of the natural suppression of the DD cross section (that is vanishing at the tree level in the limit of zero momentum transfer, see Sec. 5.3), the resonant annihilation is not necessary to reproduce the correct DM abundance. Nevertheless, to compare the models, we have found it convenient to plot the cross section in the space spanned by m DM − m 1 /2 and m DM − m 2 /2 in the vicinity of the resonance, i.e. m DM m 1 /2 and/or m DM m 2 /2. The DM and h 2 masses adopted hereafter satisfy the following constraints what implies that 57.5 GeV < m DM < 67.5 GeV and 105 GeV < m 2 < 145 GeV. 10 In Figs. 9-10 we plot maximized cross section for the Z and DM production (normalized to the SM prediction for the Zh SM production, i.e. σ SM = σ(e + e − → Zh SM )| m h SM =m 1 ) at the ILC for pGDM, VDM and FDM models, respectively. The greenish colors denote regions where the models satisfy adopted constraints showing (by color) the corresponding cross section. The cyan marks regions excluded by the SM invisible BR limit, BR(h 1 → DM) < 0.19, while the black corresponds to parameters excluded by the DD limit (see Sec. 5.3). As explained earlier, the allowed regions for the VDM and the FDM models appear in the vicinity of resonant DM annihilation. For the gray region, the expected 95% CL sensitivity limit for σ/σSM (shown in Fig. 8) is above the σ/σ SM prediction. Therefore, one can conclude that the greenish regions are those which are detectable at the ILC, assuming that the sin α is close to the value that maximizes the cross section. It turns out that usually the sin α that provides maximal cross section is just at the largest value allowed by the LHC Higgs signal measurement, sin α 0.3. The fair conclusion from inspecting Figs. 9-10 is that in the substantial part of the parameter range that was shown, the DM production can be detected at future e + e − colliders running around √ s = 250 GeV. 9 It should be remembered that in the closest vicinity of m1 = m2 one should adopt results discussed at the end of Sec. 6. However, with the resolution adopted to draw plots in this paper those effects are invisible. 10 Region in (mDM, m2) plane that corresponds to (8.1) is not a rectangle. The simplest and straightforward way to disentangle the models is to measure m 2 and m DM and then verify if the measured masses are consistent with any of the discussed models after imposing constraints. In other words, one would need to check if for the measured values of m 2 and m DM the corresponding point (m DM − m 1 /2, m DM − m 2 /2) is located in the greenish area in any of Figs. 9-10. In order to facilitate and illustrate the verification, we have made plots shown in Fig. 11. The upper-left panel shows the cyan region where the pGDM is excluded by the h 1 -invisible-BR condition, while the white region is allowed. In the upper-right panel the yellow regions denote region where the FDM model is disallowed, Figure 11. The parameter space allowed or forbidden for the discussed models. Top-left: pGDM model, top-right: FDM model, bottom-left: VDM model, bottom-right: the three models combined. Light-and dark-gray regions denote violation of perturbativity conditions (5.14). Note that the |κ| < 4π condition is not violated in any place of the considered range of parameters. while white color stands again for region that agrees with all the constraints. Similarly, the lower-left panel shows forbidden (magenta) and allowed (white) regions for the VDM model. The lower-right plot combines results for all the models; again, white denotes the region where all the models are allowed. As it is seen, there exist regions where two or even three models coexist. However, there is also, in the lower-right panel, the magenta region where only the pGDM may exist. Therefore, the very first step in an attempt to disentangle the models should be a measurement of m 2 and m DM and its verification against the results shown in Fig. 11. Even though there is a substantial region of full degeneracy (white), there exists also significant area where some valuable conclusions could be drawn. It is even conceivable that this measurement would be consistent with the spin 0 (pGDM) hypothesis only. Now, we would like to focus on estimating chances to disentangle pairs of the models by the measurement of the normalized cross section σ/σ SM . In order to verify this option, we plot in Figs. 12 and 13 differences between model predictions and compare them against the expected experimental precision given by the limit provided by Fig. 8. As previously, it turns out that the highest differences are obtained for sin α as large as allowed, i.e. sin α 0.3. More reddish color indicate parameter regions for which models that are being compared are easier to disentangle since there an absolute value of the corresponding difference of cross sections is larger. It is clear that the disentanglement is a very ambitious task. It seems that only the VDM and the pGDM could be relatively easily disentangled by the measurement of e + e − → Z + · · · cross section at future e + e − colliders operating near √ s = 250 GeV if the parameters are in the more reddish regions of Fig. 12. Figure 12. The difference between predictions for the pGDM and the VDM. The gray region denotes parameter space for which the difference is smaller than the limit of Fig. 8. The models are compared in the region where both of them are consistent with the data, see Fig. 11.

Summary
In this analysis, our goal was to investigate how could one disentangle models of dark matter of different spin at future e + e − colliders operating near √ s = 250 GeV. For that purpose, we adopted the ILC project with the CoM energy at √ s = 250 GeV. Our strategy was pragmatic and phenomenological. We considered three nearly simplest models of dark matter of spin 0, 1 and 1 /2. The models adopted here were not the "simplified" ones, discussed often in a phenomenological literature on dark matter; in contrast, they were simple but attractive, consistent and renormalizable quantum field theories. In spite of dark-matter-spin differences, the models considered here share exactly the same parameter space, so the comparison point by point was meaningful.
In order to verify if a model was testable, we had adopted expected 95% CL sensitivity for the measurement of the e + e − → Z + · · · cross section obtained for the ILC project. It has been assumed that only the Z boson is reconstructed without any other detector activity. Predictions of the models were calculated taking into account the dark matter abundance, indirect and direct detection experiments and the collider constraints on the Higgs-boson invisible branching ratio and limits on the mixing angle (present in all the models). That way, regions of the parameter space where the cross section would be measurable were obtained for each of the models. We have also discussed the possibility to disentangle the models by a measurement of the cross section. It turned out that the most optimistic case is the detection of differences between the pseudo-Goldstone dark matter (spin 0) and the vector dark matter (spin 1). Regions of the parameter space where no model is allowed or some models could coexist were also determined.

Acknowledgments
The work was partially supported by the National Science Centre (

A Higgs boson self-energies and decay widths
Here we collect results for imaginary parts of two-point functions Π ij (Q 2 ) for Higgs bosons h i,j . Each self-energy is a sum of contributions of loops with various intermediate states 11 being on-shell: where DM stands for dark matter particle A, X or ψ while q denotes SM quarks and l denotes SM leptons. These contributions are given by (see [45] for the VDM case 12 ): where V 111 ≡ 3m 2 1 sin 3 α v S + cos 3 α v , (A.8) 11 We omit tadpole and seagull diagrams. The h 1 's and h 2 's partial widths can be calculated as The widths relevant for this project are therefore given by