Multi-component dark matter: the vector and fermion case

Multi-component dark matter scenarios constitute natural extensions of standard single-component setups and offer attractive new dynamics that could be adopted to solve various puzzles of dark matter. In this work we present and illustrate properties of a minimal UV-complete vector-fermion dark matter model where two or three dark sector particles are stable. The model we consider is an extension of the Standard Model (SM) by spontaneously broken extra U(1)X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U(1)_X$$\end{document} gauge symmetry and a Dirac fermion. All terms in the Lagrangian which are consistent with the assumed symmetry are present, so the model is renormalizable and consistent. To generate mass for the dark-vector Xμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_\mu $$\end{document} the Higgs mechanism with a complex singlet S is employed in the dark sector. Dark matter candidates are the massive vector boson Xμ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_\mu $$\end{document} and two Majorana fermions ψ±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\psi _\pm $$\end{document}. All the dark sector fields are singlets under the SM gauge group. The set of three coupled Boltzmann equations has been solved numerically and discussed. We have performed scans over the parameter space of the model implementing the total relic abundance and direct detection constraints. The dynamics of the vector-fermion dark matter model is very rich and various interesting phenomena appear, in particular, when the standard annihilations of a given dark matter are suppressed then the semi-annihilations, conversions and decays within the dark sector are crucial for the evolution of relic abundance and its present value. Possibility of enhanced self-interaction has been also discussed.


Introduction
The experimental data from the WMAP [1] and more recently the Planck [2] collaborations provided an independent and a e-mail: aqeel.ahmed@fuw.edu.pl b e-mail: mateusz.duch@fuw.edu.pl c e-mail: bohdan.grzadkowski@fuw.edu.pl d e-mail: michal.iglicki@fuw.edu.pl indisputable confirmation for the presence of dark matter (DM) in the Universe. Nevertheless, in spite of a huge theoretical and experimental effort, its nature is still unknown. Till now only gravitational interactions of DM have been detected in a series of independent observations like the flatness of rotation curves of spiral galaxies [3], gravitational lensing [4], and collision of galaxy clusters with its pronounced illustration known as the Bullet cluster [5]. All attempts to detect DM non-gravitational interactions with ordinary matter have failed so far implying more and more stringent limits on DMnucleon cross-section, see e.g. LUX [6] and XENON-1T [7] results. The most popular models of DM are based on the assumption that it is composed of weakly interacting massive particles (WIMPs). Unfortunately, it turns out that the WIMP scenarios suffer from various difficulties when confronted with observations on small cosmological scales. For instance, the "too-big-to-fail" [8,9] and the "cusp-core" [10][11][12][13] problems are widely discussed in the literature. In particular, the DM densities inferred in the central regions of DM dominated galaxies are usually smaller than expected from WIMP simulations [14,15]. It turns out that an appealing solution to those problems is to assume that dark matter selfinteracts strongly [16]. The assumption of self-interacting dark matter (SIDM) implies that central (largest) DM density could be reduced. Usually, self-interacting DM scenarios require the presence of light DM and also light DM mediators.
Dark matter could also be searched for through indirect detection experiments, which assume that in regions of large DM density, its pairs are likely to annihilate. Then secondary particles released in this process, e.g. gamma rays, neutrinos, electrons, positrons, protons, and anti-protons, could be observed on Earth, which could reveal some properties of DM. Independent analysis of the Fermi-LAT data [17] by various groups have shown an excess of gamma ray in the energy range 1−3 GeV that can be interpreted as a result of DM annihilation in the Galactic Center. Besides the 1-3 GeV excess gamma rays, there exists an observation of unidentified 3.55 keV X-ray line found by [18] and [19]. As shown by several groups, this unknown X-ray line can also be explained by DM annihilation. To explain the indirect signals relatively large DM mass is needed, e.g. ∼ 50 GeV in [20].
In this work, we propose a minimal UV complete vectorfermion DM model that predicts two or three stable dark states. Our model involves an extension of the SM by an Abelian dark gauge symmetry U (1) X . The model is minimal in a sense that it contains only three new fields in the dark sector; a dark gauge boson X μ , a Dirac fermion χ , and a complex scalar S, which serves as a Higgs field in the hidden sector. They are singlets under the SM gauge group but they are all charged under the dark U (1) X gauge symmetry and therefore they interact with each other. The complex scalar S acquires a vacuum expectation value (vev) and gives mass to the dark gauge boson X μ by the dark sector Higgs mechanism. It also contributes to the mass of the Dirac fermion χ through the Yukawa coupling. Moreover, the presence of the Dirac mass for the fermion introduces a mixing of its chiral components. After diagonalization of the mass matrix, the Dirac fermion splits into two Majorana fermions ψ ± with mass eigenvalues m ± . As a result, after the dark sector symmetry breaking we have three potentially stable interacting particles: a dark vector and two Majorana fermions. Their stability is ensured by a residual Z 2 ×Z 2 discrete symmetry, which also dictates the possible form of dark sector interactions. The communication with the visible (SM) sector proceeds only via the Higgs portal κ|H 2 ||S 2 |.
Our minimal vector-fermion dark matter model has many attractive features. First of all, the very fact that in the multi-component DM literature, the vector-fermion dark matter possibility has not been studied in detail 1 speaks for itself and therefore the goal of this work is to provide an extensive analysis of the minimal vector-fermion scenario. Some of the interesting features of the model are: (1) the presence of a second scalar helps to achieve the stability of the SM Higgs potential even at the tree level, see e.g. [82][83][84], (2), a possibility of enhancing vector component self interactions, see e.g. [85], (3) a very small/large mass splitting among the dark sector states (vector and Majorana fermions) are possible without large tuning of the parameters, (4) our model is a gauged version of the model considered by Weinberg [86] for different purposes, and (5) more importantly the minimality of the model; since there is only one parameter, the dark gauge coupling coupling g x , which controls the dynamics in the dark sector, including the conversion, semi-annihilation and decay processes. In this work, we are going to illustrate the relevance of conversions, semi-annihilations, and decays in the vector-fermion DM model.
The paper is organized as follows. In Sect. 2 the vectorfermion (2-3 component) model of DM is presented. Solutions of three coupled Boltzmann equations are discussed in Sect. 3 focusing on conversion, semi-annihilation, and decay processes. There we show results of a detailed scan over the parameter space of the model satisfying the observed total relic density and direct detection constraints. In Sect. 4 we focus on the region with large self-interaction crosssection. Section 5 contains summary and conclusions. Moreover, we supplement our work with an "Appendix A", collecting details of the derivation of Boltzmann equations, and an "Appendix B", describing the method adopted to obtain constraints for a multi-component DM model by direct detection experiments.

Vector-fermion multi-component dark matter model
In this section, we explore the possibility of having a multicomponent dark matter model with a vector boson and a Dirac fermion (charged under a dark sector gauge symmetry) which may serve as dark matter candidates. As mentioned in the Introduction, we consider a minimal extension of the SM by an additional U (1) X gauge symmetry of the dark sector with a complex scalar field S and a Dirac fermion χ , both charged under the dark-sector gauge group. We employ the Higgs mechanism in the dark sector such that the vev of the complex scalar S generates a mass for the U (1) X gauge field X μ . The dark-segment fields have the following quantum numbers We assume that none of the SM fields is charged under the dark gauge symmetry U (1) X . We can write down the Lagrangian for our simplest vectorfermion MCDM model as where L SM is the SM Lagrangian, L DS is the dark-sector Lagrangian, and L portal is the Lagrangian for the Higgs portal interactions between the SM and the dark sector, Note that the portal coupling is the only communication between the SM and the dark sector. Above, in Eq. (3) X μν ≡ ∂ μ X ν − ∂ ν X μ is the field strength tensor for the X μ field and D μ is the covariant derivative defined as where g x is the coupling constant corresponding to the gauge group U (1) X , whereas q x are the U (1) X charges 1 and 1 2 for S and χ (as defined in 1), respectively. Moreover, in Eq. (3) m D is the Dirac mass, y x is the dark Yukawa coupling and C ≡ −iγ 2 γ 0 is the charge-conjugation operator, where γ 0 , γ 2 are the gamma matrices. It is important to note that the dark sector is invariant with respect to the charge conjugation symmetry C under which the dark fields transform as follows: The symmetry forbids a kinetic mixing between the weak hypercharge U (1) Y and the dark U (1) X , so that X μ can not decay into SM particles. It is useful to collect the full scalar potential of our model, Tree-level positivity/stability of the above scalar potential requires the following conditions to be satisfied: It is straightforward to find the minimization conditions for the scalar potential (7) as where H ≡ (0, v/ √ 2) and S ≡ v x / √ 2 are the vevs of the respective fields. 2 In order to have both vevs non-zero (v provides masses to the SM gauge bosons and the dark-sector scalar vev v x generates mass for the dark vector) we require κ 2 > 4λ H λ S and the values of vevs are: We expand the Higgs doublet and the singlet around their vevs as: where π 0,± and σ are the would-be Goldstone modes that disappear in the unitary gauge. Hence, only the scalars h and φ correspond to the physical states. The mass squared matrix for the scalar fluctuations (h, φ) is given by, The above mass squared matrix M 2 can be diagonalized by the orthogonal rotational matrix R, such that, and (h 1 , h 2 ) are the two physical mass-eigenstate Higgs bosons with masses (m 2 h 1 , m 2 h 2 ), defined in terms of (h, φ) as The mixing angle α could be expressed in terms of mass matrix elements as follows: We will later adopt sinα as an independent input parameter while scanning over model parameters. The masses for the two eigenstates h 1 and h 2 are In this analysis we will treat the two Higgs masses as independent parameters. It is always assumed that h 1 is the observed state with m h 1 = 125 GeV. The other mass eigenstate could be either lighter or heavier than h 1 . After the spontaneous symmetry breaking (SSB) the SM fermions acquire mass from the SM Higgs doublet, whereas the dark sector fermion χ receives mass from the dark-scalar S through the Yukawa interaction term and from the Dirac mass term. After the SSB the dark fermionic sector Lagrangian can be rewritten as, In the above Lagrangian the field χ and its charge-conjugate χ c mix through the Yukawa interactions. One can easily diagonalize the above Lagrangian in terms of the following Majorana mass-eigenstates ψ ± (= ψ c ± ) with mass eigenvalues, In the new basis we can rewrite the above dark fermionic Lagrangian as, The mass splitting between ψ + and ψ − is controlled by the Yukawa coupling and v x : Hereafter, we will assume y x > 0 and since v x could be chosen positive therefore we have m + > m − . Note that the above Lagrangian is invariant with respect to a discrete symmetry Z 2 × Z 2 , under which the SM fields are even whereas the dark sector fields transform non-trivially, as given in Table 1. It is easy to see that the above Z 2 is a direct consequence of the charge conjugation symmetry (6). It is worth to notice that since v x is real (without compromising any generality) therefore the charge conjugation remains unbroken so that X μ cannot decay into SM particles, see also [59]. On the other hand Z 2 is implied by the U (1) X gauge symmetry. Note that the Z 2 is responsible for the stability of ψ − since it is lighter than ψ + . It is interesting to notice that our model has also a discrete Z 4 symmetry 3 under which the three dark matter components are charged. The Z 4 charges are: with n = (1, ± 1 2 , 0), whereas all the SM particles are neutral under this symmetry. However, for our analysis the Z 2 × Z 2 (along with Z 2 ) completely specifies all the relevant properties of Z 4 , therefore hereafter we only consider As it can be seen from the above fermionic Lagrangian there is only one interaction term between all three components of the dark segment (the vector boson X μ and the two Majorana fermions ψ ± ): In Fig. 1, we collect all the relevant vertices that desckineticribe interactions of the dark segment of the theory. As depicted in Fig. 2, depending on the masses of dark particles, there are three cases 4 in which two or all three particles could be stable and serve as dark matter: -The first case is when m + > m − + m X , the Majorana fermion ψ + will decay into a stable vector X μ and a stable Majorana fermion ψ − . This is a 2CDM case, the white area (left) in Fig. 2. -The second case is when m X > m + + m − , the vector X μ will decay into two stable Majorana fermions ψ ± . This is a 2CDM case, the white area (right) in Fig. 2. -The third case is when m + +m − > m X > m + −m − , so that none of the three particle will decay and hence all are stable. This is a 3CDM case, shown as gray region in Fig. 2. Note that the boundaries (right/left) of the gray region correspond to the case when m X = m + ± m − .

Input parameters
Here we outline the strategy adopted to investigate the vectorfermion dark matter (VFDM) model. First of all, we would like to count the free parameters in the model and the number of constraints. There are five parameters in the scalar potential (7), namely, mass parameters μ H , μ S , quartic couplings λ H , λ S and the portal coupling κ. Additionally, there are three parameters in the dark gauge and fermionic part, i.e. the dark gauge coupling g x , the fermion Dirac mass m D , and the Yukawa coupling y x , see Eq. (3). In total there are eight parameters of the model, however there are two constraints from the SM Higgs vev v and the Higgs mass m h 1 , hence leaving six parameters free. We adopt the following set as an independent input parameters: Note that the mixing angle sinα is constrained by the SMlike Higgs couplings with the gauge bosons. We employed | sinα| ≤ 0.33, the current LHC 2σ bound [87].
Remaining parameters could be expressed in terms of the input set as follows (note the potential mass parameters μ H , μ S are already traded for v, v x , see Eq. 10) : Note that the Yukawa coupling y x = Δm/(2v x ) = Δm/(2m X ) ×g x , therefore for fixed m X the Yukawa coupling y x is proportional to g x . In other words, for fixed m X , the vev v x must vary if g x is being changed, implying a variation of y x . This remark is important hereafter, e.g. for fixed Δm and m X one has to adjust g x in order to satisfy constraints from direct detection experiments even in the case when DM is dominated by ψ ± .

Vector-fermion dark matter phenomenology
In this section we present coupled Boltzmann equations for the evolution of number density of dark matter particles (X μ , ψ + , ψ − ) in VFDM model. Figures 3, 4 and 5 contain Feynman diagrams relevant for collision terms for the annihilation, semi-annihilation and conversion processes, respectively. It is assumed that dark matter components scatter against SM particles frequently enough so that their temperatures are the same as that of the thermal bath. The Boltzmann equations for the DM components (X μ , ψ + , ψ − ), can written as: Semi-annihilation diagrams for the dark particles X, ψ ± . In these processes two of the dark sector particles annihilate to a dark sector particle and a SM particle in the thermal bath where σ i jkl v M∅l ≡ σ i jkl v M∅l is the thermally averaged cross-section for the process i j → kl as defined in Eq. (A.17). Above h i = h 1 , h 2 and φφ denote all the allowed SM particles including h 1 , h 2 . In the above Boltzmann equations (30)-(32), the first term 3Hn i is the usual term in an expanding universe with Hubble parameter H . The second term in Eqs. (30)-(32) is the standard annihilation term for each dark particle corresponding to Feynman diagrams Fig. 3, whereas, the third, fourth and fifth terms are capturing the effects of semi-annihilations shown in Feynman diagrams Fig. 4, and the sixth and seventh terms are conversion processes within the dark sector shown in Fig. 5. Note the last term in these Boltzmann equations is for the case considered in Sec. 3.1.1 when ψ + is unstable and it decays to stable particles X μ and ψ − . However, for the case discussed in Sec. 3.1.2 where X μ is unstable and it decays to stable particles ψ + and ψ − , we replace ψ + ↔ X and change signs of the last terms in Eqs. (30)- (32). Moreover, for the case studied in Sec. 3.1.3 when all three dark particles X μ , ψ ± are stable then the last terms in the above Boltzmann equations are zero. The details of the derivation of above collision terms are presented in "Appendix A".
After solving the Boltzmann equations we calculate the present relic density of the dark species as, where s 0 is the total entropy density today, ρ cr is the critical density, m i is the mass of the dark particle and Y i is the yield of the dark particle today. Total dark matter relic density is a sum of the individual relics, i.e.
The total relic density Ω tot h 2 is compared to the observed dark matter relic density Ω obs h 2 = 0.1197 ± 0.0022 by the PLANCK satellite [2]. The simple form of the derived Boltzmann equations relies on the assumption of kinetic equilibrium between visible and hidden sectors which is maintained by frequent scatterings of dark matter species on relativistic SM states. Effects of kinetic decoupling are negligible in the calculation of relic densities only if its temperature T kd is smaller or comparable to the chemical decoupling temperature, T cd , at which DM annihilation is no longer effective. T kd can be estimated by comparing the scattering rates Γ s (T ) for the processes DM SM → DM SM with the Hubble rate H (T ) [88].
where n r is the equilibrium density of relativistic states in the thermal bath, σ s v is the thermally averaged scattering cross-section and T /m DM describes momentum transfer at each collision. m DM is the mass of the dominant DM component, so either m X or m ± . Considering scatterings of ψ ± and X on SM quarks and leptons we find using the expansion of the scattering amplitudes [89] that where ξ = 1 if scattering of X dominates and ξ = 2 √ m X m ± /(m + −m − ) for the ψ ± domination. Choosing typical values that are used in further discussions g x ∼ sin α ∼ 0.1 and m DM ∼ O(100) GeV, we obtain T kd 0.5 GeV which is well below the temperature of the chemical decoupling in this case (T cd m DM /20). It becomes comparable to T cd only for degenerate Higgs masses, but we checked that for the parameters we adopt always T kd < T cd .
Even without solving the above coupled set of Boltzmann equations (30)-(32) some generic observations could be made: -Conversions are present even in the absence of the dark sector self-interactions, an existence of a mediator is the only requirement. On the other hand, the existence of semi-annihilations and decays of dark particle depend on the presence a vertex with three dark states, which have different transformation rules under the dark symmetry (such that a singlet could be formed), in our model such an interaction is in Eq. (24). -When, for a given dark matter species, a standard annihilation channel is suppressed then its abundance might be very sensitive to the presence of other ingredients of the dark segment. In this case semi-annihilation plays a major role, e.g. if ψ − ψ − → h i h i (or any SM states) is suppressed then ψ − can still disappear, for instance, through ψ − ψ + → X μ h i followed by unsuppressed annihilation of pairs of X μ , see also [45]. In other words, X μ can work as a catalyzer that enables disappearance of ψ − . In this case, it is possible that the presence of other (ψ + and/or X μ ) dark components might be crucial for the determination of the asymptotic abundance of the major DM element. Also, decays within the dark sector may play a relevant role in the determination of the final abundance. -Standard WIMPs decouple from thermal equilibrium at m/T ∼ 20−25, which implies that the heavy states decouple earlier (large T ). However, in the multi-component scenario, it might be possible that the decoupling of a heavier dark component is delayed so that it happens later than that of a lighter one. The effect is again a consequence of interactions with remaining dark matter states.
In particular, as a proof of principle, below we consider two interesting possibilities in our specific vector-fermion DM scenario: Small y x implies suppressed ψ ± ψ ± annihilation, so ψ ± dominates the dark matter abundance. Since the annihilation is slow therefore Y ψ ± is controlled by semi-annihilation which is sensitive to g x and to the presence of other dark components. In order to have semi-annihilation controlled exclusively by g x one should assume m In this case, one expects fast ψ ± ψ ± annihilation and so that X μ may dominate the dark matter abundance. If in addition sinα 1 then X X annihilation would be suppressed so that Y X shall be controlled by semi-annihilation and conversion processes which are sensitive to the gauge coupling g x and Yukawa coupling y x . In both cases, X μ would be effectively replaced by ψ ± , which then would disappear through enhanced standard annihilation.
It is worth to emphasize the importance of the gauge coupling g x between all the dark components X μ , ψ + and ψ − . This is the most relevant coupling which determines interesting aspects of dynamics of dark matter density evolution. If the gauge coupling g x was very small then the model would reduce to a simple sum of two non-interacting components (ψ + , ψ − ), i.e. a rather uninteresting scenario. Note that g x is a consequence of the presence of the (−, −) state and, as illustrated by our vector-fermion model, existence of this state is quite natural. Note that if the (−, −) state would have been absent then only two fermion dark components (ψ + , ψ − ) would have been allowed by the stabilizing symmetries. However, in this case, only annihilation and conversion processes -without semi-annihilations and decayswould have been allowed. Again not a very appealing scenario.

Multi-component cases
In the following sections we consider various interesting setups with two or three dark particles. Matrix elements squared needed for collision terms in the Boltzmann equations (30)- (32) are computed by employing the CalcHEP [90], whereas for thermal averaging and solutions of the Boltzmann equations, we adopt our dedicated C++ code. 5 5 Note that the presence of three stable dark matter components is quite generic in models with two interacting stable states. In this case, even the 2-component version of micrOMEGAs [91] is not applicable as the code assumes there are at most two DM sectors within which particles remain in equilibrium. Therefore for the case of 3-component DM, we adopt our dedicated code which employs the full set of three Boltzmann equations.

2CDM: a vector and a Majorana fermion as dark matter
In this section we show results for the scenario with m + > m − + m X , so that ψ + is unstable and can decay into lighter Majorana fermion ψ − and the vector boson X μ . Fig. 6 shows results of a scan over sinα, g x , m X for fixed m h 2 , m − and Δm = (m X +10) GeV. All the points satisfy the correct relic density (for the total abundance) observed by PLANCK at 5σ and the recent direct detection experimental bound from LUX2016 at 2σ . Note that, since the Yukawa coupling to dark fermions is proportional to Δm/v X = (Δm/m X )g X therefore for fixed Δm and m X the annihilation cross-sections for the both ψ − and X are proportional to (g X sin 2α) 2 therefore, if the remaining parameters are fixed, the requirement of correct DM abundance determines g X sin 2α. Then the mass of the DM component decides which component contributes more to the Ω DM .The left panels are for close h 1 and h 2 masses, m h 2 = 120 and m h 2 = 130 GeV which allows for cancellation between h 1 and h 2 contributions to the direct detection cross-section on nuclei and therefore for consistency with LUX2016 data even for coupling constants that are not so small (g x = 0.3 − 1.0 with majority of points located above g x ∼ 0.3). Note that also sin α does not need to be particularly small due to efficient cancellation between h 1 and h 2 contributions to the direct detection. The right panels correspond to parameter points such that around m X = 200 GeV there is a resonant enhancement of X X annihilation through the h 2 s-channel exchange. Of course, there is also a non-resonant contribution from ψ − ψ − annihilation. In the case of h 2 resonance, in order to satisfy the relic abundance constraint, the relevant coupling constants must be small, i.e. g X sin α 1, it turns out that for the region of sin α considered here the gauge coupling constants must be in the range g x ∼ 0.03 − 0.13.
Lower panels in Fig. 6 show the abundances of both components separately, Ω 2cdm X and Ω 2cdm − , for our model and corresponding Ω 1cdm X calculated in a 1-component VDM assuming the Majorana fermions are decoupled (which is achieved by m D → ∞). The abundance and direct detection cross-section in the VDM model were calculated for the same parameters as those adopted in the 2-component model, just truncated to m h 2 , m X , sinα, g x . It is worth to focus at X masses between 110 and 200 GeV in the lower-left panel, where for each given m X there is a clear shift upwards of the X abundance in the 2CDM scenario (filled triangles are above empty boxes). This is a nice illustration how the presence of the other components of the dark sector may influence the abundance of X . Similarly, in the lower-right plot, an analogous shift is observed for m X < ∼ 200 GeV, this shift is enhanced by small sinα, e.g. for sinα ∼ 0.1 the X abundance receives an extra factor ∼ 10 2 . Usually, the presence of Fig. 6 Upper panels: The plots show results of scans over the parameter space of a 2CDM case, where X μ and ψ − are stable. All the points satisfy the correct relic density observed by PLANCK at 5σ and the recent direct detection experimental bound from LUX2016 at 2σ . Lower panels: The plots show results of scans over the parameter space of a 1CDM (assuming the fermions ψ ± are decoupled from the dark sector and hence X μ is the only dark matter particle) and 2CDM (where X μ and ψ − are stable) cases, respectively. In the 2CDM case, all the points are taken from the corresponding upper panel plots. Note that for the 1CDM case, in most of the parameter space, the single dark matter X μ is under-abundant, whereas in the 2CDM the presence of second component ψ − provides the remaining relic density. The point denoted by the black star located in the upper-left panel at m X = 245 GeV corresponds to the same parameters as those adopted in the middle panel of Fig. 7 other constituents of the dark sector implies both an upward shift of Ω X , when compared to the 1CDM, and provides a necessary extra contribution by Ω 2cdm − in order to satisfy the abundance constraint.
In the plots of Fig. 7, and similar figures in the following sections, we show dark matter yields Y i (x) ≡ n i /s (s is the total entropy density and x is defined as x ≡ m X /T ) for different species i = ψ + , ψ − and X μ . Note that we plot bare values of yields, with no extra normalization adopted. Moreover, the tables in Fig. 7 and similar figures in the following section contain first two non-vanishing coefficients of thermallyaveraged cross-sections (pb) expanded in powers of x −1 , given by σ i jkl v M∅l = a N x −N + a N +1 x −(N +1) + · · · , and the decay width Γ ψ + →X ψ − (GeV). The plots in Fig. 7 illus-trate solutions of the Boltzmann equations for three selected sample points. The middle panel shows solutions for parameters that reproduce correct total DM abundance and also satisfy the direct detection LUX limits, so the corresponding point is also present in the scan results shown in the upperleft panel of Fig. 6 as a black star . In order to illustrate the relevance of the g x coupling, the left, middle and right plots are obtained for g x = 0.02, 0.2 and 1, respectively, while other parameters remain unchanged. It is clear that the dependence of the abundance for the major DM component (ψ − ) on g x is very strong. The gray dots and boxes show results obtained using the micrOMEGAs code for 2CDM [91]. As it is seen in the plots they agree very well with the solid lines which correspond to solutions obtained from our The tables contain first two non-vanishing coefficients of thermallyaveraged cross-sections (pb) expanded in powers of x −1 , given by σ i jkl v M∅l = a N x −N + a N +1 x −(N +1) + · · · , and the decay width Γ ψ+→X ψ− (GeV) dedicated C++ code that solves the set of three Boltzmann equations (30)- (32). In order to identify the most important processes for a given parameter set, in the tables below the panels in Fig. 7 we collect the first two non-vanishing coefficients in the expansion of thermally averaged cross-sections in powers of x −1 .
Let's look closer at the middle table of Fig. 7. The crosssections shown there correspond to the point in the parameter space marked by which is located in the upper left panel of Fig. 6 at m X = 245 GeV. As seen from the middle upper panel of Fig. 7 the abundance is dominated by ψ − , and X is a sub-leading component abundance of which is by nearly two orders of magnitude smaller than for ψ − , while the abundance of ψ + is absolutely negligible. Note that both ψ − and X decouple from equilibrium roughly at the same temperature, this is the first signal that there must be some correlation between annihilation mechanisms responsible for their disappearance. Since the abundance of ψ + could be neglected the only relevant processes are X X → SM, ψ − ψ − → SM and X X → ψ − ψ − . Note that the ratio of cross-sections for the first and the second process is ∼ 1.8 while their abundances differ by almost two orders of magnitude, therefore the process for additional depletion of X abundance must be X X → ψ − ψ − . This illustrates how sub-leading components may influence the abundance of a dominant component.

2CDM: two Majorana fermions as dark matter
In this section we show results for the scenario with m X > m − +m + , so X μ is unstable and can decay into the Majorana fermions ψ − and ψ + . Figure 8 shows results of a scan over sinα, g x , m − for fixed m h 2 , m X and Δm = 100 GeV (left panel) or 50 GeV (right panel). All the points satisfy the correct relic density (for the total abundance) observed by PLANCK at 5σ and the recent direct detection experimental bound from LUX2016 at 2σ . In this case the ψ − turns out to be the dominant DM component in most of the parameter space. The second Higgs boson mass was chosen to be m h 2 = 120, 125, 130 GeV and 390, 400, 410 GeV in the left and right panels, respectively. Therefore the left panel allows for partial cancellation between an exchange of h 1 and h 2 both for ψ ± annihilation diagrams and also for ψ ± -nuclei scattering process to avoid the direct detection limits even if couplings are not small. Since 50 GeV ≤ m − ≤ 200 GeV therefore one can observe both a resonance behavior at m − ∼ m h 1 /2 ∼ m h 2 /2 in ψ ± ψ ± annihilation trough s-channel h 1,2 exchange and Fig. 8 The plots show results of scans over the parameter space of a 2CDM case, for which ψ + and ψ − are stable. All the points satisfy the correct relic density observed by PLANCK at 5σ and the recent direct detection experimental bound from LUX2016 at 2σ . The point denoted by the star located in the right panel at m − = 199.25 GeV corresponds to the same parameters as those adopted in the middle panel of Fig. 9 also a threshold effect at m − ∼ m h 1 ∼ m h 2 for annihilation into h i h j final state. The vertical structure observed in Fig. 8 around m − ∼ 60 GeV corresponds to a domination of ψ ± ψ ± → h * i → V V,f f . As it is seen from the plot large values of g x are needed, this is a consequence of partial cancellation between h 1 and h 2 exchange. On the other hand, the independence on g x could be understood as a result of resonance enhancement around m − ∼ m h 1 /2 ∼ m h 2 /2: even a tiny change of m − can compensate large variation of g x . The other triple-branch structure that starts around m − ∼ 120 GeV corresponds to a threshold for the process ψ ± ψ ± → h i h j . Its initial steepness represents the opening of the h i h j final state that must be compensated by suppression of g x in order to generate correct dark matter abundance.
The right panel of Fig. 8 with its three distinct branches corresponds to a vicinity of the resonance at m − ∼ m h 2 /2 = (390, 400, 410)/2 GeV. In this case g x must be small to compensate the resonance enhancement, therefore direct detection limits are easily satisfied.
In Fig. 9, we illustrate solutions of the Boltzmann equations for three sample points in the parameter space. The middle panel corresponds to the correct abundance and is in agreement with the LUX upper limits on the DM-nucleon cross-section, which is also present in the scan results shown in the right panel of Fig. 8 as a black star . As in the previous case of stable X μ and ψ − , here we also observe relevance of semi-annihilation and conversion processes and strong g x dependence. One can see a discrepancy with micrOMEGAs results, which are substantial especially in the left panel. The reason for that is the influence of unstable component X μ , which can be properly described only with a set of 3 coupled Boltzmann equations, whereas in micrOMEGAs one has to assume it is in chemical equilibrium with one of the stable components.

3CDM: a vector and two Majorana fermions as dark matter
In this section we show results for the scenario with m + + m − > m X > m + − m − , so all the three dark components are stable. Figure 10 contains results of a scan performed adopting our dedicated code 6 that solves the set of the three Boltzmann equations (30)- (32). All the points presented in Fig. 10 satisfy the correct relic density (for the total abundance) observed by PLANCK at 5σ and the direct detection experimental bound from LUX2016 at 2σ . The scan is performed over m − , g x with fixed values of sinα = (0.05, 0.1, 0.2), m X = (200, 150, 100) GeV and m h 2 = (200, 120, 50) GeV. Note that the left and right panels of Fig. 10 are for the same data-set but for different filling style, in the left and right panel the filling corresponds to m X and m h 2 , respectively. Here we have tested sensitivity to Δm ≡ m + − m − focusing on small Δm = (0.1, 1, 10) GeV. As it is seen from the Fig. 10, for m − ∼ m + > ∼ 200 GeV relatively small U (1) X coupling is required, g x = 0.5 − 1, in order to suppress too fast ψ ± ψ ± s-channel annihilation. Note that y x = Δm g x /(2m X ) therefore this annihilation is already quite strongly suppressed by the small Yukawa coupling y x . An important final state is tt, so if m − ∼ m + < ∼ m t this annihilation channel closes so that even large g x is allowed/necessary, as observed in the figure. Figure 11 shows solutions of the Boltzmann equations for three sample points in the parameter space of 3CDM. The The tables contain first two non-vanishing coefficients of thermallyaveraged cross-sections (pb) expanded in powers of x −1 , given by σ i jkl v M∅l = a N x −N + a N +1 x −(N +1) + · · · , and the decay width Γ X →ψ+ψ− (GeV) Fig. 10 These plots show results of scan over the parameter space of a 3CDM case, where X μ and ψ ± are stable. In this scan we fix three different values of m X , m h2 and Δm and vary m − and g x as shown in the plots. The left and right plots represent the same data set however points markers are changed from m X (left) to m h2 (right), whereas the color represent the values of dark Yukawa coupling y x . All the points shown satisfy the correct relic density observed by PLANCK at 5σ and the recent direct detection experimental bound from LUX2016 at 2σ . The point denoted by the black star in the these plots at m − = 246.4 GeV corresponds to the same parameters as those adopted in the middle panel of Note that for these points the dark fermion Yukawa coupling is very small, i.e. y 1, as a result the direct annihilation processes for ψ ± are inefficient, therefore the main annihilation processes are through the semi-annihilations and conversions to X which further annihilate to SM middle panel represents the point marked as a black star in the scan results of Fig. 10 that satisfies relic abundance and LUX2016 constraints. As in the previous cases for 2CDM, here we also observe relevance of semi-annihilation and conversion processes and strong g x dependence on the yields of dark matter components.

Limiting cases
In this section we are going to discuss special regions in the parameter space of the model that result in simpler, models of DM.

The vector dark matter (VDM) model
If mass splitting between ψ + and ψ − is small comparing to m X , then, for fixed g X , the Yukawa coupling y X is suppressed as y X = Δm/(2v x ) = (Δm/m X )(g X /2). Therefore, in this limit, since Yukawa couplings become irrelevant, the model might be reduced to the 1-component VDM model (see e.g. [83,85] where the same notation as here has been adopted). Note however that even though fermionic DM decouple from the SM, nevertheless it is still present and may influence cosmological dynamics and contribute to the observed amount of DM. In order to enable efficient ψ ± annihilation it is sufficient to assume that 2m D > m X + m i and/or m X > m i (m D ≡ (m + + m − )/2) so that at present only X μ contributes to the observed DM abundance and could be successfully fitted by tuning g X . Samples of parameter sets that imply proper Ω DM and fit in the VDM limit are shown in Table. 2. In Ref. [85] the VDM has been analyzed focusing on the possibility of enhancing self-interaction, and some regions of the parameters space where elastic X X scattering is amplified and all other constraints are satisfied have been found.

The fermion dark matter (FDM) model
Another interesting limit of our model is a renormalizable model of fermionic DM, see e.g. [92][93][94][95]. Those models usually employ an extra singlet real scalar field that couples to the SM Higgs doublet via the Higgs portal and to a singlet dark Dirac fermion as well. In the fermionic DM limit of our model, we recover a model of a Majorana singlet DM that couples to a complex scalar S. Since our model is invariant under local U (1) X therefore in the limit of small gauge coupling, g X 1, and substantial mass-splitting between ψ + and ψ − (so for enhanced Yukawa coupling, y X ), effectively we obtain a renormalizable model symmetric under a global U (1) X of a single Majorana dark fermion ψ − interacting with S. The scalar S controls communication between dark sector and the SM. Examples of parameter sets that imply proper Ω DM and fit in the FDM limit are shown in Table 2. The model is slightly more restrictive than those considered earlier in [92][93][94][95] since our DM is a Majorana fermion and the scalar potential is more restricted, as being invariant under the global U (1) X , however predictions of those models are similar, see [94]. In order to obtain fermionic (Majorana) DM model one should reduce g X in order to decouple X μ , keeping in mind that certain minimal interaction strength is necessary for X μ to maintain kinetic equilibrium. Note that in order to reduce the model to a single fermionic DM model, we have to remove somehow ψ + and X μ . The easiest way to get rid of ψ + is to assume that it is the heaviest dark state, so that it will have a chance to decay quickly. If the following mass ordering, m + > m X > m − is fulfilled, then indeed the dominant DM component is the Majorana fermion ψ − , while other dark components disappear. The model contains, of course, two scalar Higgs bosons that mix in the standard manner and play a role of mediators between the dark sector and the SM. Self-interaction in the FDM model has been discussed in [96]. It turns out that for the self-interactions to be sufficiently strong, the scalar mediator, h 2 , has to be very light what implies well-known problems [97] in the early Universe if h 2 is present during the era of BBN.

The fermion dark matter (FDM) model with a stable vector mediator
Another interesting limit of our model has been considered very recently in [98]. If, in our model, Δm → 0 then the Yukawa coupling y X vanishes and masses of Majorana fermions ψ + and ψ − become degenerate. Then our model reduces to the model considered in [98], which is just a model of a Dirac fermion as a DM and a stable vector mediator. Their [98] mediator corresponds to our vector component of DM, X μ , while the dark Dirac fermion is an analog of our degenerate Majorana dark fermions ψ + and ψ − . The authors of [98] show that the model can indeed predict enhanced DM self-interaction while satisfying all existing experimental constraints if mass of the stable vector mediator is of the order of 1 MeV. This has been also confirmed in the appropriate (Δm → 0, m X ∼ O(1) MeV) region of the parameter space of our model. In this case the DM abundance is dominated by mass-degenerate ψ ± , even though formally it is a 3-component case (3CDM) if Δm < m X .

Distinguishing limiting cases
As it has been discussed in the previous section the model discussed here simplifies in various regions of the parameter space where it reduces to a single-component dark matter mode i.e. FDM or VDM. In remaining parts of the parameter space it describes a genuine 2 or 3 component dark matter. In this context it is natural to rise the question how could one disentangle those three possibilities. Some attempts to address this sort of question have already been made in the literature, see e.g. [99,100], where the authors considered a possibility to disentangle spin, 0, 1/2 or 1 dark matter at e + e − future colliders. The VDM model they considered was the same as the limiting version of our model discussed in Sect. 3.2.1, however their FDM was slightly different than ours from Sect. 3.2.2. Of course, they did not discuss multi-component scenario. An exhaustive discussion, that takes into account all existing experimental constraints within a single model that allows for 2 or 3 DM components is still missing. Such an analyzes lies beyond the scope of this paper however it shall be investigated in the near future [101]. Nevertheless few comments are here in order.
-Direct detection Contributions to DM-nucleon scattering consists of the sum of standard σ X −N and σ ψ ± −N cross-sections that are not sensitive to the presence of all the 2-3 DM components, rather this is a sum of contributions that exist in 1-component models. There exists however a more interesting inelastic scattering process which is sensitive to the multi-component nature of the model considered here, i.e. ψ + N → ψ − X N, note that all the dark particle are involved, so that this process might provide a signature of the multi-component scenario or perhaps some useful correlation with other observables. This process could be enhanced (and therefore efficiently constrained) for small m X which on the other hand helps to enhance ψ ± self-interaction.

-Indirect detection
Similarly indirect detection experiments, besides standard X X → SM and ψ ± ψ ± → SM contributions receive also more interesting one ψ + ψ − → Xh i followed by h i decays. -Colliders e + e − colliders provide a clean environment that might be used to test the multi-DM scenario considered here. Namely one can investigate the process e + e − → Z * → h i Z → χ D χ D Z with χ D = X or ψ ± . Energydistribution of Z might be adopted to gain some information on the invisible objects being produced. Initial estimation indicates that in some regions of the parameter space, for sufficiently large luminosity one should be able to disentangle 1-and 2-3 component scenarios.

Self-interacting DM
It is well known that the cosmological small-scale structure problems, such as the 'cusp vs. core' and the 'too-big-tofail' problems could be ameliorated if DM self-interaction was sufficiently strong at the dwarf galaxy scale [14,16,[102][103][104][105][106][107], the required value of the cross-section is where σ T ≡ dΩ(1 − cos θ)dσ/dΩ is the so-called momentum transfer cross section between DM particles. However, DM self-scattering cross-section as large as σ T / m DM 10 cm 2 /g turns out to be disallowed by observations at the cluster scale with the typical constraint σ T /m DM < 1 cm 2 /g [5,[108][109][110][111]. Therefore in the following we will try to find a region in the parameter space where 0.1 cm 2 /g < σ T /m DM < 1 cm 2 /g. ss A possible strategy that may generate large DM self-interaction is to introduce a mediator which is much lighter than the DM particles. In the VFDM model, there are two options, the mediator could be either h 2 or DM component X μ . As shown in [85] the choice of light h 2 implies number of severe constraints therefore here we will focus on the case of light vector DM component, which may serve as a mediator in elastic ψ + ψ − scattering. The main contribution to the amplitude for ψ + ψ − → ψ + ψ − comes from the t-channel X μ -exchange. The transfer cross-section for the two Majorana eigenstates interacting with a vector mediator was discussed in [112]. In case of the small mass splitting m + − m − m D , we can use the result obtained in a Born approximation for the Dirac fermion [113] The above perturbative result is valid if g 2 X m − /m X < ∼ 16π and (m + − m − ) m D . The substantial enhancement could In addition in case of small mass splitting the contribution from t-channel h 2 exchange is suppressed. For m − m X also the s-channel X μ exchange could be neglected. Here we consider the scenario with three stable components, therefore In the region of parameters considered here abundance of ψ ± and X μ might be comparable, however for instance domination of ψ ± could also be reached by facilitating X X annihilation by assuming m 2 < m X , then appropriate g X might be adjusted to tune the proper total DM abundance. Note that the Yukawa coupling y X remains small so that the potential relevance of h 2 mediation is therefore limited. If (m + − m − ) m D then both indirect detection of ψ ± ψ ± annihilating at present time and the cross-section for ψ ±nucleon scattering would be sufficiently suppressed. The mixing angle is as usually assumed to be small, sin θ 1, what provides additional suppression of both direct and indirect detection. Concluding, it seems that there exists the region of parameter space consistent with the data and providing substantial self-interaction of DM components and large ratio of masses (m D m X ) for DM components. This illustration of possibility for enhanced self-interaction in our model is located in a region of parameter space, which is similar to the limit considered in Sect. 3.2.3, i.e. for fermionic Dirac DM and stable vector mediator discussed very recently in [98].
In the Fig. 12, we show results of detailed scans focused on that region. The relic abundance was calculated using micrOMEGAs code [91] by placing ψ + and ψ − in one dark sector and X in another. Here we assume that both fermions are kept in equilibrium with each other by the efficient exchange processes ψ ± X ↔ ψ ± h 2 . The scan was made over masses in the range m X ∈ [1,15] MeV, m D ∈ [1, 10] GeV for fixed values of m h 2 ∈ {1, 2, 5} MeV and sin α ∈ {10 −5 , 10 −6 , 10 −7 }. Parameter g x was fitted imposing the condition that density of fermions satisfies relic abundance constraint whereas contribution of X is negligible. The latter is achieved by the effective annihilation X X → h 2 h 2 if h 2 is lighter than X . We choose the mass splitting m + − m − = 10 −5 GeV. In this way we can ensure that both states are present with nearly the same relic abundance. Moreover as it leads to the suppression of the Yukawa coupling, therefore we can avoid the indirect detection bounds. Another strong constraint comes from the limit on the invisible Higgs decay h 1 → h 2 h 2 . It results in the bound sin α ≤ 10 −5 , which on the other hand suppresses the DM-nucleon scattering cross-section to the range which is in agreement with direct detection experiments.
Similar scenario was discussed in the case of Dirac fermion in [98]. Since here we focus on the small mass splitting therefore our model effectively also contains a Dirac dark fermion and a stable vector as in [98]. Therefore, our results for σ T /m D shown in Fig. 12 indeed agree with those obtained in [98]. Note however, this accordance takes place only in this particular region of the parameter space, while in general the models are quite different, for instance, by the presence of Yukawa interactions that are allowed in our model due to specific assignments of dark charges. More comprehensive analysis of our model will be presented elsewhere.

Summary and conclusions
Multi-component dark matter scenarios are natural extensions of a simple WIMP dark matter. They predict more than one stable component in a dark sector and therefore they constitute a much richer dynamical structure. In this work we have presented a minimal UV-complete vector-fermion DM model with two or three stable particles. Its dynamical properties were discussed. Our vector-fermion DM model involves a dark sector with a U (1) X gauge symmetry. The dark matter contents are the dark gauge boson X μ , a Dirac fermion χ , and a complex scalar S, all are charged under the dark U (1) X gauge symmetry and are neutral under the SM gauge symmetry. Moreover, all the SM particles are neutral under the dark U (1) X gauge symmetry. The dark sector communicates with the visible sector (SM) through the Higgs portal κ|H | 2 |S| 2 . To generate the dark gauge boson X μ mass we have employed the Higgs mechanism in the dark sector.
After the dark sector spontaneous symmetry breaking and mass diagonalization, our vector-fermion DM scenario comprises a dark vector X μ and two dark Majorana fermions ψ ± . Out of eight free parameters of the model, the SM Higgs vev v = 246 GeV and the SM-like Higgs mass m h 1 = 125 GeV are fixed which leaves us with six independent parameters. We have chosen the physical basis where the six independent parameters are four masses m X , m ± , m h 2 , the mixing angle sinα, and the dark gauge coupling g x . To guarantee perturbativity we assumed g x ≤ 4π . We have employed sinα ≤ 0.33, which is consistent with the 2σ constraint from current measurements of the SM-like Higgs boson couplings to the SM gauge bosons at the LHC. Our VFDM model has an exact charge conjugation symmetry and the dark gauge symmetry which result in an accidental discrete Z 2 × Z 2 symmetry. The charge assignments under this Z 2 × Z 2 symmetry are: X μ (−, −), ψ + (−, +), ψ − (+, −) and h 1,2 (+, +) (also all SM gauge bosons and fermions are even under both discrete symmetries). The dynamics of the dark sector is mainly controlled by the gauge coupling g x which couples the three dark fields, i.e. X μψ+ γ μ ψ − .
In this work, we have analyzed the dynamics of the dark sector in the thermal freeze-out paradigm by solving the three coupled Boltzmann equations for the dark sector species (X μ , ψ + , ψ − ). The vector-fermion DM dynamics turns out to be different in many ways than the standard single component WIMP dark matter scenarios. In our model, depending on the masses of the dark sector particle m X , m + , m − there are the following three distinct cases where either two or all three dark sector particles are stable. To demonstrate the importance of three stable dark matter particles, we have illustrated in Fig. 11 the case with two Majorana fermions nearly degenerate in mass, i.e. y 1, hence their standard annihilations are suppressed, due to small Yukawa couplings, and the semiannihilations are most important for their annihilations.
Note that all the points presented in our scans satisfy the total relic density Ω tot h 2 at 5σ as observed by PLANCK and also the 2σ direct detection bound from LUX2016. Moreover, to understand the dark matter dynamics we have shown the evolution of the yields of dark matter components for each of the above cases for selected benchmark points, supplemented by tables containing all cross-sections for processes involved in collision terms. Also we compared our results for two-component cases with those obtained from the micrOMEGAs code [91] and satisfactory agreement has been found, see Figs. 7 and 9.
We have also discussed limiting cases of the model that are realized in appropriate regions of the parameter space. One of them corresponds to a model with Dirac fermion DM and stable vector mediator, this is an interesting scenario. A possibility of self-interacting DM has also been addressed and the region of parameter space where σ T /m DM can be substantially enhanced has been found.
To summarize, the absence of any direct, indirect or collider signatures of dark matter suggests a direction that leads beyond the single component WIMP-like dark matter. In particular, multi-component dark matter scenarios offer very rich dynamical structures which could solve current dark matter puzzles. In this work, we have presented a minimal renormalizable vector-fermion dark matter model where the presence of gauge symmetry and charge conjugation in the hidden sector implies the existence of two or three stable (vector and/or Majorana fermions) dark matter particles. The dynamics of the dark sector in our model is primarily controlled by a single parameter, the dark gauge coupling g x through the interaction ψ + γ μ ψ − −ψ − γ μ ψ + X μ which connects all dark sector states X μ , ψ + and ψ − . Such an interaction allows semi-annihilation and decay processes within the dark sector. We have explored the parameter space of our two/three-component VFDM scenarios requiring the correct total relic density and compliance with the current direct detection bounds.
We review here the derivation of the Boltzmann equation for the evolution of number density for the multi-component dark matter in the homogeneous and isotropic universe. Let us consider a generic dark sector which allows interactions with the visible sector and the dark sector, i.e., annihilations (co-annihilations), semi annihilations, conversions, and decays. We can write the Boltzmann equation for χ i particle as dn i dt where H is the Hubble expansion parameter, whereas, C i and D i are collision and decay terms for the χ i particle with the visible sector and the dark sector. In general, the collision term C i may involve 2 → 2, 3 → 2 and similar scattering processes, however, in order to keep the discussion simple, we focus only on 2 → 2 processes. We assume that all the dark components have the same temperature as the thermal bath. Similarly, the decay term D i may involve more than 2-body decays however we limit our-self to 2-body decays, the generalization is straightforward. For 2 → 2 scattering processes and the 2-body decay processes, we can write down the collision and the decay terms as follows, where the summation is over all possible interactions of χ i with the visible sector as well as dark-sector particles. The collision C i j→kl and the decay D i→ jk terms read: where the phase space integrand is, with |M i j→ jk | 2 and |M i→ jk | 2 being the matrix element squared summed over initial and final spins for the reaction i j → kl and i → jk, respectively. Above the factors of the form (1 ± f i ) are due to the spin statistics, the plus sign for bosons and the minus sign for fermions. Here f i denotes the distribution function of a given kind of particles, connected with the number density as follows: with g i being the number of spin degrees of freedom.
Hereafter it is assumed that appropriate symmetry factors for initial [114] and final 7 states are included in |M| 2 . We adopt the following assumptions: -Time reversal (T) invariance holds, so the amplitudes satisfy, where m is the mass of dark matter species, T is temperature, and μ i is the chemical potential), so that the Bose-Einstein (for bosons) and the Fermi-Dirac (for fermions) distribution functions could be approximated by the Maxwell-Boltzmann distribution functions, -In the absence of quantum degeneracies (which is assumed since the particles form a very dilute gas), in (A.3-A.4) the blocking and stimulated emission factors can be neglected, so 1 ± f i 1 will be adopted, -The initial chemical potentials are negligible, -Standard Model particles are in thermal equilibrium with the thermal bath, -Scattering processes with the thermal bath enforce kinetic equilibrium (also after decoupling and out of chemical equilibrium), so that phase-space distribution functions for particles involved in the collision satisfy [115] wheref i (E) is the thermal Maxwell-Boltzmann equilibrium distribution function for zero chemical potential and With the above assumptions we can rewrite the above collision term as, 14) The thermally averaged cross section is defined as, where the equilibrium number densityn i is defined as, After integrations over the momenta and changing variables we can rewrite the above equations as whereσ i j→kl (s) is the total cross-section averaged over initial and summed over final spins (= σ i jkl /(g i g j )) while K 1,2 are the Bessel functions of second kind and p 2 i j is a square of the incoming particle momenta in the center of mass frame, 19) which is related to the Møller velocity (A.13) by Following similar steps as above for calculating the C i jkl function (A.3), one can calculate the D i→ jk function (A.4). Assuming T invariant amplitudes for the decaying process, M i→ jk = M jk→i and using the fact that thermal equilibrium distributions satisfy the following relation, we can rewrite the decay contribution as, × |M i→ jk | 2f ī n i n i −n i n j n k n jnk . (A.22) The decay width is defined as usually where |M i→ jk | 2 is a matrix element squared summed over initial and final spins. After performing the integration over the outgoing momenta, D i→ jk = − Γ i→ jk n i −n i n j n k n jnk .
(A. 24) where thermally averaged decay rate Γ i→ jk is defined as, where K 1,2 are the Bessel functions andΓ i→ jk is the width averaged over the initial and summed over final spins, i.e. Γ i→ jk ≡ Γ i→ jk /g i .

Appendix B: Direct detection of multi-component dark matter
In the MCDM scenario, the standard direct detection bounds given by experimental groups in terms of DM-nucleon scattering cross-section and DM mass cannot be imposed, unless one of the components is responsible for nearly all recoil events [57,62,116]. Furthermore, the combined differential event rate in multi-component case may have a distinctive shape, which allows to discriminate it from single-component scenario [117]. In a general case, one has to confront the theoretical predictions with the results of experiments to put a constraint on the parameter space of the model. In our analysis, we follow [118]. The differential recoil event rate for a given DM component i can be written as [119] where ρ i = Ω i /Ω tot × 0.3GeV/cm 3 is the local density of that DM component, σ i N is its nucleus scattering crosssection, μ i N is the reduced mass of DM-nucleus system, F is the nuclear form factor, which we take as the conventional Helmi form and the function η i is a mean inverse speed of the DM particles in the local earth frame For the velocity distribution f G (v) in our Galaxy we use a truncated Maxwell-Boltzmann distribution with v esc = 550 km/s.
where v 0 = 220 km/s is the mean DM velocity relative to galaxy and N esc is the normalization factor. The distribution of DM as observed from the Earth takes into account its velocity v e relative to the galactic halo rest frame (B.29) The total d R/d E R differential recoil event rate is obtained by summing the rates for all DM components. Various DM direct detection experiments measure different kinds of detection signals, eg. prompt scintillation signal S1, ionization charge signal S2, the electron equivalent energy or energy released in photons. To put a constraints on a region of DM parameter space, one has to compute the expected experimental signal from the recoil event rate d R/d E R of multi-component DM obtained above. We focus on the predictions for the S1 signal measured by LUX experiment [6]. Following [118] we count the number of events N in the signal range S1 ∈ [S1 a , S2 b ] as described in [120] N [S1 a ,S1 b ] = E x S2 b S1 a d S1 ∞ n=1 (S1)Gauss(S1|n, σ ) where additional S2 efficiency S2 (E R ) = θ(E R − 3keV) is cutting the recoil energies from below and ν(E R ) = g 1 L y E R is the averaged expected number of photoelectrons from a given recoil event, which is calculated based on LUX gain factor g 1 = 0.0985 and photon yield L y adopted from the middle plot of Fig. 1. in [121]. The Poisson distribution gives the probability of obtaining n photoelectrons, which in the detector produce signal S1 normally distributed around n with σ = n(σ 2 PMT + g 1 ), where σ PMT is the single-photon resolution [122]. We include also the detector efficiency for events passing analysis selection criteria (S1), taken as a black curve from Fig. 2 in [6], and calculate the expected signal taking into account the total exposure E x = 4.47×10 4 kg × days.
We assume that all candidate events observed by LUX agree with the background-only model, using S1 a = 1 and S1 b = 50 we constraint the number of events in the range N [1,50]