WIMP and SIMP Dark Matter from the Spontaneous Breaking of a Global Group

We propose and study a scalar extension of the Standard Model which respects a $\mathbb{Z}_3$ symmetry remnant of the spontaneous breaking of a global $U(1)_\text{DM}$ symmetry. Consequently, this model has a natural dark matter candidate and a Goldstone boson in the physical spectrum. In addition, the Higgs boson properties are changed with respect to the Standard Model due to the mixing with a new particle. We explore regions in the parameter space taking into account bounds from the measured Higgs properties, dark matter direct detection as well as measurements of the effective number of neutrino species before recombination. The dark matter relic density is determined by three classes of processes: the usual self-annihilation, semi-annihilation and purely dark matter $3 \to 2$ processes. The latter has been subject of recent interest leading to the so-called `Strongly Interacting Massive Particle' (SIMP) scenario. We show under which conditions our model can lead to a concrete realization of such scenario and study the possibility that the dark matter self-interactions could address the small scale structure problems. In particular, we find that in order for the SIMP scenario to work, the dark matter mass must be in the range $7-115$ MeV, with the global symmetry energy breaking scale in the TeV range.


Introduction
The existence of dark matter (DM) comprising approximately 85% of the matter content in the Universe, as supported by several astrophysical and cosmological observations, is a strong evidence that the highly successful Standard Model (SM) of particle physics is incomplete: there must be new physics beyond the Standard Model (BSM physics) (for a review, see e.g. [1][2][3][4][5]).
If DM is made out of new particles, these must be electrically neutral, have a lifetime longer than the age of the Universe and be produced in the observed amount. There are several models containing particles that are suitable candidates for DM, such as supersymmetric (SUSY) extensions of the SM with a discrete R-parity symmetry that makes the lightest supersymmetric particle stable. However, at this point there is no evidence for this class of models from searches conducted at the Large Hadron Collider (LHC) and the simplest SUSY models are under stress due to naturalness issues (see, e.g. [6]).
Another class of BSM scenarios with a DM candidate consists in extending the SM with an augmented scalar sector furnished with a discrete symmetry ensuring the stability of the DM particle. If the extra scalars are postulated to be SM singlets, they can interact with SM matter only through their couplings to the Higgs doublet, the so-called Higgs portal scenarios.

JCAP04(2015)012
This class of models may lead to invisible Higgs decays, allows for strongly self-interacting DM [7,8] and may also alleviate the vacuum stability issue of the SM [9].
The simplest model in this class is the SM with an extra real scalar field respecting a discrete Z 2 symmetry [10][11][12]. Constraints on this scenario have been recently updated [13]. In fact, requiring the model to respect the observed DM abundance, direct detection bounds and the limits on the invisible Higgs width severely constrains the parameters, placing the extra scalar mass in the range between 53 and 63 GeV (see also [14]).
Another possibility is to consider scalar DM models with a Z 3 symmetry. One particular feature of these models is that they allow for semi-annihilation processes [15][16][17][18], where the annihilation of two DM particles can give rise to a DM particle in the final state, as for example XX →XΨ, where X and Ψ stand for a DM and a non-DM particle respectively. In addition, they predict other interesting number-changing processes involving only DM particles, such as XXX → XX. These inelastic processes can be relevant for the chemical equilibrium and freeze-out of DM particles and their effect must be included in the Boltzmann equation controlling the DM density. Remarkably, Z 3 models with that feature can also address the 'core versus cusp' [19][20][21][22] and 'too big to fail' [23,24] problems because they may exhibit sizable DM self-interactions. In this case, the DM in this type of scenarios has been dubbed Strongly Interacting Massive Particle (SIMP) [25]. 1 Z 3 models have been recently studied, assuming that this symmetry is either just a global symmetry as in [27] or that it is a remnant of a local U(1) DM symmetry spontaneously broken to Z 3 [28].
In this paper we explore an intermediate situation, where a global U(1) DM symmetry is spontaneously broken to a discrete Z 3 symmetry by postulating the existence of a scalar field with charge 3 under the global group and a non-vanishing vacuum expectation value [29]. 2 Due to this, in addition to a natural DM candidate there is a Goldstone boson (GB) in the physical spectrum. In our study, we include the novel processes mentioned above in the determination of the DM relic abundance and consider the role of GB in its production. In doing so we take into account bounds from the measured Higgs properties, DM direct detection as well as measurements of the effective number of neutrino species before recombination. In addition, we find under which conditions our model can be a concrete implementation of the SIMP scenario. This paper is organized as follows. In the next section we introduce the model and discuss the different mechanisms that could generate the DM relic density. In section 3 we analyze the constraints arising from the effective number of neutrino species, Higgs physics and direct detection experiments. Following that, in section 4 we focus on the self-and semiannihilating scenarios and show that in these cases DM behaves like a Weakly Interacting Massive Particle (WIMP). Subsequently in section 5, we address the possibility of having DM as a Strongly Interacting Massive Particle (SIMP). In particular we study how to implement the 3 → 2 annihilation mechanism for its production and show that this naturally leads to self-interacting cross-sections that typically address small scale structure problems in astrophysics. Finally, in section 6 we present our conclusions.

JCAP04(2015)012
2 Z 3 dark matter from the breaking of a global U(1) DM

Description of the model
We postulate a dark sector with a global U(1) DM which is spontaneously broken into a Z 3 symmetry. This is achieved by considering, in addition to the usual scalar doublet H, two complex scalar fields φ X and X which are singlets under the SM group. While the SM particles do not transform under the U(1) DM , the new scalars φ X and X have charges equal to 3 and 1 respectively. This choice leads naturally to a Z 3 symmetry after the spontaneous breaking of the continuous group [29]. In fact, if we assume that in the vacuum state X = 0 and φ X ≡ v φ = 0, then the breaking SU (2) To see this explicitly, let us consider the most general renormalizable scalar potential After symmetry breaking we can write the fields in the unitary gauge as where v H = 246 GeV is the vacuum expectation value of H, which gives mass to the SM particles via the Brout-Englert-Higgs mechanism [31]. Therefore the part of the potential involving X becomes where The potential in eq. (2.3) is manifestly invariant under the remnant Z 3 ⊂ U(1) DM that acts non-trivially only on the field X in the following way: X → e 2π 3 i X. It is precisely because of this reason that the particle associated to X can not decay and is thus identified with the DM. In other words, the Z 3 group ensures the stability of the DM in the present model. Notice also that while X andX belong to different representations of Z 3 , they have the same mass and therefore they both contribute to the total amount of DM.
Due to the λ φH coupling there is a mixing between theh andρ fields, which determines the masses for the physical states h and ρ where we identify h with the SM Higgs boson with mass m h ∼ 125 GeV and the mixing angle θ is defined by

JCAP04(2015)012
From this it follows that the masses and v φ are related by On the other hand, the field η remains massless and it is hence the GB associated to the symmetry breaking of the global U(1) DM . Consequently, the physical spectrum of this model consists of the real scalars h, ρ and the GB η, along with the complex scalar X. There are ten free parameters in the potential. Fixing m h and v H to their physical values reduces the number of independent parameters to eight. We will use the following parameters to characterize the model: two masses (m ρ and m X ), the mixing angle θ and five quartic couplings λ φ , λ X , λ 3 , λ φX and λ HX . According to the potential (2.1), these couplings are real with the exception of λ 3 , whose phase can be absorbed by a field redefinition, and therefore the dark sector does not introduce any CP -violation.
Let us point out that, with the exception of the term involving λ 3 , all the terms of eq. (2.3) contain pairs of DM fields. Hence, the λ 3 coupling characterizes the Z 3 DM phenomenology. On the other hand, λ HX controls the so-called Higgs portal because it connects the SM with the DM particles via the Higgs boson. Likewise, λ φX connects the ρ and η bosons with the DM, and we therefore call it the GB portal.
Some comments are now in order. Stability of the potential requires that it should be bounded from below for large values of the fields. One could also require perturbativity of scattering amplitudes up to high energy scales, such as the Planck scale, which typically imply that the coupling constants should not be too large at the electroweak scale. Some of these issues were discussed in [28]. For the purposes of this phenomenological analysis, we just assume that these conditions are met; however, we do require our couplings to be in the perturbative regime, |λ| < 4π at the electroweak scale.
The model was implemented in FeynRules [32,33] and the output was used both in CalcHEP [34] and MicrOMEGAs [35,36] in the phenomenological studies discussed in the following sections.

Dark matter relic abundance
Due to the Z 3 symmetry, in addition to the usual self-annihilation of X andX, there are other processes that change the number of DM particles, and hence contribute to the thermal abundance of DM. For instance, two X particles can semi-annihilate into aX and a Z 3 singlet state, or three X particles can annihilate into a XX pair. All these processes can be divided in three classes: self-annihilation [37], semi-annihilation [15,16] and 3 → 2 annihilations [25]. The first category corresponds to the usual processes where two DM particles annihilate into non-DM particles. In addition, semi-annihilation occurs when two DM particles produce a single DM particle in the final state. Finally, 3 → 2 processes involve only DM particles both in the initial and final states.
If we assume no asymmetry between X andX, the total DM number density n associated to both particles satisfies the Boltzmann equation [25,38] dn dt + 3Hn = −2 n n eq 2 − 1 γ self − n n eq n n eq − 1 γ semi − n n eq 2 n n eq − 1 γ 3→2 , (2.9) where H is the Hubble parameter and n eq is the equilibrium number density for a given DM temperature T . The thermally-averaged interaction rates γ are defined in appendix A, and  Figure 1. Freeze-out of the DM particles for the different thermal production modes, for selfannihilation, semi-annihilation and the 3 → 2 process, giving rise to the measured relic abundance (see text for details). The equilibrium density is also depicted. A DM mass m X = 1 GeV was assumed. are given in terms of the cross-sections by γ self = 1 2 n 2 eq σv self , γ semi = n 2 eq σv semi , γ 3→2 = n 2 eq σv 2→3 = n 3 eq σv 2 3→2 . (2.10) The three terms on the right-hand-side of eq. (2.9) correspond to the self-annihilation, the semi-annihilation, and the 3 → 2 annihilation processes, respectively. Here σv is the standard thermally-averaged velocity-weighted cross-section for 2 → n reactions and σv 2 is the equivalent quantity for 3 → n processes, as discussed in appendix A. Let us note that the Boltzmann equation should contain other terms corresponding to three DM particles annihilating into a two-body final state with none or one DM particle. However these interactions are subdominant (cf. appendix B) and will be neglected in the following analysis. Notice that σv 2→3 = n eq σv 2 3→2 as shown in appendix A using CP conservation and the principle of the detailed balance [39].
In this work we solve the Boltzmann equation (2.9) in order to estimate the DM relic density abundance today, which is given by where s is the entropy density of the photon plasma. That value must agree with the one measured by the Planck collaboration [40]: For the sake of illustration, figure 1 shows the evolution of the comoving number densities of the DM as a function of x ≡ m X /T obtained by solving the Boltzmann equation for the different thermal production modes considered separately: self-annihilation (red), semiannihilation (blue) and the 3 → 2 process (orange). We take m X = 1 GeV with constant σv self , σv semi and σv 2 3→2 chosen so that the final abundance equals the measured one.  Since this plot is only for illustration, we assume that all the relativistic particles have the same temperature (this might not always be the case, as shown in section 3.1). The corresponding equilibrium density is also depicted in black. As shown in the figure, although the 3 → 2 process freezes out at larger x compared to the self-and semi-annihilations, the former reaches its final DM abundance much faster.
In sections 4 and 5, we analyze further the different regions of parameter space where the self-annihilation, semi-annihilation and 3 → 2 processes are dominant and study the phenomenological consequences, taking into account the constraints discussed in the following section.

Constraints on the model
In this section we discuss some constraints arising from the production and decay of the SM-like Higgs, DM direct detection, as well as from the measurement of the effective number of neutrinos species before recombination.

Effective number of neutrino species
Being massless particles, GBs can fake the effects of neutrinos in the Cosmic Microwave Background (CMB). The existence of these new relativistic degrees of freedom, the so-called Dark Radiation, is characterized by the effective number of neutrino species N eff (for a recent discussion see, e.g. [41]). Measurements of the anisotropies in the CMB performed by the Planck collaboration in 2013 indicate that N eff = 3.30 ± 0.27 [40], consistent with the SM prediction. 3 In this model, the contribution to the value of N eff from the GB is determined by the moment in which it goes out of equilibrium with the thermal plasma. 4 If T d η and T d ν 3 After the completion of this work, new preliminary results form the Planck collaboration indicate an even better agreement with the SM prediction N eff = 3.13 ± 0.32 [42]. These results are consistent with BBN [43]. 4 The equilibrium between the SM and the dark sector is achieved in the Early Universe by a number of processes. For instance by hh → hρ, which is only suppressed by the mixing angle θ, and which could establish the equilibrium at temperatures above the TeV scale. We corroborate this for the scenarios considered in this work. are respectively the decoupling temperatures of the GBs and neutrinos from the plasma and if the GBs are not reheated after their decoupling from the SM, then where g(T ) is the number of relativistic degrees of freedom at a given temperature. We plot the resulting contribution as function of the decoupling temperature on the left panel of figure 2. From this it is clear that when the decoupling temperature is greater than the muon mass, N eff is reasonably close to the SM prediction. Consequently, we have to make sure that the interaction rates keeping the equilibrium between the GB and the SM become smaller than the Hubble expansion rate at temperatures greater than the muon mass. In fact, assuming m ρ greater than a few GeV is enough to have a diluted GB contribution to N eff . In order to show this, we consider the case T d η ≈ m µ (for a detailed analysis, see e.g., [44,45]). Then the equilibrium takes place via reactions such as η f ↔ η f and ff ↔ η η, where f denotes a SM fermion, as shown in figure 3. Since the corresponding amplitudes are proportional to the fermion mass, processes involving muons are dominant. Then, the departure of equilibrium occurs when the Hubble expansion becomes comparable to the interaction rate with muons. Exploiting the fact that the latter depends only on the three parameters m ρ , λ φ and θ (as well as on T d η ), in figure 2 we show the values of the mixing angle for which the interaction rate equals the Hubble parameter for different values of λ φ . The regions below the solid lines correspond to decoupling temperatures higher than the mass of the muon and thus to a diluted contribution to N eff . From figure 2 we conclude that such contribution is small when the ρ scalar is heavier than about 5 GeV or when the mixing angle is smaller than about 10 −5 .

Higgs sector
The enlarged scalar sector increases the number of Higgs decay channels. The new decay modes are into pairs of ρ, η and DM particles, when kinematically allowed. The last two channels (h → η η and h → XX) contribute to the invisible decay of the Higgs, whereas the first one (h → ρ ρ) only does it partially when the ρ pair itself decays into a couple of η or DM particles. The corresponding decay rates are: where one can see that its rate is proportional to λ 2 3 . Nevertheless, this channel is phase-space suppressed and therefore subdominant compared to the two-body decay modes.
Correspondingly, the decay channels for the ρ boson are In addition, due to theh −ρ mixing, the ρ boson can also decay into SM particles, with a rate give by Γ(ρ → SM SM) = Γ(h SM → SM SM) × sin 2 θ, taking m h SM → m ρ in the right-hand-side of the expression. Because of the same reason, all the decay widths of the h boson into SM particles decrease compared to the SM ones by a factor cos 2 θ.
With these decay rates, we calculate the invisible Higgs branching ratio and apply the bound reported in [46]. That analysis considers simultaneously a universal modification of the Higgs boson couplings to SM particles as well as the possibility of an invisible Higgs decay, concluding that where κ is the universal modification for the Higgs coupling. In our case we identify κ = cos θ and the invisible Higgs decay branching ratio with Although this is the most stringent constraint, we also require that the total decay width never exceeds the experimental bound [47] Γ tot Direct searches for a light neutral scalar produced in association with a Z in the OPAL detector [48] at LEP also apply. By identifying such a particle with the ρ scalar, it is possible to set an upper bound on the mixing angle.
We apply these limits for all the scenarios considered in this work. In the right panel of figure 2 we show them in the particular case when the scalars can not decay into DM particles. The dashed lines represent the upper bounds arising from the Higgs invisible decay of eq. (3.7), whereas the upper hatched region is excluded because of direct searches of the ρ scalar in OPAL.

JCAP04(2015)012
X h X q q X ρ X q q Figure 5. Diagrams responsible for DM direct detection. Figure 6. Tree level diagrams for the annihilation between X andX. Here f corresponds to any SM fermion and V to any massive gauge boson.

Dark matter direct detection
The scattering of DM particles off nuclei in direct detection experiments takes place via the t-channel exchange of h and ρ bosons as depicted in figure 5. Such process only generates a spin-independent signal with a DM-nucleon cross-section given by where m N denotes the nucleon mass and f N ≈ 0.27 is a constant that depends on the nucleon matrix element [36]. For GeV Dark Matter, the resulting expression is constrained by comparing it to the upper bound obtained by the LUX collaboration [49].

Weakly interacting dark matter
If self-or semi-annihilations dominate the production of DM, the latter is in kinetic equilibrium with the SM particles and the GB (see appendix A). As discussed in section 3.1, due to the constraint on N eff , this equilibrium can only exist at temperatures greater than the muon mass. Since the freeze-out temperature is roughly m X /25, we thus consider only selfand semi-annihilations in the GeV range. Because of this mass range and since the couplings involved in DM phenomenology are typically small, self and semi-annihilation scenarios describe Weakly Interacting Massive Particles (WIMP). We now discuss them separately.

Self-annihilating scenario
If the self-annihilation of DM dominates over the other two processes, that is if γ self is much greater than γ semi and γ 3→2 , the DM production proceeds via the familiar Lee-Weinberg scenario [37]. In this case, as shown in figure 6, the DM abundance is obtained by selfannihilation into η η, ρ ρ, h ρ and pairs of SM particles, until the reactions freeze-out at x ∼ 25.  DM annihilates into SM fermions and vector bosons via the s-channel exchange of a h or a ρ boson. For the other annihilation channels the contact term and the t-channel exchange of a X are also present. The only exception is the reaction XX → ηη where the t-channel exchange does not exist. All these diagrams are proportional to λ HX and λ φX as opposed to the ones corresponding to the semi-annihilation and the 3 → 2 mechanisms which are proportional to λ 3 . As a result, self-annihilation dominates whenever λ 3 is much smaller than λ HX and λ φX .
If we neglect the mixing angle, there are two regimes depending on the relative size of the Higgs and GB portals λ HX and λ φX . On the one hand, when the latter dominates, DM annihilates into the dark sector, that is, into GB and ρ pairs. The two channels are democratically distributed provided that m X ≫ m ρ . On the other hand, when λ HX ≫ λ φX , the annihilation produces SM particles pairs. Due to the s-channel Higgs exchange, the heaviest kinematically allowed pair tend to be produced. These two regimes are illustrated in the upper panel of figure 7, where we show the parameter space giving rise to the measured DM relic abundance,  for different combinations of m X , λ HX and λ φX . There we take m ρ = 5 GeV, sin θ = 10 −5 , λ φ = 0.1 as well as λ 3 = 0 in order to avoid semi-annihilations and 3 → 2 annihilations. We calculate the relic abundance using MicrOMEGAs-3. From the figure it is possible to see that large masses in general require large couplings. We can also see that when λ HX dominates, due to the existence of the Higgs resonance and kinematic thresholds, the dependence on the mass is more complicated than in the other case. This is exemplified in the lower panels of figure 7, where the dependence on the DM mass is shown explicitly. The Higgs boson funnel appears as well as the kinematic openings of the annihilations into W + W − , Z Z and h h. The LUX experiment rules out parts of the parameter space corresponding to DM masses O (10) GeV and sizable λ HX couplings. This corresponds to the hatched regions in figure 7.
Similarly, the left (right) panel of figure 8 shows in blue the regions of the parameter space that could give rise to the measures DM relic abundance in the plane [m X , λ HX ] ([m X , λ φX ]) after marginalizing with respect to λ φX (λ HX ); that is, each plot includes all the possible values of the corresponding quartic coupling. The empty regions, corresponding to high values for the couplings, always generate a DM underabundance. In addition, both the LUX and the Higgs invisible decay exclusions are shown in the figure as hatched regions. Again, those areas correspond to the points that could be excluded for some combinations of the parameters. The unhatched regions are always allowed.

Semi-annihilating scenario
In the previous section we assumed that λ 3 was negligible with respect to λ HX and λ φX . In general this may not be the case and therefore the semi-annihilation and the 3 → 2 process have to be taken into account. Although the latter process can still dominate, here we simply assume that it is subdominant and work under that hypothesis. In section 5.2 we will address the general case and discuss under which circumstances the 3 → 2 reaction can be safely ignored. Under these conditions, the generation of relic abundance occurs not only via the self-annihilation of DM particles but also via the semi-annihilation processes DM DM → DM S, where S could be h, ρ or η.  Figure 9. Tree level diagrams for the semi-annihilation process X X →X S where S could be h, ρ or η. tree-level are shown in figure 9. Notice that the semi-annihilation into GBs is always possible and can only be suppressed by taking λ 3 very small.
When λ 3 is much greater than λ φX and λ HX , self-annihilations can be neglected as well as the second and third diagrams of the semi-annihilation process in figure 9; in other words, only the contact interaction diagram contributes. Hence, the relic abundance scales like Ω DM h 2 ∼ 1/ σv semi ∼ (m X /λ 3 ) 2 . This behavior is illustrated in figure 10, where we depict the regions of the parameter space [m X , λ 3 ] giving rise to the observed relic density when λ φX = λ HX = 0. There we also show the equivalent regions when we relax this assumption, while still keeping λ φX = λ HX for the sake of simplicity. In this plot we choose m ρ = 5 GeV, sin θ = 10 −5 and λ φ = 0.1. In the case where λ φX = λ HX are sizable, both the Higgs boson funnel and the kinematic openings of the annihilations into W + W − , ZZ and hh are visible. Similarly, the LUX experiment rules out parts of the parameter space corresponding to DM masses O(10) GeV. These parts are shown in the figure 10 as hatched regions. As we did for the self-annihilation case, in figure 11 we show the region of the parameter space that could give rise to the measures DM relic abundance in the plane [m X , λ 3 ] after marginalizing with respect to λ HX and λ φX . The white region, corresponding to large values for λ 3 , always generates a DM underabundance. Moreover, assuming that the cross-section for the semi-annihilation into a GB is never greater than about σv semi ∼ 3 · 10 −26 cm 3 /s, this translates into the upper bound   Figure 11. Semi-annihilation scenario. Same as figure 10 but marginalizing with respect to λ φX and λ HX . We also include the constraint from the Higgs invisible decay width.
which can be seen in figure 11. In addition, we also show both the LUX and the Higgs invisible decay exclusion limits.

Dark matter self-interactions
In spite of the great success of the ΛCDM scenario in explaining the formation and evolution of cosmic structures, there are some problems at small scales (see, e.g. [50]). Simulations of this scenario result in cuspy density profile of halos towards its center and a large number of small satellite halos. These results are challenged by observations. Although the inclusion of baryonic effects can ameliorate the discrepancies, it is still unclear whether it is necessary to change the physics of the DM sector. For instance, self-interacting DM with a strength [51][52][53][54][55][56][57][58]: 0.1 cm 2 /g σ scatter m X obs 10 cm 2 /g (5.1) can help in solving the observed discrepancies. Yet, the Bullet Cluster [59][60][61] and recent analysis of the constraints from halo shapes [54][55][56], favor the following upper bound on the DM self-interacting cross section (at velocities greater than 300 km/s) σ scatter m X 1 cm 2 /g .
The corresponding elastic scattering processes in our model are shown in figure 12. From the diagrams it is possible to see that the exchange of a ρ or a h boson lighter than the DM can possibly lead to non-perturbative effects such as long range interactions. In this section we limit ourselves to the perturbative scenario and we therefore assume m ρ , m h ≫ m X . A posteriori we will find that this is also necessary for the DM production to take place via the 3 → 2 annihilation process (see section 5.2). Under this assumption the elastic scattering cross-section in the non-relativistic limit reads  Figure 12. Diagrams for the elastic scattering process X X → X X and XX → XX.  where Z ≡ λ 3 v φ /m X . The Z parameter can not be arbitrarily large. In fact, perturbative unitarity requires the s-wave amplitude of the process XX → XX to be bounded from above [62]. Considering m ρ and m h ≫ m X , this translate into which implies the following weaker bounds The process XX → XX also sets an upper bound on the same parameters, however it is never as stringent as the one of eq. (5.4). Figure 13 shows the regions of the plane [λ X , Z] fulfilling eq. (5.1) for different values of the DM mass. In addition, we depict the areas disfavored by the Bullet Cluster for those masses (i.e. eq. (5.2)) and we also hatch the regions where perturbative unitarity is lost. In this figure a limit on the DM mass is manifest, corresponding to m X 173 MeV , (5.6) -14 -

JCAP04(2015)012
which is the mass that saturates the perturbative unitary bound of (5.5) and that is still in agreement with the astrophysical lower bound of eq. (5.1).

3 → 2 dark matter scenario
Instead of the usual 2 → 2 annihilation scenarios described in previous sections, another mechanism for the reduction of the comoving DM number density is the annihilation via 3 → 2 processes. Here we discuss this case, assuming that it is the dominant process for the relic density computation, or equivalently that γ self , γ semi ≪ γ 3→2 at the epoch of the freeze-out. As explained in the appendix A, it is necessary for the DM to release the kinetic energy produced in the 3 → 2 reactions. This can be fulfilled if DM is in kinetic equilibrium with either SM particles or the GBs. The first possibility is realized by means of DM scattering off electrons. However, it requires a large Higgs portal. This can be understood from the effective Lagrangian between electrons and DM. Assuming m h , m ρ ≫ m X and θ ≪ 1, we have and one can write where m e is the electron mass. In order to keep the kinetic equilibrium between the DM and the electrons from this interaction, it is necessary that ǫ 5 × 10 −9 [25], or equivalently, that λ HX 1.5 × 10 3 (m e /m X ), which is too large if eq. (5.6) holds. As a result, the only possibility is to have kinetic equilibrium between the DM with the GBs. The corresponding rate depends on λ φX and for m h , m ρ ≫ m X is mostly driven by a contact interaction (see figure 6). Because of this, it is not suppressed by a heavy scale in contrast to the previous case. Moreover, although this scattering rate should be greater than the Hubble parameter at the epoch of the freeze-out, we must insure that the rate for the selfannihilation process into GBs -which is obtained by crossing symmetry-is negligible. This is possible because Γ kin /Γ ann ∼ n η /n DM ≫ 1 at the freeze-out, due to the sub-dominance of n DM during the radiation-dominated epoch. We conclude that λ φX must be sizable but not too large in order to avoid self-annihilations into GBs.
In section 3.1, it was implicitly assumed that the GB did not heat up after its decoupling from the SM plasma. In fact, this is not the case here because DM is not in kinetic equilibrium with the SM and transfers all its entropy exclusively to the GBs when it annihilates via the 3 → 2 mechanism. As a result, eq. (3.1) must be modified accordingly in order to account for this effect and thus T d η must be larger than the muon mass so that the contribution to N eff is negligible. This always takes place here because we consider m h , m ρ ≫ m X ∼ 100 MeV, which implies that the processes establishing the equilibrium are very suppressed for temperatures around the muon mass.
With this in mind, we can consider the 3 → 2 reactions in this model. All of them are equivalent, up to CP conjugation, to eitherXX X → X X or X X X → XX. In figures 14 and 15 we show the corresponding diagrams at tree-level. The corresponding cross-sections can be calculated analytically if we work in the non-relativistic approximation (see appendix B). In fact, the resulting expressions can be further simplified in the case γ self ≪ γ 3→2 because of the following reason. As for self-annihilation, some of the 3 → 2 diagrams are proportional to λ HX and λ φX . In fact, this happens for those in the first and Figure 14. Tree level diagrams for the annihilation X X ↔ XXX. The first two rows contain processes that are mediated via the exchange of h, ρ and X whereas the last two rows involve solely DM particles.

JCAP04(2015)012
second rows of figures 14 and 15 simply because they involve exchange of ρ or h particles. In contrast, the diagrams of the third and fourth lines depend uniquely on different couplings, namely λ 3 and λ X . Consequently, the contribution of the first set of diagrams is presumably negligible, if self-annihilations are subdominant. The fact that m ρ and m h ≫ m X reinforces this argument. However, notice that for this to be true it is not necessary that λ φX vanishes. Under these assumptions we find where Z = λ 3 v φ /m X was already introduced. Because self-and semi-annihilations into GBs are always kinematically allowed, eq. (5.9) must be compared with the cross-sections associated to these processes We already discussed the first process. Regarding the second one, although σv semi and σv 2 3→2 have different mass dimensions and can not be compared directly, it is clear that the semi-annihilations are suppressed for small λ 3 , and that the 3 → 2 is important when Z is large. Moreover, a suppression in the 3 → 2 process due to λ 3 can be compensated with a large v φ with respect to the DM mass. Under these circumstances we can make sure that γ self , γ semi ≪ γ 3→2 , while still having kinetic equilibrium between the GBs and the DM. We remark that the DM mass is not required to be in the GeV range because we assume no equilibrium between the GBs and SM plasma and consequently there is no significant contribution to N eff . In order to study the regions of the parameter space where the 3 → 2 reproduces the observed relic abundance, instead of solving numerically the Boltzmann equations directly we work in the freeze-out approximation.

The 3 → 2 freeze-out approximation
When the temperature of the plasma is much greater than the DM mass, the Boltzmann eq. (2.9) has the solution n = n eq , which corresponds to the chemical equilibrium. However, once the temperature drops below the DM mass, the rate for the process 2 → 3 is kinematically suppressed compared to the 3 → 2 reaction. As a result the DM number density deviates from its equilibrium value. In fact, it eventually becomes much greater than its equilibrium value n ≫ n eq . Hence, for these temperatures and when there are no self-and semi-annihilations, eq. (2.9) can be approximated by When the DM particles are non-relativistic, the cross-section σv 2 is independent of the temperature. In that case, using standard methods, eq. (5.11) admits the following solution where FO stands for the freeze-out, g are the number of relativistic degrees of freedom and x = m X /T . Here T is the DM temperature which equals that of the GB plasma. It is remarkable that, in contrast to the standard self-annihilation scenario, here the cross-section -17 -

JCAP04(2015)012
that matches the observed relic density depends on the DM mass. In order to estimate x FO , it is necessary to establish when the annihilation rate per particle n 2 eq σv 3→2 drops below the expansion rate of the Universe. Using eq. (5.12) it is found that this happens when freeze-out temperature satisfies x FO = 20.6 + log m X 100 MeV Because the second term of this equation depends logarithmically on the DM mass, it is a very good approximation to take x FO ≈ 20. It is somewhat interesting to calculate the value of σv 2→3 associated to eq. (5.12). To that end, notice that at the freeze-out the expansion rate is H = n eq σv 2→3 = n 2 eq σv 2 3→2 and as a result n eq = H/ σv 2 3→2 . According to eq. (A.9) we obtain σv 2→3 = H σv 2 3→2 = 1.08 × 10 −9 GeV −2 x FO g −0.5 FO ≈ 5.9 × 10 −26 cm 3 /s . (5.14) In contrast to eq. (5.12), notice that this equation is independent of the DM mass and is similar to the usual condition for the WIMP scenario.
We use this formalism in order to calculate the relic density. In figure 16 we show the regions of the plane [λ X , Z] in agreement with the observed relic density for different values of the DM mass when the 3 → 2 process dominates the DM production. They are obtained using eqs. (5.9) and (5.12), taking into account an uncertainty of 25% in order to account for the error in the freeze-out and the non-relativistic approximations. As apparent from the figure, there are two regimes. On the one hand, for small values of λ X the DM abundance is set by Z. On the other hand, when λ X is of order one both variables are relevant for the DM dynamics. In particular notice that Z can not vanish. Figure 16 also includes the region where perturbative unitarity is lost, which implies an upper limit on the DM mass as shown in the figure: m X 115 MeV .
(5.15) Let us note that this limit, obtained from the relic density, is comparable to the one in eq. (5.6) coming from demanding significant DM self-interactions.

Self-interactions and the 3 → 2 mechanism
It is interesting that the variables that control DM self-interactions also determine the 3 → 2 process, as shown in eqs. (5.3) and (5.9). This suggests that both phenomena have a common origin, which we attribute to the fact that DM is strongly interacting. Figure 17 illustrates this by showing the regions of the parameter space where the relic abundance is generated via the 3 → 2 mechanism and where the self-interaction cross-section simultaneously fulfills the astrophysical hints from small scale structures, as in eq. (5.1). In addition, we show the areas disfavored by the Bullet Cluster (i.e. eq. (5.2)). All this has been done under the assumption that perturbative unitarity is valid. The figure conclusively shows that the DM mass is bounded both from above and below 7 MeV m X 115 MeV , (5.16) which is in agreement with the findings of [25]. Similarly, for the Z parameter we find   We illustrate this scenario with three benchmark points which are defined in table 1 and shown in figure 17. They correspond to different masses in the MeV range, which reproduce simultaneously the observed relic abundance via the 3 → 2 mechanism and the hints on the self-interaction cross-section. Whereas points B and C give rise to a σ scatter /m X in agreement with the Bullet Cluster observations, point A is disfavored due to its small mass. In every case λ 3 , λ HX and λ φX are chosen small enough in order to ensure that semi-annihilations and self-annihilations are subdominant (see table 1). Nevertheless, we require that DM is in kinetic equilibrium with the GBs via the λ φX coupling, as argued before. We again find that the scale at which the U(1) DM symmetry breaks down has to be at least in the TeV range σ scatter /m X [cm 2 /g] 6.9 0.91 0.27 Table 1. Benchmark points that generate the measured DM relic abundance mainly via the 3 → 2 processes. Here sin θ = 10 −5 and λ HX = 10 −4 . Similarly, in order to keep the DM in kinetic equilibrium until its freeze-out we set λ φX = 10 −6 . m ρ and v φ are derived quantities.
so that Z is of order one. We would like to remark that all these points are in agreement with collider constraints because the ρ boson is very heavy and the invisible decay width of the Higgs is very small. In fact, the most important contribution to the latter comes from the process h → XX with a branching ratio of approximately 1.5 · 10 −5 . The other channels (ρρ, ηη, XXX andXXX) contribute less than ∼ 10 −12 .

Conclusions
We have presented a scalar extension of the Standard Model (SM) consisting of two additional complex fields with charges one and three under a global U(1) DM symmetry. The nonzero vacuum expectation value of the latter field spontaneously breaks the global symmetry down to a remnant Z 3 . This leads to the appearance of a Goldstone boson (GB), two Higgs-like particles, one of which is identified with the SM scalar, and a complex field that, due to the discrete symmetry, can not decay and is therefore a candidate for dark matter (DM). This model is constrained by different observations. In first place, the enlargement of the scalar sector could potentially modify the dynamics of the recently found SM scalar. In particular, its decay into invisible particles due to opening of new channels involving the GB and the DM. Secondly, the exchange of the Higgs-like particles can mediate elastic scattering between the DM and SM particles leading to detectable signals in direct detection experiments for DM. Lastly, GB can mimic neutrinos in the Cosmic Microwave Background, significantly altering the so-called effective number of neutrinos N eff . We have found that these constraints are always satisfied if the Higgs-like particles are heavier than ∼ 5 GeV, and if their mixing angle is sufficiently small, namely smaller than about 10 −5 .

JCAP04(2015)012
In these regions of the parameter space we have investigated the thermal production of DM. On the one hand, DM can self-annihilate either into SM particles by means of the Higgs portal, or into the non-DM additional scalars via the GB portal. Similarly, DM can semi-annihilate, that is two DM particles can be converted into another DM particle and a non-DM scalar. Because the GB is massless, both annihilation processes into GBs are always kinematically open and they therefore play an important role in the DM production. For DM in the GeV range, the typical values for the Higgs and GB portals as well as the coupling responsible for semi-annihilations are similar to the weak couplings of the SM. Because of this reason, when one of these two processes dominates, DM behaves as a weakly interacting massive particle (WIMP).
On the other hand, our model provides for a novel DM production mechanism, the so-called 3 → 2 mechanism consisting of the annihilation of three DM particles into two of them, which is possible due to the Z 3 symmetry. We find regions in the parameter where this mechanism is the dominant one. By studying the freeze-out approximation of the corresponding Boltzmann equation, and assuming the absence of non-perturbative effects, we have found that only DM masses in the MeV range can reproduce the observed relic abundance. In addition, for this mechanism to dominate over the other production modes, the Higgs and GB portals are required to be suppressed. Also, it is found that the U(1) DM symmetry must break down at least at the TeV scale. This naturally leads to relatively strong cubic DM self-interactions. Because of this, this scenario is a realization of the strongly interacting massive particle (SIMP) paradigm. Remarkably, the same mass range and similar DM cubic couplings naturally give rise to self-interacting cross-sections which could solve the discrepancies between observations and simulations of small-scale structures in the Universe, as shown in eq. (5.1). In fact, by performing a dedicated analysis of the parameter space we have found that in this scenario the DM mass must be in between 7 and 115 MeV. In order to illustrate our findings, we have studied three benchmark points (see table 1), compatible with all experimental constraints, in which the relic abundance is produced via the 3 → 2 mechanism and where the DM self-interactions are in the range suggested by astrophysics.
Direct and indirect detection experiments as well as collider searches will continue closing in on the parameter space of the WIMP scenario, due to the comparatively large couplings with the SM particles. Although this is more challenging for the SIMP scenario, astrophysics opens up a new window on the dynamics of DM self-interactions. This is the case, for instance, of the benchmark point A which is in tension with observations of the Bullet Cluster and recent analysis of the constraints of halo shapes. Consequently, new astrophysical observations can shed light on the SIMP scenario, potentially constraining further its parameter space.

A The Boltzmann equation
In this appendix we consider the Boltzmann equation for the DM particle and how to calculate the corresponding annihilation rates. Since we consider no asymmetry between X andX, we can think of the Z 3 charge as an internal degree of freedom for the DM. As a result the DM density is given by where f (E, t) is the DM phase-space distribution. The evolution of the DM density is given by the Boltzmann equation where H is the Hubble parameter, while N α and N β are the the number of DM particles in the initial and final states. In addition, the interaction rates are given bŷ where |M α→β | 2 is the square invariant amplitude including possible symmetry factors in both the initial and final state, respectively. Without loss of generality we assume N α > 1. For instance, in self-annihilation processes N α = 2 and N β = 0, in semi-annihilations N α = 2 and N β = 1, and for the 3 → 2 reactions N α = 3 and N β = 2.
We assume now that the non-DM particles are in chemical equilibrium at a common temperature T during the period of interest. Their phase-space distribution is thus given by f eq (E, T ) ≡ exp(−E/T ). Although this is not true for the DM particles, they are assumed to be in kinetic equilibrium, as discussed below. If this is the case, because of symmetry considerations in a Friedmann-Lemaître-Robertson-Walker Universe, their phase-space density f (E, t) is proportional to the chemical equilibrium distribution f eq (E, T ), with a proportionality factor independent of the three-momentum [63]. This assumption implieŝ where γ α→β corresponds to the rate of eq. (A.
3) when f i = f eq i . Furthermore, due to the principle of detailed balance and CP conservation γ α→β = γ β→α . Then, after adding over all the possible self-, semi-and 3 → 2 annihilation processes, the Boltzmann equation can be cast as dn dt + 3Hn = −2 n n eq 2 − 1 γ self − n n eq n n eq − 1 γ semi − n n eq 2 n n eq − 1 γ 3→2 (A. 5) or as dn dt + 3 H n = − σv self n 2 − n 2 eq − σv semi n 2 − n n eq − σv 2 3→2 n 3 − n 2 n eq , (A Kinetic equilibrium. When self or semi-annihilations are non-negligible, because of crossing symmetry, the transition amplitudes leading to them also induce scattering processes and these keep DM in kinetic equilibrium with its decaying products. In fact, when the DM particles become non-relativistic, the scattering processes occur at a faster rate than the annihilations because Γ kin Γ ann = n non-DM σv kin n DM σv ann ∼ n non-DM n DM ≫ 1 . (A. 10) where in the last step we take relativistic non-DM particles, which are more abundant during the radiation-dominated era. Because of this, the kinetic equilibrium is still maintained much after the DM freeze-out. In this work, when self and semi-annihilations can not be neglected, we further assume that all the non-DM particles are in equilibrium with each other. In particular, this means that the GB and the SM particles are in equilibrium. In order to avoid non-negligible contributions to N eff , as explained in section 3.1, such equilibrium can only exist at temperatures greater than the muon mass. Because of this and because the freeze-out temperature is roughly m X /25, in this work we only consider self-and semi-annihilations in the GeV range.
The situation is different in the opposite case, that is, when the 3 → 2 mechanism dominates. In this case, non-DM particles do not participate in the DM production, and therefore scattering processes between DM and GBs or SM particles do not necessarily take place. Nevertheless, the 3 → 2 processes convert a significant fraction of the mass into DM kinetic energy. If the latter is not radiated away, DM heats up altering significantly structure formation [64][65][66]. 5 We can circumvent this by allowing sufficiently large DM couplings to the SM plasma or the GBs, so that the former remains in equilibrium with one of the latter [25]. As shown in section 5, the first possibility requires huge couplings because the scattering rate is suppressed by the electron Yukawa coupling. Therefore, we assume the only possibility left, that is, that GBs and DM are in equilibrium. Because of eq. (A.10), this is not in contradiction with the fact that self-or semi-annihilations are negligible. Notice also that in this case the DM mass is not required to be in the GeV range because we assume no equilibrium between the GBs and SM plasma.

B The non-relativistic approximation
In order to calculate γ for the different processes with only DM particles in the initial state, we take f eq = exp(−E/T ) in eq. (A.3). For temperatures much smaller than the DM mass, because of the Boltzmann suppression, the integration over the phase-space of the initial state JCAP04(2015)012 Table 2. Square amplitudes in the non-relativistic limit. Here θ, λ HX and λ φX are neglected. We use the notation V = v φ mX and Z = λ 3 V.
is dominated by small three-momenta. In fact, by taking each of them equal to zero one finds where the subscript NR means that all the initial state particles are taken at rest. After performing the two-body final state phase-space integration, this implies that the cross-sections are

JCAP04(2015)012
Notice that for the self-and semi-annihilation cases, these expressions are just the s-save piece of the partial wave expansion of the thermal averaged cross-sections.
In table 2, we show all the square amplitudes that contribute to these cross-sections when λ HX , λ φX and θ are neglected. For the sake of completeness, we also include the square amplitude for other processes with three DM particles in the initial state: XXX → Xη and XXX → ηη. Notice that when the 3 → 2 mechanism is relevant -that is, when either Z or λ X is order one and λ 3 is small -those processes have a suppressed amplitude. We also show some processes with four DM particles in the initial state. Because λ X and Z are bounded from perturbative unitarity considerations as implied by eq. (5.4), those processes can not dominate over the 3 → 2 mechanism.