A multi-component SIMP model with $U(1)_X \rightarrow Z_2 \times Z_3$

Multi-component dark matter scenarios are studied in the model with $U(1)_X$ dark gauge symmetry that is broken into its product subgroup $Z_2 \times Z_3$ \'{a} la Krauss-Wilczek mechanism. In this setup, there exist two types of dark matter fields, $X$ and $Y$, distinguished by different $Z_2 \times Z_3$ charges. The real and imaginary parts of the $Z_2$-charged field, $X_R$ and $X_I$, get different masses from the $U(1)_X$ symmetry breaking. The field $Y$, which is another dark matter candidate due to the unbroken $Z_3$ symmetry, belongs to the Strongly Interacting Massive Particle (SIMP)-type dark matter. Both $X_I$ and $X_R$ may contribute to $Y$'s $3\rightarrow 2$ annihilation processes, opening a new class of SIMP models with a local dark gauge symmetry. Depending on the mass difference between $X_I$ and $X_R$, we have either two-component or three-component dark matter scenarios. In particular two- or three-component SIMP scenarios can be realised not only for small mass difference between $X$ and $Y$, but also for large mass hierarchy between them, which is a new and unique feature of the present model. We consider both theoretical and experimental constraints, and present four case studies of the multi-component dark matter scenarios.


Introduction
For a few decades or so, the Weakly Interacting Massive Particle (WIMP) paradigm has been one of the mainstreams of dark matter (DM) physics (see, for example, Refs. [1,2] for a recent review). Within the WIMP scenarios, DM mass could be of the electroweak scale and its interactions with the Standard Model (SM) particles have the weak-interaction strength, and thermal relic density could be realised by freeze-out mechanism. One would have anticipated to observe DM particle by (in)direct detections and at colliders. However the null results from the LHC searches for new particles and direct detections of DM with mass in the electroweak scale lead us to consider more seriously new paradigms beyond the WIMP paradigm.
As of now, parameter spaces for most of the single-component DM models are being tightly constrained by various experiments. Single-component DM models are simple for theoretical ideas on the model building, phenomenological analyses, and comparison with various experimental data. However, there is no compelling reason that the DM density of the current Universe is composed of only one component. For example the visible sector of the Universe contributes only ∼ 5% to the energy density of the current Universe, and it has a number of stable or long-lived particles; electron, proton, photon, and three species of active light neutrinos. Presently the energy density of the dark Universe in the form of matter is about five times more than the visible Universe. Therefore it would be an oversimplification to assume that the dark sector of the Universe consists of a single-component DM. Even if the dark sector of the current Universe is dominated by a single-component DM, there could be more DM species in the earlier Universe that would modify the evolution of the early Universe. They may leave some footprints that can be observed in the current Universe, e.g. on Cosmic Microwave Background or matter power spectra.
One important feature that has to be taken into account seriously in DM physics is that DM particle should be absolutely stable or its lifetime should be much longer than the age of the Universe. Usually the DM stability is supposed to be guaranteed by some dark symmetries such as Z 2 . If this dark symmetry is a global symmetry, it will be violated by gravity effects. For example, let us consider a real scalar DM model with dark Z 2 symmetry under which the real scalar DM field S is Z 2 -odd, namely transforming like S → −S, whereas all the SM fields are Z 2 -even. Then the DM S will be stable at renormalizable level. However there will be Z 2 -breaking dimension-5 operators such as Then S will decay too fast to be a good cold-DM candidate even for Λ ∼ M Pl where M Pl is the Planck mass, unless its mass is very light m S O(1) eV. The Z 2 breaking from dimension-6 operators will not be harmful at all, since one can make DM lifetime long enough for Λ M Pl , like the proton lifetime in the SM with B-and L-violating dimension-6 operators. One can evade this problem of dimension-5 operators that violate dark global symmetry by promoting the dark global symmetry to a local dark gauge symmetry. In quantum field theory, these could be achieved by assuming that DM carries a dark charge that is either exactly conserved or approximately conserved. Finally the symmetry associated with the exactly conserved charge is implemented into a local gauge symmetry [22] (see also Sec. IV B and C of Ref. [23] for more discussion on this issue). Alternatively, the approximately conserved dark charge may be an accidental symmetry of an underlying dark gauge theory, like the U (1) baryon number in the SM being an accidental symmetry of the SM and being broken only by dimension-6 operators. Thus, for either absolutely stable or long-lived DM, we end up with a local gauge theory 1 . This strategy was the cornerstone of the SM of particle physics, and it could be extended to DM model building. In DM models with a dark gauge symmetry, there are dark gauge bosons and dark Higgs bosons in addition to the DM particles with various dark charges. Furthermore, there could be more than one generation of DM particles, as in the SM. There could be more than one stable or long-lived particles in the dark sector in general.
One of the most interesting multi-component scenarios proposed so far would be the dark QCD (DQCD) with a light Nambu-Goldstone boson DM (dark pion DM). One can consider both WIMP [44][45][46][47][48][49][50][51] and SIMP scenarios [4,8,14,20] in DQCD models . In the multi-component SIMP models proposed so far in the literature, various DM components are related with each other. For example, dark pion DM with different dark flavours in DQCD models are related with each other by a dark flavour symmetry. They are not independent with each other when we consider their origin, and their masses cannot be arbitrarily different. In the case of the dark pion DM in DQCD, for example, one assumes that the dark pions are pseudo-Nambu-Goldstone bosons. This picture requires that m q i + m q j Λ dc , where Λ dc is the dark confinement scale, which is nothing but the good old partiallyconserved-axialvector-current (PCAC) condition, in order that dark pion DMs are good Nambu-Goldstone bosons. On the other hand, dark pion DM mass is given by m 2 π ∼ m q Λ dc , which is bounded from above by Λ 2 dc . Therefore, one cannot imagine a wide separation in the dark pion DM masses with different flavours. They are bounded by ∆m 2 π ∼ Λ 2 dc . In other words, we can say that they are not genuine multi-component SIMP models where arbitrary masses are allowed, including the large mass hierarchy. For genuine multi-component DM models, different DM components should have different masses (and different lifetimes in the case of a decaying DM), different (dark) charges, and carry even different spins.
In this paper we introduce two independent fields X and Y with different dark charges (and even different spins) from the beginning, and show that two-or three-component SIMP scenarios can be realised. Therefore the model we present in the following could be taken as a real multi-component SIMP scenario. Our study is not based on the effective field theory approach. The model constructed in this work is renormalizable and gauge invariant under the SM and the dark gauge symmetry groups, and thus a UV-complete one. Our conclusion can thus be taken as reliable in the entire energy range considered.
At this point, it is worthwhile to remind ourselves that WIMP and SIMP are about mechanisms to reproduce observed DM relic density in a given DM model, and not about the underlying particle physics models for DM. Within the same DM model, different mechanisms can be realised depending on the parameter values and mass scales. Likewise, either single-component or multi-component WIMP or SIMP scenarios can be realised within the same DM model, depending on the parameter values and mass scales of various particles in the model. The model discussed in the following shows how rich the dark sector can be, depending on the parameter spaces (couplings and mass scales), and how DM particles are stabilised (U (1) dark charge assignments in our case). The multi-component SIMP scenarios are realised only in a particular corner of the parameter space. In a generic parameter space, we expect that single-component SIMP scenarios or WIMP scenarios would be realised, but these two options shall not be discussed in detail, since they are not qualitatively different from the literature of the SIMP/WIMP scenarios, and thus are not the main interest of this work.
In this work we consider a multi-component DM model with U (1) X dark gauge symmetry broken into its product subgroup Z 2 × Z 3 á la Krauss-Wilczek mechanism [52]. There exist two types of DM fields, X and Y , which are distinguished by different U (1) X charges with the following symmetry breaking pattern: U (1) X → Z 2 × Z 3 . The stability of X is guaranteed by the unbroken Z 2 symmetry, while the unbroken Z 3 symmetry ensures the stability of Y . Then we can realise a multi-component SIMP scenario through various DM number-changing scattering processes such as Y Y Y → XX. Note that X can be either scalar or fermion, whereas Y should be scalar. The relic densities of different DM fields are determined through interactions between the dark fields. For the Z 2 -charged field X, number-changing 2 → 2 processes are responsible for X's relic density. On the other hand, the relic density of the Z 3 -charged field Y is given by number-changing 3 → 2 processes. As the processes responsible for the relic densities of the DM fields occur within the dark sector and both X and Y participate in the number-changing processes, we dub our model a multi-component SIMP DM scenario. The field X receives different masses for its real and imaginary parts from the vacuum expectation value (VEV) of the dark Higgs, breaking the U (1) X symmetry. Due to the interactions between X and Y , both X I and X R may contribute to Y 's 3 → 2 processes. It opens a new class of SIMP models. Depending on the mass gap between X I and X R , two-component or three-component DM scenarios can be realised. One of the most distinctive features of this model is that the different DM fields X and Y may have vastly different masses, while significantly contributing to the total DM relic density today. This is in sharp contrast to most of the studies on multicomponent SIMP where each component could coexist when the mass difference is small enough [9,15,20]. For completeness we consider both cases with large mass hierarchy and small mass difference.
It should be emphasised that the present work, which considers multi-component SIMP scenarios with a large mass hierarchy, is mostly motivated by theoretical reasons, as a kind of existence proof of such a case within a theoretically well-defined DM model. Still there is an interesting phenomenological consequence. The coupling used to maintain the kinetic equilibrium between the DM sector and the SM sector may also be used to detect the boosted Y (from X annihilation) in DM and neutrino detectors. The boosted DM signal is one of the most important predictions of the multi-component DM scenarios with large mass hierarchy.
The paper is organised as follows. We first describe our model setup in Section 2, where we set the notations, explain the U (1) X symmetry breaking, and identify DM candidates. In Section 3, we study phenomenology of the multi-component DM scenarios. The evolution of the DM system is described by four coupled Boltzmann equations, each for Y , X I , X R , and Z . We numerically solve the Boltzmann equations and present two benchmark cases for the two-component scenario and two benchmark cases for the three-component scenario. In Section 3.2, we consider theoretical constraints, such as perturbativity and unitarity, as well as experimental constraints, including LHC constraints for the invisible Higgs decay, the DM direct detection, and Z searches. One of the crucial constraints for a SIMP-type DM is the sufficient release of its kinetic energy. We address this issue in detail in Section 3.2. Taking into account all of the constraints we present our results in Sections 3.3 and 3.4. Finally we conclude in Section 4.

Model
We consider a dark U (1) X gauge symmetry with three different types of scalars in the dark sector: two DM candidates X and Y plus the dark Higgs field φ. We assign the U (1) charges of these fields as follows 2 : (q X , q Y , q φ ) = (1/2, 1/3, 1) . (2.1) Then, the gauge-invariant renormalizable Lagrangian of this model is given by where L SM is the SM Lagrangian without the Higgs sector, H is the SM Higgs field, and is the U (1) Y -U (1) X gauge kinetic mixing parameter. The covariant derivative of the dark sector is defined as D µ ≡ ∂ µ − iq i g XXµ , with g X being the dark gauge coupling and the U (1) X -charge assignments of dark fields q i (i = X, Y, φ) are defined in Eq. (2.1). The scalar potential is given by where v φ and v are VEVs of the dark Higgs and the SM Higgs fields. After U (1) X symmetry breaking by a nonzero VEV φ = v φ / √ 2, the last line of Eq. (2.3) still has Z 2 × Z 3 discrete gauge symmetries 3 : We note that µ Xφ has the mass dimension one. Expanding the dark Higgs and the SM Higgs around their VEVs, v φ and v, and writing X = (X R + iX I )/ √ 2, the real and the 2 We can also consider another possibility with spin-1/2 dark fermion ψ with q ψ = 1/2, instead of dark scalar X. We shall not pursue this case in detail in this paper, since the qualitative features of multicomponent SIMP scenarios would be the same. 3 See also Refs. [53][54][55] and Refs. [23,[56][57][58] for scalar DM models with U (1) dark gauge symmetry broken into Z2 and Z3, respectively. imaginary parts of X get different masses from the U (1) X symmetry breaking. The mass splitting between X R and X I is generated by the µ Xφ term [53]: Taking µ Xφ > 0 without loss of generality, we see that X I is always lighter than X R . The mass-squared matrix of the Higgs fields is given by The mass eigenvalues are then obtained as follows: The mass eigenstates, which we denote by h 1,2 , are related to the interaction eigenstates as with the mixing angle (2.10) With the U (1) Y -U (1) X gauge kinetic mixing, the mass matrix of the gauge boson can be diagonalised as follows. The kinetic mixing can be removed in terms ofB µ andX µ defined as [59], where we definedW µ for notational consistency. We then diagonalise the mass matrix for where c w ≡ cos θ w , s w ≡ sin θ w , c ζ ≡ cos ζ, s ζ ≡ sin ζ, and θ w is the Weinberg angle. Here ζ is defined as , (2.13) where mẐ = (g 2 1 + g 2 2 )v/2, mX = g X v φ , c ≡ cos , and s ≡ sin . Subsequently, we obtain the following relations between the original fields and the mass eigenstates: with t ≡ tan (similarly for ζ below). The masses of the Z and Z gauge bosons are given by In the following we assume small mixing angles, i.e. α 1 and 1, hence h 2 ≈ h and h 1 ≈ h , and m 2 Z ≈ m 2 Z and m 2 Z ≈ m 2X . In our model, due to the unbroken Z 2 × Z 3 symmetry, fields with U (1) X charges of 1/2 and 1/3 will be stable, becoming DM candidates. The X I field with q X = 1/2 is absolutely stable and makes one of the DM components. The heavier one X R decays into X I , Z ( * ) , and could be another good DM candidate, if its lifetime is long enough compared with the age of the Universe. When the mass splitting is so small, m X R − m X I min{2m Y , m Z }, that the X R decay channels such as X R → X I , Z and X R → X I , Y, Y * are closed, the X R field again becomes a good DM candidate (see Appendix A for the three-body decay of X R ). The field Y with q Y = 1/3 is another DM candidate because of the unbroken Z 3 symmetry. Then for example could be possible through the dark Higgs and/or dark gauge boson exchanges, which is a new class of SIMP models based on spontaneous breaking of dark We are primarily interested in the setup where the field Y is the lightest among the DM candidates. To implement the SIMP scenario for Y , we assume that the dark Higgs and the dark gauge boson are heavier than Y . Otherwise, the DM will be pair-annihilated into dark Higgs/dark gauge boson, which is not of interest in this work. Similarly, if either the mixing between the SM Higgs boson and dark Higgs or the kinetic mixing between gauge bosons of U (1) Y and U (1) X is sizeable, the DM pair-annihilating into SM particles will be important once kinematically open. We shall therefore focus on the parameter spaces that have small mixing between SM Higgs and dark Higgs, as well as small kinetic mixing in the gauge sector. We note that this setup is also useful to suppress the DM-nucleon scattering cross section, which is tightly constrained by the DM direct detection experiments. On the other hand, it is necessary for Y to lose its kinetic energy through elastic scattering with SM particles [3,7]. For this, a somewhat light Z is needed for the SIMP-SM kinetic equilibrium condition as we discuss in Section 3.2. We focus on the case where m Z > 2m Y ; thus Z is not stable and will decay. As long as the decay is fast enough, Z cannot be a DM candidate. For the dark Higgs, we assume that it is much heavier than Y , X I,R , and Z such that the dark Higgs does not play any significant role in the dynamics. The cross section for the process X, X → Y, Y * features a four-point vertex 5 which, in the large-m h limit, is proportional to λ 2 XY . The relic density of the X field will be diluted through this process if λ XY is large. When the dark Higgs is somewhat lighter, couplings λ Xφ and λ Y φ need to be set appropriately such that the X, X → Y, Y * is suppressed via the destructive interferences between the contact process and dark Higgs mediated process. We first assume that λ XY is small enough such that this contact X-annihilation process is not efficient in the large-m h limit. We then comment on the case of destructive interferences, where λ XY may take a larger value, presenting parameter sets with a lighter dark Higgs and larger dark Higgs-portal couplings.
In summary the DM candidates are X R , X I , and Y in our setup. Our interest is to investigate whether two-or three-component DM scenarios can be realised for the following four cases: The cases i) and ii) will be the two-component DM scenarios, while iii) and iv) correspond to the three-component scenarios. We do not consider the inverted mass hierarchy, i.e. m X I m X R m Y . On one hand, in order for the Z 3 -charged field Y to serve as a SIMP DM, the mass should be less than ∼50 MeV, in the large-m h limit, as number-changing 3 → 2 processes are dominant. On the other hand, a DM lighter than ∼5 MeV is heavily constrained by cosmology such as Big Bang Nucleosynthesis and Cosmic Microwave Background, as explored in Ref. [60]. Therefore, we do not consider the inverted mass hierarchy in this work, and focus on the aforementioned cases i)-iv).

Generalities
Under the assumptions described in the previous section, we have four coupled Boltzmann equations for describing the evolution of the system consisting of {Y, X I , X R , Z } 6 . 5 We note that contact interactions between Z2-charged field and Z3-charged field are absent at renormalizable level in the fermionic Z2 DM, where the scalar field X is replaced by a spin-1/2 dark fermion ψ. 6 Since the dark Higgs is assumed to be sufficiently heavier than the rest of the dark sector fields, it does not play any significant role in the dynamics.
Schematically the Boltzmann equation for the species i is given bẏ i n j n k n eq j n eq k − Γ i→j,k,l n i − n eq i n j n k n l n eq j n eq k n eq l − σv i,j→l,m n i n j − n eq i n eq j n l n m n eq l n eq m − σv 2 i,j,k→l,m n i n j n k − n eq i n eq j n eq k n l n m n eq l n eq m , where n i (n eq i ) is the (equilibrium) number density of the species i, H the Hubble parameter, and˙≡ d/dt with t being the cosmic time.
Here σv and σv 2 are the thermally-averaged cross sections, and Γ is the thermallyaveraged decay rate. The thermally-averaged cross section σv for the i, j → l, m process is given by where F and S are the spin and symmetry factors, E i,j the energy of incoming particle, and p f the 3-momentum of final-state particle. Similarly, the thermally-averaged cross section σv 2 for the i, j, k → l, m process is given by For the decay rate, see Appendix A. We utilise FeynRules [61] and CalcHEP [62] to compute the scattering amplitudes. We work in the nonrelativistic framework and take only s-wave contributions. Instead of n i and t, it is convenient to use the yield Y i and the dimensionless time variable x, defined as where s(T ) = (2π 2 /45)g * s T 3 is the entropy density with T being the temperature of the heat bath. Note that, during the radiation-dominated era, The full Boltzmann equations are summarised in Appendix B. We choose m i = m Y . The DM relic density is then given by where the subscript "0" indicates the present value, ρ c is the critical energy density, and h ≈ 0.67 is the scaling factor for the Hubble parameter.
The input parameters for the model under consideration are chosen as In the following, we fix m Z = 200 MeV, m h = 30 GeV, λ X = 0.025, and = 2 × 10 −4 . We assume that all the quartic and portal couplings are non-negative. Furthermore we assume that λ Y φ is positive. As explained in the previous section, λ XY is required to be small in order to implement multi-component DM scenarios in the large-m h limit. We first perform the analysis with λ XY = 0 and give a maximum value of λ XY that do not jeopardise the multi-component scenarios realised by our chosen benchmark points. We then present another sets of parameter values with a somewhat smaller m h together with large dark Higgs-portal couplings, λ Xφ and λ Y φ . In this case the destructive interference occurs and λ XY may take a larger value. The mixing angle between the SM Higgs and the dark Higgs is chosen to be small, α = 10 −2 , to evade the invisible Higgs decay constraint. Consequently, couplings to the SM Higgs field are not important. Before presenting our results for the multi-component DM scenarios, let us consider theoretical and experimental constraints.

Constraints
We take into account the following theoretical and experimental constraints.

• Unitarity
The unitarity bounds are given by for i, j, l, m = {X I , X R , Y }. The unitarity bound is shown in Fig. 2 in the (Ω y h 2 , g X ) plane.
• Kinetic equilibrium If the 3 → 2 processes are dominant, the kinetic energies of the DM candidates keep piling up. Thus, the kinetic energies need to be released. This can be achieved when the DM particles have a way to talk to the SM thermal bath, so that elastic scattering processes release the kinetic energy. See Refs. [17,63,64] for details. Following Ref. [63], we compute the momentum relaxation rate γ of a DM particle scattered off the relativistic SM thermal bath species, namely the electron. We then require the relaxation rate to be larger than or compatible with the kinetic energy production from the 3 → 2 processes, i.e. γ Hx 2 KD , until the kinetic decoupling time x KD ≡ m Y /T KD . We set the kinetic decoupling time to be similar to the freezeout time, x f . In our study, we use the Z -portal interaction to communicate with the relativistic SM heat bath. Since the relativistic SM degree of freedom is electron at T ∼ O(MeV), we consider the DM-e elastic scattering process. The most dominant Thus we focus on the kinetic energy release of the field Y . From the kinetic equilibrium condition we see that the Z mass is bounded as follows: The kinetic decoupling region is sketched in Fig. 1.
• DM relic density The current DM relic density is given by [65] Ω DM h 2 = 0.1200 ± 0.0012 . (3.11) • Invisible Higgs decay The Higgs invisible decay is bounded [66] as at 90% C.L., where BR h is branching ratio of the Higgs invisible decay. We choose a small enough mixing angle α = 10 −2 to satisfy the bound.
• Direct detection For the DM direct detection bounds, we consider Xenon10 [67], Xenon1T [68,69], and expected SENSEI-100 bound [70]. For all the benchmark scenarios that we consider (Table 1), the DM-electron cross sections are tiny, and thus safe from the direct detection constraints. The DM direct detection bounds are depicted in Fig. 2. The DM-quark scattering bounds are unimportant since the DM masses are less than 1 GeV in our benchmark scenarios.
• DM self-scattering cross section The SIMP-type DM models naturally have a sizeable self-scattering cross section that helps to solve the small-scale problems mentioned in Section 1. We require [21,[81][82][83] 0.1 cm 2 /g < σ self m DM < 10 cm 2 /g . The Bullet Cluster [84,85] imposes a stricter constraint, σ self /m DM < 1 cm 2 /g. See also Ref. [86]. A similar bound is obtained from cosmological simulations with selfinteracting DM on the scales of galaxies and galaxy clusters [87,88].
For the m DM , we take an effective DM mass. For the cases i) and iii) the dominant contribution to σ self comes from the Y self-scattering cross section, as X I and X R are much heavier than Y . We thus choose m DM = m Y . For the case of ii) both X I and Y may contribute; thus we choose m DM = (m Y + m X I )/2. The case iv) is more complicated. In this case X R , X I , and Y all contribute. We thus consider all the contributions with the effective masses. All the contributions are weighted by the fractions of the relic density. In summary, (3.16) In Fig. 2, σ self /m DM = 0.1 cm 2 /g, 1 cm 2 /g, and σ self /m DM > 10 cm 2 /g are shown.
Taking all of the listed constraints into account we now present four benchmark results for the multi-component DM scenarios.

Two-component Scenario
When the mass of X R is sufficiently larger than m X I + m Z , X R cannot be a DM candidate as it decays into X I and Z . Thus we have the two-component scenario. We consider two We have the same plot for ii), iii), and iv) since the same m Y , m Z , and are used. Furthermore, we do not see a clearly visible difference between the left and right plots since the difference in m Y is small. The blue points represent our benchmark cases shown in Table 1. cases The input parameter values are shown in Table 1.
We solve the Boltzmann equations given in Appendix B for the given parameters, and the solutions are shown in Fig. 3. We see that, in the case of m X I m Y , X I freezes out first at x ≈ 5, followed by the Y freeze-out at x ≈ 20. When m X I ≈ m Y , both X I and Y freeze out at almost the same time x ≈ 20. In the case i), we find that X I (Y ) constitutes about 23.7% (76.3%) of the total relic density, while in the case ii), the fractions of X I and Y relics are respectively 19.7% and 80.3%. Thus we identify these scenarios as twocomponent DM scenarios. The self-scattering cross sections are given by 1.49 cm 2 /g and 3.83 cm 2 /g for cases i ) and ii ), respectively. We summarise the results in Table 2. For each case, the unitarity bound (grey), direct detection bound (orange), and large self-scattering cross section bound σ self /m DM > 10 cm 2 /g (red) are shown for different values of the dark gauge coupling g X and the fraction of the Y relic density Ω y h 2 /0.12. The black lines represent DM-electron scattering cross section values, 10 −40 cm 2 (dashed), 10 −41 cm 2 (dot-dashed), and 10 −42 cm 2 (dotted) from top to bottom, respectively. The red lines correspond to two different values of the self-scattering cross section, 0.1 cm 2 /g (dotted) and 1 cm 2 /g (dashed). For details of the used constraints, see the main text in Section 3.2. The blue points represent our benchmark cases shown in Table 1. Here we vary only g X and the fraction of the Y relic density, with the condition (Ω y + Ω X I + Ω X R )h 2 = 0.12. For the two-component DM cases i) and ii), Ω X R = 0 is chosen, while for the three-component DM cases iii) and iv), Ω X R = Ω X I is assumed, except for the self-scattering cross section for which we follow the method outlined in Section 3.2. The other input parameters are fixed (see Table 1). We note that the correct present-day DM relic density constraint is not strictly imposed, except at the benchmark points. Thus, it is important to note that not all the white region is allowed. On the other hand, this approach allows one to comprehensively understand the various constraints and where our benchmark cases lie. The full picture requires an intensive parameter scan by numerically solving the coupled Boltzmann equations which is beyond the scope of the current work.

Three-component Scenario
When the mass gap between X R and X I is small enough, X R may become a good DM candidate. We take the mass splitting to be small enough, m X R − m X I min{2m Y , m Z }, such that the X R decay channels such as X R → X I , Z and X R → X I , Y, Y * are kinematically closed 7 . Thus, together with Y and X I , we have a three-component scenario. We consider two cases: iii) m X R ≈ m X I m Y and iv) m X R ≈ m X I ≈ m Y . The input parameter values are shown in Table 1. We chose a large dark gauge coupling, g X ∼ O(1).
While the decay channels of the field X R are kinematically forbidden, it is still possible that the field X R may be converted into the field X I via the conversion process X R , Y → X I , Y . We numerically solve the full Boltzmann equation and found that the mass gap should be less than ∼10 keV, i.e. m X R − m X I 10 keV to suppress such a process, thereby achieving three-component scenarios.
The solutions of the Boltzmann equations are shown in Fig. 3. Since m X R ≈ m X I we do not see any clear deviations between the X R evolution and X I evolution. Similar to the cases i) and ii), we see that, in the case of m X I m Y , X I freezes out first at x ≈ 5, followed by the Y freeze-out at x ≈ 20. When m X I ≈ m Y , both X I and Y freeze out at almost the same time x ≈ 20. In the case iii), we find that X I and X R (Y ) constitute about 20% each (60%) of the total relic density, while in the case iv), the fractions of X I and X R relics are about 9% each, with 82% of Y relic. Thus we identify these scenarios as three-component DM scenarios. The self-scattering cross sections are given by 1.66 cm 2 /g and 4.04 cm 2 /g for cases iii ) and iv ), respectively. We summarise the results in Table 2.
Let us comment on the parameter λ XY before we conclude. In our benchmark cases i)-iv) we chose λ XY to be zero so that the contact interaction between X and Y gets suppressed, helping realizations of the multi-component DM scenarios in the large-m h limit. For the chosen parameter values of the benchmark cases, we present the upper bound on λ XY to achieve such multi-component scenarios. We first require that Ω Y /Ω DM 90% be the condition for a multi-component DM scenario. For the case i), we obtained the maximum value of λ XY = 1.5 × 10 −4 with λ Y = 5.45, in which case the relic of the Y (X I ) field consists 90% (10%) of the total relic density. For the case ii), the maximum value is given by λ XY = 4.7 × 10 −5 with λ Y = 4.62. The Y (X I ) relic consists 90% (10%) of the total relic density. For the case iii) (case iv)), λ XY = 4.0 × 10 −4 (3.8 × 10 −4 ) is the maximum value, with λ Y = 3.71 (4.95) that gives Ω Y /Ω DM = 90% and Ω X I,R /Ω DM = 5%. Note that the value of λ Y is modified accordingly in order to satisfy the correct relic density On the other hand, in the cases iii) and iv) all the DM candidates freeze out, and thus we have three-component scenarios. We note that there is no visible difference between the X I relic and X R relic in the cases iii) and iv) due to the small mass gap. A large gauge coupling g X ∼ O(1) is chosen for our benchmark points. In this case, once the mass gap becomes larger than ∼10 keV, a significant amount of the relic of X R is converted into the X I relic, mainly through X R , Y → X I , Y , becoming a two-component scenario. In all cases we see that Z follows its equilibrium state.
of Ω DM h 2 = 0.12. The rest of the input parameters are kept to be the same as in Table 1. As such, the value of λ XY needs to be small to realise multi-component DM scenarios. This is not a general statement, however, since the dark Higgs field is assumed to be massive compared to the other dark fields. If the dark Higgs is chosen to be somewhat lighter and/or the couplings λ Xφ and λ Y φ are somewhat larger, the value of λ XY may become larger, due to different signs in terms between the contact interaction and the dark Higgs mediated processes, causing cancellation, i.e. destructive interferences. In Table 3 we present such cases; we choose m Y = 37.5 MeV, m Z = 200 MeV, m h = 10 GeV, g X = 1, λ X = 1, λ Xφ = λ Y φ = 10, and λ Y = 5 for all the four cases.
Case Ω X R /Ω DM Ω X I /Ω DM Ω y /Ω DM σ self /m DM (cm 2 /g)  Table 2. The fractions of relic density for each DM candidate field and the self-scattering cross sections are shown for the four benchmark cases. The input parameters are summarised in Table  1. The total DM relic density is Ω DM h 2 = 0.12. As expected from the mass gap between X I and X R , cases i) and ii) give rise to two-component DM scenarios, while three-component DM scenarios are realised in the cases iii) and iv). The self-scattering cross sections are somewhat larger than 1 cm 2 /g imposed by the Bullet Cluster constraint [84][85][86] (see also Refs. [87,88] where the similar bound is obtained from cosmological simulations with self-interacting DM), but well within the bound, 10 cm 2 /g.  Table 3. Input parameters with non-zero λ XY and outcomes of DM relics and the self-scattering cross sections. The rest of the input parameters are chosen as follows: m Y = 37.5 MeV, m Z = 200 MeV, m h = 10 GeV, g X = 1, λ X = 0.1, λ Xφ = λ Y φ = 10, and λ Y = 5. All the cases are free from the constraints listed in Section 3.2. We see that the multi-component DM scenarios are still realised with sizeable λ XY values. This is due to the destructive interference between the X-Y contact interaction and the dark Higgs mediated processes.

Conclusions
In this paper, we constructed multi-component dark matter scenarios with U (1) X dark gauge symmetry broken into Z 2 × Z 3 . Two independent types of dark matter fields, namely X and Y , distinguished by their different charges under Z 2 ×Z 3 , emerge as separately stable dark matter candidates. The relic densities of different dark matter fields are governed by interactions between the dark fields within the dark sector through number-changing 2 → 2 and 3 → 2 processes. Unlike many other SIMP models in the literature, X and Y are completely independent of each other, with different masses and spins in principle. Due to this reason we call our model a genuine multi-component SIMP dark matter scenario. A nonzero vacuum expectation value of the dark Higgs field gives mass splitting to the real and imaginary parts of the X field. When the mass gap between the real part X R and imaginary part X I is large, only the lighter field X I becomes a stable dark matter candidate which, together with Y , gives rise to a two-component dark matter scenario. On the other hand, once the mass gap is small enough, both X R and X I may contribute to the dark matter relic density; we thus have a three-component scenario. In our analysis, the dark gauge coupling g X is chosen to be O(1), and we found that the mass gap should be less than ∼10 keV, i.e. m X R − m X I 10 keV, in order to achieve three-component scenarios.
The unbroken Z 3 symmetry ensures the stability of the Y field, which becomes another dark matter candidate of the SIMP-type. Due to the interactions between X I,R and Y , both X I and X R may contribute to Y 's 3 → 2 annihilation processes through the Z and the dark Higgs. These kinds of processes open a new class of SIMP models.
We presented four choices of the parameter set for the multi-component dark matter scenarios, two for a two-component dark matter scenario, i) m X R m X I m Y and ii) m X R m X I ≈ m Y , and two for a three-component dark matter scenario, iii) m X R ≈ m X I m Y and iv) m X R ≈ m X I ≈ m Y , by solving the full coupled Boltzmann equations and considering both theoretical and experimental constraints. Figures 1 and 2 show a set of constraints. The parameter choices are given in Table 1 and the solutions of the Boltzmann equations are shown in Fig. 3. For each case we computed fractions of the relic density. The results, together with values of the self-scattering cross section, are summarised in Table 2. We found that, in the large-m h limit, the λ XY coupling, which dilutes the relic density of the X field via the X-Y contact interaction, needs to take a small value to realise a multi-component dark matter scenario. On the other hand, the λ XY coupling may take a larger value once the dark Higgs field becomes lighter and the dark Higgs portal couplings, λ Xφ and λ Y φ , become larger. In this case destructive interferences occur and the dilution of the X relic is suppressed even with a larger value of λ XY . In Table 3 we present four parameter sets of this case.
Our findings of the multi-component SIMP-type dark matter scenarios are new and interesting. One of the most distinctive features of our model is that the different dark matter fields X and Y may have very different mass scales, while significantly contributing to the total dark matter relic density today, as demonstrated in the cases i) and iii). In the current work a scalar Z 2 -charged field is considered. It is certainly possible to have a fermionic Z 2 -charged field 8 as studied e.g. in Refs. [54,90,91]. Furthermore, we investigated four characteristic show-cases for the multi-component dark matter scenario. In order to fully understand the model and its phenomenology, an intensive parameter scan is better to be performed. We shall tackle them in future works.
8 Note that the Z3-charged dark matter cannot be fermionic.
The evolution of the number density of particle 0 is then governed by the following Boltzmann equation: and we assumed 1 ± f i ≈ 1 and CP conservation, i.e. |M| 2 0→1,2,3 = |M| 2 1,2,3→0 . Now, using f i = (n i /n eq i )f eq i and f eq 0 = f eq 1 f eq 2 f eq 3 , we obtain In terms of the yields, the Boltzmann equation can be rewritten as follows: where x ≡ m/T with m being a mass variable. For the decay process X R → X I + Y + Y * , we have The decay rate is given by

B The Full Boltzmann Equations
In this appendix, we present the full Boltzmann equations for Y , X R , X I , and Z . The Boltzmann equation for Y is, using Y Y = Y Y * = Y y /2, given by The Boltzmann equations for X I and X R are Finally, the Boltzmann equation for Z is given by