On thermal production of self-interacting dark matter

We consider thermal production mechanisms of self-interacting dark matter in models with gauged $Z_3$ symmetry. A complex scalar dark matter is stabilized by the $Z_3$, that is the remnant of a local dark $U(1)_d$. Light dark matter with large self-interaction can be produced from thermal freeze-out in the presence of SM-annihilation, SIMP and/or forbidden channels. We show that dark photon and/or dark Higgs should be relatively light for unitarity and then assist the thermal freeze-out. We identify the constraints on the parameter space of dark matter self-interaction and mass in cases that one or some of the channels are important in determining the relic density.


Introduction
Dark matter(DM) is the dominant component of matter in the Universe, and evidences for dark matter from galaxy rotation curves, gravitational lensing, and Cosmic Microwave Background have been getting more diverse and precise. For instance, the averaged relic density of dark matter is inferred by Planck data to be Ω DM h 2 = 0.1198 ± 0.0015 [1]. Weakly Interacting Massive Particles (WIMP) have been a well-motivated candidate for dark matter carrying weak interaction and weak-scale mass. The freeze-out mechanism for producing dark matter in the early Universe sets the annihilation cross section of dark matter to σv ann ∼ pb · c, which enables us to take different approaches to test the WIMP scenario. Therefore, there have been a lot of complementary efforts 1 for discovering WIMP dark matter, from direct detection, indirect detection and collider searches, but there have been null results until now. In particular, various direct detection experiments such as XENON100 [3] as well as LUX [4,5] and PandaX-II [6] have quite much constrained the WIMP-nucleon scattering cross section at the order of sub-zepto barn (10 −46 cm 2 ).
On the other hand, dark matter is assumed to be collisionless, namely, carry no selfinteractions, in Standard Cosmology, the so called ΛCDM. But, the numerical simulation with collisionless dark matter would lead to cuspy DM profiles that are not consistent with observed galaxies (core-cusp problem), as well as too many sub-halos (missing satellite problem) and too large masses for sub-halos (too-big-to-fail problem). These are the so called small-scale problems at galaxy scales [7,8]. Although the inclusion of baryons and supernova feedback in simulations might resolve such tensions in massive galaxies [9], the small-scale problems persist in the lowest mass galaxies where the discrepancy exists [10]. Therefore, small-scale problems may call for strong dark matter self-interactions, leading to σ self /m DM = 0.1 − 10 cm 2 /g or σ self /m DM ∼ barn for m DM ∼ 1 GeV.
Strongly Interacting Massive Particles(SIMP) [11] have recently drawn attention, due to the fact that the thermal freeze-out with 3 → 2 annihilation [12] allows for a large self-scattering of light dark matter. It is the Boltzmann suppression factor associated with an extra dark matter particle in the 3 → 2 process that naturally generates a hierarchy between the thermal annihilation cross section of about pico-barn and the self-scattering cross section of about barn [15]. However, the Boltzmann suppression factor is more or less fixed to e −x f with x f = m DM /T f at freeze-out temperature T f . Therefore, the relic density condition needs a relatively large self-interaction of SIMP dark matter, that is on the verge of violating unitarity or perturbativity and is in a tension with the bounds from Bullet cluster and halo shapes in most of the parameter space. There have been quite a few works in the literature for proposing concrete models to realize the SIMP [13][14][15][16][17][18][19] and its variations [20][21][22]. There have been growing interests in detecting the light dark matter of sub-GeV scale from direct detection [23][24][25][26][27] and cosmic rays [28].
In this article, we consider a complex scalar dark matter with gauged Z 3 symmetry in light of self-interacting dark matter [16]. The Z 3 coming from the spontaneous breaking of a local U (1) d stabilizes dark matter while the resultant dark photon and dark Higgs can contribute to the determination of the relic density. In particular, the semi-annihilation of dark matter [29] into heavier dark photon or dark Higgs, the so called forbidden channels [30,31], can be suppressed by a Boltzmann factor where m DM is the DM mass and m i is the mass of dark photon or dark Higgs. Then, taking ∆ i 1, it is possible to accommodate a smaller self-interaction of dark matter being compatible with the relic density, as compared to the SIMP case. Interestingly, the cubic self-coupling of dark matter in our model contributes to both SIMP and forbidden channels.
Furthermore, there exists a standard 2 → 2 annihilation of light dark matter into a pair of SM particles, in the presence of a Z portal coupling [32]. In this case, the smallness of the standard 2 → 2 annihilation (SM-annihilation) could be attributed to the smallness of the gauge kinetic mixing between U (1) d and hypercharge gauge group. We provide the general discussion on the thermal production of self-interacting dark matter in our model, in cases that one or some of SM-annihilation, SIMP and forbidden channels are relevant.
The paper is organized as follows. We begin with a review on the model with gauged Z 3 symmetry for dark matter and discuss the resultant mass spectrum and vacuum stability of the model. Then, we make general discussion on self-scattering and Boltzmann equation and kinetic equilibrium condition in the model. Next we study three mechanisms for thermal production, namely, SM-annihilation, SIMP and forbidden channels, in the presence of non-decoupled dark photon and dark Higgs, and discuss the constraints on the model from the relic density, self-scattering and various collider searches for light dark matter. Then, conclusions are drawn. There is one appendix dealing with the 2 → 2 forbidden channels involving dark photon and dark Higgs in our model.

Model for self-interacting dark matter
We consider dark matter as a complex scalar χ having a charge q χ = +1 under the dark local U (1) d symmetry, which is spontaneously broken to Z 3 by the VEV of another complex scalar φ with charge q φ = +3. Thus, the remaining discrete Z 3 symmetry 2 ensures the stability of scalar dark matter χ [16,33].
The Lagrangian for SM singlet scalars, χ, φ, and the SM Higgs doublet H, is given [16,33] by where the field strength tensor for dark photon is H, and the gauge kinetic mixing between dark photon V µ and hypercharge gauge boson B µ is introduced by sin ξ. Then, the dark photon communicates between dark matter and the SM particles through the gauge kinetic mixing. Here, the scalar potential is We note that the presence of a dark Higgs φ allows a triple coupling for χ after the U (1) d is spontaneously broken. Therefore, the corresponding κ coupling leads to SIMP processes as well as (forbidden) semi-annihilation processes of dark matter, which will be relevant for the later discussion.

Mass spectrum
For a nonzero VEV of dark Higgs field with φ = 1 √ 2 v d , the U (1) d symmetry is broken to a discrete subgroup Z 3 and dark photon gets massive and can mix with photon and Zboson. After expanding the dark Higgs as φ = 1 √ 2 (v d +h d ) and taking the SM Higgs doublet to be H T = 1 √ 2 (0, v ew + h), the dark Higgs can mix with the SM Higgs by Higgs-portal interaction, λ φH . Then, the SM and dark Higgs bosons are mixed [16] by where h 1 , h 2 are mass eigenstates. The mass eigenvalues of Higgs-like states are and the mixing angle is On the other hand, the effective mass of dark matter is given by m 2 χ,eff = m 2 χ + 1 2 λ λχ v 2 d + 1 2 λ χH v 2 ew , but we can absorb the contributions from symmetry breaking into the bare mass of dark matter. The details of interaction terms for dark/SM Higgses and dark photon can be found in Ref. [16].
Moreover, the mass eigenvalues of Z-boson and dark photon are where m 2 Z = 1 4 (g 2 + g 2 )v 2 and m 2 V = 9g 2 d v 2 d , and the mixing angle between Z-boson and dark photon is given by .
Then, taking ζ −s W ξ for m V m Z , we obtain the current interactions for dark photon as where ε ≡ c W ξ, and J µ EM and J µ d are electromagnetic, neutral and dark currents, respectively. In this case, we get m 2 ≈ 3g d v d ≡ m Z . See the appendix A of Ref. [16] for the details.

Vacuum stability
The absolute vacuum stability requires the potential to be bounded from below, meaning that V > 0 for large field values away from the local minimum with V = 0. In this section, for simplicity, we focus on the vacuum stability in the hidden sector with dark Higgs and dark matter scalars only. Although the mixing quartic couplings with the SM Higgs can affect our discussion too, they can be safely ignored, when they take positive or small values as compared to couplings in the hidden sector.
Taking φ = 1 √ 2 α and χ = 1 √ 2 β e iγ for large field values, the vacuum stability is determined by the quartic couplings in V DM , which becomes After minimizing the potential for γ along any field values with α = 0 and β = 0, the above hidden sector potential becomes Therefore, the vacuum stability conditions are given by and f (X min ) > 0 with where X min is the global minimum satisfying f (X min ) = 0. Then, solving f (X min ) = 0, we obtain the third condition (15) as with where 24λ φ and Q ≡ λ φχ 6λ φ . For instance, for κ = 0, the third vacuum stability condition (15) becomes trivial for λ φχ > 0, while it is given by 4λ φ λ χ − λ 2 φχ > 0 for λ φχ < 0, which is the standard result for two scalar fields with a mixing quartic coupling. On the other hand, for λ φχ = 0, the third vacuum stability condition (15) The vacuum stability condition equivalent to eq. (15) can be also derived by the condition that there is no real solution to the quartic polynominal f (X), in a more explicit form [35], and Then, the conditions, (15), (17) and (18), turn out to be equivalent.
For negative Higgs mixing quartic couplings with λ φH < 0 and λ χH < 0, there are corresponding vacuum stability conditions for them too. But, in the later analysis, we assume λ φH , λ χH to be positive if they are nonzero, so there is no extra conditions for vacuum stability. The general discussion on the vacuum stability conditions with arbitrary λ φH and λ χH are given in Ref. [35].
For the later sections, we will impose the vacuum stability conditions, (12) and (18), for the consistency of the vacuum breaking the U (1) d .

Dynamics of self-interacting dark matter
We discuss the self-scattering of dark matter and resulting unitarity bounds and present the general Boltzmann equation for the early Universe in our model. Then, we comment on the kinetic equilibrium condition for dark matter and the elastic scattering between light dark matter and electron.

Dark matter self-scattering and unitarity bounds
The squared amplitude for the χχ → χχ self-scattering is given [16] by On the other hand, the squared amplitude for the χχ * → χχ * self-scattering is given [16] by Therefore, in the non-relativistic limit for dark matter, the effective scattering cross section, The perturbativity and unitarity bounds on the DM couplings are given as follows, In the later sections, we will impose the above unitarity and perturbativity conditions for the consistency of the model.

General Boltzmann equation
Assuming CP conservation in the dark sector, we obtain the general Boltzmann equation for dark matter number density in our model, n DM = n χ + n χ * , with n χ = n χ * , as In principle, three annihilation processes, SM-annihilation, SIMP and forbidden channels can contribute equally in determining the number density of dark matter. In the next sections, we discuss the cases where one or some of annihilation processes become dominant.
In particular, in order to make the model unitarity up to relatively large masses for selfinteracting dark matter, it is necessary to introduce relatively light dark photon and/or dark Higgs so forbidden channels can be important too.

Kinetic equilibrium and DM detection
We assume that dark matter keeps in kinetic equilibrium during the freeze-out process, meaning that n SM σv χ,SM > H, where n SM is the equilibrium number density of the SM particles and σv χ,SM is the scattering cross section between dark matter and the SM particles in thermal bath. Then, we require a nonzero coupling between dark matter and the SM particles. To that purpose, Higgs portal or Z portal interactions in our model would be appropriate. It turns out that Higgs portal could not be used for kinetic equilibrium of sub-GeV light dark matter, because of small Yukawa couplings. We note that there are other possibilities that can be also consistent with observations, if dark matter is in kinetic equilibrium with dark radiation, namely, dark photon in our case. If dark matter were decoupled from both the SM and dark radiation, the dark sector could undergo an epoch of heating [12] so it would be unacceptable for structure formation.
In the later discussion on the SM-annihilating dark matter, a minimum value of the gauge kinetic mixing is needed for the correct relic density. Then, one has to take into account the bounds from direct detection as well as Z searches at colliders.
For Z portal interaction, the kinetic scattering cross section for χf → χf with f being the SM leptons is, in the early Universe, when leptons carries about the DM momentum, given [16] by Due to cross symmetry, a nonzero kinetic scattering cross section leads to the annihilation of dark matter into ff . For sub-GeV dark matter annihilating into leptons, the X-ray and gamma-ray searches can impose strong bounds on the corresponding annihilation cross section [28]. But, in our case, as will be shown in the next section, the annihilation cross section is velocity-suppressed, so there is no limit from indirect detection [16].
Similarly, for m e , m χ , m Z p m χ v DM at present, the DM-electron elastic scattering cross section with Z -portal interaction, that is relevant for direct detection, is given [16] by where µ ≡ m e m χ /(m e + m χ ) is the reduced mass of the DM-electron system. In the later section, we will show the parameter space that could be accessible by direct detection with semi-conductor or superconductor detectors [24,25]. The region for m χ vs ε, that is consistent with the SIMP dark matter, has been also shown to be constrained by direct detection and Z searches [16].

Thermal freeze-out from allowed channels
We discuss the thermal production of light dark matter from the 2 → 2 annihilations into a pair of SM particles and the 3 → 2 annihilations due to DM self-interactions.

SM-annihilating dark matter
For m h 1 , m Z m χ , the 2 → 2 annihilation channels are kinematically forbidden. Furthermore, for small self-couplings of dark matter, the 3 → 2 annihilation processes are also suppressed, namely, n 2 DM σv 2 3→2 < n DM σv 2→2 or n 2 DM σv 2 3→2 < H. In this case, dark matter annihilates dominantly into a pair of the SM particles.
As a result, the Boltzmann equation (23) is approximated to where σv 2→2 ≡ 1 2 σv χχ * →f f , which is given [16], before thermal average, by with v ew = 246 GeV, and We note that the Z -portal contribution in the first line of (27) is p-wave suppressed while the Higgs-portal contribution in the second line of (27) is suppressed by lepton Yukawa couplings. Thus, the model is not constrained by gamma-ray searches from the galactic center [28] or CMB constraints at recombination [1]. The SM-annihilating process with Z -portal interaction is still relevant for producing a right relic density from freeze-out.  Consequently, for 1 2 (σv) χχ * →ff = a + bv 2 , we get the relic density as This is the standard formula for the relic density in the case of SM annihilation, except that DM mass is taken to be sub-GeV.
In Fig. 1, we have shown the parameter space for m Z vs ε in (red) solid lines for dark matter m χ = 150(300) MeV on left (right) and dark gauge coupling g d = 1, 5, 10, being consistent with the relic density. Electron g − 2 limit and muon g − 2 favored region are shown in yellow and orange colors while the bound from EWPT is given in green. Monophoton + MET bounds from BaBar(improved) and Belle2(expected) [36] are shown in blue region and pink, light-blue and black dashed lines. The contour with elastic scattering cross section between dark matter and electron being given by σ χe = 10 −40 cm 2 are also shown in dotted lines. We find that the region that is consistent with the relic density can be probed by semi-conductor or superconductor detectors [24,25].
Since light dark matter annihilates into light fermions such as muons and electrons, Higgs portal interactions are Yukawa-suppressed, so they give negligible contributions to the DM annihilation. Nonetheless, non-negligible mixing quartic couplings λ χH and λ φH , or Higgs mixing angle, would lead to additional Higgs decay modes with decay rates given by, Then, additional Higgs couplings are bounded by Higgs data of signal strengths and/or searches for Higgs invisible decays at the LHC. The combined VBF, ZH and gluon fusion productions of Higgs boson at CMS lead to the bound, BR(h 2 → χχ * ) < 0.24 at 95% CL [37], while the bounds from the VBF [38] and ZH [39] Higgs productions at ATLAS are BR(h 2 → χχ * ) < 0.29 and BR(h 2 → χχ * ) < 0.75, respectively. As a result, the bound on the Higgs invisible decay leads to |y h 2 χ * χ |/v 0.010. On the other hand, the Higgs signal strength is bounded to µ > 0.81 at 95% CL from ATLAS/CMS data combined [40]. Thus, the Higgs mixing angle is bounded as sin θ < 0.44, which satisfied in our case, because sin θ λ φH v ew v d /m 2

SIMP dark matter
For m h 1 , m Z m χ and small couplings between messenger fields and the SM particles, all the 2 → 2 annihilation channels are kinematically forbidden or small. Then, the 3 → 2 annihilation process for dark matter becomes dominant, namely, n 2 DM σv 2 3→2 > n DM σv 2→2 or H > n DM σv 2→2 .
Consequently, ignoring the 2 → 2 annihilation processes, the Boltzmann equation (23) is approximated to The squared amplitude for χχχ * → χ * χ * scattering is, in the non-relativistic limit, given [16] by Likewise, the squared amplitude for χχχ → χχ * scattering is, in the non-relativistic limit, given [16] by Then, the effective 3-to-2 annihilation cross section appearing in the above Boltzmann equation is obtained as Therefore, the correct relic density fixes the ratio, m χ /α eff , which directly predicts the self-scattering cross section, σ self ∼ α 2 eff /m 2 χ . In Fig. 2, we have solved the relic density condition for λ χ and identified the parameter space of m χ vs R (DM cubic coupling), that is excluded by unitarity (in blue), perturbativity (in red) and vacuum stability (in green). Contours with self-scattering cross section with σ self /m χ = 0.1, 1, 10 cm 2 /g are shown in black dotted, dashed and dot-dashed lines, respectively. We have set g d = 0.1, λ φχ = 0.4, ∆ Z = 4 and ∆ h 1 = 0.5 where ∆ i ≡ (m i − m χ )/m χ . We find that the newly included vacuum stability bound is less severe than unitarity bound. The dark matter masses are bounded to be smaller than 150 MeV, due to perturbativity and unitarity.
On the other hand, in Fig. 3, we also drew the parameter space of ∆ Z vs ∆ h 1 on left and ∆ Z vs g d on right, that are excluded by unitarity (in blue) and vacuum stability Figure 4: Feynmann diagrams for forbidden channels with Z .
(in green). The black dotted, dashed and dot-dashed lines correspond to contours with self-scattering cross section as in Fig. 2, except that the dot-dashed line on right is for σ self /m χ = 2 cm 2 /g. As a consequence, in the allowed parameter space, dark Higgs mass is close to dark matter mass while dark photon mass can be much heavier than dark matter mass. Thus, for relatively light dark Higgs, the forbidden channels, χχ * → h 1 h 1 and χχ → h 1 χ * , can be also important in determining the relic density, as will be shown in the later sections. From the right plot in Fig. 2, the self-scattering cross section can be large, being insensitive to the choice of dark gauge coupling, as far as dark photon mass is large as well.

Thermal freeze-out from forbidden channels
When dark photon and/or dark Higgs boson masses are close to dark matter mass, they can contribute to the relic density through the forbidden channels, provided that the corresponding 2 → 2 cross sections are large enough. Therefore, we still allow for a large self-scattering of dark matter. We study the parameter space that is consistent with the relic density, first in the case with light dark photon, then the case with light dark Higgs and finally the case where both dark photon and dark Higgs are light. Here, we assume that dark photon and/or dark Higgs boson are in kinetic equilibrium during the freeze-out process.

The case with
In the case where m χ < m Z and dark Higgs is much heavier than the other particles, the forbidden channels involving Z as shown in Fig. 4 contribute to the Boltzmann equation.  Then, the Boltzmann equation (23) is approximated to The detailed balance conditions at high temperature are and σv χχ→Z χ * = 2n eq Z n eq DM σv Z χ * →χχ Then, we can rewrite the Boltzmann equation by using the detailed balance conditions, (39) and (40), as follows, Relic, Δ Z' =0.5  Figure 6: Parameter space of R vs m χ for forbidden channels with Z . The red lines satisfy the relic density and the blue region is excluded by unitarity. Dotted, dashed and dotdashed lines correspond to self-scattering cross sections, σ self /m χ = 0.001, 0.01, 0.1 cm 2 /g. We took ∆ Z = 0.4 or 0.5 and the value of λ χ saturates the vacuum stability bound, and g d = 0.1, λ φχ = 0.4 and ∆ h 1 = 10.
where λ ≡ s(m χ )/H(m χ ) with s(m χ ) = 2π 2 45 g * s m 3 χ and 1/H(m χ ) = 3.02g . As a result, the approximate solution to the Boltzmann equation (41) is given by Setting (σv) Z Z →χχ * = a and (σv) Z χ * →χχ = bv 2 , which leads to σv Z Z →χχ * = a and σv Z χ * →χχ = 6b/x, where the detailed expressions for a and b are given in eqs. (A.1) and (A.2), we get the DM abundance as with Consequently, the relic density is determined to be Ω DM h 2 = 5.20 × 10 −10 GeV −2 g * 10.75 Then, the 2 → 2 annihiation cross sections can be large, due to the inverse of the Boltzmann suppression factor, e ∆ Z x f , appearing in the relic density. Therefore, the self-scattering cross section of dark matter can be large enough. But, for a small ∆ Z , dark matter selfinteraction can be smaller than in the SIMP case, being compatible with the relic density.
In Fig. 5, we depicted the relic density as a function of ∆ Z , only with forbidden channels involving Z on left and with both SIMP and forbidden channels on right, varying the self-interaction and mass of dark matter, (R, m χ ), between R = 0.1 − 1 and m χ = 10 MeV − 1 GeV. For both plots, we took g d = 0.1, λ φχ = 0.4 and ∆ h 1 = 10. On the right plot of Fig. 5, as ∆ Z gets larger than about 0.5 the relic density approaches a certain fixed value, that is determined by the SIMP processes dominantly. However, for ∆ Z 0.5, the relic density become sensitive to the value of ∆ Z , as well as to R and m χ .
In Fig. 6, we also show the parameter space of m χ vs R, satisfying the relic density in red lines, for ∆ Z = 0.4 and 0.5, from bottom to top. We took dark matter masses to be larger than 150 MeV to cover beyond the maximal value allowed by unitarity in the SIMP case. Furthermore, we chose the value of λ χ such that the vacuum stability bound is saturated and set g d = 0.1, λ φχ = 0.4 and ∆ h 1 = 10. The blue region is excluded by unitarity and the contours with self-scattering cross section, σ self /m χ = 0.001, 0.01, 0.1 cm 2 /g, are shown in dotted, dashed and dot-dashed lines, respectively. For relatively heavy DM masses, the self-scattering cross section is smaller than in the SIMP case, being consistent with the relic density and unitarity. Therefore, the forbidden channels are crucial to keep the model perturbative for the wide range of masses for light dark matter.

The case with
In the case where m χ < m h 1 and dark photon is much heavier than the other particles, the forbidden channels involving h 1 as shown in Fig. 7 contribute to determining the relic density. In this case, the Boltzmann equation (23) is approximated to Similarly to the case with light Z , the detailed balance conditions at high temperature are Therefore, we can rewrite the Boltzmann equation by using the detailed balance conditions, (47) and (48), as follows, As in the case with Z channels, the approximate solution to the above Boltzmann equation is then given by with Consequently, the relic density is determined to be Ω DM h 2 = 5.20 × 10 −10 GeV −2 g * 10.75 In Fig. 8, we depicted the relic density as a function of ∆ h 1 , only with forbidden channels involving h 1 on left and with both SIMP an forbidden channels on right, varying (R, m χ ) such that R = 0.1 − 1 and m χ = 10 MeV − 1 GeV. For both plots, we took g d = 0.5, λ φχ = 0.4 and ∆ Z = 10. Inclusion of the SIMP processes on the right plot clearly shows a resonance behavior at ∆ h 1 ∼ 1 or m h 1 ∼ 2m χ , drastically changing the relic density to much smaller values. But, the same resonance also appears in the self-scattering of dark matter as can be seen in eq. (20), so it would be in a tension with the bound from Bullet cluster. Below the resonance region on the right plot, between ∆ h 1 ∼ 0.5−1, there appears a similar plateau with a fixed relic density, that is dominantly determined by the SIMP processes.
In Fig. 9, we also show the parameter space of m χ vs R, satisfying the relic density in red lines, for ∆ h 1 = 0.75 and 0.85, from bottom to top. We chose the value of λ χ such that the vacuum stability bound is saturated and set g d = 0.5, λ φχ = 0.4 and ∆ Z = 10. The blue region is excluded by unitarity and the contours with self-scattering cross section are shown similarly to those in Fig. 6.

5.3
The case with m χ < m Z ∼ m h 1 When Z and dark Higgs are comparably light, they both can contribute comparably to the forbidden channels at the same time. In this case, from eqs. (41) and (49), we obtain the approximate Boltzmann equation (23) as  Therefore, the DM relic abundance becomes In this case, the relic density is given by Ω DM h 2 = 5.20 × 10 −10 GeV −2 g * 10.75 In Fig. 10, we depicted the relic density as a function of ∆ Z = ∆ h 1 , when dark photon and dark Higgs are degenerate in mass, for varying (R, m χ ) between R = 0.1 − 1 and m χ = 10 MeV − 1 GeV. For both plots, we took g d = 0.3 and λ φχ = 0.4. As in the case with light Z or h 1 in the previous subsections, there is a similar dependence on ∆ Z as well as (R, m χ ).
In Fig. 11, we showed the parameter space for m χ vs R, that explains the observed relic density in red lines, for ∆ Z = ∆ h 1 = 0.8 and 0.9, from bottom to top. We chose the value of λ χ such that the vacuum stability bound is saturated and set g d = 0.3 and λ φχ = 0.4. The blue region is excluded by unitarity and the contours with self-scattering cross section are shown similarly to those in Fig. 6.
Relic, Δ Z' =Δ h1 =0.9  Figure 11: Parameter space of R vs m χ for forbidden channels with Z and h 1 . The red lines satisfy the relic density and the blue region is excluded by unitarity. Dotted, dashed and dot-dashed lines correspond to self-scattering cross sections, σ self /m χ = 0.001, 0.01, 0.1 cm 2 /g. We took ∆ Z = ∆ h 1 = 0.8 or 0.9 and the value of λ χ saturates the vacuum stability bound, and g d = 0.3, λ φχ = 0.4.

Conclusions
We have considered the thermal production of self-interacting dark matter in models with Z 3 gauged symmetry. We showed that standard 2 → 2 annihilation and hidden sector annihilations (3 → 2 annihilation and forbidden channels) can contribute equally in determining the relic density. In particular, dark photon and dark Higgs in the model must be kept light for unitarity, so they both can contribute to the processes of dark matter annihilation. In particular, we found that forbidden channels with semi-annihilation such as χχ → χ * Z or χχ → χ * h 1 assist a thermal production of light dark matter with larger masses than in the SIMP case, but keeping a sizable self-scattering of dark matter. Depending on the value of the self-scattering cross section favored by small-scale problems, we can identify the relevant thermal production mechanisms for self-interacting dark matter, in the same model.