Multi-component dark matter from a hidden gauged SU(3)

We study Dark Matter (DM) phenomenology with multiple DM species consisting of both scalar and vector DM particles. More specifically, we study the Hidden Gauged SU(3) model of Arcadi {\it et al}. Before proceeding to the Hidden Gauged SU(3) model, we study the relic abundances of simplified multi-species DM scenarios to gain some insights when multiple species and interactions are included. In the Hidden Gauged SU(3) model, because of the large parameter space, we restrict ourselves to three representative benchmark points, each with multiple DM species. The relic densities for the benchmark points were found using a program developed to solve the coupled Boltzmann equations for an arbitrary number of interacting DM species with two particles in the final state. For each case, we varied the mass of the DM particles and then found the value of the dark SU(3) gauge coupling that gives the correct relic density. We found that in some regions of the parameter space, the DM would be difficult to observe in direct detection experiments while easier to observe in indirect detection experiments and vice versa, so that complementary measurements could help pinpoint the details of the Hidden Gauged SU(3) model. Important to this, is that even for moderate changes in input parameter values, the relative relic density of each species can change significantly resulting in large changes in the observability of multi-species DM by direct or indirect detection.


I. INTRODUCTION
The nature of Dark Matter (DM) is one of the biggest puzzles in particle physics. The evidence for DM in the Universe comes from a wide variety of astrophysical and cosmological observations and there has been considerable theoretical effort to understand its nature (for recent reviews, see for example Ref. [1][2][3][4][5]). Up until now, much of the work has focused on simple DM sectors, typically with only one or possibly two DM candidates. These scenarios are increasingly constrained by experimental measurements. Therefore, workers in the field are exploring ever more complex DM scenarios including models with multiple DM species (see for example, ). In this context, there is a growing body of work studying details of "freezeout" for multiple interacting DM species, e.g. [9, 10, 13, 14, 19, 25-27, 36, 37, 39, 40]. To expand on this theme, we developed a computer program capable of calculating the relic densities of an arbitrary number of interacting DM species and, in this work, use it to explore the behavior for several scenarios and the implications for DM direct and indirect detection experiments.
As a specific example of a model with multiple DM species, we re-examine the Hidden Gauged SU(3) model of Arcadi et al. [31]. This model has the potential to have numerous DM species depending on the choice of parameter values from the breaking of the new SU(3). The paper of Arcadi et al. [31] chose specific values of the parameters that reduced the multi-species DM problem to effectively two DM species to enable them to implement the model in the computer program MicrOMEGAs * apoulin@physics.carleton.ca † godfrey@physics.carleton.ca [42]. It turns out that different values of the parameters can lead to much more complex and interesting scenarios with multiple DM species to explore. Thus, a main focus of this paper is to explore different values of the parameters of the Hidden Gauged SU(3) model that lead to a more complex DM particle spectrum. We found that some scenarios with multiple DM species give the correct relic density, but because the dominant components have a small direct detection cross section, the DM will be difficult to observe in direct detection experiments. However, in some cases when the DM is not observable by direct detection experiments, it could be observable in indirect detection measurements and vice versa. These points were also made by Arcadi et al. [31] and by many others for other DM scenarios. A small sampling of examples is given by Ref. [1-5, 15, 20, 35]. The Hidden Gauged SU(3) model has a very rich phenomenology with details of the phenomenology very much dependent on the choice of parameter values. Some of our results convey a warning that approximations can lead to erroneous conclusions. For example, in a scenario with multiple DM species, a moderate change in the input parameter values can lead to major differences in the relative relic densities of the stable DM species. This would have a major impact on whether DM could be observed by direct or indirect detection experiments. An important lesson from our study is that models with multiple species of DM particles can lead to a very rich phenomenology and that one needs to be careful when studying relic densities in multi-species DM models.
The parameter space of the Hidden Gauged SU(3) model is large. Integrating the coupled Boltzmann equations to obtain the relic density for multiple DM species is computationally intensive so that it is impractical to perform a complete scan of the parameter space. Instead, we choose a small number of representative benchmark points in the parameter space and examine their phenomenology. For two of the benchmark points, there are 4 stable DM species, while in the third, there are 7. These cases will be examined in detail in Section V.
Section II gives the details of the coupled Boltzmann equations we use to calculate the relic abundance for multi-species DM scenarios [40]. This is the core ingredient of the paper needed to calculate the subsequent results for the Hidden Gauged SU(3) model. Before proceeding to the Hidden Gauged SU(3) model, we explore two simple multi-species DM scenarios in Section III. In that section, we assume three or four scalar DM species with a range of masses and assumptions for the interactions between them. We do so to gain some insights before proceeding to the more complicated Hidden Gauged SU(3) model. In Section IV, we give the details of the Hidden Gauged SU(3) model. In that section, we start by writing down the Lagrangian, then compute the mass eigenstates for the scalar and vector sectors in terms of the Lagrangian parameters, the vacuum expectation values, and the new gauge coupling. We then write down expressions for the theoretical constraints on these parameters. In Section V, we explore the phenomenology of three benchmark points for the Hidden Gauged SU(3) model that satisfy the theoretical constraints and use the relic density to constrain the parameters of the model which are then used to calculate direct detection and indirect detection cross sections. Finally, in Section VI, we summarize our results.

II. THE BOLTZMANN EQUATIONS FOR MULTI-SPECIES DM
Our purpose is to study models of DM that include multiple DM species. The first step in such a study is to calculate the relic density for a DM sector with an arbitrary number of interacting DM species. To calculate the DM relic density, we assume that the dark sector was in thermal and chemical equilibrium at some early time which allows us to use the coupled Boltzmann equations in the form written down by Dienes, Huang, and Thomas [40] which is a generalization of the Boltzmann equation given in Kolb and Turner [43] (see also Ref. [44]). These themes have been explored previously by Edsjo and Gondolo [39] and others, for example [9,19,25,26,30,36,41].
In what follows, we simply reproduce the coupled Boltzmann equations as given by Dienes, Huang, and Thomas [40] which we will use to calculate each of the individual relic densities and refer the interested reader to their work for details [40]. We only include the details needed to reproduce our results including some of the equations and definitions that we used from Kolb and Turner [43] and refer the interested reader to that monograph for more complete details.
We wish to calculate the relic density of each of the DM species, φ i , which is given by where m i , n i , and ρ i = m i n i are the mass, number density, and energy density of particle φ i , ρ C is the critical density, G is Newton's constant, and H is the Hubble parameter. For species φ i in kinetic equilibrium, the number density is given by where g i is the number of internal degrees of freedom and the phase space distribution, f i ( p i ), is given by where T is the temperature of the photon bath, p i , E i , and µ i are the three-momentum, energy, and chemical potential associated with species φ i , and the plus and minus signs are for fermionic and bosonic species respectively. We are working in the comoving frame in which v = 0 for all φ i . In a frame in which v = 0, the f i ( v) only depend on the magnitude of v and not its direction. In general, the Boltzmann equations can include an arbitrary number of species participating in reactions and decays [40]. We will restrict ourselves to simplified cases where: the Boltzmann equations for all relevant species can be expressed in terms of number densities, n i , rather than the phase space distributions, f i ; we assume CP is conserved so that the matrix elements satisfy |M(ij → kl)| 2 = |M(kl → ij)| 2 for cross sections and |M(i → jk)| 2 = |M(jk → i)| 2 for decays; and that all visible sector particles and all dark sector particles are in thermal equilibrium due to rapid interactions among themselves until long after the lightest DM species have frozen out. To minimize the complexity of the problem, we restrict ourselves to 1 → 2 and 2 → 2 type interactions. With these simplifications, the Boltzmann equations become [40]: , where in Eq. 4, the indices i, j, k, l run over the darksector particles and a, b run over visible sector particles. The first term on the right-hand side is the dilution due to the expansion of the Universe and H is the Hubble parameter. σ xy→zw v , which will be given below, represents the thermally averaged annihilation cross section of two initial state particles, xy, to two final state particles, zw, where xyzw represent either visible or dark sector particles. Likewise, Γ x→yz represents the partial decay width of particle x into particles yz. Eq. 4 assumes that all particles are in thermal equilibrium. In Eq. 4, the thermally averaged cross sections are given by: where g i is the number of internal degrees of freedom of particle i, σ xy→zw is the cross section for the process is the magnitude of the relative velocity of φ x and φ y in the comoving frame. In the absence of a degenerate Fermi species or a Bose condensate, the phase space distributions of the DM species reduce to the classical Maxwell-Boltzmann distribution in the non-relativistic limit [43]. For sufficiently rapid interactions, the particles approach kinetic and chemical equilibrium depending on the nature of the interaction. This latter situation places constraints on the chemical potentials of the species involved which depend on the reaction. For example, the reaction φ i + φ j ↔ φ k + φ l would imply the constraint µ i + µ j = µ k + µ l . For the visible particles, the chemical potential µ is small or zero which results in the chemical potential of the DM species being small or zero in the absence of an asymmetry. We will therefore treat all equilibrium chemical potentials as zero. Under these conditions we approximate f i (v i ) by the equilibrium phase space distribution: where v i is the speed of the particle. With these approximations, the equilibrium number density, n eq i , is given by where K 2 (x) is the modified Bessel function of the second kind. The standard approach for solving Eq. 4 is to scale out the effect of the expansion of the Universe by considering the evolution of the number density in a comoving volume. This is done by using the entropy density, s, as a fiducial quantity by defining [43]: where the entropy density in the comoving volume in any cosmological epoch is defined by T is the photon temperature and ρ and p are the total energy density and the pressure, respectively, expressed in terms of the photon temperature T [43]. g * S represents the number of interacting degrees of freedom in the thermal bath which can be found by rewriting Eq. 11 as g * S = 45 2π 2 T 4 (ρ + p).
Y i is the actual number per comoving volume and Y eq i is the equilibrium number per comoving volume. In addition, we rescale the Hubble constant H(T ) = H(m)/x 2 with x = m/T where m is any convenient mass scale [43] which we will take to be the mass of the heaviest stable DM particle. In the radiation dominated epoch, H(m) is given by [43] H(m) = 2π 3/2 √ 45 where m P l is the Planck mass and g * represents the total number of effectively massless degrees of freedom and is given by where ρ R is the total energy density of all species in equilibrium. Expressions for ρ R , g * S , and g * can be found in Kolb and Turner [43]. In computing g * and g * S we used the Fermi-Dirac or Bose-Einstein distributions as appropriate in the sum. For the region around the QCD phase transition, the correct degrees of freedom becomes ambiguous and requires lattice QCD to compute. We use the values found in [45] for this region. With these substitutions, Eq. 4 is recast as:

III. SIMPLE MULTI-COMPONENT DM SCENARIOS
Before proceeding to the Hidden Gauged SU(3) DM model, we start by exploring the relic abundance of two simplified multi-species DM scenarios. Models with multiple species and different types of interactions can be complicated. Deconstructing how the different terms in Eq. 15 affect the number densities can give us insights into the coupled Boltzmann equations and help us understand which terms are important and which ones can be neglected. Both of the scenarios initially consist of multiple DM species freezing out by only self-annihilating into Standard Model (SM) particles. We then add additional interactions between the DM species and observe the effects. Finally, we include all interactions that would be present due to crossing symmetry to obtain the final result. We assume that the thermally averaged cross sections for all interactions are s-wave, and thus independent of speed, and use the approximation σv = 0.1 pb for all the interactions. In a particular model, these cross sec-tions would not necessarily be equal, but this assumption is reasonable for the purpose of exploring the behavior of Eq. 15 for simplified scenarios. In all cases, we assume that the SM particles are in thermal equilibrium with the photon bath and treat them as relativistic with the appropriate degrees of freedom. Plots of the DM number per comoving volume versus x = m/T are shown in Figures 1 and 2 where we take m to be the mass of the heaviest stable DM particle.
In all cases, the solid curves represent the equilibrium number densities, the dashed lines represent the scenarios with only self-annihilation to SM particles which we label as the "base", the lines with a combination of dashes and dots include the base interactions and a single additional interaction, ignoring the contributions from crossing symmetry, and finally the dotted lines show the results when all the additional interactions are included. We label the species as φ i with the index i increasing with increasing species mass.
In Fig. 1, we consider three DM species, φ 1 (green), φ 2 (blue), and φ 3 (red), with masses of 100 GeV, 120 FIG. 1. DM number densities as a function of x = m/T for three DM species. The solid curves represent the equilibrium density, the dashed curves represent the scenario with only self-annihilation to SM particles which we label as the "base", the curves with both dots and dashes represent scenarios with the base interactions plus a single additional interaction as labeled in the plot, and the dotted curve represent the result including all the additional interactions of the previous curves. The green, blue, and red curves represent the number density of φ1, φ2, and φ3 with masses of 100 GeV, 120 GeV, and 150 GeV, respectively. Note that for the green and blue curves, the various curves sit very closely on top of each other.
FIG. 2. DM number densities as a function of x = m/T for four DM species. The line labelling is as in Fig. 1. The green, blue, red, and purple curves represent the number density of φ1, φ2, φ3, and φ4 with masses of 100 GeV, 120 GeV, 150 GeV, and 180 GeV, respectively. Note that for the green, blue, and red curves, the various curves sit very closely on top of each other.
GeV, and 150 GeV, respectively. We will start with the most general application of Eq. 15 that is applicable to this case and examine each of the contributions in turn. For this case, Eq. 15 reduces to the following set of three coupled equations: However it turns out that for our choice of masses the this term also has a neglible contribution compared to the first line in Eq. 20. The net result is that the interactions added to the base case do not significantly alter the φ 2 relic abundance.
Finally, we consider the temperature where φ 3 freezes out; Y 3 = cY eq 3 . In this case, Y 1 ≈ Y eq 1 and Y 2 ≈ Y eq 2 . This gives: In this case the largest reaction rate is for φ 1 + φ 3 ↔ φ 1 +φ 2 , corresponding to the third line of Eq. 21, followed by φ 2 + φ 3 ↔ φ 1 + φ 1 , corresponding to the second line of Eq. 21, because of the hierarchy of the number densities, . This is what we see in Fig. 1 where the curves with the additional interactions stay in thermal equilibrium longer, resulting in a smaller relic density for φ 3 . So for the case of annihilation of φ 3 with φ 1 , there is a small amount of φ 3 which sees a large amount of φ 1 so that the interaction φ 1 + φ 3 ↔ φ 1 + φ 2 depletes φ 3 significantly but barely affects the φ 1 density which is what is seen in Fig. 1. Because φ 2 has a smaller relic density than φ 1 , it doesn't deplete φ 3 as much, which is also observed in Fig. 1.
In Fig. 2 we repeat the exercise for four DM species, φ 1 (green), φ 2 (blue), φ 3 (red), and φ 4 (purple), with masses of 100 GeV, 120 GeV, 150 GeV, and 180 GeV, respectively. For this case Eq. 15 leads to a set of four coupled equations. φ 1 , φ 2 , and φ 3 behave similarly to the previous case so to avoid repetition we will focus on the behavior of the φ 4 number density. The differential equation for Y 4 is given by: To obtain some insights into the behavior of the φ 4 number density we start by considering the temperature where φ 4 freezes out, Y 4 = cY eq 4 with c 1, and assume that in this region These expressions can explain why the curves, but fails to explain other details of this plot. As pointed out above, the number density of a DM species depends on its mass so that there is less of φ 3 to annihilate with than φ 1 and φ 2 so that φ 4 is depleted less by its iteractions with φ 3 than with the other DM species. Next, we consider the behavior of the φ 4 number density for the curves when the φ 4 + φ 2 ↔ φ 1 + φ 3 and φ 4 +φ 1 ↔ φ 2 +φ 3 interactions are present in the temperature region around φ 3 and φ 4 freezeout where Y 4 = c 1 Y eq 4 and Y 3 = c 2 Y eq 3 with c 1 c 2 1: In this temperature region the factors, (1 − c 2 c −1 1 ), in the third and fourth equations in Eq. 24, are small so the φ 4 number density is dominated by the first equation in Eq. 24. The 2nd equation is not relevant here as it applies to the φ 4 + φ 3 ↔ φ 1 + φ 2 case which we are not considering here. The result is that the curves for values of Y 4 shortly after φ 4 departs from the equilibrium curve.
Lastly, we examine the region where φ 2 freezes out corresponding to the temperature where Y 2 = cY eq 2 with c 1 as before. We assume that at the temperature range we are considering. As the temperature decreases, c, which is a measure of the difference between Y and Y eq , increases, so that will become negative and the reverse reaction will dominate. At the same time, remains positive and continues to deplete φ 4 . This is why in Fig. 2 The dotted line shows the result when all processes are included. It can be seen that the reverse reaction is important up until φ 1 starts to freeze out, at which point the φ 1 + φ 4 ↔ φ 2 + φ 3 interaction becomes dominant and starts depleting φ 4 .
These simple examples show that even for a simple multi-species DM sector the behavior of the relic abundance can be quite complicated and the resulting relic abundances can vary substantially depending on the details of the underlying model. It is not always a priori obvious how the details of the model will work out and one needs to understand how each of the interactions contribute to the final result. All things being equal, heavier DM species will be more affected by their interactions with lighter DM species than vice versa. The interactions between multi-species SM can lead to experimental consequences. Say, for example, that a model has many DM species but only one interacts with a portal to the visible universe. One could imagine a scenario where the model could account for the measured DM relic density, but the component that interacts with the visible universe has a very small relic density, making it unlikely that it could be observed in direct or indirect detection experiments.
In the next sections, we will explore in detail a specific multi-species DM model, the Hidden Gauged SU(3) model. Because the parameter space is large, we will restrict ourselves to several representative benchmark points with interesting phenomenology.

IV. THE HIDDEN GAUGED SU(3) MODEL
In the remainder of this paper, we study the phenomenology of the Hidden Gauged SU(3) model of Arcadi et al. [31] which consists of spin-1 and spin-0 states. In their paper, they examined two representative limiting cases that made their numerical analysis tractable. We extend their studies to consider different points in the parameter space that leads to more complex DM scenarios. We mention that others have also studied hidden SU(3) models [26,28,33].

A. Details of the Model
For completeness and clarity, we start by reproducing the details of the Hidden Gauged SU(3) model of Arcadi et al. [31]. The model consists of a gauged SU(3) which is fully broken by two complex scalar triplets, Φ 1 and Φ 2 , so that all the new gauge bosons acquire a mass. These new scalars are not charged under the SM gauge groups and can only interact with the SM through the SM Higgs doublet. To simplify the Lagrangian and insure additional stable states, a Z 2 symmetry is also imposed such that which has the effect of only including even powers of the scalar triplets in the Lagrangian. Following Ref. [31], and imposing this additional Z 2 symmetry, we write the Lagrangian as: where Here, G a µν = ∂ µ A a ν − ∂ ν A a µ +gf abc A b µ A c ν is the field strength tensor for the hidden SU(3) whereg is the gauge coupling and f abc are the SU(3) structure constants. The covariant derivative is given by D µ = ∂ µ − igA a µ t a where t a are the Gell-Mann matrices (see for example Ref. [46] where the t a = 1 2 λ a ).
In general, the scalar doublet and triplets will have the form: where the v i are the vacuum expectation values (vevs), the superscript indicates to which of the Φ i the field belongs, and the subscripts indicate the position in the triplet and whether it is the real or imaginary component. When all three new vevs are taken to be non-zero, we have some freedom in defining the Goldstone bosons of the hidden SU(3). One convenient way to define them is which leaves the physical states We will see in section IV B 3 that in order to have a consistent vacuum state, we must have v 1 v 3 (λ 4 +λ 5 ) = 0. Since we require v 1 = 0 to fully break the SU(3), we will study the cases where we take v 3 = 0 and allow λ 4 = −λ 5 . For completeness, we give the details of the λ 4 = −λ 5 case in the appendix for the interested reader who may wish to study this possibility.
In the v 3 = 0 case, the resulting Goldstone bosons are given by where we define s β = sin β = . This leads to the convenient definition for the physical states given by and the resulting form of the triplets in the unitary gauge is given by Spontaneous symmetry breaking in the SM happens in the usual way so that the Higgs doublet in the unitary gauge is given by: where v is the SM vev.
In the unitary gauge, we write down the scalar mass matrix in the form L ⊃ 1 2 Φ T m 2 scalar Φ where Φ T = (h, ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 ) are the scalars in the model and where We see that because the second row of Φ 1 and the third row of Φ 2 do not have a vev, ϕ 3 and ϕ 4 cannot mix with any of the other scalars. In addition, because we do not have CP-violation in the scalar sector, ϕ 3 and ϕ 4 cannot mix with each other so both are mass eigenstates. To differentiate the mass states from the gauge states, we relabel the mass eigenstates as ϕ 3 = H and ϕ 4 = χ. We label the remaining mass eigenstates as h 1 , h 2 , and h 3 and parameterize the mixing between them using three angles θ 1 , θ 2 , and θ 3 as follows: Here, m hi is the mass of h i , s i = sin θ i , c i = cos θ i , for i = 1, 2, 3. The advantage of this parameterization is that, in the limit of small θ 2 , we can define an effective mixing angle θ = θ 1 + θ 3 which mixes h and ϕ 2 . This is convenient since the mixing between the SM Higgs and the dark scalars is constrained to be small to be consistent with precision electroweak measurements.
Ref. [31,47] finds that sin θ 0.3 − 0.4 depending on the mass of the second scalar. For all our benchmark points, we will take sin θ = 0.1.
In the vector boson sector, there is only mixing between A 3 and A 8 through the following mass matrix We define the mass eigenstates as: where c α = cos α, s α = sin α, and We note that α ∈ 0, π 3 and that t α = tan α > 0. We label all the mass eigenstates with a prime to distinguish between the mass and gauge eigenstates. The resulting vector boson masses are given by We can also write m A 3 and m A 8 explicitly in terms of the vevs as From this form, we can see that A 3 µ is the lightest vector boson while A 8 µ is the heaviest, or one of the heaviest. When v 1 = v 2 , A 8 µ becomes degenerate with A 6 µ and A 7 µ . This model has a global Z 2 × Z 2 × Z 2 symmetry remaining from the breaking of the SU(3) with Z 2 being an additional Z 2 to the model studied in Ref. [31].
The Z 2 symmetry is a subgroup of two remaining U(1) symmetries, one for the vectors and one for the scalars. The U(1) for the vectors is characterized by the matrix Under the transformation A a µ t a → U A a µ t a U † , setting δ 1 = π takes A 1,2,4,5 µ → −A 1,2,4,5 µ and does not change the other vector fields. The U(1) for the scalars is described by taking Φ 1,2 → e iδ2 Φ 1,2 . The Z 2 subgroup corresponds to choosing δ 1 = π and δ 2 = π. The Z 2 symmetry is associated with complex conjugation which is an outer automorphism of SU(3). The Z 2 symmetry is one that does not appear in Ref. [31]. This symmetry is a result of having added the extra Z 2 in equation 26. This added symmetry results in the scalar potential being invariant under rephrasing φ 3 and φ 4 . Independently, a U(1) subgroup of SU(3) characterized by the matrix leaves the the field strength tensor invariant. In order for the scalar kinetic terms to be invariant, we must make a specific choice for the rephasing of φ 3 and φ 4 , as well as for δ 3 . This results in the following transformations: We summarize the charges of the physical states under the Z 2 × Z 2 × Z 2 symmetries in Table I. A minus sign in the table represents multiplying the field by −1 under the symmetry.
From these charges, we can identify the expected DM species in this model. From the charges in Table I, we expect a minimum of 3 stable DM species and a maximum of 7 DM species. Because the mass degeneracies, m A 1 = m A 2 and m A 4 = m A 5 , put restrictions on which decays are kinematically allowed, the actual minimum number of stable DM species is 4.

B. Theoretical Constraints on Lagrangian Parameters
The Hidden Gauged SU(3) model has numerous parameters. Much of the parameter space can be excluded by theoretical consistency of the scalar potential and by experimental measurements. In the following subsections, we lay out the details of the theoretical constraints.

Unitarity
The scalar couplings in the potential can be bounded by perturbative unitarity of the 2 → 2 scalar field and vector boson scattering amplitudes.
The partial wave amplitudes, a J , are related to the matrix element of the process, M, by: where J is the (orbital) angular momentum and P J (cos θ) are the Legendre polynomials. Perturbative unitarity requires that the zeroth partial wave amplitude, a 0 , satisfies |a 0 | ≤ 1 or |Re a 0 | ≤ 1 2 . Because the 2 → 2 scalar field scattering amplitudes are real at tree level, we adopt the second, more stringent constraint. We will use this to constrain the magnitudes of the scalar quartic couplings.
We work in the high energy limit where only four-point diagrams contribute to 2 → 2 scalar scattering since all diagrams involving scalar propagators are suppressed by the square of the collision energy. Thus, the dimensionful couplings m 2 H , m 11 , and m 22 are not constrained directly by perturbative unitarity. In the high energy limit, it is  Under these conditions, only the zeroth partial wave amplitude contributes to M resulting in the constraint |Re a 0 | < 1 2 being equivalent to |M| < 8π. To obtain a minimal list of constraints, we must ensure that all the eigenvalues of the coupled-channel scattering matrix M, which includes each possible combination of two scalar fields in the initial and final states, satisfy this constraint. We include a symmetry factor of 1/ √ 2 for each pair of identical particles in the initial or final states.
We use the fields in the form used in equation 33. This leads to a large 120 × 120 dimension scattering matrix. However, the part of the scalar potential that leads to the quartic interactions obeys nine Z 2 discrete symmetries which bring the matrix M into block diagonal form which greatly simplifies the problem. In a 2 → 2 scattering process, the initial and final states must have the same charges under these Z 2 symmetries. The charges we assign for the scalar fields that leave the potential invariant under these symmetries are given in Table II. We can understand the origins of these charges by examining the various terms in the potential. The first four symmetries arise because the doublet, H, only appears as ) so changing the sign of any one of the fields does not change the term. The fifth through seventh symmetries can be understood by looking at the quadratic terms involving the triplets which appear in the potential, namely |Φ 1 | 2 , |Φ 2 | 2 , and Φ † 1 Φ 2 . From this we see that there are never any terms which multiply fields from different scalar triplet rows. The eighth symmetry is related to the extra Z 2 symmetry given in Eq. 26 that we imposed on the potential. Finally, the last symmetry comes from taking the complex conjugate of the potential. The resulting constraints are given by: where This last polynomial comes from finding the eigenvalues of the 3 by 3 matrix: The bounds from Eq. 70 can be translated into the following bounds using the same technique used in Ref. [48]:

Scalar Potential Bounded from below
We next study the constraints on the scalar couplings by requiring that the potential be bounded from below. The constraints that must be satisfied at tree level for the scalar potential to be bounded from below can be determined by only considering the quartic terms of the potential because these dominate at large field values. We follow the approach of Ref. [48] and [49]. First, we make the following definitions: The parameters ζ and ω are bounded by Making these substitutions, we can write the quartic part of the potential as and The fraction in Eq. 80 is always positive and grows with the overall field excursion r. The x T Ay term in Eq. 80 can be positive or negative; we require it to be positive to ensure that the potential is bounded from below. This term can be expressed as a biquadratic in tan γ 2 with coefficients being other biquadratics in tan γ 1 . A biquadratic of the form a + bz 2 + cz 4 will be positive for all values of z if the following conditions are satisfied: This provides us with the following constraints This last bound can be split up into the different cases for λ 4 > 0, λ 5 < 0, −λ 4 − λ 5 for λ 4 < 0, λ 5 < 0. (90)

Avoiding alternative minima
In this section, we discuss the structure of the vacuum of the model to confirm that the choice of vevs corresponds to the absolute minimum of the potential. First, we write the potential in terms of the vevs: III. Input parameters for the benchmark points. mχ, mH, and m A 1 were allowed to vary for benchmark point A, B, and C, respectively, along with the gauge coupling which was varied to obtain the correct relic density. In all cases, we take the mixing in the scalar sector to only be between h1 and h3 with a mixing angle of sin θ = 0.1. We also set m h 1 equal to the SM Higgs mass.
Minimizing the potential leads to the following constraints: Note that v 2 ∂V ∂v3 = v 3 ∂V ∂v2 + (λ 4 + λ 5 )v 2 1 v 3 . From this, we conclude that to have a consistent vacuum configuration requires that (λ 4 + λ 5 )v 2 1 v 3 = 0. Since we need v 1 = 0 to fully break the hidden SU(3) and setting λ 4 = −λ 5 is an unmotivated tuned choice, we will assume that v 3 = 0. Given a collection of input parameters, we wish to ensure that we are in the deepest minima. From the inputs, we can obtain values for the scalar couplings, the vevs, and g, but not m 2 H , m 2 11 , and m 2 22 . We use Eq. 92 to 95 to solve for these parameters which ensures that our choice of parameters is one of the extrema of the potential. We then find all possible values for (v, v 1 , v 2 , v 3 ) which satisfy Eq. 92 to 95 and check to see which values lead to the smallest value of the potential. If the values of the vevs with which we started do not lead to the deepest minima, we must reject that point.

V. RESULTS AND DISCUSSION FOR THE HIDDEN GAUGED SU(3) MODEL
The Hidden Gauged SU(3) model has too many parameters to make a full scan of the parameter space practical. Therefore, we will explore the implications of three rep-resentative benchmark points. The parameter values for benchmark point A were chosen to explore the region where one of the DM species could annihilate to SM particles through a Higgs resonance. Benchmark point B was chosen to explore a region where no such resonant effects would occur. Finally, benchmark point C was chosen to explore the region where v 1 and v 2 have similar values and where there are many more stable DM particles.
We will vary the mass of one of the DM species and find the gauge couplingg which results in the correct relic density. We vary the mass of χ, H, and A 1 for benchmark point A, B, and C, respectively. Some of the hidden sector particles are stable because symmetry respecting decays are kinematically forbidden. Thus, benchmark points A and B have 4 stable DM species, while benchmark point C has 7 DM species. Our purpose is to study some of the features that arise from having many stable DM species. The benchmark points we consider are defined in Table III. For all the benchmark points we take θ 2 = 0 and combine θ 1 and θ 3 into a single angle θ, which we set to sin θ = 0.1 which is consistent with experimental constraints coming from the SM Higgs couplings [31,47]. We also ignore any unstable dark sector species by assuming that they would have decayed out long before freezeout occurred.

A. DM Relic Density from the Hidden Gauged SU(3) Model
When varying the mass of one of the DM species for a benchmark point, we used a binary search to find a value of the gauge coupling which would result in the relic density lying in the measured range of Ω CDM h 2 = 0.1186 ± 0.002 [45] . Figure 3 shows the resulting values for the coupling as a function of the DM mass for each benchmark point. If the required value of the coupling fell outside the unitarity bound of √ 4π, we did not consider it. This was the case for benchmark point A for a range of m χ masses above the threshold for annihilation via the Higgs boson resonance as can be seen by the lack of points between half the Higgs mass and about 90 GeV.
Another insightful variable is the ratio of the relic density of vector particles to the total relic density. This is defined by This variable shows how the composition of the DM changes as we vary parameters such as the gauge coupling. It is very useful when looking at direct and indirect detection since the scalar and vector particles interact very differently with detectors. Figure 3 shows how f A varies for each benchmark points.

B. DM Direct Detection
We now examine the constraints obtained from direct detection measurements using the latest results from the XENON1T experiment [50,51]. Because there are multiple DM species that can each interact differently with the detector, we cannot simply calculate a cross section and compare it with the reported experimental bound. Instead, we must find the predicted number of events that the experiment would have seen, and compare that to what the experiment observed. The theoretical rate of events is given by where f i = Ωi Ωtot and [15] dR Here, E R is the nuclear recoil energy in the detector, σ 0 iN is the DM nucleus cross section at zero momentum transfer, ρ loc i is the local DM density for the ith species, µ iN = m i m N /(m i + m N ) is the reduced mass of the nucleus DM system, F i (E R ) is a nuclear form factor, and I i (E R ) is the mean of the inverse of the speed, v −1 i , of the ith species in the DM halo for a given E R . Explicit formulas for I i (E R ) and F i (E R ) can be found in Ref. [15]. Once the differential rate is calculated, we integrate it over the energy range appropriate for the experiment and multiply the result by an efficiency factor which could in principle be different for each species. In the case for XENON1T, the energy range for single nuclear recoil events is 4.9 to 40.9 keV, and we assume a detection efficiency of 89%, which for simplicity we took to be the same for all DM species and recoil energies.
For the parameter values to be allowed, the calculated theoretical rate should fall below the experimental limit. At a confidence level of 1 − β, a background rate ν b , and a total of n obs observed events, we can find the upper limit on the signal rate ν up s by numerically solving: For XENON1T [50,51], we use ν b = 735 and n obs = 739.
To compare the experimental rate to the theoretical prediction, the experimental rate needs to be given in terms of events per unit time per unit mass which is achieved by dividing the experimental values by the exposure. For XENON1T, this is 1.0 t×yr. At the 95% confidence level, the resulting limit on the rate is 1.58 × 10 −4 kg −1 day −1 .
[51]. For the cross section, our benchmark points have two mediating particles between the dark and visible sectors, namely h 1 and h 3 . The cross section for DM species i and a nucleon is thus given by: where c θ and s θ are cos θ and sin θ, respectively, f N is the SM Higgs effective coupling to nucleons which is approximately 0.30 ± 0.03 [52], m n is the nucleon mass which we take to be 0.93895 GeV, m i is the mass of the DM species in question, v is the SM vev, and g iih1 and g iih3 are the couplings between the DM and h 1 and h 3 , respectively. In the case of scalar DM, we take the Feynman rule to be of the form −ig, and in the case of vector DM, the Feynman rule is of the form −igg µν . To transform Eq. 100 to a cross section with the entire nucleus, we use [15]: where A is the atomic number of the nucleus, µ iN is the reduced mass between the nucleus and the DM particle, and µ in is the reduced mass between the nucleon and the DM particle. The scalar DM particles (χ and H) interact weakly with the nuclear target. This is because the scalar-h 1 coupling is weak due to our choice of λ H11 = 0 and λ H22 to be small and h 3 interacts with nucleons via its SM Higgs component which is small (recall sin θ = 0.1). As a consequence, when DM is mainly scalar it is weakly constrained by direct detection. The results for each benchmark point are shown in Fig 4. Again, the gap between half the Higgs mass and roughly 90 GeV for benchmark point A is because the value of the coupling that gives the correct relic abundance is not allowed because it violates the unitarity constraint. In  Fig 4(a), we see that there is an allowed region at high mass and in the Higgs resonance region. In Fig 4(b), the gauge coupling is quite small and has a large fraction of scalar DM so direct detection does not impose any significant constraints on benchmark point B. In Fig 4(c), we see that direct detection excludes most of the mass region due to the large vector fraction and large gauge coupling. There is however the tip of a resonance which allows for efficient annihilation of vectors around 340 GeV which remains allowed. This corresponds to when the mass of A 3 is just below 300 GeV and where there is significant annihilation of A 3 through the 600 GeV h 2 .
We see that f A is non-trivial in each case, that is to say, it is not close to being exactly 1 or exactly 0. This non-trivial behavior tells us that it is important to consider all the relevant particles when doing the freeze-out calculation, not just the lightest or most weakly coupled. In the cases where a resonance effect is present, it is particularly important to be precise, as these regions require quite different couplings constants to obtain the observed relic density.

C. DM Indirect Detection
Dwarf spheroidal satellite galaxies (dSphs) are typically DM dominated so are a good place to study DM [53,54]. The Fermi collaboration [55] has acquired 6 years of data from observing 15 dSphs and have released bounds for WIMP DM annihilation based on the observed gamma ray flux. They considered the following representative final states for the DM annihilation: e + e − , µ + µ − , τ + τ − , uū, bb, and W + W − [55].
We translate these bounds into constraints on our model by considering the cross sections of pairs of identical DM particles into the above final states scaled by their fraction squared and their branching ratio. Explicitly, we will compare: where σ ii→SM +SM is the cross section from two identical DM particles to one of the SM final states considered by the Fermi collaboration and σ tot ii is the sum of the cross sections to all possible final states. Of all the final states considered by Fermi, only the W + W − final state has the potential to provide useful constraints as the scaled cross sections to the other final states are orders of magnitude below the Fermi constraints. This is because the cross sections to these fermions are proportional to the mass squared of the fermion which are small compared to the W mass squared. Fig. 5 shows the results for each of the benchmark points. In all cases, the scalar DM ends up having a significantly larger scaled cross section to W + W − than any vector DM because the vector DM an-nihilates preferentially to other DM particles. We see that although the constraint is close, indirect detection does not impose any constraints. For benchmark point A the experimental limits are quite close to the theoretical predictions while for benchmark points B and C the experimental limits are about an order of magnitude above the theoretical predictions. Thus, in all three cases there is the potential to rule out large regions of parameter space with moderate improvements to the experimental bounds. It should also be noted that a more sophisticated analysis that considers the resulting spectrum from multiple DM species co-annihilating to the same final state could result in more stringent constraints.
These results lead to the interesting observation that if the DM is dominated by vector species, then it will be more easily detected in direct detection experiments while if the DM is dominated by scalar species, then it would more likely to be detected by indirect detection experiments. This could lead to an observation in one type of experiment, but no signal in the other and emphasizes the importance of the complementarity of multiple types of searches. Arcadi et al. [31] and many others have also pointed this out, e.g. [1-5, 15, 20, 35].
To end this section, we note that the constraints used are conservative. We only considered constraints from the annihilation of one species at a time. This bound could be improved with a more detailed study of the gamma-ray spectrum [56,57] or cosmic-ray spectra [58] resulting from multiple species annihilating or co-annihilating and comparing this to the observed spectrum from dSphs [23].   [55] while the solid lines are the predictions of the model. Note that the purple and green dashed line are almost superimposed because they have very similar masses leading to similar constraints. We also note that the solid orange line does not appear because it is too small.

VI. SUMMARY AND CONCLUSIONS
This paper reported on a study of multi-species DM, in particular, the Hidden Gauged SU(3) model. We calculated the DM relic density to constrain the parameters of the model. Using the resulting allowed values of the parameters we studied the prospects for observing multispecies DM in direct detection and indirect detection measurements. Before examining the Hidden Gauged SU(3) model, we studied the relic abundance of simplified scenarios to gain some insight into how adding different types of interactions between the DM species would affect the relic abundance. We then examined a specific model, the Hidden Gauged SU(3) model of Arcadi et al. [31]. With many DM particles, it was not possible to implement a detailed scan of the parameter space. Instead, we studied three representative benchmark points. In all cases, we fixed the gauge coupling to give the observed relic density. When the relic density was dominated by vector particles, we found that it could more easily be constrained by direct detection experiments while if the relic abundance was dominated by scalars, indirect detection has the potential to be more constraining in the future, even with conservative bounds. This reinforces the need to search for DM with a broad experimental program. An important result was that for the different benchmark points, the required value for the gauge coupling could vary quite significantly based on the choice of parameter values. This can lead to significantly different relic abundances among particle species which in turn, changes the detection prospects of a particular parameter point. We also saw that for each benchmark point, the ratios of vectors to scalars was never dominated by either scalars or vectors and that the DM consisted of some non-trivial mixture of both. This shows the importance of not only including the lightest stable particle in the freeze-out calculation but to include all relevant particles.
This case also results in more mixing between the gauge bosons. They are described by the following mass