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$ 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_\mu$ the Higgs mechanism with a complex singlet $S$ is employed in the dark sector. Dark matter candidates are the massive vector boson $X_\mu$ and two Majorana fermions $\psi_\pm$. 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 confirmed the presence of dark matter (DM) in the Universe. 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 DM-nucleon 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, this is the renown WIMP paradigm. Unfortunately, it turns out that the WIMP scenario suffers from various difficulties when confronted with observations on small cosmological scales. For instance "too-big-to-fail" [8,9] and the "cuspcore" [10][11][12][13] problems are widely discussed in the literature. In particular, the DM densities inferred in 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 self-interacts 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.
DM could also be searched for in 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 ray, neutrinos, electrons, positrons, protons and antiprotons, could be observed on Earth and 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 ray, 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].
Since very different DM masses are needed to solve the small-scale problems (through self-interaction) and to interpret the excess of gamma rays, therefore in order to accommodate both observations, a multi-component DM seems to be a natural option.
Motivated by the above arguments we will consider first a generic multi-component dark matter (MCDM) scenario, and then a simple renormalizable model that also predicts 2 or 3 DM components. Of course, the dynamics of multi-component DM is much richer than that of a simple WIMP, and therefore attractive to study by itself, even without any phenomenological direct application. Various multi-component DM scenarios have been proposed and studied so far . Models of multi-component DM were also considered and adopted in astrophysical simulations, e.g. by [73][74][75][76].
The paper is organized as follows. In Sec. 2 a generic 2-3 component scenario is sketched. Based on the generic setup, in Sec. 3 a vector-fermion (2-3 component) model of DM is presented and solutions of Boltzmann equations are discussed. Section 4 contains summary and conclusions. Appendix A collects details of the derivation of Boltzmann equations, while the appendix B describes the method adopted to obtain constraints for a multi-component DM model by direct detection experiments.

A generic 2-3 component dark matter
This section is supposed to convince readers that setups with two or three dark components offer physics much richer than simple WIMP scenarios which only involve annihilation processes to SM. Let us assume that there exists a dark sector composed of particles with possibly different spins, either stable or unstable depending on their quantum numbers and masses, which are singlets of the SM gauge group. Stability of DM candidates will be ensured by a discrete Z 2 ×Z 2 symmetry 1 . We introduce two states, χ andχ odd with respect to Z 2 and Z 2 , respectively, any of them may serve as DM candidate. In fact, it is natural to introduce also a state odd under both symmetries, here it is denoted byφ. All the SM states are even under Z 2 ×Z 2 and for convenience we represent them by φ. Note that this setup allows for a non-trivial coupling between fields that belong to the dark sector: χχφ is invariant under Z 2 ×Z 2 . This scenario is summarized in Tab. 1.
The following primary vertices for interactions involving members of the dark sector Table 1. Quantum numbers under stabilizing discrete symmetry Z 2 ×Z 2 . The fields χ,χ andφ belong to the dark sector, whereas φ denotes mediator/SM fields.
are permitted by the symmetry: Note that we have not listed quartic vertices above. In case of the scalar particles, there could be dimension-4 primary quartic vertices, however we assume their relevance is small. In Fig. 1, we show all possible Feynman diagrams built from vertices presented in (2.1), which are relevant for collision terms in Boltzmann equations describing evolution of number densities of dark particles. There appear the following four classes of processes contributing collision terms which are relevant for the dark matter abundance: where φ represents the mediator/SM particles.
Depending on their masses, the setup discussed here contains 2 or 3 stable dark particles. The vertex g χχφ given in Eq. (2.1) could allow for decays of any of the three involved particles if its mass is larger than the sum of masses of other two particles. As a result there are three distinct cases where one of the dark sector particles will be unstable and hence the remaining two particles would be stable providing a 2-component dark matter scenario. Following are these 2-component dark matter cases: In the cases when the above conditions are not met, then all of the three particles would be stable and hence constitute a 3-component dark matter scenario. A schematic diagram showing all the above 2-and 3-component cases in the mχ − mφ plane is shown in Fig. 2. Even without solving Boltzmann equations for the setup considered here some generic observations could be made:  terms for corresponding Boltzmann equations are much more interesting, e.g. semiannihilation and conversion between dark elements are possible.
• 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 χχ → φφ is suppressed then χ can still disappear, for instance, through χφ →χφ followed by unsuppressed annihilation of pairs ofχ. In other words,φ can work as a catalyzer that enables disappearance of χ. In this case, it is possible that the presence of other (χ orφ) dark components might be crucial for the determination of the asymptotic abundance of the major DM element. Also, semi-decays may play a relevant role in the determination of the final abundance.
• Standard WIMPs decouple from thermal equilibrium at m/T ∼ 20 − 25, therefore heavy states decouple earlier. 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 the following part of this paper, we are going to illustrate the above generic observations within a well-defined renormalizable 2-3 component DM model.

Vector-fermion multi-component dark matter
After gaining some intuition from the model-independent discussion in the previous section, here we are going to construct a specific model of 2-3 component dark matter, discuss its properties, and investigate solutions of corresponding Boltzmann equations.
We explore the possibility of having a multi-component 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. We consider the simplest 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 vacuum expectation value (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 under the gauge group 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 vector-fermion 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.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 3.1), respectively. Moreover, in Eq. (3.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. Note that our model is an extension of the model presented in [77,78] and is a gauged version of the model recently considered by Weinberg [79]. 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 (3.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, (3.12) The above mass squared matrix M 2 can be diagonalized by the orthogonal rotational matrix R, such that, where (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 (3.14) 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-eingenstates ψ ± (= ψ c ± ) with mass eigenvalues, m ± = m D ± yv x . (3.20) 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 Tab. 2. It is easy to see that the above Z 2 is a direct consequence of the charge conjugation symmetry (3.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 [80]. 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 under which the three dark matter components are charged. The Z 4 charges are: 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 Z 2 ×Z 2 .
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 ψ ± ): As depicted in Fig. 3, depending on the masses of dark particles, there are three cases 3 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 Majorna fermion ψ + will decay into a stable vector X µ and a stable Majorana fermion ψ − . This is a 2CDM case, the white area (left) in Fig. 3. We consider m + to have a fixed value and the gray region represent parameter space where the all three dark sector particles (X µ , ψ + , ψ − ) are stable, whereas the white regions represent the 2-component scenarios where ψ + and X µ are unstable, respectively.
• 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. 3.
• 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. 3. 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 vector-fermion 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 (3.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.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 SM-like Higgs couplings with the gauge bosons. We employed | sinα| ≤ 0.33, the current LHC 2σ bound [81]. 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. 3.10) : (3.27) Note that the Yukawa coupling y x ∝ 1/v x = g x /m X , therefore for fixed m X the Yukawa coupling y x starts to be g x dependent. In other words, for fixed m X , the vev v x must vary if g x is being changed, implying variation of y x .

Boltzmann equations for the vector-fermion dark matter
In this subsection we present coupled Boltzmann equations for the evolution of number density of dark matter particles (X µ , ψ + , ψ − ) in VFDM model. In Fig. 4, we collect all the relevant vertices that describe interactions of the dark segment of the theory. Figures 5, 6 and 7 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:  Figure 6. Semi-annihilation diagrams for the dark particles X, ψ ± . These processes annihilate two of the dark sector particles to a dark sector particle and a SM particle (h 1,2 further decays to other SM particles). Figure 7. Dark matter conversion processes involving X and ψ ± . These processes are important to keep thermal equilibrium within the dark sector.
where σ ijkl v Møl ≡ σ ijkl v Møl is the thermal averaged cross-section for the process ij → 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 . 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]. In order to calculate matrix elements needed for collision terms in the Boltzmann equations, we are using CalcHEP [82]. For thermal averaging and solutions of the Boltzmann equations, we adopt a dedicated C++ code.
Note that there is a clear correspondence between dark matter components considered in the generic discussion contained in Sec. 2 and those that are present in our vector-fermion Following our discussion from Sec. 2, there are two interesting possibilities: 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 XX 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 (in the generic setup of Sec. 2, g x corresponds to g χχφ ). 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 for 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 semi-decay -would have been allowed. Again not a very appealing scenario.

2CDM: a vector and a Majorana fermion as dark matter
In this subsection 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. 8 • Ω 1cdm x Ω 2cdm x Ω 2cdm Ω 1cdm x Ω 2cdm x Ω 2cdm  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 underabundant, 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. 9.
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σ. 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). The right panels correspond to parameter points such that around m X = 200 GeV there is a resonant enhancement of XX annihilation through the h 2 s-channel exchange. In this case, in order to satisfy the relic abundance constraint relevant coupling constants must be small, i.e. g x ∼ 0.03 − 0.13. Note that small g x allows satisfying also the direct detection constraint. Lower panels in Fig. 8 show the abundances of both components separately, Ω 2cdm X and Ω 2cdm − , for our model and corresponding Ω 1cdm X calculated in a 1-component vector dark matter (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 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. 9, 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 (m X , ∆m) = (500, 100) GeV  Figure 10. 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. 11. normalization adopted. Moreover, the tables in Fig. 9 and similar figures in the following section contain first two non-vanishing coefficients of thermally-averaged cross-sections [pb] expanded in powers of x −1 , given by σ ijkl v Møl = a N x −N + a N +1 x −(N +1) + · · · , and the decay width Γ ψ + →Xψ − [GeV]. The plots in Fig. 9 illustrate solutions of the Boltzmann equations for three selected sample points. The middle panel show solutions for parameters that reproduce correct total DM abundance and also satisfies the direct detection LUX limits, so the corresponding point is also present in the scan results shown in the upper-left panel of Fig. 8 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 [83]. As it is seen in the plots they agree very well with the solid lines which correspond to solutions obtained from our dedicated C++ code that solves the set of three Boltzmann equations (3.29 -3.31). In order to identify the most important processes for a given parameter set, in the tables below the panels in Fig. 9 we collect the first two non-vanishing coefficients in the expansion of thermally averaged cross-sections in powers of x −1 . As one can read, for instance, from the middle table, in this case the semi-annihilations (ψ + h 2 → Xψ − , Xψ + → ψ − h 2 ) and conversions (ψ + ψ + → ψ − ψ − ) are very relevant.

2CDM: two Majorana fermions as dark matter
In this subsection we show results for the scenario with m X > m − + m + , so X µ is unstable and can decay into the Majorana fermions ψ − and ψ + . Figure 10 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 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 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. 10 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. 10 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. 11, we illustrate solutions of the Boltzmann equations for three sample points in the parameter space. The middle panel corresponds to correct abundance and 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. 10 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 subsection we show results for the scenario with m + + m − > m X > m + − m − , so all the three dark components are stable. Figure 12 contains results of a scan performed adopting our dedicated code 4 that solves the set of the three Boltzmann equations (3.29 -3.31). All the points presented in Fig. 12 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. 12 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. 12, 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 Note that for these points the dark fermion Yukawa coupling is very small, i.e. y 1 as a result the direct annihilation processed for ψ ± are inefficient therefore the main annihilation processes are through the semi-annihilation processes and conversion to X with further annihilation to SM. 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 13 shows solutions of the Boltzmann equations for three sample points in the parameter space of 3CDM. The middle panel represents the point marked as a black star in the scan results of Fig. 12 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.

Summary and conclusions
It is very plausible that the dark sector of our Universe has much richer structure and dynamics than that of a simple WIMP. In particular, it is possible that the dark sector consists of more than one dark species. Motivated by such an option, we have presented a case for a multi-component dark matter (MCDM) scenario and analyzed its dynamical properties. In order to discuss the properties of a MCDM scenario we started our discussion in Sec. 2 with a generic setup involving three dark sector states χ(−, +),χ(+, −) and φ(−, −), charged under a discrete symmetry Z 2 ×Z 2 . Depending on masses and interactions within the dark sector there are two or three stable particles. The presence of more than one stable particle in the dark sector implies a very rich dynamics which allows -apart from the standard annihilations -semi-annihilation, conversion, and decay processes. In particular, if the standard annihilation processes are suppressed for a given dark particle, the primary source of its disappearance could be semi-annihilation or conversions processes which involve other dark matter components. Hence the presence of other dark matter components is crucial to forbid the over-abundance of dark matter. On the other hand, the presence of more than one dark state can also be handy in models where single dark matter state gives under-abundance such that the relic from the other components can add up to give the correct observed relic density of the Universe.
In order to illustrate the outlined generic properties of MCDM scenarios, in Sec. 3 we have considered a minimal vector-fermion UV complete model which 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 employed the Higgs mechanism in the dark sector. In summary, after the SSB and mass diagonalizations, the dark sector comprises of a dark vector X µ and two dark Majorana fermions ψ ± which interact with the SM via the Higgs portal. Out of eight free parameters of the theory, 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 . Out of these six independent parameters, except mixing angle sinα, all are free parameters (g x ≤ 4π due to perturbativity). We employed sinα ≤ 0.33, which is consistent with the 2σ from the current measurements of the SM-like Higgs couplings to gauge bosons at the LHC. Our vector-fermion dark matter model has an exact charge conjugation symmetry and the dark gauge symmetry which result in the discrete Z 2 ×Z 2 symmetry, as in the generic case of Sec. 2. The charge assignments are: X µ (−, −), ψ + (−, +), ψ − (+, −) and h 1,2 (+, +) (also all SM gauge bosons and fermions are even under this discrete symmetry). The dynamics of the dark sector is mainly controlled by the gauge coupling g x which couples the three dark fields, i.e. X µψ+ γ µ ψ − .
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 MCDM dynamics turns out to be different in many ways than the standard single component WIMP dark matter scenarios. In our VFDM 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. (iii) m + + m − > m X > m + − m − : a 3CDM scenario where all three dark sector particles (X µ , ψ + , ψ − ) are stable. As in the 2CDM cases we have performed scan over m − , g x for different choices of m X , ∆m, m h 2 and sinα.
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 evolution we have shown the yields of dark matter components for each of the above cases, which are supplemented by tables containing all the crosssections of the processes involved in the collision terms. We conclude this work by emphasizing some of the most interesting features of our MCDM model: • The presence of more than single dark matter state allows semi-annihilation, conversion and decay processes, which could play a vital role in the disappearance of a dark matter state for which the standard annihilations are suppressed or not allowed.
Hence an MCDM scenario can help to suppress an over-abundance in single dark matter cases.
• MCDM scenarios can also help to understand the absence of any direct or indirect signatures of the dark matter. For instance, if the dominant dark matter component has negligible interactions with the visible sector (i.e. standard annihilation processes are suppressed or absent) then, as noticed above, its disappearance is possible via semiannihilations or conversions such that it makes the dominant component undetectable by either direct or indirect experiment.
• The presence of multi-component dark matter states means that their relic densities add up to reproduce the observed density of the Universe. In other words, if a single component dark matter model gives under-abundance then the presence of other components can contribute to the relic density to give correct observed value.
To summarize, the absence of any direct, indirect or collider signatures of dark matter suggests a direction that leads beyond the single component dark matter. In particular, multi-component dark matter scenarios offer very rich dynamical structures which could help to solve the dark matter puzzles and improve our understanding of the dark sector.

A Boltzmann equations for the multi-component dark matter
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 th particle as follows, where H is the Hubble expansion parameter, whereas, C i and D i are collision and decay terms for the χ i th 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 on only 2 → 2 process. 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 ij→kl and the decay D i→jk terms read where the phase space integrand is, with |M ij→jk | 2 and |M i→jk | 2 being the matrix element squared summed over initial and final spins for the reaction ij → 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 [84] and final 5 states are included in |M| 2 . We adopt the following assumptions: • Time reversal (T) invariance holds, so the amplitudes satisfy, M ij→kl = M kl→ij and M i→jk = M jk→i , • m T , (m i − µ i )/T 1 (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 [85] f 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, (A.10) The thermal equilibrium distributions satisfy the following relation due to the conservation of energy,f After performing the integration over the outgoing momenta, the collision term can be written as, where v Møl is the Møller velocity 13) and the total cross-section summed over initial and final spins (A.14) The thermally average cross section is defined as, where the equilibrium number densityn i are defined as, After integrations over the momenta and changing variables we can rewrite the above equations as whereσ ij→kl (s) is the total cross-section averaged over initial and summed over final spins (= σ ijkl /(g i g j )) while K 1,2 are the Bessel functions of second kind and p 2 ij is a square of the incoming particle momenta in the center of mass frame, 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 average 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 .

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 [32,86,87]. In general case, one has to confront the theoretical predictions with the null results of experiments to put a constraint on the parameter space of the model. In our analysis, we follow [88]. The differential recoil event rate for a given DM component i can be written as [89] dR where ρ i = Ω i /Ω tot × 0.3GeV/cm 3 is the local density of that DM component, σ iN is its nucleus scattering cross-section, µ iN 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 The total dR/dE 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 dR/dE R of multicomponent DM obtained above. We focus on the predictions for the S1 signal measured by LUX experiment [6]. Following [88] we count the number of events N in the signal range S1 ∈ [S1 a , S2 b ] as described in [90] N [S1a,S1 b ] = Ex S2 b S1a dS1 ∞ 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 [91]. 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 [92]. 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 Ex = 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]