Cosmic abundances of SIMP dark matter

TL;DR

This paper refines the calculation of thermal averages for 3-to-2 dark matter annihilation processes, crucial for accurately predicting relic densities of SIMP dark matter, especially near resonance poles.

Contribution

It provides a general method for thermal averaging of 3-to-2 and higher-order interactions, emphasizing the importance of correct averaging in models with velocity dependence and resonance effects.

Findings

Correct thermal averages significantly impact relic density predictions.

Resonance poles can alter the relic abundance calculations.

Application to SIMP models demonstrates the importance of precise thermal averaging.

Abstract

Thermal production of light dark matter with sub-GeV scale mass can be attributed to $3\rightarrow 2$ self-annihilation processes. We consider the thermal average for annihilation cross sections of dark matter at $3\rightarrow 2$ and general higher-order interactions. A correct thermal average for initial dark matter particles is important, in particular, for annihilation cross sections with overall velocity dependence and/or resonance poles. We apply our general results to benchmark models for SIMP dark matter and discuss the effects of the resonance pole in determining the relic density.

Tables

Table 1. Table 1: Coefficients of thermal averaged velocity terms.

$(n,m,l)$	$(1,0,0)$	$(1,1,0)$	$(2,0,0)$	$(1,1,1)$	$(2,1,0)$	$(3,0,0)$
$c_{nml}$	$2$	$\frac{14}{3}$	$\frac{20}{3}$	$\frac{100}{9}$	$\frac{160}{9}$	$\frac{280}{9}$

Equations

$$⟨\sigma v⟩=\frac{\int d^3{v}_1{d}^3{v}_2{\delta }^3({\vec{v}}_1+{\vec{v}}_2)(\sigma v)e^{-\frac{1}{2}x(v_1^2+v_2^2)}}{\int d^3{v}_1{d}^3{v}_2{\delta }^3({\vec{v}}_1+{\vec{v}}_2)e^{-\frac{1}{2}x(v_1^2+v_2^2)}},$$

$$⟨\sigma v⟩=\frac{x^{3/2}}{2\sqrt{\pi }}{\int }_0^{\infty }dv{v}^2(\sigma v)e^{-\frac{1}{4}x{v}^2}.$$

$$(\sigma v)=\sum _{l=0}^{\infty }\frac{a_l}{l!}(v^2{)}^l.$$

$$⟨\sigma v⟩$$

$$(\sigma v{)}_R$$

$$⟨\sigma v{⟩}_R=2x^{3/2}\sqrt{\pi }{\gamma }_R\sum _{l=0}^{\infty }\frac{b_l}{l!}{F}_l(z_R;x)$$

$$F_l(z_R;x)$$

$$F_0(z_R;x)$$

$$\mathrm{Erfc}(a)\equiv \frac{2}{\sqrt{\pi }}{\int }_a^{\infty }{e}^{-t^2}dt.$$

$$⟨\sigma v{⟩}_R\approx 2\sqrt{\pi }{\gamma }_R{ϵ}_R^{1/2}{x}^{3/2}{e}^{-x{ϵ}_R}\theta ({ϵ}_R)\sum _{l=0}^{\infty }\frac{b_l}{l!}{ϵ}_R^l.$$

$$⟨\sigma v^2⟩=\frac{\int d^3{v}_1{d}^3{v}_2{d}^3{v}_3{\delta }^3({\vec{v}}_1+{\vec{v}}_2+{\vec{v}}_3)(\sigma v^2)e^{-\frac{1}{2}x(v_1^2+v_2^2+v_3^2)}}{\int d^3{v}_1{d}^3{v}_2{d}^3{v}_3{\delta }^3({\vec{v}}_1+{\vec{v}}_2+{\vec{v}}_3)e^{-\frac{1}{2}x(v_1^2+v_2^2+v_3^2)}}.$$

$$(\sigma v^2)=a_0+a_1(v_1^2+v_2^2+v_3^2)+a_2^{(1)}(v_1^2+v_2^2+v_3^2{)}^2+a_2^{(2)}(v_1^4+v_2^4+v_3^4)+\cdots .$$

$$(\sigma v^2)=\sum _{l=0}^{\infty }\frac{a_l}{l!}{\eta }^l$$

$$⟨\sigma v^2⟩$$

$$⟨(v_1^2{)}^n(v_2^2{)}^m(v_3^2{)}^l⟩$$

$$\displaystyle\langle(v^{2}_{1})^{n}(v^{2}_{2})^{m}\rangle=\frac{3\sqrt{3}x^{2}}{\pi}\,\int^{\infty}_{0}dv_{1}\,v^{2n+1}e^{-\frac{3}{4}xv^{2}_{1}}\int^{\infty}_{-\infty}dv^{\prime}_{2}\Big{(}v^{\prime}_{2}+\frac{1}{2}v_{1}\Big{)}^{2m+1}e^{-xv^{\prime 2}_{2}}.$$

$$\displaystyle\langle(v^{2}_{1})^{n}\rangle=\Big{(}\frac{4}{3}\Big{)}^{n}\frac{\Gamma(n+\frac{3}{2})}{\Gamma(\frac{3}{2})}\,x^{-n}=\langle(v^{2}_{2})^{n}\rangle=\langle(v^{2}_{3})^{n}\rangle.$$

$$⟨\sigma v^2{⟩}_{d-\mathrm{wave}}$$

$$(\sigma v^2{)}_R$$

$$⟨\sigma v^2{⟩}_R=\frac{3}{4}\pi x^3\sum _{l=0}^{\infty }\frac{b_R^{(l)}}{l!}{G}_l(z_R;x),$$

$$G_l(z_R;x)$$

$$G_0(z_R;x)$$

$$\Gamma (0,a)\equiv {\int }_a^{\infty }\frac{e^{-t}}{t}dt.$$

$$\displaystyle\langle\sigma v^{2}\rangle_{R}\approx\frac{27}{16}\pi\epsilon^{2}_{R}x^{3}e^{-\frac{3}{2}x\epsilon_{R}}\theta(\epsilon_{R})\sum_{l=0}^{\infty}\frac{b^{(l)}_{R}}{l!}\,\Big{(}\frac{3}{2}\Big{)}^{l}\epsilon^{l}_{R}.$$

$$\frac{d{n}_{\mathrm{DM}}}{dt}+3H{n}_{\mathrm{DM}}=-⟨\sigma v⟩(n_{\mathrm{DM}}^2-(n_{\mathrm{DM}}^{\mathrm{eq}}{)}^2).$$

$$\displaystyle\frac{dY_{\rm DM}}{dx}=-\lambda x^{-2}\langle\sigma v\rangle\Big{(}Y^{2}_{\rm DM}-(Y^{\rm eq}_{\rm DM})^{2}\Big{)}$$

$$\displaystyle Y_{\rm DM}(\infty)\approx\bigg{(}\lambda J(x_{f})\bigg{)}^{-1}.$$

$$J(x_f)\equiv {\int }_{x_f}^{\infty }dx{x}^{-2}⟨\sigma v⟩.$$

$$J(x_f)=\frac{a_l}{l+1}{x}_f^{-l-1}.$$

$${\Omega }_{\mathrm{WIMP}}{h}^2$$

Full text

{centering}

Cosmic abundances of SIMP dark matter

Soo-Min Choi†, Hyun Min Lee† and Min-Seok Seo‡

†Department of Physics, Chung-Ang University, Seoul 06974, Korea.

‡Center for Theoretical Physics of the Universe, Institute for Basic Science, 34051 Daejeon, Korea.

Thermal production of light dark matter with sub-GeV scale mass can be attributed to $3\rightarrow 2$ self-annihilation processes. We consider the thermal average for annihilation cross sections of dark matter at $3\rightarrow 2$ and general higher-order interactions. A correct thermal average for initial dark matter particles is important, in particular, for annihilation cross sections with overall velocity dependence and/or resonance poles. We apply our general results to benchmark models for SIMP dark matter and discuss the effects of the resonance pole in determining the relic density.

1 Introduction

Thermal dark matter, that was once in chemical equilibrium and decoupled from thermal plasma in the Universe, has been one of the plausible candidates for dark matter, in particular, under the name of Weakly Interacting Massive Particles (WIMP). Chemical equilibrium of dark matter usually requires a standard $2\rightarrow 2$ annihilation of dark matter into the SM particles, so it has provided an interesting interplay between the relic density, direct and indirect detection of dark matter at terrestrial and satellite experiments. Recently, a new mechanism for freezing out the density of dark matter from the $3\rightarrow 2$ annihilation process, coined the Strongly Interacting Massive Particles (SIMP) [1], has recently drawn special attention, due to the fact that there is no need of a large coupling between dark matter and the SM particles in this case.

Dark matter in the early Universe has once had a Maxwell-Boltzmann velocity distribution in the non-relativistic limit for which the DM annihilates only. Thus, there is a need of making a thermal average for the annihilation cross section of dark matter in order to incorporate it in the Boltzmann equation for the DM relic density. In particular, when the annihilation cross section depends strongly on the DM velocity, for instance, due to dominance of higher partial waves or resonance poles. In the case of WIMP dark matter that is based on the $2\rightarrow 2$ annihilation, it is enough to do the thermal average for the velocity of a single DM particle or the relative velocity in the center of mass frame. On the other hand, in the case of SIMP dark matter that is based on the $3\rightarrow 2$ annihilation, we need to do the thermal averages for two (relative) velocities of dark matter in the initial states. Given that the velocity dependence of the $3\rightarrow 2$ annihilation depends on the properties of dark matter [2, 3, 4, 5, 6, 8, 7] and the existence of resonance poles [8], it is worthwhile to make a systematic study of the thermal averages for $3\rightarrow 2$ and higher-order annihilation processes in general.

In this article, we present a general discussion on the thermal average of the $3\rightarrow 2$ annihilation cross section in the perturbative regime where the velocity expansion is valid and near the resonance pole that mediates between three particles in the initial state and two particles in the final state. We discuss the effects of the resonance pole on the thermal-averaged cross section as well as the relic density and compare the results to the WIMP case. Representative examples for SIMP dark matter, such as models with $Z_{n}$ discrete symmetries and dark mesons, are discussed in light of the thermal average of the $3\rightarrow 2$ annihilation cross section without or with a resonance pole.

The paper is organized as follows. We begin with a review on the thermal average of the $2\rightarrow 2$ annihilation cross section and then discuss a counterpart of the $3\rightarrow 2$ annihilation cross section without or with a resonance. Next we incorporate the thermal-averaged cross sections in the Boltzmann equations for WIMP and SIMP cases and apply our general results for known models for SIMP dark matter. We continue to generalize our discussion to the $3\rightarrow 2$ coannihilation between particles with different masses and higher-order annihilation processes. Finally, conclusions are drawn.

2 Thermal average for $2\rightarrow 2$ DM annihilations

To warm up and compare to our later discussion on $3\rightarrow 2$ processes, we first give a review on the thermal average of the standard $2\rightarrow 2$ annihilation cross section without or with a resonance. Assuming that two DM particles in the initial states have the same masses, $m_{1}=m_{2}\equiv m_{\rm DM}$, the thermal averaged $2\rightarrow 2$ cross section is given by

[TABLE]

where the momentum conservation is included as a delta function in the center of mass frame and $x\equiv\frac{m_{\rm DM}}{T}$ with $T$ being the DM temperature that is equal to the background temperature in kinetic equilibrium. In this case, the thermal average is simplified to the integral for relative velocity, $|{\vec{v}}_{1}-{\vec{v}}_{2}|\equiv v$, as follows,

[TABLE]

Suppose to take the velocity expansion of the $2\rightarrow 2$ cross section as

[TABLE]

Due to the absence of a resonance, we get the thermal average simply as

[TABLE]

Thus, we have recovered the well known results for the thermal-averaged $2\rightarrow 2$ annihilation cross section [10].

On the other hand, in the presence of a resonance $R$, the $2\rightarrow 2$ annihilation cross section for ${\chi\chi\rightarrow R\rightarrow f{\bar{f}}}$ takes the following Breit-Wigner form,

[TABLE]

where $\beta_{\chi}$ is the DM velocity, and $\eta\equiv\frac{1}{4}v^{2}$, $\epsilon_{R}\equiv\frac{m^{2}_{R}-4m^{2}_{\rm DM}}{4m^{2}_{\rm DM}}$ and $\gamma_{R}\equiv\frac{m_{R}\Gamma_{R}}{4m^{2}_{\rm DM}}$, with $m_{R},\Gamma_{R}$ being the mass and width of the resonance. Then, we obtain the general result for the thermal average with a resonance as follows,

[TABLE]

where $z_{R}\equiv\epsilon_{R}+i\gamma_{R}$ and

[TABLE]

Here, the generating integral is given by

[TABLE]

with the complementary error function being given by

[TABLE]

In particular, in the narrow width approximation with $\gamma_{R}\ll 1$, we get $F_{0}(z_{R};x)\approx\epsilon^{1/2}_{R}e^{-x\epsilon_{R}}\theta(\epsilon_{R})$ with $\theta(\epsilon_{R})=1$ for $\epsilon_{R}>0$; $\theta(\epsilon_{R})=0$ for $\epsilon_{R}<0$ and the thermal averaged cross section becomes

[TABLE]

Thus, the averaged annihilation cross section becomes a step function in the narrow width approximation, being sensitive to the resonance mass [10, 11].

In Fig. 1, we show the exact results for the averaged annihilation cross section with $s$-wave overall factor in arbitrary unit as a function of $\epsilon_{R}$ for a fixed $\gamma_{R}$ and temperature, $T=\frac{m_{\rm DM}}{15}$. In the limit of a narrow width, the averaged annihilation cross section is shown to be step-wise as in our approximate formula in eq. (10).

3 Thermal average for $3\rightarrow 2$ DM annihilations

Assuming that three DM particles in the initial states have the same masses, $m_{1}=m_{2}=m_{3}\equiv m_{\rm DM}$, the thermal averaged $3\rightarrow 2$ cross section is given by

[TABLE]

We assumed that the spins of dark matter are averaged and summed over initial and final states in $3\rightarrow 2$ processes. Then, the resulting velocity expansion of the $3\rightarrow 2$ cross section depends on the spin and parity of dark matter. For instance, in the case of fermionic SIMP, the initial states in the $3\rightarrow 2$ process can be all fermions as discussed in Ref. [15] while the case of vector SIMP was discussed [6] or will be published elsewhere [7].

In the non-relativistic limit of dark matter, taking into account the Galilean symmetry and permutation symmetry between three initial DM particles, we can take the velocity expansion of the $3\rightarrow 2$ cross section as follows,

[TABLE]

There appear more combinations of squared velocities at higher orders. We note that at the fourth order in velocities, an alternative basis can be choosen with $v^{2}_{1}v^{2}_{2}+v^{2}_{2}v^{2}_{3}+v^{2}_{3}v^{2}_{1}$ or $({\vec{v}}_{1}\cdot{\vec{v}}_{2})^{2}+({\vec{v}}_{2}\cdot{\vec{v}}_{3})^{2}+(\vec{v}_{3}\cdot{\vec{v}}_{1})^{2}$, instead of $v^{4}_{1}+v^{4}_{2}+v^{4}_{3}$, whenever it is more convenient for thermal average 111We note the following identities, $v^{2}_{1}v^{2}_{2}+v^{2}_{2}v^{2}_{3}+v^{2}_{3}v^{2}_{1}=\frac{1}{2}(v^{2}_{1}+v^{2}_{2}+v^{2}_{3})^{2}-\frac{1}{2}(v^{4}_{1}+v^{4}_{2}+v^{4}_{3})$, and $({\vec{v}}_{1}\cdot{\vec{v}}_{2})^{2}+({\vec{v}}_{2}\cdot{\vec{v}}_{3})^{2}+(\vec{v}_{3}\cdot{\vec{v}}_{1})^{2}=v^{4}_{1}+v^{4}_{2}+v^{4}_{3}-\frac{1}{4}(v^{2}_{1}+v^{2}_{2}+v^{2}_{3})^{2}$ due to ${\vec{v}}_{1}+{\vec{v}}_{2}+{\vec{v}}_{3}=0$..

3.1 Non-resonance

The thermal average of velocity terms, given by a function of $v_{1}^{2}+v^{2}_{2}+v^{2}_{3}$, namely in an $SO(9)$ symmetric form, can be easily computed in a closed form as below. Thus, we first treat them separately and next consider general terms of the form, $(v^{2}_{1})^{n}(v^{2}_{2})^{m}(v^{2}_{3})^{l}$.

First, we take the velocity expansion of the $3\rightarrow 2$ cross section in the following form with $SO(9)$ invariance,

[TABLE]

with $\eta\equiv\frac{1}{2}(v^{2}_{1}+v^{2}_{2}+v^{2}_{3})$. Then, the corresponding thermal average is given by

[TABLE]

In most cases, the most important terms appear up to $p$-wave terms that are $SO(9)$ invariant, so the above result gives rise to a good approximation for the full average. But, if the $3\rightarrow 2$ cross section is velocity-suppressed, we need to take into account the precise form of higher order terms in the velocity expansion.

There are cases where the leading terms in the velocity expansion are higher than $p$-wave, such as in the case with SIMP mesons which have leading $d$-wave terms. Thus, for more general velocity terms, we need to do the velocity integrations as

[TABLE]

where $c_{nml}$ are constant coefficients depending on $(n,m,l)$. In the case with $l=0$, the above integration can be simplified to

[TABLE]

These integrals can be calculated numerically and some of them with low $n,m,l$ are shown in Table 1. Other combinations with a fixed value of $n+m+l$ are not shown because they are the same as the one shown in Table 1 due to permutation symmetry between dark matter particles.

Instead, taking $m=l=0$, we can perform the integral in a closed form as

[TABLE]

In particular, using eqs. (14) and (17), we get the thermal average of $d$-wave terms as follows,

[TABLE]

In most of examples for $3\rightarrow 2$ processes such as SIMP mesons, it would be sufficient to consider at most the $d$-wave terms for thermal average.

3.2 Resonance

In the presence of resonances near the center of mass energy of three initial DM particles, more care is needed in the process of thermal average. In the non-relativistic limit of dark matter, the $3\rightarrow 2$ cross section for $\chi\chi\chi\rightarrow R\rightarrow\chi\chi$, before thermal average, takes a generalized Breit-Wigner form,

[TABLE]

where $\beta_{\chi}$ is the DM velocity in the two-body decay of the resonance, namely, $\beta_{\chi}\equiv\sqrt{1-4m^{2}_{\chi}/m^{2}_{R}}$, $\Phi_{3}$ is the phase space integral for the three-body decay of the resonance, $R\rightarrow\chi\chi\chi$, and $\epsilon_{R},\gamma_{R}$ are the counterparts for the $3\rightarrow 2$ resonance, given by $\epsilon_{R}\equiv\frac{m^{2}_{R}-9m^{2}_{\rm DM}}{9m^{2}_{\rm DM}}$ and $\gamma_{R}\equiv\frac{m_{R}\Gamma_{R}}{9m^{2}_{\rm DM}}$, with $m_{R},\Gamma_{R}$ being the mass and width of the resonance. We note that the three-body phase space integral $\Phi_{3}$ is proportional to $\epsilon^{2}_{R}m^{2}_{R}$ near resonance, so the three-body decay rate of the resonance is suppressed as compared to the two-body decay rate.

First, when the overall factor of the $3\rightarrow 2$ cross section is taken as a function of $\eta$ as $b_{R}=\sum_{l=0}^{\infty}\frac{b^{(l)}_{R}}{l!}\,\eta^{l}$, the resulting thermal average is given by

[TABLE]

where

[TABLE]

with $z_{R}\equiv\epsilon_{R}+i\gamma_{R}$. Here, the generating integral $G_{0}(z_{R};x)$ can be written in a closed form as follows,

[TABLE]

where the incomplete gamma function being is given by

[TABLE]

For narrow width approximation with $\gamma_{R}\ll 1$, we get $G_{0}(z_{R};x)\approx\frac{9}{4}\epsilon^{2}_{R}e^{-\frac{3}{2}x\epsilon_{R}}\theta(\epsilon_{R})$ and the thermal averaged cross section becomes

[TABLE]

We find that the averaged cross section in the SIMP case is more sensitive to the resonance mass through $\epsilon^{2}_{R}$ than in the WIMP case where the averaged cross section is proportional to $\epsilon^{1/2}_{R}$ in eq. (10). This is due to the fact that the phase space in the velocity average for three initial DM particles takes a higher power in DM velocity so it becomes more sensitive to the pole of the resonance.

In Fig. 2, we depict the analytic results for thermal-averaged $3\rightarrow 2$ annihilation cross section with $s$-wave overall factor in arbitrary unit as a function of $\epsilon_{R}$ for a fixed $\gamma_{R}$ and temperature, $T=\frac{m_{\rm DM}}{15}$. Similarly to the WIMP case, the result is sensitive to the mass of the resonance and it becomes step-wise in the limit of a narrow width.

4 Boltzmann equations for dark matter

We use the general results on thermal averages in the previous section to solve the Boltzmann equations for the relic density of WIMP or SIMP dark matter.

4.1 Boltzmann equation for WIMP

The Boltzmann equation for WIMP dark matter is given by

[TABLE]

Then, the above equation can be rewritten in terms of the relic abundance of dark matter, $Y_{\rm DM}=n_{\rm DM}/s$, as follows,

[TABLE]

where $\lambda\equiv s(m_{\rm DM})/H(m_{\rm DM})$ with $s(m_{\rm DM})=\frac{2\pi^{2}}{45}g_{*s}m^{3}_{\rm DM}$ and $1/H(m_{\rm DM})=3.02g^{-1/2}_{*}\frac{M_{P}}{m^{2}_{\rm DM}}$. Therefore, we obtain the solution to the Boltzmann equation as

[TABLE]

with

[TABLE]

Here, $x_{f}=m_{\rm DM}/T_{f}$ with $T_{f}$ being the freeze-out temperature. In the case without a resonance, when $\langle\sigma v\rangle=a_{l}x^{-l}$ from eq. (4), the $J$ factor becomes

[TABLE]

As a result, the relic density of WIMP dark matter is given by

[TABLE]

In the case with a resonance having a narrow width with $\epsilon_{R}>0$, when $\langle\sigma v\rangle=b_{l}\gamma_{R}\,\epsilon^{l+\frac{1}{2}}_{R}x^{3/2}e^{-x\epsilon_{R}}$ from eq. (10), the $J$ factor becomes

[TABLE]

In Fig. 3, we draw the ratio of $J$-factors for the $2\rightarrow 2$ annihilation cross section with $s$-wave overall factor at on- and off-resonance as a function of $\epsilon_{R}$. Thus, the large enhancement of the thermal-averaged cross section stands out in the $J$-factors, helping reducing the relic density to a right value without a large coupling. We note that the ratio of $J$-factors changes by order of magnitude, depending on $\epsilon_{R}$ below $0.1$.

4.2 Boltzmann equation for SIMP

The Boltzmann equation for SIMP dark matter is given by

[TABLE]

Similarly as in the WIMP case, we rewrite the above equation for the relic abundance of dark matter, $Y_{\rm DM}=n_{\rm DM}/s$, as follows,

[TABLE]

where $\rho\equiv s^{2}(m_{\rm DM})/H(m_{\rm DM})$. Therefore, we obtain the solution to the Boltzmann equation as

[TABLE]

with

[TABLE]

As a result, the relic density of SIMP dark matter is given by

[TABLE]

In the case without a resonance, when $\langle\sigma v^{2}\rangle=a_{l}x^{-l}$ from eq. (14), the $K$ factor becomes

[TABLE]

In the case with a resonance having a narrow width with $\epsilon_{R}>0$, when $\langle\sigma v^{2}\rangle=b_{l}\,\epsilon^{l+2}_{R}x^{3}e^{-\frac{3}{2}x\epsilon_{R}}$ from eq. (24), the $K$ factor becomes

[TABLE]

As a result, we find that the $K$-factor has a different dependence on $\epsilon_{R}$ from the one of the $J$-factor in the previous section, due to the fact that the phase space in the velocity average for the SIMP case is more sensitive to $\epsilon_{R}$ than for the WIMP case. In Fig. 4, we depict the ratio of $K$-factors for the $3\rightarrow 2$ annihilation cross section with $s$-wave overall factor at off- and on-resonance as a function of $\epsilon_{R}$. Thus, we find that the $K$-factor becomes suppressed at small $\epsilon_{R}$ unlike the WIMP case while there is an optimal value of $\epsilon_{R}$ for which the $K$-factor is maximized.

5 Benchmark models for SIMP dark matter

In this section, we discuss some benchmark models for SIMP dark matter, with or without a resonance. We first consider a complex scalar dark matter in models with discrete gauge symmetries and then dark mesons in models with hidden non-abelian gauge symmetries.

5.1 SIMP dark matter with discrete gauge symmetries

We consider discrete symmetries as remnants of a dark local $U(1)$ after it is spontaneously broken by a Higgs mechanism. Then, the $3\rightarrow 2$ processes appear with dark Higgs resonance $h^{\prime}$ for the $Z_{3}$ case [5] and with extra scalar resonance $S$ for the $Z_{5}$ case [8]. Dark matter is a complex scalar $\chi$ with $q_{\chi}=+1$ in both cases or another complex scalar $S$ with $q_{S}=+3$ in the $Z_{5}$ case. In both cases, the $3\rightarrow 2$ processes are $s$-wave so our previous discussion in Section 2.2 for the thermal average of the $SO(9)$ invariant velocity expansion applies.

After a dark local $U(1)$ is broken into a discrete symmetry $Z_{n}$ due to a VEV of a charged scalar $\phi$ with $q_{\phi}=n$, the relevant interaction terms for SIMP dark matter in the dark sector are given as follows [5, 8],

[TABLE]

Here, $v^{\prime}$ is the VEV of a dark Higgs, which is expanded as $\phi=(v^{\prime}+h^{\prime})/\sqrt{2}$. Moreover, the dark photon $Z^{\prime}$ gets mass of $m_{Z^{\prime}}=3g_{D}v^{\prime}$ or $5g_{D}v^{\prime}$ in the $Z_{3}$ or $Z_{5}$ cases. The resonance poles for $3\rightarrow 2$ processes appear at $m_{h^{\prime}}=3m_{\chi}$ in the $Z_{3}$ case and $m_{S}=3m_{\chi}$ or $m_{\chi}=3m_{S}$ in the $Z_{5}$ case. For the $3\rightarrow 2$ dominance, we need to suppress the $2\rightarrow 2$ annihilations in the dark sector, requiring that $m_{Z^{\prime}},m_{h^{\prime}}>m_{\chi}$.

In the non-relativistic limit of dark matter, the $3\rightarrow 2$ annihilation cross sections with a resonance for discrete gauge symmetries take the form,

[TABLE]

where $\epsilon_{n}=\frac{m^{2}_{n}-9m^{2}_{\chi}}{9m^{2}_{\chi}}$ and $\gamma_{n}=\frac{m_{n}\Gamma_{n}}{9m^{2}_{\chi}}$ with $m_{3}=m_{h^{\prime}}$ and $m_{5}=m_{S}$ or $m_{\chi}$, and $C_{n}$ is given by

[TABLE]

Here, $C^{\chi,S}_{5}$ denote the coefficients for $\chi$ and $S$ SIMP dark matters in $Z_{5}$ models, respectively, and ${\beta}_{\chi}\equiv\sqrt{1-4m^{2}_{\chi}/m^{2}_{h^{\prime}}}$, ${\beta}^{\prime}_{\chi}\equiv\sqrt{1-4m^{2}_{\chi}/m^{2}_{S}}$, and $\beta_{S}\equiv\sqrt{1-4m^{2}_{S}/m^{2}_{\chi}}$. The width of the resonance is approximated by the partial decay width of the two-body decay mode, $h^{\prime}\rightarrow\chi\chi^{*}$ in the $Z_{3}$ case and $S\rightarrow\chi^{*}\chi^{*}$ or $\chi\rightarrow SS$ in the $Z_{5}$ case, as follows,

[TABLE]

We note that there are also three-body decay modes of the resonance in both $Z_{3}$ and $Z_{5}$ models, but the corresponding decay rates are suppressed by extra phase space, roughly by $\epsilon^{2}_{n}/(4\pi^{2})$ for a constant squared decay amplitude, as compared to the two-body decay rates. Therefore, near the resonance with $\epsilon_{n}\lesssim 0.1$, the three-body decay contributions to the total decay rate of the resonance can be ignored.

Then, since the $3\rightarrow 2$ processes are $s$-wave in all the cases above, using the result in eq. (20), we obtain the thermal average as

[TABLE]

with $z_{n}\equiv\epsilon_{n}+i\gamma_{n}$. In the narrow width approximation, the above result becomes

[TABLE]

In Fig. 5, we show the relic density $\Omega h^{2}$ as a function of $\epsilon_{n}$ ($C_{n}$) in the left (right) panel for a fixed $C_{n}$ ($\epsilon_{n}$). These results are for the resonance cases with $s$-wave annihilation, which are applicable to $Z_{n}$ models. DM mass is chosen to $100(200)\,{\rm MeV}$ in the upper (lower) panel. Moreover, in Fig. 6, we show the parameter space for DM cubic coupling and mass satisfying the relic density measured by Planck, depending on the value of $\epsilon_{R}=0.01,0.02,0.06$ from top to bottom. The DM cubic coupling is given by $\kappa$ for the $Z_{3}$ model and $\lambda_{3}$ for the $Z_{5}$ model. Here, the narrow width approximation is assumed. As a consequence, we find that the required value of $\kappa$ for the relic density varies by a factor of $3-5$, depending on $\epsilon_{n}$. We note that we kept only the resonant channels in $Z_{n}$ models to show the dependence on the resonance pole but extra non-resonant channels to the same $3\rightarrow 2$ process can allow for a smaller $\kappa$ coupling [5, 8, 9]. Furthermore, other couplings such as $\lambda_{\chi}$ make the model consistent with the bound on the self-scattering cross section of dark matter [5, 8, 9].

5.2 Dark mesons

We consider non-abelian gauge symmetries with flavor groups in the dark sector, such as $SU(N_{c})$ gauge symmetry and $SU(N_{f})\times SU(N_{f})/SU(N_{f})$ coset space for flavor group. The Wess-Zumino-Witten(WZW) terms [12, 13] are responsible for $3\rightarrow 2$ processes for dark mesons [2, 14]. When dark quarks are charged under a dark local $U(1)$, the dark gauge boson $Z^{\prime}$ has vector-like couplings to dark quarks, resulting in dark meson couplings such as $Z^{\prime}-\pi_{i}-\pi_{j}-\pi_{k}$ and $Z^{\prime}-\pi_{i}-\pi_{j}$ [13, 4]. In this case, a gauge kinetic mixing between dark photon and SM hypercharge gauge boson allows for dark matter to be in kinetic equilibrium until freeze-out, and the extra $2\rightarrow 2$ (semi-)annihilation channels, $\pi\pi\rightarrow Z^{\prime}Z^{\prime}(\pi)$, is kinematically forbidden222We note that the forbidden channels can be still important for determining the relic density if $m_{\pi}<m_{Z^{\prime}}\lesssim 2(\frac{3}{2})m_{\pi}$ [9]. for $m_{Z^{\prime}}>m_{\chi}$. Furthermore, the $3\rightarrow 2$ process for dark mesons can have a resonance at $m_{Z^{\prime}}=3m_{\pi}$.

The effective Lagrangian for dark mesons including WZW terms is the following,

[TABLE]

where $F$ is the decay constant of dark mesons, $\pi\equiv 2T^{a}\pi^{a}$ with $T^{a}$ satisfying $[T^{a},T^{b}]=if_{abc}T^{c}$ and belonging to $SU(N_{f})\times SU(N_{f})/SU(N_{f})$ (e.g. $\lambda^{a}=2T^{a}$ being Gell-Mann matrices for $N_{f}=3$), and the covariant derivative for dark mesons is given by $D_{\mu}\pi=\partial_{\mu}\pi+ig_{D}Z^{\prime}_{\mu}[Q_{D},\pi]$. Here, $Q_{D}$ is the dark charge operator which is chosen to be ${\rm Tr}\,Q_{D}\neq 0$ and $Q^{2}_{D}=1$ for the absence of chiral anomalies [4, 14]. For the $3\rightarrow 2$ dominance, we need to suppress $\pi\pi\rightarrow Z^{\prime}Z^{\prime}(\pi)$, requiring $m_{Z^{\prime}}\gtrsim 2(\frac{3}{2})m_{\pi}$.

First, the WZW terms for dark mesons lead to the $d$-wave suppressed $3\rightarrow 2$ processes for dark mesons and the corresponding annihilation cross section takes the following form in the velocity expansion,

[TABLE]

Thus, as there is no resonance, we can make use of eq. (18) to get the thermal average as

[TABLE]

where $C_{WZW}$ depends on group factors. The result agrees with Ref. [2].

The gauged WZW terms for dark mesons lead to additional $3\rightarrow 2$ processes for dark mesons with a resonance. After the dark photon is integrated out, the resulting effective interaction is

[TABLE]

For the resonance case, we only have to replace $1/m^{2}_{Z^{\prime}}$ by $-1/(s-m^{2}_{Z^{\prime}})$ where $s$ is the center of mass energy for $3\rightarrow 2$ processes. As the gauged WZW terms lead to the effective 5-point interactions of the same form as the one of the ungauged WZW terms, the corresponding $3\rightarrow 2$ annihilation cross section is given by

[TABLE]

where $\epsilon_{Z^{\prime}}=\frac{m^{2}_{Z^{\prime}}-9m^{2}_{\pi}}{9m^{2}_{\pi}}$ and $\gamma_{Z^{\prime}}=\frac{m_{Z^{\prime}}\Gamma_{Z^{\prime}}}{9m^{2}_{\pi}}$ and $C_{gWZW}$ depends on group factors as well as the dark charge operator $Q_{D}$. Here, the decay rate of the dark photon is approximated by the two-body decay to be

[TABLE]

Then, in the narrow width approximation, using the result in eq. (24) and doing an explicit integration for the thermal average of the terms with $v^{2}_{i}v^{2}_{j}$, $i\neq j$, we get the thermal average of the additional $3\rightarrow 2$ annihilation cross section as

[TABLE]

In this case, the resulting averaged cross section has a higher power dependence on $\epsilon_{R}$ near resonance, due to the overall $d$-wave suppression of the $3\rightarrow 2$ annihilation cross section.

6 Generalizations

In this section, we generalize our previous discussion on the thermal average to the cases with non-degenerate masses in the initial states or the $n+2\rightarrow 2$ annihilation processes.

6.1 $3\rightarrow 2$ co-annihilations

The results on thermal average can be generalized to the case with non-degenerate masses in the initial states of the $3\rightarrow 2$ process [15], namely, the co-annihilation between multiple components of dark matter. In this case, we consider the momenta $p_{i}(i=1,2,3)$ instead of velocities $v_{i}(i=1,2,3$) in the integration and the velocity expansion of the $3\rightarrow 2$ annihilation cross section.

For simplicity, we take the $3\rightarrow 2$ annihilation cross section as a function of the total kinetic energy, namely, $K=\frac{p^{2}_{1}}{2m_{1}}+\frac{p^{2}_{2}}{2m_{2}}+\frac{p^{2}_{3}}{2m_{3}}$, in the non-relativistic limit. Then, the thermal average for the case with non-degenerate masses can be simply given by the one for the case with degenerate masses where $m_{\rm DM}$ is replaced by $(m_{1}+m_{2}+m_{3})/3$ in eqs. (14) or (20), depending on whether the process is non-resonant or resonant. This result is particularly useful for the $s$-wave $3\rightarrow 2$ process with non-degenerate masses. But, if the $3\rightarrow 2$ co-annihilation process is velocity-suppressed, one needs to take care of the thermal average of all the individual velocity terms, that are not necessarily $SO(9)$ invariant due to mass differences.

6.2 Higher-order DM annihilations

We can generalize our previous discussion to the thermal average for $n+2\rightarrow 2$ annihilation processes [1, 16] with initial particles having the same masses. We denote the corresponding annihilation cross section by $(\sigma v^{n+1})$ and the corresponding thermal average is given by

[TABLE]

Then, in the case of the $SO(3(n+2))$ invariant velocity expansion, namely, $(\sigma v^{n+1})=\sum_{l=0}^{\infty}\frac{a_{l}}{l!}\,\eta^{l}$ with $\eta=\frac{1}{2}(v^{2}_{1}+\cdots+v^{2}_{n+2})$, we obtain the thermal average in a simple matter as

[TABLE]

Likewise in the case of $3\rightarrow 2$ processes, in most cases, the most important terms appear up to $p$-wave terms that are $SO(3(n+2))$ invariant, so the above result gives rise to a good approximation for the full average of $n+2\rightarrow 2$ processes.

7 Conclusions

We have presented general results on the thermal average of $3\rightarrow 2$ annihilation cross sections of dark matter. The results can be important to improve the calculation of the dark matter abundances in the case with strong velocity-dependence and resonance poles. We have shown some examples on SIMP dark matter where the obtained results can be applied and have extended our discussion to the case with the $3\rightarrow 2$ co-annihilation and even higher-order annihilation processes.

Acknowledgments

We would like to thank the CERN Theory group for hospitality and participants in the CERN Theory Institute on “New Physics at Intensity Frontiers” for discussion during the final stage of the project. We would like to thank Myeonghun Park for the discussion on three-body decays. The work is supported in part by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF-2016R1A2B4008759). The work of SMC is supported in part by TJ Park Science Fellowship of POSCO TJ Park Foundation. MS is supported by IBS under the project code, IBS-R018-D1.