The $Z_7$ model of three-component scalar dark matter

We investigate, for the first time, a scenario where the dark matter consists of three complex scalar fields that are stabilized by a single $Z_7$ symmetry. As an extension of the well-known scalar Higgs-portal, this $Z_7$ model is also subject to important restrictions arising from the relic density constraint and from direct detection experiments. Our goal in this paper is to find and characterize the viable regions of this model, and to analyze its detection prospects in future experiments. First, the processes that affect the relic densities are identified (they include semiannihilations and conversions) and then incorporated into the Boltzmann equations for the dark matter abundances, which are numerically solved with micrOMEGAs. By means of random scans of the parameter space, the regions consistent with current data, including the recent direct detection limit from the LZ experiment, are selected. Our results reveal that the $Z_7$ model is indeed viable over a wide range of dark matter masses and that both conversions and semiannihilations play an important role in determining the relic densities. Remarkably, we find that in many cases all three of the dark matter particles give rise to observable signals in future direct detection experiments, providing a suitable way to test this scenario.

In the literature it is possible to find different scenarios with three dark matter particles, but most of them rely on three independent Z 2 symmetries to stabilize them. Among them, let us mention the following alternatives for the dark matter particles: one vector and 2 fermions [32], one vector or a scalar and at least two scalars [36,62], one vector, one fermion and one scalar [63], two fermions and one scalar [33], and three scalars [43]. The main novelty of our model is thus the much simpler assumption of a single Z 7 symmetry, which may be a remnant of a U (1) gauge symmetry that has been spontaneously broken at higher energies. Such a symmetry entails novel conversion and semiannihilation processes among the dark matter particles that affect their relic densities, opening up new viable regions of parameter space. Finally, let us mention that this study provides the first detailed application of a new version of micrOMEGAs for N-component dark matter that will be released shortly.
The rest of the paper is organized as follows: the model is introduced in the next section while its dark matter phenomenology is described in section III; our main results are presented in section IV, where the results of different scans of the parameter space are analyzed; in section V our main findings are contrasted against those from other Z N models; finally our conclusions are drawn in section VI.
These new fields-the dark matter particles-are assumed to be singlets of the SM gauge group whereas the SM particles are taken to be singlets under the Z 7 . The most general Z 7 -invariant scalar potential is then given by µ 2 i |φ i | 2 + λ 4i |φ i | 4 + λ Si |H| 2 |φ i | 2 + λ 412 |φ 1 | 2 |φ 2 | 2 + λ 413 |φ 1 | 2 |φ 3 | 2 + λ 423 |φ 2 | 2 |φ 3 | 2 + 1 2 µ S1 φ 2 1 φ * 2 + µ S2 φ 2 2 φ 3 + µ S3 φ 2 3 φ 1 + 2µ S4 φ 1 φ 2 φ * 3 + H.c. where H is the SM Higgs doublet. Notice that this model can describe scenarios with one, two, or three dark matter particles, depending on the relations among the masses of the new scalar fields. Figure 1 displays, in a ternary plot, the regions consistent with these different possibilities. The case of a single dark matter particle is rather similar to the scalar complex singlet model [7][8][9] whereas the two-component scenario shares most of the features of the Z 5 model [42]. In this work, we are interested only in the case where a three-component dark matter scenario is realized -the central red region in figure 1. Therefore, we assume in the following that φ 1,2,3 do not acquire a vacuum expectation value and that their masses satisfy These conditions are obtained directly from equation (2) and guarantee the stability of φ 1 , φ 2 and φ 3 .
In total, this model contains 22 new parameters (7 dimensionful and 15 dimensionless), but three of them -λ 4i -are irrelevant for the dark matter phenomenology and can be ignored in our analysis. The parameters µ 2 i (i = 1, 2, 3), on the other hand, can be conveniently traded for the physical masses M i of the scalar fields, so that the free parameters of the model may be taken to be M i (3), λ Si (3), λ 4ij (3), µ Si (4), and λ 3i (6). Without loss of generality, we will restrict our numerical analysis to the hierarchy M 1 < M 2 < M 3 while leaving the couplings free. For simplicity, Our goal in this paper is to study the dark matter phenomenology of this model. In particular, we will identify and characterize the regions, within its multi-dimensional parameter space, that lead to scenarios consistent with all present data, including the relic density and the limits from dark matter experiments. In addition, we will show that current and future direct detection experiments offer a promising way to test this scenario.

III. DARK MATTER PHENOMENOLOGY
In this model, the dark matter particles, φ 1,2,3 , couple to the SM fields only via the Higgsportal -see the λ Si terms in 2. In addition, they couple among themselves through the different trilinear and quartic terms allowed by the Z 7 symmetry. As a consequence, a large number of 2 → 2 processes can contribute to their relic densities -see figure 2 for some examples. Table I shows the complete list of processes for each of the dark matter particles, classified according to their type. They include, besides the usual annihilations into SM particles, several examples of semiannihilations and of dark matter conversions, which are typical of Z N models [39]. Half of the different types of processes, in fact, had not been studied before, as they involve three different dark matter particles. Notice that semiannihilations always involve a combination of a trilinear interaction with the higgs-portal (top diagram in figure 2), whereas dark matter conversions can be obtained either via a single quartic interaction or two trilinear ones (bottom diagrams in figure   2).
All these 2 → 2 processes must be taken into account in the set of Boltzmann equations that determine the evolution of the dark matter abundances and, ultimately, the relic densities. Let us now present the equations for the evolution of DM abundances in generic multi-component scenarios, written in a way that allows a simple numerical implementation.
We assume that all DM particles are split into sectors, labelled by letters of the Greek alphabet, and that chemical equilibrium is maintained within each sector. The Z 7 model is particularly simple because it consists of three dark matter sectors, each containing just one dark matter particle. In this paper we assume kinetic equilibrium of all DM particles with SM bath particles, so the momentum distribution of all DM particles is described by the bath temperature T , but the number density can deviate from the thermal equilibrium case and this deviation is described by a common factor in each sector.
As usual, we write equations for the abundances where n α is a particle number density and s(T ) is the entropy density. In case of thermal equilibrium in a given DM sector Here and below we use Greek indexes for thermal sets and Latin ones for particles. m i and g i are, respectively, the mass and number of degrees of freedom of particle i. We will ignore effects of Fermi and Bose statistics. Thus, we can calculate the rate of reaction for thermal equilibrium distributions and then rescale it for the case of real DM density. The number of events of a, b → c, d collision in a space-time unit for equilibrium densities is where C is a combinatoric factor: C ab = 1/2 if a = b and 1 otherwise. The thermal balance equation readsN a,b→c,d =N c,d→a,b We introduce the function In this notation, the equations for the evolution of DM abundances caused by 2 → 2 reactions where H(T ) is the Hubble expansion rate. Taking into account entropy conservation ds dt = −3Hs we get By solving these equations, the relic densities, Ω 1,2,3 , can be calculated. We rely on mi-crOMEGAs for numerically solving this set of coupled equations. Since its version 4.0, mi-crOMEGAs incorporated two-component dark matter scenarios, and has been used in the analysis of different Z N models [42,45,57]. For this work, the code has been extended to models with three dark matter particles, a version that will become public shortly [64].
Since we are dealing with a three-component dark matter scenario, it is the sum of the relic densities of each dark matter particle, Ω 1 + Ω 2 + Ω 3 , that must be compared against the observed value, Ω DM . Here Ω DM is the dark matter abundance as reported by the PLANCK [1] collaboration with an estimated theoretical uncertainty of order 10%, As we will see, this dark matter constraint imposes severe restrictions on the viable parameter space of this model.
To estimate the relevance of the different kinds of processes-annihilations, semi-annihilations, and conversions-that contribute to the relic density of φ 1 , it is convenient to use a function in micrOMEGAs that allows to omit a set of processes from the relic density calculation. For example, ExcludedFor2DM="2211" will remove all conversion processes involving pairs of φ 2 and φ 1 [65]. In particular in the following section we will compare for each component the relic density obtained from including all processes (Ω i ) with the ones taking into account only processes involving DM pair annihilation into SM particles (Ω i ).
Another relevant quantity that must be taken into account when studying the feasibility of the model is the dark matter direct detection cross section. The elastic scattering of the dark matter particles off nuclei are possible thanks to the Higgs-portal interactions, driven by λ Si . The expression for the spin-independent (SI) cross-section then reads where µ R is the reduced mass, m p the proton mass and f p ≈ 0.3 is the quark content of the proton.
Given that we have multiple dark matter particles, the quantity to be compared against the direct detection limits provided by the experimental collaborations is not the cross section itself but rather the products ξ i σ SI φ i with ξ i ≡ Ω i /Ω DM . To impose current limits from the LZ experiment [66], we use micrOMEGAs which calculates the 90% confidence level exclusion in DD experiments using the recoil energy distribution of nuclei [64,67]. For multicomponent DM model, the total energy distribution is calculated where the contribution of each component is weighted by the factor ξ i . Note that we have explicitly checked that in all cases where the DM components have similar masses this procedure gives a similar result than the one obtained by simply comparing i ξ i σ SI φ i with the 90% exclusion limit. In this work we also consider, for each of the dark matter particles, the detection prospects in future experiments, in particular DARWIN [68].

IV. THE VIABLE REGIONS
In this section we present our main results: the viable regions of the Z 7 model are found and characterized, and its detection prospects are analyzed. To that end, we randomly scan the parameter space of the model in different ways, so as to identify the most relevant effects and the parameters that determine them. From each of these scans, we obtain a large sample of viable models. A model is considered viable if it is consistent with all known experimental constraints, including the relic density and the LZ direct detection limit. The viable models are then projected onto different planes for each random scan, and their implications separately discussed.
Among the free parameters of this model, the dark matter-Higgs couplings (λ Si ) are the most relevant because they link the SM and the DM sectors, allowing the dark matter particles to reach thermal equilibrium in the early Universe -a prerequisite for freeze-out. In addition, these same couplings determine the dark matter direct detection cross sections, which provide the main avenue to test this scenario. Consequently, in all our scans the couplings λ Si will be taken to be different from zero. In fact they vary within the range: All other couplings could in principle vanish. If they all vanished simultaneously, however, the model would become equivalent to three copies of the complex singlet scalar model, which is known to be viable only for heavy DM. Indeed an upper limit on the coupling λ Si is derived from direct detection and that coupling cannot provide efficient DM annihilation in the early Universe unless the DM mass is around 6 TeV. That case is of no interest to us. What we want to study is precisely if and how the new couplings present in the Z 7 model allow to satisfy the cosmological and experimental bounds for masses smaller than those in the singlet scalar case. Since φ 1 is the lightest dark matter particle (M 1 < M 2 < M 3 ), these new couplings must facilitate the annihilation of φ 1 in the early Universe. And we already know, from the study of previous Z N models, that the trilinear couplings are more effective than the quartic ones in opening up new regions of parameters space.
Taking these facts into consideration, we have done a series of scans in which only a subset of the free parameters of the model are different from zero at a time. In this way we expect to identify the most relevant couplings as well as the new processes that allow to satisfy the experimental constraints. The four different scans we will analyze in detail are: where, in each case, the rest of free parameters has been set to zero. For instance, in scan C) the parameters µ S4 , λ 4ij , λ 3i were all set to zero. In all the scans, the parameters that are different from zero are allowed to randomly vary within the following ranges: 100 GeV ≤µ Si ≤ 10 TeV.
Next, the results for each of the scans A-D are separately presented. Since M i , λ Si = 0 is common to all, we will label them according to the other couplings that are not zero.
This scenario presents semiannihilation processes that involve the three dark matter particles, with φ k (φ j ) acting as the mediator in the s (t) channel. These processes have no counterpart within two-component dark matter scenarios and their impact on the relic densities has not been investigated before. Notice that the number density of the two incoming particles is reduced by one unit while the number density of the outgoing particle is increased by one unit. In this way, the three relic densities can be affected by a single trilinear interaction.     The detection prospects for this scenario are illustrated in figure 4, which compares the elastic scattering cross section for each of the dark matter particles against the current limit from LZ (solid line) and the expected sensitivity of future experiments (dotted lines). Remarkably, we find that all the dark matter particles might give rise to observable signals in future experiments. That is, even φ 2 and φ 3 , whose contributions to the dark matter density are below 1% (see previous figure), could be detected. What explains this counterintuitive result is the fact that the smaller dark matter fraction is compensated by larger Higgs-portal couplings -see bottom-left panel in figure 3. Note, in particular, that a significant fraction of the viable models in our scan lies within the expected sensitivity of DARWIN. Direct detection, therefore, constitutes a promising way of probing this scenario.
B. µ S1 = 0, µ S2 = 0, µ S3 = 0 For our second scan we allow the trilinear terms µ S1 , µ S2 and µ S3 to simultaneously vary. The main novelty with respect to the previous scenario is the appearance of dark matter conversion processes involving two different trilinear couplings such as the s-channel process φ 1 + φ 3 → φ 2 + φ 2 mediated by φ 3 -another example of a process involving three different dark matter particles.
Conversion processes of the type 1122 are also relevant and they may all affect the relic densities of the three dark matter particles.
The set of viable models for this case is displayed in figure 5. Regarding the relic density, we learn from the top-left panel that it is always dominated by φ 1 , with φ 2 giving at most a 10% contribution, and a much smaller one from φ 3 . As in the previous case, conversion and semiannihilation processes provide the dominant contribution to the relic density of the three components, as shown in the top-right panel. For example for Ω 1 , conversion processes involving either φ 2 and φ 3 play a role while as before, φ 1 φ 1 → φ 2 h is the relevant semiannihilation process, it reduces the density of φ 1 when M 2 + m h < 2M 1 . This process is more important for large values of µ S1 . The next row demonstrates that, for our sample of viable models, the dark matter masses vary pretty much over their entire allowed range. Thus, in contrast with the previous scan, no degeneracies are required in this case. The next two rows display the trilinear and the Higgs-portal couplings. Notice that µ S1 tends to increase with M 1 whereas µ S2 and µ S3 vary within their allowed range without a clear pattern. From the bottom-right panel, we see that, as a result of the direct detection constraint, λ S1 (blue points) is suppressed, with a maximum value below 0.1, while λ S2 (red points) tends to be large, foreseeing a significant direct detection rate. |µ S1 | (GeV) FIG. 6: Detection prospects in the scan with µ S1 , µ S2 , µ S3 = 0 (B). The solid line is the current LZ limit whereas the dash-dotted and dotted lines correspond to the expected sensitivities of LZ and DARWIN respectively . one dark matter particle, but it would be more challenging to establish that it corresponds to a three-component dark matter scenario. In any case, direct detection provides a robust way to test this possibility.
C. µ S1 = 0, λ 4ij = 0 For the third scan, we allow the quartic couplings λ 4ij as well as the trilinear term µ S1 to randomly vary. Since the µ S1 interaction does not affect φ 3 , its semiannihilation processes are absent and it is the conversions induced by λ 4ij (φ 3 φ † 3 → φ i φ † i ) that play a key role in reducing Ω 3 when λ S3 is small. As we will see, this fact allows φ 3 , the heaviest dark matter particle, to give a significant contribution to the total dark matter density, even the dominant one. The parameter µ S1 will induce semi-annihilation processes φ 1 φ 1 → φ 2 h that will reduce the abundance of φ 1 .
The viable parameter space for this case is illustrated in figure 7. An important novelty with respect to the two previous scans is that |µ S1 | (GeV) FIG. 8: Detection prospects in the scan with µ S1 , λ 4ij = 0 (scan C). The solid line is the current LZ limit whereas the dotted lines correspond to the expected sensitivities of LZ and DARWIN.
semiannihilations. In particular for φ 1 , the coupling µ S1 will induce semi-annihilation processes can vary over several orders of magnitude, with λ 413 and λ 423 clustered toward higher values. These couplings are required to efficiently deplete the φ 3 relic density via the conversions φ 3 + φ † 3 → φ 2,1 + φ † 2,1 . Finally, from the bottom-center panel, we notice that, due to the direct detection constraints, λ S1 and λ S3 are both suppressed by an order of magnitude with respect to λ S2 .
The detection prospects for this case are shown in figure 8. Note that all three particles may give rise to observable signals in DARWIN. Since the dark matter masses are in general non-degenerate, one might be able to differentiate the three signals, excluding all vanilla models of dark matter and confirming the multi-component nature of the dark sector.
In the previous scans, the dominant contribution to the dark matter density was due to φ 1 or φ 3 , so a natural question is whether φ 2 might be responsible for the bulk of the dark matter in certain cases. With this scan we answer that question in the affirmative and we also incorporate into the discussion one of the λ 3i couplings, which are specific to the Z 7 symmetry. Besides λ 34 ,  the couplings that are different from zero are µ S3 and λ 4ij . Figure 9 displays the results for this scan, which, unlike the previous ones, does not feature dark matter masses above 400 GeV or so. Thus, the dark matter masses are severely restricted in this case. Note that indeed the dark matter density is dominated by φ 1 or φ 2 , but never by φ 3 , whose contribution is always below 1% -see top-left panel. Once again dark matter conversions and semiannihilations play the dominant role in setting the relic density of the three components, we find that Ω i < Ω i except when M 1 ≈ m h /2 -see top-right panel. From the second row, we learn that the dark matter particles tend to be degenerate, with the ratios M 2 /M 1 and M 3 /M 2 rarely exceeding 1.2 outside the Higgs-resonance, and decreasing rapidly with M 1 . The trilinear coupling µ S3 increases with M 1 , reaching 10 TeV (the maximum value allowed in the scan) at M 1 ∼ 400 GeV, explaining why there are no viable models at higher masses. This trilinear coupling enters the semiannhilation processes φ 1 φ 3 → φ † 3 h and φ 3 φ 3 → φ † 1 h which contribute to the relic density of φ 1 and φ 3 respectively. Conversion processes between the first and third sector, for example those involving the quartic coupling λ 34 , also impact the relic density of φ 3 . All these processes lead to a very small value for Ω 3 , hence direct detection constraints can be evaded despite large values for λ S3 -see bottom right panel. The quartic coupling λ 412 plays an important role as it induces the conversion process φ 2 φ † 2 → φ 1 φ † 1 that reduce the density of φ 2 when the annihilation processes φ 2 φ † 2 → hh, W W... are not efficient because of the small value of λ S2 . The direct detection prospects for this case are presented in figure 10. There, we see that all three dark matter particles may give rise to observable signals in future direct detection experiments.
Given the degeneracy of the masses, however, it is unlikely that their signals could be disentangled from one another, so as to establish the multi-component nature of the dark matter.
As a final remark, the results obtained for this scan hold even in the case when all the quartic couplings λ 3i are allowed to vary in the same range. points (in red) cover most of that area, distributing themselves more or less uniformly over it. In contrast, in scans A (top-left) and D (bottom-right), the viable points tend to gather toward the central point, leaving some unpopulated areas. These results are consistent with the fact, already mentioned, that the latter require some degeneracy among the dark matter masses, whereas the former do not. Taking all the scans into consideration, we would conclude, then, that the viable points do span the full available region.
All the scans studied in the previous section present a common characteristic, namely at least one trilinear interaction term involving φ 1 is allowed, so as to facilitate the depletion of its relic density. The heavier dark matter particles can, instead, reduce their relic density via conversions into the lighter ones, without a kinematic suppression. Consequently, scans having only quartic interactions of the type λ 4ij or λ 3i do not lead to viable points unless a highly degenerate mass spectrum is invoked. That is the reason why we did not consider such possibilities.
We did carry out additional scans, which were not discussed in detail. For example, a scan where the only trilinear coupling varied is µ 2 results in a mass degeneracy between φ 1 and φ 2 , independently of the rest of the couplings. In fact, since the µ 2 interaction gives rise to neither semiannihilations nor conversion processes reducing the number density of φ 1 , the conversion processes mediated by the quartic couplings λ 412 , λ 31 and λ 36 (those affecting the number density of φ 1 and φ 2 ) must be efficient, implying a degeneracy. In another scan, we considered the possibility of variations of the trilinear couplings µ S1 and µ S2 . In this case, the relic densities can be significantly decreased, leading to a viable and interesting scenario. These results, however, are rather similar to those for the scan B, where µ S3 is also varied. In summary, the scans we have considered in the previous section show the main phenomenological aspects of Z 7 the model, illustrating its versatility and capability as a viable scenario with three dark matter particles.
We have not examined in detail the possible signals in indirect detection experiments since it is unlikely they will further restrict the viable parameter space. Such signals are expected to be similar to those found in the Z N two-component scenarios [42,45], where the semi-annihilation process φ 1 + φ 1 → φ 2 + h turns out to be the most relevant one, featuring a cross section that can reach values as high as 10 −25 cm 3 /s. If that value is extrapolated to the Z 7 model for the processes Now we proceed to briefly compare our results for the Z 7 model with those previously obtained for the Z 5 model of two-component dark matter [42]. This model comprises the trilinear scalar interactions From this, two simple scenarios may arise, depending on the free parameter that is neglected.
Whenμ S1 = 0 (whose interaction term is analogous to µ S1 in the Z 7 model), it was found that the ratio M 2 /M 1 varies over a wide range, and that the semi-annihilations entail that φ 1 always gives the dominant contribution. These findings are in utter agreement, as concerns φ 1 and φ 2 , with those obtained here in the scan C (µ S1 = 0, λ 4ij = 0). Whenμ S2 = 0, it was found instead that: M 1 becomes more restricted; the dark matter particles have to be degenerate; and the conversion process φ 1 + φ 1 → φ 2 + φ 2 induced byμ S2 plays a main role in decreasing Ω 1 . Since the analogous interaction in the Z 7 model turns to be µ S3 with the identification φ 2 → φ 3 , the results of the scan D (µ S3 = 0, λ 4ij = 0, λ 34 = 0) match with those deduced from the Z 5 model. All these observations illustrate the compatibility of our results with those for Z N two-component scenarios.
Let us finally comment about the applicability of our results to other Z N scenarios with three dark matter particles. The lowest Z N allowing simultaneous stability for three scalar fields is a Z 6 with the dark matter particles transforming as (φ 1 ∼ w 6 , φ 2 ∼ w 2 6 , φ 3 ∼ w 3 6 = −1), thus becoming a scenario with two complex fields and one real field. In this case the new trilinear and quartic interactions are given by [39] It follows that this model shares the µ S1 , µ S4 and λ 34 interaction terms with the Z 7 model. This means that the results of the scans A (µ S4 = 0) and C (µ S1 = 0, λ 4ij = 0) apply to the Z 6 model. Following the same line of thought, this conclusion extends to arbitrary Z N models [39] with dark matter fields [φ 1 , φ 2 , φ 3 ], since the µ S1 and µ S4 trilinear interactions (along with the λ 4ij quartic interactions) are always Z N invariant. Furthermore, since the Z 8 , Z 9 and Z 10 models with a different set of dark matter fields can lead to scenarios having interaction terms of the type µ S1 , µ S4 or both, they are also viable scenarios in light of our results for the Z 7 model.

VI. CONCLUSIONS
We proposed and investigated a novel and simple scenario that can explain the dark matter by extending the Standard Model with just two ingredients: a Z 7 discrete symmetry, and three complex scalar fields (the dark matter particles) that have different charges under it. This threecomponent dark matter model intrinsically differs from the two-component scenarios considered so far due to the semiannihilation and conversion processes involving three different dark matter species, which had not been studied before. After introducing the model and qualitatively discussing its dark matter phenomenology, we numerically studied its parameter space so as to determine the regions that are consistent with current data. To that end, a series of random scans were implemented, and the viable points were selected for further analysis. Our results clearly indicate that this model is indeed viable over a wide range of dark matter masses. Semiannihilations and dark matter conversions turn out to both play important roles in determining the relic densities, and, depending on the parameters, any of the scalar fields could give the dominant contribution to the observed dark matter density. We found that, for many viable points, the three dark matter particles may be observed at future direct detection experiments, providing an avenue to differentiate this scenario from the conventional one, and to test the multi-component nature of the dark matter. Finally, we argued that many of our results can be generalized to scenarios with a different Z N symmetry. All in all, the Z 7 model of scalar dark matter was found to be a viable and testable scenario for physics beyond the Standard Model that can account for the dark matter of the Universe.