Perturbative Unitarity of Strongly Interacting Massive Particle Models

Ayuki Kamadaa,b, Shin Kobayashic, and Takumi Kuwaharad

a Institute of Theoretical Physics, Faculty of Physics, University of Warsaw,
ul. Pasteura 5, PL-02-093 Warsaw, Poland

b Kavli Institute for the Physics and Mathematics of the Universe (WPI), The
University of Tokyo Institutes for Advanced Study, The University of Tokyo,
Kashiwa 277-8583, Japan

c ICRR, University of Tokyo, Kashiwa, Chiba 277-8582, Japan

d Center for High Energy Physics, Peking University, Beijing 100871, China

Dark pion is a promising candidate for the strongly interacting massive
particle dark matter. A large pion self-coupling $m_{\pi}/f_{\pi}$ tends to be
required for correct relic abundance, and hence the partial-wave amplitudes
can violate the perturbative unitarity even for the coupling within naïve
perturbative regime. We improve the partial-wave amplitudes in order to
satisfy the optical theorem. We demonstrate that the improvement is relevant
only for semi-relativistic pions, and thus this does not affect the self-
scattering cross section at the cosmic structures. We also discuss the impact
of the improvement of the $\pi\pi\pi\to\pi\pi$ scattering process, and we find
that there is an upper bound on $m_{\pi}$ at which the correct relic abundance
is never achieved even at large $m_{\pi}/f_{\pi}$ due to the optical theorem.

##  1 Introduction

Little is known about particle nature of dark matter (DM) even though the
existence of DM has been firmly confirmed by the astrophysical observations.
Strongly interacting massive particles (SIMPs) [1] are an interesting
framework for the thermal relic of sub-GeV DM: the thermal relic DM is
determined by the freeze-out of $3\to 2$ processes. Self-interactions of DM
are generically sizable to get the correct relic abundance in this framework,
and it leads to a large $2\to 2$ self-scattering. The small-scale structure of
the Universe may indicate the sizable self-scattering of DM [2] (see Ref. [3]
a review).

Dark sector is a hypothetical sector where DM resides, and has its own gauge
dynamics. As a consequence of dark strong dynamics, the dark sector would
consist of composite particles (dark hadrons) at the low-energy scale as with
the SM hadrons [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39,
40, 41, 42, 43, 44, 45, 46, 47, 48, 49] (see Ref. [50] for a review). Ref.
[51] proposed a model of the SIMP framework where the dark pion is identified
as DM. The dark pions arise as the pseudo-Nambu-Goldstone boson (pNGB) from
the strong dynamics in the dark sector. The dark pions have the $2\to 2$ self-
interaction and the $3\to 2$ number-changing process induced by the Wess-
Zumino-Witten (WZW) term [52, 53].

The chiral perturbation theory ($\chi$PT) describes the interactions among
dark pions. The pion self-coupling, which is defined by the ratio of the pion
mass and decay constant $m_{\pi}/f_{\pi}$, determines the size of the pionic
scattering processes. The $2\to 2$ self-scattering interaction arises at the
leading order terms of $\mathcal{O}(m_{\pi}^{2}/f_{\pi}^{2})$, while the $3\to
2$ number-changing interaction via the WZW term appears as an
$\mathcal{O}(m_{\pi}^{5}/f_{\pi}^{5})$ term of $\chi$PT. The pion self-
coupling is found to be larger than unity in order to explain the relic
abundance and to evade constraints on the self-scattering cross section [51].
Meanwhile, we cannot validate the perturbative expansion of $\chi$PT unless
$m_{\pi}/f_{\pi}\lesssim 4\pi$. We would not be able to ignore contributions
from resonances and the higher-order terms of the chiral Lagrangian for the
pion self-coupling near the naïve perturbative bound. For the consistent
treatment of the chiral expansion, Ref. [54] has discussed the impact of the
higher order of the chiral Lagrangian on SIMP scenarios.

We may encounter the other bound on the self-scattering cross section of dark
pions in the SIMP models even when the perturbative expansion of $\chi$PT is
valid. The unitarity of the $S$-matrix imposes constraints on the partial-wave
amplitude. The partial-wave amplitude for the $2\to 2$ self-scattering is
$T\simeq m_{\pi}^{2}/(32\pi f_{\pi}^{2})$ at the tree-level, and the
perturbative unitarity places an upper bound $1/2\beta$ on $\mathrm{Re}T$ with
$\beta$ corresponding to velocity of pions. The perturbative unitarity bound
$m_{\pi}/f_{\pi}\lesssim 4\sqrt{\pi}/\sqrt{\beta}$ gets weaker at the cosmic
structures due to the DM velocity at maximum of $\sim 10^{-2}$, while this
bound can be important for the annihilation mechanism determining the current
relic abundance since dark pions are semi-relativistic. In other words, there
are two bounds on the pion self-couplings: one originates from the limitation
of the perturbative expansion and another is the perturbative unitarity of the
scattering processes. In the chiral limit ($m_{\pi}\to 0$), the perturbative
unitarity violation will be cured by resumming multiple rescattering
processes, which is known as the “self-healing” mechanism [55]. In this paper,
we will propose the improvement of $2\to 2$ and $3\to 2$ partial-wave
amplitudes in a similar way to the “self-healing” mechanism in the non-chiral
limit since we focus on the dark pion DM.

This paper is organized as follows. We will show the chiral Lagrangian for the
SIMP models in Section 2. We discuss the perturbative unitarity of
$\pi\pi\to\pi\pi$ scattering cross section and thermally-averaged
$\pi\pi\pi\to\pi\pi$ cross section in Section 3, and we will see these cross
sections would violate the perturbative unitarity at large $m_{\pi}/f_{\pi}$.
In Section 4, we propose the improved amplitude, which satisfies the optical
theorem automatically, and then we will apply the procedure to the SIMP
models. Section 5 is devoted to conclusions of our work.

##  2 Chiral Lagrangian for SIMP

We discuss the SIMP model that is realized by the confining gauge dynamics of
a gauge group $G_{\mathrm{local}}$. We consider $N_{f}$-flavor quarks with the
mass below the dynamical scale of $G_{\mathrm{local}}$ in the UV description
of this model. This model possesses an approximate global symmetry $G$ among
quarks that is broken by the mass terms. It is believed that this model leads
to the chiral symmetry breaking, and that chiral condensation breaks the
global symmetry into the subgroup $H$. In the following, we consider three
classes of the models: (i) $G_{\mathrm{local}}=SU(N_{c})$, $G=SU(N_{f})\times
SU(N_{f})$, and $H=SU(N_{f})$; (ii) $G_{\mathrm{local}}=SO(N_{c})$,
$G=SU(N_{f})$, and $H=SO(N_{f})$; and (iii) $G_{\mathrm{local}}=USp(N_{c})$,
$G=SU(N_{f})$, and $H=USp(N_{f})$ (with even integers $N_{c}$ and $N_{f}$ [56,
57, 58, 59]).

The dark pions, which are the pNGBs of the chiral symmetry breaking, are the
fundamental degrees of freedom in the low-energy effective theory of this
model. We are interested in the dark pions for realizing SIMP framework, and
hence we focus on the chiral expansion with momentum $p^{2}\simeq m_{\pi}^{2}$
in the non-relativistic limit. The coset space $G/H$ is parameterized by
$N_{\pi}$ pNGB fields $\pi^{a}$ [60, 61], which corresponds to the broken
generators $T^{a}$ with $a=1,\cdots,N_{\pi}$. The parameterization and the
normalization of the dark pion fields are the following.

$\begin{split}\Sigma=\exp\left(\frac{2i\Pi}{f_{\pi}}\right)\,,\qquad\Pi\equiv\pi^{a}T^{a}\,,\qquad\mathrm{Tr}(T^{a}T^{b})=\frac{1}{2}\delta^{ab}\,,\end{split}$
(1)

and only for the case (iii), the non-linear sigma model field $\Sigma$ is

$\begin{split}\Sigma=\exp\left(\frac{2i\Pi}{f_{\pi}}\right)J\,,\end{split}$
(2)

where $J$ is the symplectic metric that satisfies
$J^{T}=-J\,,J^{2}=-\mathbf{1}$. Here, $f_{\pi}$ denotes the pion decay
constant. Under the residual symmetry $H$, the dark pion fields transform as
(i) adjoint representation, (ii) rank-2 symmetric tensor representation, and
(iii) rank-2 antisymmetric tensor representation for each different symmetry.
The relevant Lagrangian of the dark pions is given by

$\begin{split}\mathcal{L}=\frac{f_{\pi}^{2}}{4}\mathrm{Tr}\left[(D_{\mu}\Sigma)^{\dagger}D^{\mu}\Sigma\right]+\frac{Bf_{\pi}^{2}}{4}\mathrm{Tr}\left[M^{\dagger}\Sigma+M\Sigma^{\dagger}\right]+\mathcal{L}_{\mathrm{WZW}}\,.\end{split}$
(3)

Here, the first-two terms are the leading order (LO) terms in the $\chi$PT and
the soft chiral symmetry breaking term, and the third term is the Wess-Zumino-
Witten (WZW) term [52, 53]. We take the quark mass matrix $M$ to be invariant
under $H$, and hence the pion mass is universal $m_{\pi}^{2}=Bm_{q}$ with
$m_{q}$ being the quark mass.

The relevant Lagrangian contains the four-point interactions among the dark
pions, which induce the self scattering of the dark pions. We obtain the four-
point interaction terms from the kinetic and mass terms by expanding $\Sigma$.

$\begin{split}\mathcal{L}\supset\frac{r_{abcd}}{4f_{\pi}^{2}}\pi^{a}\pi^{b}\partial_{\mu}\pi^{c}\partial^{\mu}\pi^{d}+\frac{1}{24}\frac{m_{\pi}^{2}}{f_{\pi}^{2}}c_{abcd}\pi^{a}\pi^{b}\pi^{c}\pi^{d}\,,\end{split}$
(4)

where $r_{abcd}$ and $c_{abcd}$ are the coefficients defined by group-
theoretical constants, which is discussed in Appendix C. We will discuss the
self-scattering cross section in detail in the next section. Meanwhile, the
five-point interaction among the dark pions arises from the WZW term. The
$\pi\pi\pi\to\pi\pi$ scattering process arising from the WZW term determines
the relic abundance of the dark pions. The WZW term is given by111 Our
conventions are different from some literature. The normalizations of both
generators and $f_{\pi}$ differ factor two from Refs [53, 51], and the
resultant coefficients are same as ours. Meanwhile, the normalization of
generators differs, but the same convention is used for $f_{\pi}$ as ours in
[62].

$\begin{split}\mathcal{L}_{\mathrm{WZW}}&=\frac{2k}{15\pi^{2}f_{\pi}^{5}}\epsilon^{\mu\nu\rho\sigma}\mathrm{Tr}\left[\Pi\partial_{\mu}\Pi\partial_{\nu}\Pi\partial_{\rho}\Pi\partial_{\sigma}\Pi\right]\\\
&\supset\frac{2k}{15\pi^{2}f_{\pi}^{5}}T_{[abcde]}\epsilon^{\mu\nu\rho\sigma}\pi^{a}\partial_{\mu}\pi^{b}\partial_{\nu}\pi^{c}\partial_{\rho}\pi^{d}\partial_{\sigma}\pi^{e}\,.\end{split}$
(5)

Here, $k=N_{C}$ for $SU(N_{C})$ and $2k=N_{C}$ for $SO(N_{C})$ and
$USp(N_{C})$ [63]. $T_{[abcde]}$ denotes the anti-symmetrization of five
broken generators:

$\begin{split}T_{[abcde]}=\frac{1}{5!}\sum\mathrm{Tr}(T^{[a}T^{b}T^{c}T^{d}T^{e]})\,.\end{split}$
(6)

The four-point interactions arise at the order of $m_{\pi}^{2}/f_{\pi}^{2}$,
while the five-point interaction is at the order of $m_{\pi}^{5}/f_{\pi}^{5}$.
In accordance with the standard order-counting of the $\chi$PT using $p^{2}$
expansion in the chiral limit, the former is the LO contribution and the
latter is the next-to-leading order (NLO) contribution in the non-chiral
limit.

A large coupling tends to be required for correct relic abundance in the SIMP
scenarios. In the dark-pion realization, the number-changing process arises
from the WZW term, which is the higher order of the chiral expansion, with the
velocity suppression. Meanwhile, the $\chi$PT breaks down
$m_{\pi}\simeq\Lambda_{\chi\mathrm{SB}}$, and the cutoff scale is expected to
be

$\begin{split}\Lambda_{\chi\mathrm{SB}}\simeq\mathrm{min}\left(\frac{4\pi
f_{\pi}}{\sqrt{N_{c}}},\frac{4\pi f_{\pi}}{\sqrt{N_{f}}}\right)\,.\end{split}$
(7)

This is known as the scale estimated by naïve dimensional analysis (NDA) [64,
65] with taking into account the large-$N_{c}$ scaling and the large-$N_{f}$
scaling. The large-$N_{c}$ scaling of $f_{\pi}$ and
$\Lambda_{\chi\mathrm{SB}}$ is known to be
$f_{\pi}^{2}\simeq\mathcal{O}(N_{c})$ and
$\Lambda_{\chi\mathrm{SB}}\simeq\mathcal{O}(1)$ [66, 67, 68, 69, 70, 71]. It
is also known that some of the partial-wave amplitudes for the scattering
process of $\pi\pi\to\pi\pi$ is proportional to $N_{f}$ in the large $N_{f}$
limit [72, 73]. 222 This discussion is based on the chiral symmetry breaking
of $SU(N_{f})\times SU(N_{f})\to SU(N_{f})$ in the original literature. We
find the similar $N_{f}$ dependence of the amplitude with the isospin singlet
even for $SU(N_{f})/SO(N_{f})$ and $SU(N_{f})/USp(N_{f})$, and hence we expect
that we get the same cutoff scale up to constant of order unity.  Since the
pion self-coupling is required to be close to its naïve perturbative bound
$m_{\pi}/f_{\pi}\lesssim\Lambda_{\chi\mathrm{SB}}/f_{\pi}$ in the SIMP
scenario, the higher-order contributions of the $\chi$PT may affect the
predictions of the relic abundance and the self-scattering cross section. Ref.
[54] has discussed the impacts of the higher-order contributions of the
$\chi$PT in the context of the SIMP: in particular, the NLO and the next-to-
next-to-leading order (NNLO) contributions.

##  3 Perturbative Unitarity

We discuss the perturbative unitarity of $\pi\pi\to\pi\pi$ and
$\pi\pi\pi\to\pi\pi$ scattering processes. The partial-wave decomposition of
the invariant amplitude for elastic scattering process is given by333 We use
$32\pi$ as a normalization instead of $16\pi$ due to the identical particles.

$\begin{split}\mathcal{M}^{ab;cd}_{2\to
2}&=32\pi\sum_{R}\sum_{\ell}(2\ell+1)P_{\ell}(\cos\theta)T^{R}_{\ell}(s)P_{R}^{ab;cd}\,.\end{split}$
(8)

Here, $P_{\ell}$ is the Legendre polynomial with $P_{\ell}(1)=1$, $s$ denotes
the collision energy, and $\theta$ is the scattering angle of the final state
pions with respect to the collision axis. $P_{R}^{ab;cd}$ denotes the
projection operators: the product of two pions is projected into the
irreducible representation $R$ of residual global symmetry $H$.
$T^{R}_{\ell}(s)$ is the partial-wave elastic amplitude for the channel of a
representation $R$. The projection operators satisfy

$\begin{split}\sum_{c^{\prime},d^{\prime}}P_{R}^{ab;c^{\prime}d^{\prime}}P_{R^{\prime}}^{c^{\prime}d^{\prime};cd}&=\delta_{RR^{\prime}}P_{R}^{ab;cd}\,,\quad(P_{R}^{ab;cd})^{\ast}=P_{R}^{cd;ab}\,,\quad\sum_{a,b}P_{R}^{ab;ab}=d_{R}\,.\end{split}$
(9)

Here, $d_{R}$ denotes the dimension of the representation $R$. With this
decomposition, the total cross section for $\pi^{a}\pi^{b}\to\pi^{c}\pi^{d}$
takes the form

$\begin{split}\sigma^{ab;cd}_{2\to
2}=\frac{32\pi}{s}\sum_{R}\sum_{\ell}(2\ell+1)|T^{R}_{\ell}(s)|^{2}P_{R}^{ab;cd}\,.\end{split}$
(10)

The unitarity of the $S$-matrix imposes the optical theorem that relates the
imaginary part of the invariant amplitude $\mathcal{M}_{\mathrm{2\to 2}}$ for
the forward scattering to the total cross section $\sigma_{\mathrm{total}}$.
Below the four-pion threshold, the optical theorem is given by

$\begin{split}\mathrm{Im}\mathcal{M}_{2\to 2}=2\sqrt{s}p(\sigma_{2\to
2}+\sigma_{2\to 3})\,.\end{split}$ (11)

Here, $p=\sqrt{s}\beta/2$ is the momentum of incoming particles with
$\beta=(1-4m_{\pi}^{2}/s)^{1/2}$ denoting the velocity of the particles in the
center of mass frame. We keep $\beta$ in the following since we focus on a
non-relativistic limit. We include the inelastic cross section $\sigma_{2\to
3}$ in the right-hand side since it is also sizable in the SIMP models, and it
vanishes below the inelastic threshold.

Figure 1:  Unitarity circle on a complex plane of the partial-wave amplitude
for the elastic scattering with a fixed partial-wave amplitude for the
inelastic scattering. Unitarity circle without the inelastic scattering is
shown as the gray-dashed line.

We define the partial-wave amplitude for the $\pi\pi\to\pi\pi\pi$ process in
the non-relativistic limit. It is challenging to give the explicit form of the
decomposition of the three-body final state into its irreducible
representations, in other words to construct the explicit form of the
projection operator $P_{R}^{abc;de}$. Meanwhile, the cross section can be
decomposed into the irreducible representations of the initial state. The non-
relativistic cross section for the $\pi\pi\to\pi\pi\pi$ process may be written
in a similar way to the $\pi\pi\to\pi\pi$ process. In the non-relativistic
limit, the elastic cross section and the inelastic cross section averaged over
all initial states are given by

$\begin{split}\sigma_{2\to
2}&\equiv\frac{1}{N_{\pi}^{2}}\sum_{a,b}\sigma^{ab;ab}_{2\to
2}=\frac{8\pi}{m_{\pi}^{2}}\frac{1}{N_{\pi}^{2}}\sum_{R}\sum_{\ell}(2\ell+1)|T^{R}_{\ell}|^{2}d_{R}\,,\\\
\sigma_{2\to
3}&=\frac{8\pi}{m_{\pi}^{2}}\frac{1}{N_{\pi}^{2}}\sum_{R}\sum_{\ell}(2\ell+1)|T^{R}_{\mathrm{in},\ell}|^{2}d_{R}\,.\end{split}$
(12)

Here, we define the inelastic amplitude $|T^{R}_{\mathrm{in},\ell}|^{2}$. When
we include the $\pi\pi\to\pi\pi\pi$ inelastic scattering process, the optical
theorem for each partial-wave amplitude takes the form,

$\begin{split}\mathrm{Im}(T^{R}_{\ell})=\beta(|T^{R}_{\ell}|^{2}+|T^{R}_{\mathrm{in},\ell}|^{2})\,.\end{split}$
(13)

Let us apply this condition to perturbative analysis in the case that
inelastic scattering is negligible. The optical theorem relates a fixed-order
amplitude to the imaginary part of the higher-order amplitude: the tree-level
partial-wave amplitude for the elastic scattering $T_{\ell}^{\mathrm{tree}}$
is real, and is related to the imaginary part of the one-loop amplitude. We
give the perturbative analysis of $\lambda\phi^{4}$ theory in Appendix A, and
discuss this point explicitly. In other words, the fixed-order amplitude does
not satisfy Eq. 13. Hence, it is important to unitarize the amplitude in order
to satisfy the optical theorem automatically as the pion self-coupling is
large.

The optical theorem places upper bounds on the partial-wave amplitudes. Fig. 1
shows the unitarity circle on a complex plane of the partial-wave amplitude
for the elastic scattering. The gray dashed line depict the unitarity circle
without the inelastic scattering; its center is located at $(0,1/2\beta)$ and
its radius is $1/2\beta$. The radius shrinks in the presence of the inelastic
scattering as depicted as the solid-line circle, and the single point
$(0,1/2\beta)$ is allowed once the inelastic scattering saturates the
unitarity bound ($|T^{R}_{\mathrm{in},\ell}|^{2}=1/4\beta^{2}$). The upper
bound on the partial-wave amplitude for the elastic scattering is given by
$|T^{R}_{\ell}|\leq 1/\beta$ as shown in the figure. The equality holds when
$\mathrm{Re}T^{R}_{\ell}=0$, $|T^{R}_{\mathrm{in},\ell}|^{2}=0$, and
$\mathrm{Im}T^{R}_{\ell}=1/\beta$, corresponding to the top of the gray-dashed
circle. The perturbative unitarity is often imposed as
$|\mathrm{Re}T^{R}_{\ell}|\leq 1/2\beta$. The equality holds when the
imaginary part satisfies $\mathrm{Im}T^{R}_{\ell}=1/2\beta$ and the inelastic
scattering is negligible.

We discuss the unitarity bound on the pion self-coupling based on the partial-
wave unitarity of tree-level
$\pi^{a}(p_{1})\pi^{b}(p_{2})\to\pi^{c}(p_{3})\pi^{d}(p_{4})$ self-scattering
cross section. At the LO terms of $\chi$PT, the four-point interaction is
given by Eq. 4, and the tree-level amplitude of the process is obtained as
follows:

$\begin{split}\mathcal{M}^{ab;cd}_{2\to
2}(p_{1},p_{2};p_{3},p_{4})&=\frac{m_{\pi}^{2}}{f_{\pi}^{2}}c_{abcd}-\frac{s}{f_{\pi}^{2}}r_{abcd}-\frac{t}{f_{\pi}^{2}}r_{acbd}-\frac{u}{f_{\pi}^{2}}r_{adbc}\,.\end{split}$
(14)

Here, the Mandelstam variables are defined by
$s=(p_{1}+p_{2})^{2}\,,t=(p_{1}-p_{3})^{2}$ , and $u=(p_{1}-p_{4})^{2}$. The
singlet channel for $\pi\pi\to\pi\pi$ has the larger coefficient compared to
the others in an $SU(N_{f})\times SU(N_{f})/SU(N_{f})$ with a large $N_{f}$
limit, and the $s$-wave amplitude for this channel takes the form,

$\begin{split}T_{\ell=0}^{R=1}=\frac{N_{f}}{8\pi}\frac{m_{\pi}^{2}}{f_{\pi}^{2}}\,,\end{split}$
(15)

where we take the non-relativistic limit
$s=4m_{\pi}^{2}+4m_{\pi}^{2}\beta^{2}$. Once the pion self-coupling is in the
range of $m_{\pi}/f_{\pi}\gtrsim\sqrt{4\pi/N_{f}\beta}$, the partial-wave
amplitude goes beyond the unitarity bound. The elastic scattering cross
section can exceed the perturbative unitarity bound when the pion self-
coupling within the following range:

$\begin{split}\sqrt{\frac{4\pi}{\beta
N_{f}}}\lesssim\frac{m_{\pi}}{f_{\pi}}\lesssim\mathrm{min}\left(\frac{4\pi}{\sqrt{N_{c}}},\frac{4\pi}{\sqrt{N_{f}}}\right)\,.\end{split}$
(16)

This upper bound originates from the condition where tree and one-loop
contributions are compatible to each other.

We also discuss the unitarity bound on the partial-wave amplitude for the
inelastic scattering,
$\pi^{a}(p_{1})\pi^{b}(p_{2})\to\pi^{c}(p_{3})\pi^{d}(p_{4})\pi^{e}(p_{5})$.
The WZW term (5) provides the five-pion interaction that induces the inelastic
scattering. The invariant amplitude at the tree-level is given by

$\begin{split}\mathcal{M}_{\mathrm{tree}}^{abcde}(p_{1},p_{2},p_{3},p_{4},p_{5})=\frac{16k}{\pi^{2}f_{\pi}^{5}}T_{[abcde]}\epsilon^{\mu\nu\rho\sigma}p_{2\mu}p_{3\nu}p_{4\rho}p_{5\sigma}\,.\end{split}$
(17)

Here, $T_{[abcde]}$ is defined in Eq. 6. We compute the $\pi\pi\to\pi\pi\pi$
cross section in the non-relativistic limit, and we extract the partial-wave
amplitude for the inelastic scattering by comparing the cross section with Eq.
12. Since the amplitude has only $p$-wave component at the tree level even for
the relativistic limit, only antisymmetric representations in the irreducible
decomposition of initial two-pion states contribute to the amplitude. We find
that the partial-wave amplitude takes the form,

$\begin{split}(T^{R}_{\mathrm{in},\ell=1})^{2}=\left(\frac{1}{4\sqrt{3}\pi}\right)^{7}4k^{2}\beta
v^{8}\left(\frac{m_{\pi}}{f_{\pi}}\right)^{10}\frac{t_{R}^{2}}{d_{R}}\,,\end{split}$
(18)

where we take the non-relativistic limit
$s=9m_{\pi}^{2}+\sqrt{3}m_{\pi}^{2}v^{2}$. $t_{R}^{2}$ denotes a group
theoretical factor, which scales as $N_{f}^{5}$. $\beta$ and $v$ correspond to
the relative velocities of the two-body initial state and of the three-body
final state, respectively. In general, the inelastic partial-wave amplitude
has the velocity dependence of
$(T^{R}_{\mathrm{in},\ell})^{2}\propto\beta^{2\ell-1}v^{4+2\lambda}$ with
$\lambda$ denoting hyperangular quantum number, and the $2\to 3$ cross section
may scale as $v^{6}$ ($\lambda=1$) at the leading order of the $p$-wave cross
section. In the dark-pion realization, the partial-wave amplitude is the next-
leading order of the $p$-wave amplitude ($\lambda=2$).

The unitarity bound on the partial-wave amplitude, $1/2\beta$, places an upper
bound on the pion self-coupling $m_{\pi}/f_{\pi}$. The smallest representation
in dimension provides the largest partial-wave amplitude. In an
$SU(N_{f})\times SU(N_{f})/SU(N_{f})$ with a large $N_{f}$ limit, $d_{R}$ for
the adjoint representation scales as $N_{f}^{2}$, while $d_{R}$ for the larger
representations scale as $N_{f}^{4}$. The partial-wave amplitude for the
adjoint representation gives stronger bound on $m_{\pi}/f_{\pi}$ in a large
$N_{f}$ limit, and it takes the form,

$\begin{split}\left(\frac{m_{\pi}}{f_{\pi}}\right)^{10}\lesssim\frac{3(8\sqrt{3}\pi)^{7}}{4}\frac{1}{N_{f}^{3}N_{c}^{2}}\frac{1}{\beta^{3}v^{8}}\,.\end{split}$
(19)

This bound on $m_{\pi}/f_{\pi}$ is numerically weaker than the perturbative
unitarity bound obtained from the $2\to 2$ scattering for the velocity
$\beta\simeq v\simeq 1$.

Next, we discuss the partial-wave unitarity of the thermally-averaged cross
section for the $\pi\pi\pi\to\pi\pi$ process. We will give detailed
computations in Appendix B. When the invariant amplitude for the
$\pi\pi\to\pi\pi\pi$ has time-reversal invariance444 This is ensured by the
$CPT$ invariance in the dark pion case.  , the cross section for
$\pi\pi\pi\to\pi\pi$ in the non-relativistic limit
($s=9m_{\pi}^{2}+\sqrt{3}m_{\pi}^{2}v^{2}$) is also written by the use of the
partial-wave amplitude $|T^{R}_{\mathrm{in},\ell}|^{2}$, and hence the cross
section takes the form:

$\begin{split}\sigma_{3\to
2}v^{2}(\bar{k})=\frac{192\sqrt{3}\pi^{2}}{\bar{k}^{4}m_{\pi}}\frac{1}{N_{\pi}^{3}}\sum_{R}\sum_{\ell}(2\ell+1)4\beta^{2}|T^{R}_{\mathrm{in},\ell}(9m_{\pi}^{2})|^{2}d_{R}\,,\end{split}$
(20)

with the Jacobi momentum $\bar{k}$ characterizing the kinematic variables of
the initial state pions. The partial-wave unitarity
$|T^{R}_{\mathrm{in},\ell=1}|^{2}\leq 1/4\beta^{2}$ places an upper bound on
the cross section.

$\begin{split}\sigma_{3\to
2:\mathrm{uni}}v^{2}(\bar{k})=3\times\frac{192\sqrt{3}\pi^{2}}{\bar{k}^{4}m_{\pi}}\frac{1}{N_{\pi}^{3}}\sum_{R}d_{R}\,.\end{split}$
(21)

Here, a factor of 3 arises from the fact that $\pi\pi\pi\to\pi\pi$ scattering
involves only $p$-wave. The thermally-averaged cross section is obtained by
integrating the cross section over $\overline{k}$ with the Maxwell-Boltzmann
distribution.

$\begin{split}\langle\sigma_{3\to
2}v^{2}\rangle_{\mathrm{uni}}=3\times\frac{12\sqrt{3}\pi^{2}x^{2}}{m_{\pi}^{5}}\frac{1}{N_{\pi}^{3}}\sum_{R}d_{R}\,.\end{split}$
(22)

Here, $x=m_{\pi}/T$ with $T$ denoting temperature. This bound differs by a
factor of $1/4$ from what is obtained in Refs. [74, 75]. We use the upper
bound on the partial-wave amplitude for the inelastic scattering given by
$1/2\beta$ in the presence of the elastic scattering, while the upper bound
for the elastic scattering ($1/\beta$) has been used for constraining the
thermally-averaged cross section in Refs. [74, 75].555 We obtain the upper
bounds by the use of an independent way from Refs. [74, 75, 76]. As shown in
Appendix B, we directly compute the upper bound on the hyperangle averaged
cross section in the non-relativistic limit, and then we take the thermal
average of it. Meanwhile, the authors in Refs. [74, 75, 76] obtain the upper
bound on the thermally-averaged $n\to 2$ cross section by using the detailed
balance condition. Our result is consistent with Refs. [74, 75] up to a factor
of $1/4$ due to maximizing the partial-wave amplitude in the presence of the
elastic amplitude, but our result disagrees with the results in Ref [76].  The
summation runs only over the antisymmetric representations since the inelastic
amplitude is nonzero only for the antisymmetric representations.

Let us now consider the $\pi\pi\pi\to\pi\pi$ scattering at the leading
contribution in the chiral Lagrangian. We have given the partial-wave
amplitude for the $\pi\pi\to\pi\pi\pi$ cross section at the leading
contribution for the process in Eq. 18. The thermally-averaged cross section
for the $\pi\pi\pi\to\pi\pi$ process is666 Our result differs from the
original literature [51] by a factor of $1/3$. The difference arises from
mismatching frames of reference; the reference frame for computing the cross
section $\sigma_{3\to 2}$ is generally different from the frame for the
thermal averaging. In particular, the cross section explicitly depends on
momenta of particles in the dark pion SIMP. As for the thermal averaging, the
cross section should be written in the Galilean invariant way. We discuss this
point in Appendix B in detail.

$\begin{split}\langle\sigma_{3\to
2}v^{2}\rangle=\frac{5\sqrt{5}k^{2}m_{\pi}^{5}}{6\pi^{5}x^{2}f_{\pi}^{10}}\frac{t^{2}}{N_{\pi}^{3}}\,.\end{split}$
(23)

where $t^{2}\equiv\sum(T_{[abcde]})^{2}$. We show the explicit forms of
$t^{2}$ in different symmetry structures in Appendix C. $t^{2}/N_{\pi}^{3}$
approaches to $1/N_{f}$ at the large-$N_{f}$ limit.

The thermally-averaged cross section will exceed the unitarity bound given by
Eq. 22 in the large $m_{\pi}/f_{\pi}$ limit. Once we take the naïve
perturbative bound of the pion self-coupling, $m_{\pi}/f_{\pi}\simeq 4\pi$,
the cross section takes the form

$\begin{split}\langle\sigma_{3\to
2}v^{2}\rangle\simeq\frac{5\sqrt{5}(4\pi)^{5}}{6x^{2}m_{\pi}^{5}}\frac{4}{3N_{f}}\,,\end{split}$
(24)

in an $SU(N_{f})\times SU(N_{f})/SU(N_{f})$ with a large $N_{f}$ limit. The
annihilation cross section is safe from the unitarity violation at large $x$,
namely at low temperature, since the cross section is suppressed $x^{-2}$
while the unitarity bound is proportional to $x^{2}$. The cross section,
however, will be beyond the unitarity bound at small $x$ even for the pion
self-coupling below the naïve perturbative bound.

##  4 Improved Cross Section

At the large pion self-coupling, both of $\pi\pi\to\pi\pi$ and
$\pi\pi\pi\to\pi\pi$ scattering processes at tree-level can be beyond the
perturbative unitarity even for the coupling within the naïve perturbative
range. In this section, we propose the improvement of the partial-wave
amplitude to take into account the partial-wave unitarity of the cross
section. We may use perturbative analysis for the pion self-coupling within
the naïve perturbative range given by Eq. 16, but we have to modify the
amplitude to satisfy the optical theorem. In particular, we need to
incorporate the imaginary part of the amplitude at the LO into an improved
amplitude. We define the improved amplitude whose real and imaginary parts
appropriately reproduce the original amplitude at the LO as follows:

$\begin{split}T^{\mathrm{imp},R}_{\ell}&=\frac{T^{\mathrm{tree},R}_{\ell}+i\beta(T^{\mathrm{tree},R}_{\mathrm{in},\ell})^{2}}{1-i\beta[T^{\mathrm{tree},R}_{\ell}+i\beta(T^{\mathrm{tree},R}_{\mathrm{in},\ell})^{2}]}\,,\\\
|T^{\mathrm{imp},R}_{\mathrm{in},\ell}|^{2}&=\frac{(T^{\mathrm{tree},R}_{\mathrm{in},\ell})^{2}}{[1+(\beta
T^{\mathrm{tree},R}_{\mathrm{in},\ell})^{2}]^{2}+(\beta
T^{\mathrm{tree},R}_{\ell})^{2}}\,.\end{split}$ (25)

This improved partial-wave amplitude automatically satisfies Eq. 13 at any
$m_{\pi}/f_{\pi}$: namely,
$\mathrm{Im}(T^{\mathrm{imp}})=\beta(|T^{\mathrm{imp}}|^{2}+|T^{\mathrm{imp}}_{\mathrm{in}}|^{2})$
for each $R$ and $\ell$. In the absence of the inelastic scattering amplitude
$(T^{\mathrm{tree}}_{\mathrm{in}})^{2}$, this improvement of the partial-wave
amplitudes coincides with the “self-healing” mechanism in which the
perturbative unitarity violation is cured by resumming multiple rescattering
processes [55].

Figure 2:  Diagrammatic interpretation of our improvement in $2\to 2$ and
$2\to 3$ amplitudes. (Left): an effective vertex (denoted by a meshed blob)
includes the tree-level four-point vertex and the one-loop vertex induced by
the tree-level five-point vertex. (Right): the improved amplitude for $2\to 2$
and $2\to 3$ processes.

Our improvement amplitudes are not either unique or complete to remedy the
situation violating the unitarity condition, but are associated with resumming
the multiple rescattering processes that arise not only from the four-point
vertex but also from the tree-level five-point vertex. Fig. 2 shows
diagrammatic interpretation of invariant amplitudes with our improved partial-
wave amplitudes. A blob in this figure denotes the effective four-point vertex
including the tree-level four-point vertex and the one-loop contributions with
cuts and the insertion of the tree-level five-point vertex. The vertical
dotted-lines depict the cutting of intermediate propagators with replacement:

$\begin{split}\frac{i}{k^{2}-m_{\pi}^{2}+i\varepsilon}\to-2\pi
i\delta(k^{2}-m_{\pi}^{2})\,,\end{split}$ (26)

where $k$ denotes the loop momentum. These improved amplitudes correspond to a
geometric sum of multiple self-rescattering processes, in other words a
geometric sum of diagrams where all loops are cut as shown in Fig. 2. This
improved amplitude consists only of the tree-level amplitude and cut diagrams,
and hence this procedure is not appropriate when full one-loop corrections are
important.

###  4.1 Analytic Properties

Before showing the impact of the improved amplitude on the SIMP phenomenology,
we discuss the analytic properties of the improved amplitudes, such as the
scattering phase shift and the poles. The phase shift
$e^{2i\delta_{\ell}^{R}}\equiv 1+2i\beta T^{R}_{\ell}$ may be written as

$\begin{split}\cot\delta_{\ell}^{R}=\frac{1+i\beta T^{R}_{\ell}}{\beta
T^{R}_{\ell}}\,.\end{split}$ (27)

Meanwhile, the phase shift in the low-energy limit can be expanded

$\begin{split}p^{2\ell+1}\cot\delta_{\ell}^{R}\simeq-\frac{1}{a_{\ell}^{2\ell+1}}+\frac{p^{2}}{2r_{e,\ell}^{2\ell-1}}\,,\end{split}$
(28)

where $p=\sqrt{s}\beta/2$ denotes the relative momentum of the system.
$a_{\ell}$ and $r_{e,\ell}$ are known as the scattering length and the
effective range, respectively. The phase shift is generically complex in the
presence of the inelastic scattering, such as DM annihilation [77, 78].

Using our improved partial-wave amplitude, the phase shift is given by the
tree-level partial-wave amplitude as follows.

$\begin{split}\cot\delta_{\ell}^{\mathrm{imp},R}=\frac{1}{\beta
T^{\mathrm{tree},R}_{\ell}+i\beta^{2}(T^{\mathrm{tree},R}_{\ell,\mathrm{in}})^{2}}\,.\end{split}$
(29)

The phase shift with our improved amplitude is indeed complex since the
partial-wave amplitudes are real at tree level. The imaginary part of the
phase shift leads to shrinking the radius of the unitarity circle of elastic
scattering as shown in Fig. 1. As for $s$-wave amplitudes, the elastic
amplitude is proportional to $\beta^{0}$ in the non-relativistic limit and
there is no inelastic scattering. On the other hand, the $p$-wave elastic
amplitude is proportional to $\beta^{2}$. Therefore, we confirm that
$\beta^{2\ell+1}\cot\delta_{\ell}^{\mathrm{imp},R}$ is analytic at $\beta=0$
even in our improvement procedure for the application to the dark pion SIMP
models. The inelastic amplitude vanishes $\beta\lesssim\sqrt{5}/3$ due to the
kinematics, while the $\beta$-dependence of the square of the inelastic
amplitude is compatible with the expression of the effective range theory. We
note that, assuming that the imaginary part of the phase shift is negligible
compared to the real part, the phase shift is given by the tree-level partial-
wave amplitude as $\mathrm{Re}(\cot\delta_{\ell}^{\mathrm{imp},R})\simeq(\beta
T_{\ell}^{\mathrm{tree},R})^{-1}$ and
$\mathrm{Im}(\cot\delta_{\ell}^{\mathrm{imp},R})\simeq-(T^{\mathrm{tree},R}_{\ell,\mathrm{in}}/T_{\ell}^{\mathrm{tree},R})^{2}$.
Under the Born approximation, the phase shift based on the tree-level
amplitude is given by $\delta_{\ell}\simeq\beta T^{\mathrm{tree}}_{\ell}$:
i.e., $\cot\delta_{\ell}\simeq(\beta T^{\mathrm{tree}}_{\ell})^{-1}$, which is
reproduced by the expression of $\delta_{\ell}^{\mathrm{imp},R}$.

Our improved amplitudes would involve a pole. There is no inelastic channel
for $s$-wave amplitudes, and the amplitude is proportional to $\beta^{0}$ at
the LO. Therefore, the pole would appear at

$\begin{split}\beta_{\mathrm{pole}}=-\frac{i}{T^{\mathrm{tree},R}_{\ell=0}}\,.\end{split}$
(30)

It is indicated that $|T^{\mathrm{tree},R}_{\ell=0}|\gtrsim 1$, i.e.,
$|\beta|\lesssim 1$, in the non-relativistic limit. The bound state
($\mathrm{Im}\beta_{\mathrm{pole}}>0$) or the virtual level
($\mathrm{Im}\beta_{\mathrm{pole}}<0$) would appear depending on the sign of
the partial-wave amplitudes.

As for the $p$-wave amplitudes, the amplitude for the elastic scattering is
proportional to $\beta^{2}$ at the LO: we set the $\beta$-independent part of
the amplitude as
$T^{\mathrm{tree},R}_{\ell=1}=\beta^{2}t^{\mathrm{tree},R}_{\ell=1}$. The
partial-wave amplitude for inelastic scattering
$(T^{\mathrm{tree},R}_{\ell})^{2}$ is proportional to $\beta$, and hence we
set the $\beta$-independent part as
$(T^{\mathrm{tree},R}_{\ell})^{2}=\beta(t^{\mathrm{tree},R}_{\ell})^{2}$. As
far as the inelastic channels is negligible, the pole would appear at

$\begin{split}\beta_{\mathrm{pole}}^{3}=-\frac{i}{t^{\mathrm{tree},R}_{\ell=1}}\,.\end{split}$
(31)

We can find a pole on the imaginary axis of $\beta$, and hence the pole
corresponds to the bound state ($t^{\mathrm{tree},R}_{\ell=1}>0$) or the
virtual level ($t^{\mathrm{tree},R}_{\ell=1}<0$). Meanwhile, the tree-level
$2\to 2$ scattering amplitude accidentally vanishes for some channel: for
example, higher antisymmetric representation in $SU(N_{f})$. In such a case,
the pole would appear at

$\begin{split}\beta_{\mathrm{pole}}^{3}=-\frac{1}{(t^{\mathrm{tree},R}_{\mathrm{in},\ell=1})^{2}}\,,\end{split}$
(32)

which is never positioned along the imaginary axis of $\beta$. Therefore, the
pole corresponds to neither bound state nor virtual level.

###  4.2 Application to dark pion SIMP

Now, we discuss the impact of the improvement on the SIMP models. Using the
tree-level amplitudes given by Eqs. 4 and 5, we obtain the improved amplitudes
for $\pi\pi\to\pi\pi$ and $\pi\pi\pi\to\pi\pi$ processes. The improved
amplitudes coincide with the tree-level amplitudes in the non-relativistic
limit, $\beta\to 0$. The improvement has a low impact on the pion $2\to 2$
self-scattering cross section in cosmic structures due to the DM velocity of
$\sim 10^{-2}$ at maximum. Meanwhile, this improvement procedure will be
efficient when the pions are semi-relativistic. The tree-level partial-wave
amplitudes scale as $T^{\mathrm{tree}}\propto(m_{\pi}/f_{\pi})^{2}$ and
$(T^{\mathrm{tree}}_{\mathrm{in}})^{2}\propto(m_{\pi}/f_{\pi})^{10}$, and
hence the inelastic scattering more significantly affects the improvement at
the large $m_{\pi}/f_{\pi}$.

Figure 3:  Trajectories of the improved partial-wave amplitude on the complex
plane of the amplitude for the elastic scattering as $m_{\pi}/f_{\pi}$
increased. The gray-dashed line corresponds to the unitarity circle (without
the inelastic scattering), whose center is at $(0,1/2\beta)$ and radius is
$1/2\beta$. We assume $SU(N_{f})\times SU(N_{f})/SU(N_{f})$ with $N_{f}=4$ and
$N_{c}=4$, and the adjoint representation for $R$. Each line corresponds to
fixed $v=\overline{k}/\overline{\mu}$: $v=0.1$ (brown), $v=0.7$ (cyan), and
$v=2.0$ (red). There are two markers on each line, which correspond to
$m_{\pi}/f_{\pi}=4\pi/\sqrt{N_{c}}$ and $4\pi$.

We plot the improved $p$-wave amplitude for the elastic scattering,
$T^{\mathrm{imp}}$, on a complex plane in Fig. 3 for the adjoint channel of
$SU(4)\times SU(4)/SU(4)$. Both elastic and inelastic scattering amplitudes
exist at the tree-level in this channel. In the figure, we take different
velocities $v=\overline{k}/\overline{\mu}$: $v=0.1$ (brown), $v=0.7$ (cyan),
and $v=2.0$ (red).777 One can estimate the maximum relative velocity of
$n$-body identical-particle collision as $v\leq n^{n/2(n-1)}$ by assuming that
each particle comes from respective vertex of a regular $n$-sided polygon with
the same velocity. In particular, the relative velocity can be $v\simeq 2.3$
in the three-body collision.  Once we specify the symmetry structure and the
channel, the improved amplitude is a function depending only on
$m_{\pi}/f_{\pi}$ and velocity $v$. As $m_{\pi}/f_{\pi}$ increases from zero
to infinity, the amplitude for fixed $v$ moves from bottom to top in the
figure. We put marks on each lines at $m_{\pi}/f_{\pi}=4\pi/\sqrt{N_{c}}$
($\blacktriangle$) and $4\pi$ ($\blacksquare$). The improved amplitude
$T^{\mathrm{imp}}$ is on the unitarity circle (with the radius of $1/2\beta$)
when the inelastic scattering is negligible. For small $m_{\pi}/f_{\pi}$,
therefore, each line follows the unitarity circle since the tree-level
partial-wave amplitude for the inelastic scattering is suppressed by a higher
power of $m_{\pi}/f_{\pi}$ than that for the elastic scattering.
$T^{\mathrm{imp}}_{\mathrm{in}}$ can be comparable to $T^{\mathrm{imp}}$ at a
certain point of $m_{\pi}/f_{\pi}$, and thus the radius from the center of the
circle shrinks as we see in Fig. 1. Since the tree-level partial-wave
amplitude for the inelastic scattering, $(T_{\mathrm{in,},\ell=1}^{R})^{2}$,
is proportional to $v^{8}$ as shown in Eq. 18, the path goes inside the circle
for smaller $m_{\pi}/f_{\pi}$ as the velocity $v$ gets larger. The improved
partial-wave amplitudes approach to $(0,1/\beta)$ for sufficiently large
$m_{\pi}/f_{\pi}$. As shown in Eq. 25, $T^{\mathrm{imp}}\to i/\beta$ and
$T^{\mathrm{imp}}_{\mathrm{in}}\to 0$ when the tree-level inelastic scattering
gets increased.

Figure 4:  Comparison of the $\pi\pi\pi\to\pi\pi$ (thermally-averaged) cross
sections. The symmetry structure is assumed to be $SU(4)\times SU(4)/SU(4)$,
and we take $m_{\pi}=0.1\,\mathrm{GeV}$. We take $m_{\pi}/f_{\pi}=4\pi$ in top
panels and $m_{\pi}/f_{\pi}=4\pi/\sqrt{N_{c}}$ in bottom panels. (Left): Shown
is the cross sections $\sigma v^{2}$ with improvement (cyan), the cross
section without improvement (orange), and unitarity bound (black) as a
function of velocity $v=\overline{k}/\overline{\mu}$. (Right): Shown is the
thermally-averaged cross section $\langle\sigma v^{2}\rangle$ with improvement
(cyan), without improvement (orange), and unitarity bound (black) as a
function of $x=m_{\pi}/T$.

We compare the (thermally-averaged) cross section with improvement (cyan),
that without improvement (orange), and the unitarity bound (black) in Fig. 4.
We choose different $m_{\pi}/f_{\pi}$ for the top panels and for the bottom
panels. We use the improved partial-wave amplitude
$T^{\mathrm{imp}}_{\mathrm{in}}$ in the $\pi\pi\pi\to\pi\pi$ cross section
given by Eq. 20. In the left panels, we show the cross section $\sigma_{3\to
2}v^{2}$ as the function of velocity $v=\overline{k}/\overline{\mu}$. Since
the tree-level partial-wave amplitude $(T^{\mathrm{tree}}_{\mathrm{in}})^{2}$
is proportional to $v^{8}$, the improved cross section follows the original
cross section at small $v$. The improved cross section gets larger as $v$
increases, and the cross section almost saturates the unitarity bound at a
certain $v$. The partial-wave amplitude for each representation does not
simultaneously saturate the unitarity bound at a certain $v$. As we see in
small windows in Fig. 4, the improved cross section does not fully coincide
with the unitarity bound, but approaches to the bound. For a larger $v$, the
improved partial-wave amplitude $(T_{\mathrm{in}})^{2}$ follows the scaling of
$v^{-8}$, and hence the improved cross section is suppressed compared to the
unitarity-saturated cross section, whose partial-wave amplitude follows the
scaling of $v^{0}$.

The thermally-averaged cross section with the improvement hardly saturates the
unitarity bound as shown in the right panels of Fig. 4 since $\sigma v^{2}$
with the improvement saturates the unitarity bound just in a narrow range of
$v$. A larger $v$ is required for saturating the unitarity bound for a smaller
$m_{\pi}/f_{\pi}$, and hence the improvement hardly affects the thermally-
averaged cross section for a smaller $m_{\pi}/f_{\pi}$ due to the Boltzmann
suppression factor. In particular, the tree-level inelastic amplitude is quite
suppressed at small $v$ in the dark pion realizations, which is proportional
to $v^{8}$, and hence the improvement has little impact on the cross section
for the small $v$. This is the main reasons why the improvement is less
effective even for the coupling close to the naïve perturbative bound
$m_{\pi}/f_{\pi}\simeq 4\pi/\sqrt{N_{c}}$. We note that there is a small
difference between (thermally-averaged) cross sections with and without
improvement even at small $v$ (large $x$), which originates from the
improvement by the $2\to 2$ amplitude not by the $2\to 3$ amplitude.

Figure 5:  Shown are $m_{\pi}$-$m_{\pi}/f_{\pi}$ plots. We assume different
symmetry structures for each panel: (left) $SU(N_{f})\times
SU(N_{f})/SU(N_{f})$, (middle) $SU(N_{f})/SO(N_{f})$, and (right)
$SU(N_{f})/USp(N_{f})$. We assume $N_{f}=4$ and $N_{c}=4$. The thick-cyan
lines show required $m_{\pi}/f_{\pi}$ for given $m_{\pi}$ (values shown on the
left-vertical axis): the thick-solid lines show the result using the improved
cross section, while the thick-dashed lines show the result without
improvement. The thin-red lines show the self-scattering cross section with
$m_{\pi}/f_{\pi}$ on thick lines (values shown on the right-vertical axis),
and each line-type corresponds to that of thick lines. The red-shaded bands
correspond to $0.1\,\mathrm{cm^{2}/g}\leq\sigma_{\mathrm{scatter}}/m_{\pi}\leq
1\,\mathrm{cm^{2}/g}$. The horizontal solid-lines show the unitary bound.

We illustrate the impact of the improvement of the $\pi\pi\pi\to\pi\pi$ cross
section in Fig. 5. We compute the relic abundance of the SIMP DM by the use of
an approximate formula, which ignores the temperature change of the effective
degrees of freedom $g_{\ast}$ during freeze-out, given by Ref. [79]:

$\begin{split}\Omega h^{2}\simeq\frac{1.05\times
10^{-10}\,\mathrm{GeV}^{-2}}{g_{\ast}^{3/4}(T_{f})m_{\pi}[K(x_{f})/M_{\mathrm{Pl}}]^{2}}\,,\qquad
K(x_{f})=\int_{x_{f}}^{\infty}\frac{dx}{x^{5}}\langle\sigma_{3\to
2}v^{2}\rangle\,.\end{split}$ (33)

Here, $M_{\mathrm{Pl}}=2.4\times 10^{18}\,\mathrm{GeV}$ and we assume
$x_{f}=m_{\pi}/T_{f}=20$. We use the effective degrees of freedom $g_{\ast}$
at the freeze-out temperature given by Ref. [80].

In Fig. 5, we compare the required $m_{\pi}/f_{\pi}$ for the collect relic
abundance $\Omega h^{2}\simeq 0.12$ with the thermally-averaged cross section
with and without the improvement (cyan lines): solid-thick lines for the
improved cross section, while dashed-thick lines for the cross section without
improvement. The self-scattering cross section
$\sigma_{\mathrm{scatter}}/m_{\pi}$, obtained along with $m_{\pi}/f_{\pi}$
required for the correct relic abundance, is shown as thin-red lines; these
line-types correspond to that of cyan lines, namely, solid (dashed) lines for
the cross section with (without) improvement. We take the number of flavor and
colors to be fixed as $N_{f}=4$ and $N_{c}=4$ in all panels. Each panel shows
the different symmetry structure. The horizontal solid-line shows the
unitarity bounds $m_{\pi}/f_{\pi}=4\pi$ and
$m_{\pi}/f_{\pi}=4\pi/\sqrt{N_{c}}$ given by Eq. 7, and the red-shaded band
indicates the self-scattering cross section in a range of
$0.1\,\mathrm{cm^{2}/g}\leq\sigma_{\mathrm{scatter}}/m_{\pi}\leq
1\,\mathrm{cm^{2}/g}$.

The required $m_{\pi}/f_{\pi}$ gets larger in order to explain the correct
relic abundance as the pion mass gets increased. The number-changing cross
section is suppressed by the improvement procedure compared to the tree-level
cross section for a large $m_{\pi}/f_{\pi}$, and hence the required
$m_{\pi}/f_{\pi}$ is larger than the case without the improvement. We find
that any large value of $m_{\pi}/f_{\pi}$ cannot explain the correct relic
abundance above $m_{\pi}\simeq 700\,\mathrm{MeV}$. However, we have to
incorporate the higher-order terms of the $\chi$PT when the pion self-coupling
gets close to the naïve perturbative bound.

We note that there seems to be no room for the allowed range of $m_{\pi}$ when
we take the conservative upper bound on the pion self-coupling
$m_{\pi}/f_{\pi}$ given by Eq. 7, which is suppressed by $\sqrt{N_{c}}$. The
improvement is less effective for the coupling close to the naïve perturbative
bound $m_{\pi}/f_{\pi}\simeq 4\pi/\sqrt{N_{c}}$. Meanwhile, this bound gets
mild when we take a weaker bound, e.g. $m_{\pi}/f_{\pi}\lesssim 4\pi$, and
then the allowed range of $m_{\pi}$ can be enlarged. As approaching to
$m_{\pi}/f_{\pi}\simeq 4\pi$, the improvement numerically modifies
$m_{\pi}/f_{\pi}$ at the ten-percent level.

##  5 Conclusion and Discussions

The dark pion is a promising realization of the SIMP models. The pion self-
coupling $m_{\pi}/f_{\pi}$ tends to be large in order that the SIMP framework
works, and hence, it is expected that scattering cross sections easily violate
the perturbative unitarity. In this study, we have proposed the improvement of
partial-wave amplitudes, which automatically hold the optical theorem. We have
found that the impact of the improvement on phenomenology is less significant
when the DM velocity is small. Therefore, the improvement of the cross section
does not change the self-interaction of dark pions at the cosmic structures.
Meanwhile, we have also examined the impact of the improvement on
$\pi\pi\pi\to\pi\pi$ processes, whose final states are semi-relativistic. The
pion self-coupling required to achieve the correct relic abundance is hardly
changed below the naïve perturbative bound $m_{\pi}/f_{\pi}\simeq
4\pi/\sqrt{N_{c}}$ by the improvement. The change will be more efficient near
the original NDA limit $m_{\pi}/f_{\pi}\simeq 4\pi$ [64], and it can be up to
$\mathcal{O}(10)\,\%$. The impact of the improvement will be comparable to
that of the NNLO correction of the $\chi$PT studied in Ref. [54], which is of
order of $10\,\%$. $m_{\pi}/f_{\pi}$ required for the correct relic abundance
drastically changes for $m_{\pi}\gtrsim 300\,\mathrm{MeV}$, and we cannot find
$m_{\pi}/f_{\pi}$ explaining the abundance for $m_{\pi}\gtrsim
700\,\mathrm{MeV}$. Since $m_{\pi}/f_{\pi}$ is beyond the naïve perturbative
bound in the mass range, there are uncertainties due to the resonances and the
higher-order corrections of the $\chi$PT.

Ref. [54] has discussed the impacts of the NLO and NNLO contributions of
chiral Lagrangian. They include only the real part of the amplitude since the
imaginary part provides only the higher-order contribution to the cross
section. On the other hand, their amplitude does not satisfy the optical
theorem. Our improved amplitudes incorporate the imaginary part, but consists
of the LO amplitudes of the elastic and inelastic scattering processes. We may
include the higher-order contributions systematically by using the inverse
amplitude method, which automatically satisfies the optical theorem, in the
absence of the inelastic channel [81]. However, it is not straightforward to
include inelastic channels in the method. This is beyond scope of this study.

In this study, we have focused only on the dark pion realization of the SIMP
scenario. The amplitude for the number-changing process is suppressed in the
following ways in the dark pion realization. The WZW term arises as the NLO
contributions with a one-loop suppression factor in the chiral perturbation.
Furthermore, the inelastic amplitude is suppressed by $v^{8}$ in the dark pion
realization, and it is a sub-leading contribution even in the $p$-wave
contributions. In other realizations, such as the dark glueball realization
[82, 83] and the vector SIMP model [84], the inelastic amplitude can be
leading-order contribution and is not suppressed by velocity. There,
unitarizing the amplitude will be more important within the naïve perturbative
range in other realizations.

### Acknowledgement

A. K. acknowledges partial support from Grant-in-Aid for Scientific Research
from the Ministry of Education, Culture, Sports, Science, and Technology
(MEXT), Japan, 18K13535 and 19H04609; from World Premier International
Research Center Initiative (WPI), MEXT, Japan; from Norwegian Financial
Mechanism for years 2014-2021, grant nr 2019/34/H/ST2/00707; and from National
Science Centre, Poland, grant 2017/26/E/ST2/00135 and DEC-2018/31/B/ST2/02283.
This work is also supported by the Advanced Leading Graduate Course for Photon
Science and the JSPS Research Fellowships for Young Scientists (S.K.).

##  Appendix A One-loop amplitude in $\lambda\phi^{4}$ theory

Let us consider a $\lambda\phi^{4}$ theory of a scalar field $\phi$ with a
mass of $m$. The tree-level amplitude for $\phi\phi\to\phi\phi$ is given by
$i\mathcal{M}^{\text{tree}}=-i\lambda$, while the one-loop amplitude is given
by

$\begin{split}i\mathcal{M}^{\text{1-loop}}&=-\frac{i\lambda^{2}}{32\pi^{2}}\int^{1}_{0}dx\left[\ln\frac{m^{2}-x(1-x)s-i\epsilon}{\mu^{2}}\right.\\\
&\qquad\left.+\ln\frac{m^{2}-x(1-x)t-i\epsilon}{\mu^{2}}+\ln\frac{m^{2}-x(1-x)u-i\epsilon}{\mu^{2}}\right]\,.\end{split}$
(34)

Here, $x$ denotes the Feynman parameter and $\mu$ denotes the renormalization
scale. The Mandelstam variables $s\,,t\,,$ and $u$ are given in the center-of-
mass frame by

$\begin{split}s=4m^{2}(1+\beta^{2}\gamma^{2})\,,\quad
t=-2m^{2}\beta^{2}\gamma^{2}(1-\cos\theta)\,,\quad
u=-2m^{2}\beta^{2}\gamma^{2}(1+\cos\theta)\,,\end{split}$ (35)

where $\beta^{2}=1-4m^{2}/s$ and $\gamma^{-2}=1-\beta^{2}$. A branch cut
appears along $\mathrm{Re}(s)>4m^{2}$.

We obtain the partial wave amplitudes by multiplying the Legendre polynomial
and integrating over the angle $\theta$. We focus only on the $s$-wave partial
amplitude. The tree-level amplitude is given by

$\begin{split}iT_{0}^{\text{tree}}\equiv\frac{1}{64\pi}\int^{1}_{-1}d\cos\theta(i\mathcal{M}^{\text{tree}})=-\frac{i\lambda}{32\pi}\,,\end{split}$
(36)

while the one-loop amplitude is

$\begin{split}iT_{0}^{\text{1-loop}}&\equiv\frac{1}{64\pi}\int^{1}_{-1}d\cos\theta(i\mathcal{M}^{\text{1-loop}})\\\
&=\frac{i\lambda^{2}}{512\pi^{3}}-\frac{i\lambda^{2}}{1024\pi^{3}}\int^{1}_{0}dx\left[\ln\frac{m^{2}-x(1-x)s-i\epsilon}{\mu^{2}}\right.\\\
&\qquad\left.+2\ln\frac{m^{2}}{\mu^{2}}+\frac{1+4x(1-x)\beta^{2}\gamma^{2}}{2x(1-x)\beta^{2}\gamma^{2}}\ln[1+4x(1-x)\beta^{2}\gamma^{2}]\right]\,.\end{split}$
(37)

Here, we follow the definition of the partial-wave amplitude Eq. 8. The one-
loop corrected amplitude Eq. 34 is totally symmetric under the exchange of
$s\,,t\,,$ and $u$. Therefore, the partial-wave amplitude has another branch
cut, which corresponds to branch cuts along $\mathrm{Re}(t)>4m^{2}$ or
$\mathrm{Re}(u)>4m^{2}$. In this region, $\mathrm{Re}(s)<0$, and hence
$\beta^{2}\gamma^{2}<0$. The another branch cut along $\mathrm{Re}(s)<0$
originates from the last term.

As far as we consider the physical region of the scattering process,
$\mathrm{Re}(s)>0$, the imaginary part of $T_{0}^{\text{1-loop}}$ appears only
from the first logarithmic term. The argument inside the logarithmic function
is negative only when $(1-\beta)/2<x<(1+\beta)/2$, and hence the imaginary
part of the one-loop partial-wave amplitude is

$\begin{split}\mathrm{Im}(T_{0}^{\text{1-loop}})&=-\frac{\lambda^{2}}{1024\pi^{3}}\int^{1}_{0}dx\mathrm{Im}\left[\ln\frac{m^{2}-x(1-x)s-i\epsilon}{\mu^{2}}\right]\\\
&=\frac{\lambda^{2}}{1024\pi^{2}}\sqrt{1-\frac{4m^{2}}{s}}=\beta(T_{0}^{\text{tree}})^{2}\,.\end{split}$
(38)

Here, we confirm the optical theorem Eq. 13 without the inelastic scattering
at the LO of the coupling $\lambda$ in the $\lambda\phi^{4}$ theory.

Finally, we discuss the coupling expansion of the total cross section and the
importance of the imaginary part. We find the total cross section up to the
one-loop level as follows:

$\begin{split}\sigma_{\mathrm{total}}&=\frac{32\pi}{s}|T_{0}^{\text{tree}}+T_{0}^{\text{1-loop}}|^{2}\\\
&\simeq\frac{\lambda^{2}}{32\pi
s}\left[1+\frac{2\lambda}{32\pi^{2}}+\beta^{2}\left(\frac{\lambda}{32\pi}\right)^{2}+\left(\frac{\lambda}{32\pi^{2}}\right)^{2}\right]\,.\end{split}$
(39)

Here, we ignore the logarithmic factors of the real part of the one-loop
correction in this expression. In the square bracket, the first term
corresponds to the tree-level contribution, the second is the interference
between the tree-level and the real part of the one-loop correction, the third
shows the imaginary part of the one-loop correction, and the last is the
contribution from the real part of the one-loop amplitude. We note that the
last-two terms must be affected by the interference between the tree-level and
2-loop amplitudes. All loop corrections similarly contribute to the cross
section when $16\pi^{2}\lesssim\lambda$, and hence the coupling expansion is
no longer validated. The improved amplitude in the text incorporates only the
imaginary part, and hence the improvement is expected to be important as the
coupling is just below the naïve perturbative bound.

##  Appendix B N-body Cross Sections and Thermal Averaging

In this appendix, we discuss the $N$-body cross section and its thermal
averaging.

###  B.1 Cross sections in the non-relativistic limit

We consider the $N$-body system, which is composed of $N$ particles with
momentum $\vec{p}_{i}$ and mass $m_{i}$ ($i=1,\cdots,N$). Jacobi coordinates
of the system are often useful to simplify the $N$-body system. When the
particle $i$ is positioned at $\vec{x}_{i}$, the Jacobi coordinates of the
system are defined by

$\begin{split}\vec{r}_{i}&=\frac{1}{m_{1}+\cdots+m_{i}}\left[\sum_{j=1}^{i}m_{j}\vec{x}_{j}-\sum_{j=1}^{i}m_{j}\vec{x}_{i+1}\right]\,,\;(i=1,\cdots,N-1)\,,\\\
\vec{r}_{N}&=\frac{1}{m_{1}+\cdots+m_{N}}\sum_{j=1}^{N}m_{j}\vec{x}_{j}\,.\end{split}$
(40)

Here, $\vec{r}_{i}\,(i=1,\cdots,N-1)$ represents the relative coordinates
between the position of the particle $i+1$ and the center-of-mass position for
the sub-system of $i$ particles. $\vec{r}_{N}$ denotes the center-of-mass
position of the total system. The total angular momentum We introduce Jacobi
coordinates of the momentum space, which we call the Jacobi momenta
$\vec{k}_{i}$ $(i=1,\cdots,N-1)$ and the center-of-mass momentum
$\vec{k}_{N}$.

$\begin{split}\vec{k}_{i}&=\frac{1}{m_{1}+\cdots+m_{i+1}}\left[m_{i+1}\sum_{j=1}^{i}\vec{p}_{j}-\sum_{j=1}^{i}m_{j}\vec{p}_{i+1}\right]\,,\;(i=1,\cdots,N-1)\,,\\\
\vec{k}_{N}&=\vec{p}_{1}+\cdots+\vec{p}_{N}\,.\end{split}$ (41)

The Jacobi momenta corresponds to the relative momentum of a particle $i+1$
and the center-of-mass of subsystem with $i$ particles. The total angular
momentum is written as the vector product of the Jacobi coordinates and Jacobi
momenta.

$\begin{split}\vec{L}\equiv\sum_{i=1}^{N}\vec{x}_{i}\times\vec{p}_{i}=\sum_{i=1}^{N}\vec{r}_{i}\times\vec{k}_{i}\,.\end{split}$
(42)

In terms of $\vec{k}_{i}$ $(i=1,\cdots,N)$, the total energy of the system is

$\begin{split}\sum_{i}\frac{\vec{p}_{i}^{2}}{2m_{i}}=\sum_{i}\frac{\vec{k}_{i}^{2}}{2\mu_{i}}\,.\end{split}$
(43)

Here, $\mu_{i}\,(i=1,\cdots,N-1)$ denotes the reduced mass of a particle $i+1$
and the center-of-mass of subsystem with $i$ particles,

$\begin{split}\frac{1}{\mu_{i}}=\frac{1}{m_{1}+\cdots+m_{i}}+\frac{1}{m_{i+1}}\,,\quad(i=1,\cdots,N-1)\,,\end{split}$
(44)

and $\mu_{N}$ denotes the total mass of all particles,
$\mu_{N}=m_{1}+\cdots+m_{N}$.

We compute the phase space integral by the use of the Jacobi momenta. The
Jacobian of this transformation is unity, and hence the $N$-body phase space
integral is computed as follows.

$\begin{split}d\Phi_{N}&\equiv\prod_{i=1}^{N}\frac{d^{3}p_{i}}{(2\pi)^{3}2m_{i}}(2\pi)^{4}\delta^{(4)}(P-\sum
p_{i})\\\
&=\prod_{i=1}^{N}\frac{d^{3}k_{i}}{(2\pi)^{3}2m_{i}}(2\pi)^{4}\delta(E-\sum\frac{k_{i}^{2}}{2\mu_{i}})\delta^{(3)}(\vec{P}-\vec{k}_{N})\,.\end{split}$
(45)

Here, $P$ denotes the total 4-momentum $P=(E,\vec{P})$. One of advantages in
using the Jacobi momenta is the simplification of the momentum integral with
the delta function. The delta function of 3-momentum is given as a function of
a single variable $\vec{k}_{N}$ by the use of the Jacobi momenta, and hence we
can easily carry out the $k_{N}$-integral. We take the center-of-mass frame of
$\vec{p}_{i}$, with $\vec{P}=\vec{0}$, in the following. We can simplify the
phase space integral by replacing the variables
$\vec{\kappa}_{i}=\sqrt{\overline{\mu}_{N}/\mu_{i}}\vec{k}_{i}$ with

$\begin{split}(\overline{\mu}_{N})^{N-1}=\mu_{1}\cdots\mu_{N-1}=\frac{m_{1}\cdots
m_{N}}{m_{1}+\cdots+m_{N}}\,.\end{split}$ (46)

The phase space integral takes the following form:

$\begin{split}d\Phi_{N}&=\prod_{i=1}^{N}\left(\frac{1}{2m_{i}}\right)\frac{d\overline{k}_{N}d\Omega_{3N-3}}{(2\pi)^{3N-4}}\overline{k}_{N}^{3N-4}\delta\left(E-\frac{\overline{k}_{N}^{2}}{2\overline{\mu}_{N}}\right)\,.\end{split}$
(47)

Here, $d\Omega_{d}$ denotes the integration over the solid angle in $d$
dimensions, and
$\overline{k}_{N}^{2}\equiv\kappa_{1}^{2}+\cdots+\kappa_{N-1}^{2}$. The solid
angle in $d$ dimensions is given by $\Omega_{d}=2\pi^{d/2}[\Gamma(d/2)]^{-1}$.
We note that the subscripts for $\overline{\mu}_{N}$ and $\overline{k}_{N}$
just denote the $N$-body phase space, not $N$-th variables.

Now, we consider 2-body elastic and inelastic cross sections of identical
scalar particles $\chi_{i}$ in the non-relativistic limit. The scattering
amplitudes for these processes are generally the functions of solid angles for
the initial states and final states, even though these are independent of the
angles for the 2-body scattering. Integrating the amplitude over the angle of
the final state, and averaging over the angle of the initial state, we obtain
the averaged cross sections [85]. For our purpose, we consider the $2\to 2$
elastic scattering and the $2\to N$ inelastic scattering. These cross sections
take the form

$\displaystyle\sigma_{2\to 2}$
$\displaystyle=\frac{2!}{\Omega_{3}}\left(\frac{2\pi}{\overline{k}_{2}}\right)^{2}\frac{1}{N_{\chi}^{2}}\sum_{R}\sum_{\ell}(2\ell+1)\times
4\beta^{2}|T^{R}_{\ell}|^{2}d_{R}\,,$ (48) $\displaystyle\sigma_{2\to N}$
$\displaystyle=\frac{2!}{\Omega_{3}}\left(\frac{2\pi}{\overline{k}_{2}}\right)^{2}\frac{1}{N_{\chi}^{2}}\sum_{R}\sum_{\ell}(2\ell+1)\times
4\beta^{2}|T_{\mathrm{in},\ell}^{R}|^{2}d_{R}\,.$ (49)

Here, $T^{R}_{\ell}$ and $|T_{\mathrm{in},\ell}^{R}|^{2}$ are the partial-wave
amplitude integrated over the phase space of the initial state. $N_{\chi}$
denotes the number of degrees of freedom of $\chi_{i}$.
$\beta=(1-4m^{2}/s)^{2}$ for the identical (up to flavor indices) particles
with the mass $m$, where $\sqrt{s}$ is the collision energy. The prefactor
originates from the two-body phase space integral including a factor for the
identical particles, $2!$. A factor of $4\beta^{2}$ in front of the partial-
wave amplitude originates from the definition of the $S$-matrix elements
$S_{\ell}$ for the partial-wave $\ell$, $S_{\ell}=1+2i\beta T_{\ell}$.

Once the system has the time-reversal invariance, the partial-wave amplitudes
for the $2\to N$ scattering and for the $N\to 2$ scattering are given by
$|T_{\mathrm{in},\ell}^{R}|^{2}$. In a similar way to the 2-body scattering,
we define the initial-angle averaged scattering cross section for $N\to 2$
process in the non-relativistic limit.

$\begin{split}\sigma_{N\to
2}v^{N-2}&=\frac{N!}{\Omega_{3N-3}}\left(\frac{2\pi}{\overline{k}_{N}}\right)^{3N-4}\frac{1}{N_{\chi}^{N}}\sum_{R}\sum_{\ell}(2\ell+1)\times
4\beta^{2}|T_{\mathrm{in},\ell}^{R}|^{2}d_{R}\,.\end{split}$ (50)

Here, the prefactor is due to averaging the $N$ initial states. We may define
the velocity for the $N$-body states as follows:
$v\equiv\overline{k}_{N}/\overline{\mu}_{N}$.

Now, we compute the thermally-averaged cross section. Another advantage of the
use of the Jacobi momenta is that the Jacobi momenta are invariant under the
Galilei transformation except for $\vec{k}_{N}$. Indeed, the momentum
variables of the thermal averaging do not necessarily coincide with those of a
scattering system, and hence the difference of the momentum variables does
affect the thermal averaging unless we consider only $s$-wave process. The
Jacobi momenta describe the relative motion of the subsystem except for
$\vec{k}_{N}$ describing the motion of the center-of-mass. The cross section
must be expressed only by the Jacobi momenta (except for $\vec{k}_{N}$) since
the cross section depends only on the relative motions of scattering
particles. The thermally-averaged cross section is defined as

$\begin{split}\langle\sigma_{N\to
2}v^{N-1}\rangle&\equiv\prod_{i=1}^{N}\frac{1}{n_{i}}\int\frac{d^{3}p_{i}}{(2\pi)^{3}2m_{i}}e^{-\frac{p_{i}^{2}}{2m_{i}T}}(\sigma_{N\to
2}v^{N-1})\,,\\\
n_{i}&\equiv\int\frac{d^{3}p_{i}}{(2\pi)^{3}2m_{i}}e^{-\frac{p_{i}^{2}}{2m_{i}T}}\,.\end{split}$
(51)

We substitute the integral variables $\vec{p}_{i}\to\vec{k}_{i}$, and the
Jacobian is unity as stated before. The integral of the cross section with
respect to $\vec{k}_{N}$ is trivially carried out since the cross section does
not depend on $\vec{k}_{N}$. We may simplify the integral again by replacing
the variables $\vec{k}_{i}\,(i=1,\cdots,N-1)$ to $\overline{k}_{N}$ and
hyperangles as shown before, and the phase space integral of the initial state
may be rewritten as

$\begin{split}\langle\sigma_{N\to
2}v^{N-1}\rangle&=\frac{1}{n}\prod_{i=1}^{N}\left(\frac{1}{2m_{i}}\right)\int\frac{d\overline{k}_{N}d\Omega_{3N-3}}{(2\pi)^{3N-4}}\overline{k}_{N}^{3N-4}e^{-\frac{\overline{k}_{N}^{2}}{2\overline{\mu}_{N}T}}(\sigma_{N\to
2}v^{N-1})\,,\\\
n&\equiv\prod_{i=1}^{N}\left(\frac{1}{2m_{i}}\right)\int\frac{d\overline{k}_{N}d\Omega_{3N-3}}{(2\pi)^{3N-4}}\overline{k}_{N}^{3N-4}e^{-\frac{\overline{k}_{N}^{2}}{2\overline{\mu}_{N}T}}\,.\end{split}$
(52)

Here, averaging over the angle variables gives the initial-angle averaged
cross section (50). Taking the partial-wave amplitudes to be the maximum value
$|T_{\mathrm{in},\ell}^{R}|^{2}\leq(1/2\beta)^{2}$, we obtain the upper bound
of the thermally-averaged cross section for $N\to 2$ process.

$\begin{split}\langle\sigma_{N\to
2}v^{N-1}\rangle_{\mathrm{uni}}=\frac{\left(2\pi
x\right)^{\frac{3N-5}{2}}}{m_{\chi}^{3N-4}}\frac{N\sqrt{N}}{N_{\chi}^{N}}\sum_{R}\sum_{\ell}(2\ell+1)d_{R}\,.\end{split}$
(53)

Here, $R$ and $\ell$ run only over representations whose amplitudes are non-
zero, $x=m_{\chi}/T$, and we take the identical mass to be $m_{\chi}$. We
again note that we independently obtain the upper bounds from Refs. [74, 75,
76]. As we have stated in the text, our result is consistent with Refs. [74,
75] up to a factor of $1/4$ due to maximizing the partial-wave amplitude, but
our result disagrees with the results in Ref [76] by a factor of
$N\sqrt{N}/2\sqrt{2}$.

###  B.2 Application to the SIMP models

The tree-level five-pion amplitude is given by Eq. 17 in the text. We compute
the thermally-averaged cross section of the $3\to 2$ process. First, we
compute the $2\to 3$ cross section in the non-relativistic limit, and then we
obtain the $3\to 2$ cross section by the use of the relation between cross
sections Eqs. 49 and 50 (with $N=3$). We take the center-of-mass frame of the
initial pions of the $2\to 3$ process, $\vec{p}_{4}+\vec{p}_{5}=\vec{0}$, and
the collision energy is assumed to be $\sqrt{s}$, and thus the 4-momentum of
the initial pion is given by

$\begin{split}p_{4}=\frac{\sqrt{s}}{2}(1,\vec{\beta})\,,\qquad
p_{5}=\frac{\sqrt{s}}{2}(1,-\vec{\beta})\,,\end{split}$ (54)

with $|\vec{\beta}|=(1-4m_{\pi}^{2}/s)^{1/2}$. The amplitude takes the form

$\begin{split}\mathcal{M}=-\frac{16k}{\pi^{2}f_{\pi}^{5}}T_{[abcde]}\frac{s}{2}\vec{\beta}\cdot(\vec{p}_{1}\times\vec{p}_{2})\,.\end{split}$
(55)

We carry out the three-body phase-space integral for the final-state pions to
obtain the $2\to 3$ cross section. We introduce a practical technique to
compute the six-dimensional hyperspherical integral.
$\overline{k}_{3}^{3}=\kappa_{1}^{2}+\kappa_{2}^{2}$ since we take the center-
of-mass frame. We may parametrize each momentum as follows.

$\begin{split}\vec{\kappa}_{1}&=\overline{k}_{3}\cos\theta(\sin\theta_{1}\cos\phi_{1}\,,\sin\theta_{1}\sin\phi_{1}\,,\cos\theta_{1})\,,\\\
\vec{\kappa}_{2}&=\overline{k}_{3}\sin\theta(\sin\theta_{2}\cos\phi_{2}\,,\sin\theta_{2}\sin\phi_{2}\,,\cos\theta_{2})\,.\end{split}$
(56)

Here, $(\theta_{i}\,,\phi_{i})$ denotes the solid angle for the momentum
$\vec{\kappa}_{i}$, and $\theta$ distributes $\overline{k}_{3}$ to
$\kappa_{1}$ and $\kappa_{2}$. This assignment of angles is called as the
Delves coordinate [86, 87]. The integral measure is

$\begin{split}\int
d\Omega_{6}=\int^{\pi/2}_{0}d\theta\sin^{2}\theta\cos^{2}\theta
d\Omega^{1}_{3}d\Omega^{2}_{3}\,,\end{split}$ (57)

where $d\Omega^{i}_{3}$ denotes the solid-angle integral for the momentum
$\vec{\kappa}_{i}$. The three-body phase space integral in the non-
relativistic limit takes the form:

$\begin{split}d\Phi_{3}&=\prod_{i=1}^{3}\left(\frac{1}{2m_{i}}\right)\frac{1}{(2\pi)^{5}}\overline{\mu}_{3}\overline{k}_{3}^{4}d\theta\sin^{2}\theta\cos^{2}\theta
d\Omega^{1}_{3}d\Omega^{2}_{3}\,.\end{split}$ (58)

Here, $\overline{k}_{3}^{2}=2\overline{\mu}_{3}E$. Meanwhile, the relativistic
form of the three-body phase space integral is given by

$\begin{split}d\Phi_{3}&=\frac{ds_{23}}{16}\frac{1}{(2\pi)^{5}}d\Omega^{1}_{3}d\Omega^{2}_{3}\overline{\beta}\left(\frac{m_{1}^{2}}{s}\,,\frac{s_{23}}{s}\right)\overline{\beta}\left(\frac{m_{2}^{2}}{s_{23}}\,,\frac{m_{3}^{2}}{s_{23}}\right)\,,\\\
\overline{\beta}\left(a,b\right)&\equiv\frac{1}{2}\sqrt{1-2(a+b)+(a-b)^{2}}\,.\end{split}$
(59)

where $s_{23}=(p_{2}+p_{3})^{2}$. The relativistic counterpart of the Delves
angle $\theta$ is the Dalitz variable $s_{23}$. $s_{23}$ takes a value within
the range of $(m_{2}+m_{3})^{2}\leq s_{23}\leq(\sqrt{s}-m_{1})^{2}$. We may
consider the non-relativistic limit of the integral by parameterizing as
$\sqrt{s_{23}}=m_{2}+m_{3}+E\cos^{2}\theta$ with
$E=\sqrt{s}-(m_{1}+m_{2}+m_{3})$ being the center-of-mass energy in the limit.
We expand $\overline{\beta}$ and $ds_{23}$ in small $E$ limit as follows.

$\begin{split}\overline{\beta}\left(\frac{m_{1}^{2}}{s}\,,\frac{s_{23}}{s}\right)\overline{\beta}\left(\frac{m_{2}^{2}}{s_{23}}\,,\frac{m_{3}^{2}}{s_{23}}\right)&\simeq\sin\theta\cos\theta\frac{2\overline{\mu}_{3}^{3}E}{m_{1}m_{2}m_{3}(m_{2}+m_{3})}\,,\\\
ds_{23}&\simeq 4E(m_{2}+m_{3})d\theta\sin\theta\cos\theta\,,\end{split}$ (60)

and therefore, we obtain the non-relativistic form of the phase-space integral
given by Eq. 58 by replacing $E=(\overline{k}_{3}^{2}/2\overline{\mu}_{3})$.
The cross section for the $2\to 3$ in the non-relativistic limit is

$\begin{split}\sigma_{2\to
3}&=\frac{1}{3!}\frac{1}{2s\beta}\frac{1}{N_{\pi}^{2}}\int
d\Phi_{3}|\mathcal{M}|^{2}=\frac{1}{3!}\frac{1}{(2\pi)^{5}}\frac{k^{2}}{6\pi}\frac{\overline{k}_{3}^{8}\overline{\mu}_{3}s\beta}{m_{\pi}^{3}f_{\pi}^{10}}\frac{1}{N_{\pi}^{2}}\sum_{R}t_{R}^{2}\,.\end{split}$
(61)

Here, $t_{R}^{2}$ is the decomposition of $t^{2}$ into each representations.
We carry out the three-body phase-space integral by following Eq. 58. By
comparing the cross section with Eq. 49 (with $N=3$), we find the partial-wave
amplitude (18) for the inelastic scattering in the text.

The thermal averaging of the $3\to 2$ annihilation is defined by Eq. 52 (with
$N=3$). Since the partial-wave amplitude is averaged over hyperangles, the
hyperangle integral is trivially carried out and only
$\overline{k}_{3}$-integral remains. We find the thermally-averaged cross
section as follows.

$\begin{split}\langle\sigma_{3\to
2}v^{2}\rangle=\frac{5\sqrt{5}k^{2}m_{\pi}^{5}}{6\pi^{5}x^{2}f_{\pi}^{10}}\frac{1}{N_{\pi}^{3}}\sum_{R}t_{R}^{2}\,.\end{split}$
(62)

Our result is smaller than the cross section in the original paper [51] by a
factor of $1/3$. The difference could arise from the Galilei transformation of
the invariant amplitude Eq. 55, which is computed in the center-of-mass frame.
$\vec{p}_{1}\times\vec{p}_{2}$ is not Galilean invariant, and in fact it
transforms as
$\vec{p}_{1}\times\vec{p}_{2}\to\vec{p}_{1}\times\vec{p}_{2}+m_{2}\vec{p}_{1}\times\vec{v}+m_{1}\vec{v}\times\vec{p}_{2}$
under $\vec{p}_{i}\to\vec{p}_{i}+m_{i}\vec{v}$ with $\vec{v}$ being a constant
velocity. The corresponding part may be rewritten in terms of the Jacobi
momentum as follows.

$\begin{split}\vec{\beta}\cdot(\vec{p}_{1}\times\vec{p}_{2})=\vec{\beta}\cdot(\vec{k}_{1}\times\vec{k}_{2})+\frac{2}{3}\vec{\beta}\cdot(\vec{k}_{1}\times\vec{k}_{3})\,.\end{split}$
(63)

As we mentioned, $\vec{k}_{3}=0$ in the center-of-mass frame. Hence, the
second term vanishes since the amplitude (55) is computed in the center-of-
mass frame, and then we obtain Eq. 62. Meanwhile, we reproduce the result in
Ref. [51] by keeping the second term artificially and by keeping $k_{3}$
integral with the Maxwell distribution (with $\mu_{3}$). We note again that it
is important to rewrite the cross section in terms of the Jacobi momenta,
namely the Galilean-invariant way, especially for computing thermally-averaged
cross section.

##  Appendix C Group-theoretical Computation

We consider the irreducible decomposition of products of broken generators of
$G/H$. It is convenient to use the $SU(N_{f})$ generators for describing the
broken generators even for $G/H=SU(N_{f})/SO(N_{f})$ and
$SU(N_{f})/USp(N_{f})$. The fundamental generators are denoted by $T^{a}$,
whose normalization is defined by

$\begin{split}\mathrm{tr}(T^{a}T^{b})=\frac{1}{2}\delta^{ab}\,.\end{split}$
(64)

The Fierz identity for the generator is

$\begin{split}(T^{a})_{i}^{~{}j}(T^{a})_{k}^{~{}l}=\frac{1}{2}\left(\delta_{i}^{l}\delta_{k}^{j}-\frac{1}{N_{f}}\delta_{i}^{j}\delta_{k}^{l}\right)\,.\end{split}$
(65)

The commutation and anticommutation of the generators gives the group-
theoretical constants, $d$ and $f$ as

$\begin{split}[T^{a},T^{b}]=if^{abc}T^{c}\,,\qquad\\{T^{a},T^{b}\\}=d^{abc}T^{c}+\frac{1}{N_{f}}\delta^{ab}\,,\end{split}$
(66)

and we use these constants for the short-hand notation.

###  C.1 $SU(N_{f})\times SU(N_{f})/SU(N_{f})$

First, let us consider $G=SU(N_{f})\times SU(N_{f})$ and $H=SU(N_{f})$. The
broken generator is the adjoint representation of $SU(N_{f})$, and the direct
product of the generator is symbolically given by

$\begin{split}\mathrm{Ad}\otimes\mathrm{Ad}=1_{S}\oplus\mathrm{Ad}_{A}\oplus\mathrm{Ad}_{S}\oplus
T_{A}\oplus\overline{T}_{A}\oplus Y_{S}\oplus S_{S}\,.\end{split}$ (67)

Here, $\mathrm{Ad}$ denotes the adjoint representation of $SU(N_{f})$. The
subscripts $S$ and $A$ denote symmetrization and anti-symmetrization of
indices of adjoint representations in the left-hand side. Other
representations are summarized in Table 1. As seen in the table, several
representations vanish for $N_{f}\leq 3$. For $N_{f}=2$, $\mathrm{Ad}_{S}$,
$T_{A}\oplus\overline{T}_{A}$, and $Y_{S}$ do not appear and the decomposition
is given by

$\begin{split}\mathbf{3}\otimes\mathbf{3}=\mathbf{1}_{S}\oplus\mathbf{3}_{A}\oplus\mathbf{5}_{S},\end{split}$
(68)

and for $N_{f}=3$, $Y_{S}$ does not appear and the decomposition is given by

$\begin{split}\mathbf{8}\otimes\mathbf{8}=\mathbf{1}_{S}\oplus\mathbf{8}_{A}\oplus\mathbf{8}_{S}\oplus\mathbf{10}_{A}\oplus\overline{\mathbf{10}}_{A}\oplus\mathbf{27}_{S}\,.\end{split}$
(69)

Table 1: Representations in $SU(N_{f})$ and their dimensions  Reps. | Dims.
---|---
$\mathrm{Ad}$ | $N_{f}^{2}-1$
$T(\overline{T})$ | $\dfrac{1}{4}(N_{f}^{2}-1)(N_{f}^{2}-4)$
$Y$ | $\dfrac{1}{4}(N_{f}+1)N_{f}^{2}(N_{f}-3)$
$S$ | $\dfrac{1}{4}N_{f}^{2}(N_{f}+3)(N_{f}-1)$

The projection operators with adjoint indices are defined by

$\begin{split}P_{R}^{ab;cd}=4(T^{a})_{j}^{~{}i}(T^{b})_{l}^{~{}k}(P_{R})^{jl;i^{\prime}k^{\prime}}_{ik;j^{\prime}l^{\prime}}(T^{c})_{i^{\prime}}^{~{}j^{\prime}}(T^{d})_{k^{\prime}}^{~{}l^{\prime}}\,.\end{split}$
(70)

Here, the subscript $R$ stands for the representations, and a factor of 4 is
for the Dynkin normalization.
$(P_{R})^{jl;i^{\prime}k^{\prime}}_{ik;j^{\prime}l^{\prime}}$ denotes the
projection operators for the fundamental indices, which reflect the symmetries
of corresponding Young tableaux. Now, we decompose the product
$(T^{a})_{j}^{~{}i}(T^{b})_{l}^{~{}k}$ into its irreducible representations,
which corresponds to
$(T^{a})_{j}^{~{}i}(T^{b})_{l}^{~{}k}(P_{R})^{jl;i^{\prime}k^{\prime}}_{ik;j^{\prime}l^{\prime}}$.
We use the short-hand notation $a_{j}^{~{}i}\equiv(T^{a})_{j}^{~{}i}$. We
define the traceless part of the product $a_{j}^{~{}i}b_{l}^{~{}k}$ as
follows.

$\begin{split}\tilde{a}_{j}^{~{}i}\tilde{b}_{l}^{~{}k}&\equiv
a_{j}^{~{}i}b_{l}^{~{}k}-\frac{1}{(N_{f}^{2}-4)(N_{f}^{2}-1)}\left[3N_{f}\delta_{j}^{~{}i}\delta_{l}^{~{}k}-(N_{f}^{2}+2)\delta_{j}^{~{}k}\delta_{l}^{~{}i}\right]\mathrm{tr}(ab)\\\
&+\frac{1}{(N_{f}^{2}-4)}\left[\delta_{j}^{~{}i}(a_{l}^{~{}m}b_{m}^{~{}k}+a_{m}^{~{}k}b_{l}^{~{}m})+\delta_{l}^{~{}k}(a_{j}^{~{}m}b_{m}^{~{}i}+a_{m}^{~{}i}b_{j}^{~{}m})\right]\\\
&-\frac{2}{N_{f}(N_{f}^{2}-4)}(\delta_{l}^{~{}i}a_{m}^{~{}k}b_{k}^{~{}m}+\delta_{j}^{~{}k}a_{l}^{~{}m}b_{m}^{~{}i})-\frac{N_{f}^{2}-2}{N_{f}(N_{f}^{2}-4)}(\delta_{l}^{~{}i}a_{j}^{~{}m}b_{m}^{~{}k}+\delta_{j}^{~{}k}a_{m}^{~{}i}b_{l}^{~{}m})\,.\end{split}$
(71)

Here, the trace is given by $\mathrm{tr}(ab)=a_{j}^{~{}i}b_{i}^{~{}j}$. Using
this quantity, we can decompose $a_{j}^{~{}i}b_{l}^{~{}k}$ into their
irreducible parts. We introduce $(P_{R}^{ab})_{jl}^{ik}$ as a projection of
$a_{j}^{~{}i}b_{l}^{~{}k}$ into its irreduceble representation $R$

$\begin{split}a_{j}^{~{}i}b_{l}^{~{}k}=\sum_{R}(P_{R}^{ab})_{jl}^{ik}\,,\end{split}$
(72)

and each representation takes the form

$\displaystyle(P_{S}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{4}(\tilde{a}_{j}^{~{}i}\tilde{b}_{l}^{~{}k}+\tilde{a}_{j}^{~{}k}\tilde{b}_{l}^{~{}i}+\tilde{a}_{l}^{~{}i}\tilde{b}_{k}^{~{}k}+\tilde{a}_{l}^{~{}k}\tilde{b}_{j}^{~{}i})\,,$
(73) $\displaystyle(P_{Y}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{4}(\tilde{a}_{j}^{~{}i}\tilde{b}_{l}^{~{}k}-\tilde{a}_{j}^{~{}k}\tilde{b}_{l}^{~{}i}-\tilde{a}_{l}^{~{}i}\tilde{b}_{k}^{~{}k}+\tilde{a}_{l}^{~{}k}\tilde{b}_{j}^{~{}i})\,,$
(74) $\displaystyle(P_{T}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{4}(\tilde{a}_{j}^{~{}i}\tilde{b}_{l}^{~{}k}+\tilde{a}_{j}^{~{}k}\tilde{b}_{l}^{~{}i}-\tilde{a}_{l}^{~{}i}\tilde{b}_{k}^{~{}k}-\tilde{a}_{l}^{~{}k}\tilde{b}_{j}^{~{}i})\,,$
(75) $\displaystyle(P_{\overline{T}}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{4}(\tilde{a}_{j}^{~{}i}\tilde{b}_{l}^{~{}k}-\tilde{a}_{j}^{~{}k}\tilde{b}_{l}^{~{}i}+\tilde{a}_{l}^{~{}i}\tilde{b}_{k}^{~{}k}-\tilde{a}_{l}^{~{}k}\tilde{b}_{j}^{~{}i})\,,$
(76) $\displaystyle(P_{\mathrm{Ad}_{S}}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{2(N_{f}^{2}-4)}\biggl{[}N_{f}\delta_{l}^{~{}i}\\{a,b\\}_{j}^{~{}k}+N_{f}\delta_{j}^{~{}k}\\{a,b\\}_{l}^{~{}i}-\delta_{j}^{~{}i}\\{a,b\\}_{l}^{~{}k}-\delta_{l}^{~{}k}\\{a,b\\}_{j}^{~{}i}$
$\displaystyle\qquad\left.-4\delta_{l}^{~{}i}\delta_{j}^{~{}k}\mathrm{tr}(ab)+\frac{4}{N_{f}}\delta_{j}^{~{}i}\delta_{l}^{~{}k}\mathrm{tr}(ab)\right]\,,$
(77) $\displaystyle(P_{\mathrm{Ad}_{A}}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{2N_{f}}\left(\delta_{l}^{~{}i}[a,b]_{j}^{~{}k}-\delta_{j}^{~{}k}[a,b]_{l}^{~{}i}\right)\,,$
(78) $\displaystyle(P_{1}^{ab})_{jl}^{ik}$
$\displaystyle=\frac{1}{N_{f}^{1}-1}\left(\delta_{l}^{~{}i}\delta_{j}^{~{}k}-\frac{1}{N_{f}}\delta_{j}^{~{}i}\delta_{l}^{~{}k}\right)\mathrm{tr}(ab)\,.$
(79)

Here, $\\{a,b\\}_{j}^{~{}i}\equiv
a_{j}^{~{}m}b_{m}^{~{}i}+b_{j}^{~{}m}a_{m}^{~{}i}$ and $[a,b]_{j}^{~{}i}\equiv
a_{j}^{~{}m}b_{m}^{~{}i}-b_{j}^{~{}m}a_{m}^{~{}i}$. By multiplying
$4c_{i}^{~{}j}d_{k}^{~{}l}$, we obtain the projection operators with indices
of pions, $P_{R}^{ab;cd}$. The projection operators are

$\displaystyle P_{1}^{ab;cd}$
$\displaystyle=\frac{1}{N_{f}^{2}-1}\delta^{ab}\delta^{cd}\,,$ (80)
$\displaystyle P_{\mathrm{Ad}_{A}}^{ab;cd}$
$\displaystyle=\frac{1}{N_{f}}f^{abe}f^{cde}\,,$ (81) $\displaystyle
P_{\mathrm{Ad}_{S}}^{ab;cd}$
$\displaystyle=\frac{N_{f}}{N_{f}^{2}-4}d^{abe}d^{cde}\,,$ (82) $\displaystyle
P_{T}^{ab;cd}$
$\displaystyle=\frac{1}{4}(\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc})-\frac{1}{2N_{f}}f^{abe}f^{cde}$
$\displaystyle\qquad+\frac{i}{8}(f^{ace}d^{bde}+d^{ace}f^{bde}-f^{ade}d^{bce}-d^{ade}f^{bce})\,,$
(83) $\displaystyle P_{\overline{T}}^{ab;cd}$
$\displaystyle=\frac{1}{4}(\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc})-\frac{1}{2N_{f}}f^{abe}f^{cde}$
$\displaystyle\qquad-\frac{i}{8}(f^{ace}d^{bde}+d^{ace}f^{bde}-f^{ade}d^{bce}-d^{ade}f^{bce})\,,$
(84) $\displaystyle P_{Y}^{ab;cd}$
$\displaystyle=\frac{N_{f}-2}{4N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})+\frac{N_{f}-2}{2N_{f}(N_{f}-1)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad-\frac{1}{4}(d^{ace}d^{bde}+d^{ade}d^{bce})+\frac{N_{f}-4}{4(N_{f}-2)}d^{abe}d^{cde}\,,$
(85) $\displaystyle P_{S}^{ab;cd}$
$\displaystyle=\frac{N_{f}+2}{4N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})-\frac{N_{f}+2}{2N_{f}(N_{f}+1)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad+\frac{1}{4}(d^{ace}d^{bde}+d^{ade}d^{bce})-\frac{N_{f}+4}{4(N_{f}+2)}d^{abe}d^{cde}\,.$
(86)

These projection operators satisfy

$\begin{split}\sum_{R}P_{R}^{ab;cd}=\delta^{ac}\delta^{bd}\,,\qquad\sum_{c^{\prime},d^{\prime}}P_{R}^{ab;c^{\prime}d^{\prime}}P_{R^{\prime}}^{c^{\prime}d^{\prime};cd}=\delta_{RR^{\prime}}P_{R}^{ab;cd}\,.\end{split}$
(87)

To verify the second equality, we use several identities with $d^{abc}$ and
$f^{abc}$, such as

$\begin{split}f^{acd}f^{bcd}&=-N_{f}\delta^{ab}\,,\qquad
d^{acd}d^{bcd}=\frac{N_{f}^{2}-4}{N_{f}}\delta^{ab}\,,\qquad\\\
f^{abe}f^{cde}&=d^{ace}d^{bde}-d^{ade}d^{bce}+\frac{2}{N_{f}}(\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc})\,,\\\
f^{abe}d^{cde}&=d^{ace}f^{bde}+d^{ade}f^{bce}\,,\\\
f^{abe}f^{cde}&=-f^{ace}f^{dbe}-f^{ade}f^{bce}\,,\end{split}$ (88)

and

$\begin{split}\mathrm{Tr}(F^{a}F^{b}F^{c}F^{d})&=\delta^{ab}\delta^{cd}+\delta^{ad}\delta^{bc}+\frac{N_{f}}{4}(d^{abe}d^{cde}-d^{ace}d^{bde}+d^{ade}d^{bce})\,,\\\
\mathrm{Tr}(F^{a}F^{b}D^{c}D^{d})&=\frac{N_{f}^{2}-4}{N_{f}^{2}}(\delta^{ab}\delta^{cd}-\delta^{ac}\delta^{bd})+\frac{N_{f}^{2}-8}{4N_{f}}(d^{abe}d^{cde}-d^{ace}d^{bde})+\frac{N_{f}}{4}d^{ade}d^{bce}\,,\\\
\mathrm{Tr}(D^{a}D^{b}D^{c}D^{d})&=\frac{N_{f}^{2}-4}{N_{f}^{2}}(\delta^{ab}\delta^{cd}+\delta^{ad}\delta^{bc})+\frac{N_{f}^{2}-16}{4N_{f}}(d^{abe}d^{cde}+d^{ade}d^{bce})+\frac{N_{f}}{4}d^{ace}d^{bde}\,,\end{split}$
(89)

where $(F^{a})^{bc}=-if^{abc}$ and $(D^{a})^{bc}=d^{abc}$, and $\mathrm{Tr}$
denotes a trace of adjoint indices. The first equation in Eq. 89 is derived
from computing $\mathrm{tr}([[[[T^{e},T^{a}],T^{b}],T^{c}],T^{d}]T^{e})$ in
two different ways: one way is to use the commutation relation, and another
way is to reduce the number of generators by use of the Fierz identity and
directly compute the traces of four generators. More complicated relations
would be required for the decomposition of three products of generators, which
also requires the dedicated analysis for the resummation of $3\to 2$ initial
states. The trace of projection operators gives the dimension of the
representation,

$\begin{split}\sum_{a,b,c,d}P_{R}^{ab;cd}P_{R}^{cd;ab}=d_{R}\,,\end{split}$
(90)

where $d_{R}$ refers to the dimensions of the representation $R$ (see Table 1
for the dimensions).

###  C.2 $SU(N_{f})/SO(N_{f})$

Let us consider the symmetry breaking $SU(N_{f})\to SO(N_{f})$. The number of
unbroken generators is $N_{f}(N_{f}-1)/2$, and hence the number of broken
generators is $(N_{f}+2)(N_{f}-1)/2$. The broken generator satisfies
$T^{a}=(T^{a})^{T}$, and hence the pions transform as a second rank symmetric
tensor field. The product of the symmetric tensors is decomposed as follows.

$\begin{split}s\otimes s=1_{S}\oplus\mathrm{Ad}_{A}\oplus s_{S}\oplus
T_{A}\oplus Y_{S}\oplus S_{S}\,.\end{split}$ (91)

Here, the subscripts $S$ and $A$ denote symmetrization and anti-symmetrization
of indices of symmetric representations in the left-hand side. Table 2 shows
the representations and their dimensions. As seen in the table, several
representations vanish for $N_{f}\leq 3$. For $N_{f}=3$, $Y_{S}$ does not
appear and the decomposition is given by

$\begin{split}\mathbf{5}\otimes\mathbf{5}=\mathbf{1}_{S}\oplus\mathbf{3}_{A}\oplus\mathbf{5}_{S}\oplus\mathbf{7}_{S}\oplus\mathbf{9}_{S}\,.\end{split}$
(92)

Table 2:  Representations in $SO(N_{f})$ and their dimensions  Reps. | Dims.
---|---
$s$ | $\dfrac{1}{2}(N_{f}-1)(N_{f}+2)$
$\mathrm{Ad}$ | $\dfrac{1}{2}N_{f}(N_{f}-1)$
$Y$ | $\dfrac{1}{12}(N_{f}-3)N_{f}(N_{f}+1)(N_{f}+2)$
$S$ | $\dfrac{1}{24}(N_{f}^{2}-1)N_{f}(N_{f}+6)$
$T$ | $\dfrac{1}{8}(N_{f}^{2}-1)(N_{f}+4)(N_{f}-2)$

It is convenient to make the upper indices of the broken generator lower since
the broken generators are symmetric: $(T^{a})_{ij}\equiv(T^{a})_{i}^{~{}j}$.
The Fierz identity for the broken generator is given by

$\begin{split}(T^{a})_{ij}(T^{a})_{kl}=\frac{1}{4}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}-\frac{2}{N_{f}}\delta_{ij}\delta_{kl}\right)\,.\end{split}$
(93)

Here, the summation over $a$ runs only over the broken generators. We
construct the projection operators for the product of broken generators as
with the $SU(N_{f})$. The projection operators with the labels of the broken
generators are defined by

$\begin{split}P_{R}^{ab;cd}=4(T^{a})_{ij}(T^{b})_{kl}(P_{R})_{ijkl;i^{\prime}j^{\prime}k^{\prime}l^{\prime}}(T^{c})_{i^{\prime}j^{\prime}}(T^{d})_{k^{\prime}l^{\prime}}\,.\end{split}$
(94)

Now, we decompose the product $(T^{a})_{ij}(T^{b})_{kl}$ into its irreducible
representations, which corresponds to
$(T^{a})_{ij}(T^{b})_{kl}(P_{R})_{ijkl;i^{\prime}j^{\prime}k^{\prime}l^{\prime}}$.
We use the short-hand notation $a_{ij}\equiv(T^{a})_{ij}$. We define the
traceless part of the product $a_{ij}b_{kl}$ as follows.

$\begin{split}&\tilde{a}_{ij}\tilde{b}_{kl}\equiv a_{ij}b_{kl}\\\
&\quad-\frac{1}{(N_{f}^{2}-4)(N_{f}-1)(N_{f}+4)}\left[2(3N_{f}+2)\delta_{ij}\delta_{kl}-(N_{f}^{2}+4)(\delta_{il}\delta_{jk}+\delta_{ik}\delta_{jl})\right]\mathrm{tr}(ab)\\\
&\quad+\frac{2}{(N_{f}+4)(N_{f}-2)}(\delta_{ij}a_{km}b_{ml}+\delta_{ij}a_{lm}b_{mk}+\delta_{kl}a_{im}b_{mj}+\delta_{kl}a_{jm}b_{mi})\\\
&\quad-\frac{4}{(N_{f}^{2}-4)(N_{f}+4)}(\delta_{il}a_{km}b_{mj}+\delta_{jl}a_{km}b_{mi}+\delta_{ik}a_{lm}b_{mj}+\delta_{jk}a_{lm}b_{mi})\\\
&\quad-\frac{N_{f}^{2}+2N_{f}-4}{(N_{f}^{2}-4)(N_{f}+4)}(\delta_{il}a_{jm}b_{mk}+\delta_{jl}a_{im}b_{mk}+\delta_{ik}a_{jm}b_{ml}+\delta_{jk}a_{im}b_{ml})\,.\end{split}$
(95)

The trace is defined by $\mathrm{tr}(ab)\equiv a_{ij}b_{ji}$. Using this
quantity, we can decompose $a_{ij}b_{kl}$ into their irreducible parts. We
introduce $(P_{R}^{ab})_{ijkl}$ as a projection of $a_{ij}b_{kl}$ into its
irreducible representation $R$

$\begin{split}a_{ij}b_{kl}=\sum_{R}(P_{R}^{ab})_{ijkl}\,,\end{split}$ (96)

and each representation takes the form

$\displaystyle(P_{S}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{6}(\tilde{a}_{ij}\tilde{b}_{kl}+\tilde{a}_{kj}\tilde{b}_{il}+\tilde{a}_{ki}\tilde{b}_{jl}+\tilde{a}_{lj}\tilde{b}_{ik}+\tilde{a}_{il}\tilde{b}_{kj}+\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(97) $\displaystyle(P_{Y}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{6}(2\tilde{a}_{ij}\tilde{b}_{kl}-\tilde{a}_{kj}\tilde{b}_{il}-\tilde{a}_{ki}\tilde{b}_{jl}-\tilde{a}_{lj}\tilde{b}_{ik}-\tilde{a}_{il}\tilde{b}_{kj}+2\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(98) $\displaystyle(P_{T}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{2}(\tilde{a}_{ij}\tilde{b}_{kl}-\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(99) $\displaystyle(P_{s}^{ab})_{ijkl}$
$\displaystyle=-\frac{2}{(N_{f}+4)(N_{f}-2)}\biggl{[}\delta_{ij}\\{a,b\\}_{kl}+\delta_{kl}\\{a,b\\}_{ij}$
$\displaystyle\qquad-\frac{N_{f}}{4}\left(\delta_{il}\\{a,b\\}_{jk}+\delta_{ik}\\{a,b\\}_{jl}+\delta_{jl}\\{a,b\\}_{ik}+\delta_{jk}\\{a,b\\}_{il}\right)$
$\displaystyle\qquad\left.-\frac{4}{N_{f}}\delta_{ij}\delta_{kl}\mathrm{tr}(ab)+(\delta_{il}\delta_{jk}+\delta_{ik}\delta_{jl})\mathrm{tr}(ab)\right]\,,$
(100) $\displaystyle(P_{\mathrm{Ad}}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{2(N_{f}+2)}\left(\delta_{il}[a,b]_{jk}+\delta_{ik}[a,b]_{jl}+\delta_{jl}[a,b]_{ik}+\delta_{jk}[a,b]_{il}\right)\,,$
(101) $\displaystyle(P_{1}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{(N_{f}-1)(N_{f}+2)}\left[\delta_{il}\delta_{jk}+\delta_{ik}\delta_{jl}-\frac{2}{N_{f}}\delta_{ij}\delta_{kl}\right]\mathrm{tr}(ab)\,.$
(102)

By multiplying $4c_{ij}d_{kl}$, we obtain the projection operators with
indices of pions, $P_{R}^{ab;cd}$. Since the broken generators are parts of
$SU(N_{f})$ generators, we may write the projection operators in terms of the
group-theoretical constants of $SU(N_{f})$, $f_{abc}$ and $d_{abc}$, as in the
previous subsection. Hence, the commutation and anticommutation relations of
the broken generators are given by

$\begin{split}[T^{a},T^{b}]=if^{ab\bar{c}}T^{\bar{c}}\,,\qquad\\{T^{a},T^{b}\\}=d^{ab\bar{c}}T^{\bar{c}}+\frac{1}{N_{f}}\delta^{ab}\,,\end{split}$
(103)

where the indices with bars run over all generators of $SU(N_{f})$, but those
without bars denote the broken generators. We find $f^{ab\bar{c}}$ is non-zero
only when $\bar{c}$ runs over the unbroken generators, while $d^{ab\bar{c}}$
is non-zero only when $\bar{c}$ runs over the broken generators, which are
known as Cartan decomposition.

We find the projection operators $P_{R}^{ab;cd}$ in $SU(N_{f})/SO(N_{f})$ as
follows.

$\displaystyle P_{1}^{ab;cd}$
$\displaystyle=\frac{2}{(N_{f}-1)(N_{f}+2)}\delta^{ab}\delta^{cd}\,,$ (104)
$\displaystyle P_{\mathrm{Ad}}^{ab;cd}$
$\displaystyle=\frac{2}{N_{f}+2}f^{ab\bar{e}}f^{cd\bar{e}}\,,$ (105)
$\displaystyle P_{s}^{ab;cd}$
$\displaystyle=\frac{2N_{f}}{(N_{f}+4)(N_{f}-2)}d^{ab\bar{e}}d^{cd\bar{e}}\,,$
(106) $\displaystyle P_{T}^{ab;cd}$
$\displaystyle=\frac{1}{2}(\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc})-\frac{2}{N_{f}+2}f^{ab\bar{e}}f^{cd\bar{e}}\,,$
(107) $\displaystyle P_{Y}^{ab;cd}$
$\displaystyle=\frac{N_{f}-2}{3N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})+\frac{2(N_{f}-2)}{3N_{f}(N_{f}-1)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad-\frac{1}{3}(d^{ac\bar{e}}d^{bd\bar{e}}+d^{ad\bar{e}}d^{bc\bar{e}})+\frac{N_{f}-4}{3(N_{f}-2)}d^{ab\bar{e}}d^{cd\bar{e}}\,,$
(108) $\displaystyle P_{S}^{ab;cd}$
$\displaystyle=\frac{N_{f}+4}{6N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})-\frac{2(N_{f}+4)}{3N_{f}(N_{f}+2)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad+\frac{1}{3}(d^{ac\bar{e}}d^{bd\bar{e}}+d^{ad\bar{e}}d^{bc\bar{e}})-\frac{N_{f}+8}{3(N_{f}+4)}d^{ab\bar{e}}d^{cd\bar{e}}\,.$
(109)

These are consistent with the projection operators given in Ref. [88].

The trace of projection operators gives the dimension of the representation as
with the previous subsection,

$\begin{split}\sum_{a,b,c,d}P_{R}^{ab;cd}P_{R}^{cd;ab}=d_{R}\,,\end{split}$
(110)

where $d_{R}$ refers to the dimensions of the representation $R$ (see Table 2
for the dimensions). We confirm properties of the projection operators by
taking the trace of the projection operators:

$\begin{split}\sum_{R}\sum_{a,b}P_{R}^{ab;ab}=\frac{1}{4}(N_{f}-1)^{2}(N_{f}+2)^{2}\,,\qquad\sum_{a,b,c,d}P_{R}^{ab;cd}P_{R^{\prime}}^{ab;cd}=\delta_{RR^{\prime}}d_{R}\,.\end{split}$
(111)

Here, we use the Fierz identity for the broken generators [see Eq. 93] instead
of Eq. 88 modified to apply it for the broken generators.

###  C.3 $SU(N_{f})/USp(N_{f})$

Let us consider the symmetry breaking $SU(N_{f})\to USp(N_{f})$. The number of
unbroken generators is $N_{f}(N_{f}+1)/2$, and hence the number of broken
generators is $(N_{f}-2)(N_{f}+1)/2$. The broken generator satisfies
$T^{a}J=J(T^{a})^{T}=-(T^{a}J)^{T}$ with $J$ being the symplectic metric, and
hence the pions transform as a second rank antisymmetric tensor field. The
product of the antisymmetric tensors is decomposed as follows.

$\begin{split}a\otimes a=1_{S}\oplus a_{S}\oplus\mathrm{Ad}_{A}\oplus
T_{A}\oplus Y_{S}\oplus S_{S}\,.\end{split}$ (112)

Here, the subscripts $S$ and $A$ denote symmetrization and anti-symmetrization
of indices of pions (or $\pi^{a}T^{a}$) in the left-hand side. Table 3 shows
the representations and their dimensions. As seen in the table, several
representations vanish for $N_{f}\leq 6$. For $N_{f}=4$, $a_{S}$, $T_{A}$, and
$S_{S}$ do not appear and the decomposition is given by

$\begin{split}\mathbf{5}\otimes\mathbf{5}=\mathbf{1}_{S}\oplus\mathbf{10}_{A}\oplus\mathbf{14}_{S}\,,\end{split}$
(113)

and for $N_{f}=6$, $S_{S}$ does not appear and the decomposition is given by

$\begin{split}\mathbf{14}\otimes\mathbf{14}=\mathbf{1}_{S}\oplus\mathbf{14}_{S}\oplus\mathbf{21}_{A}\oplus\mathbf{70}_{A}\oplus\mathbf{90}_{S}\,.\end{split}$
(114)

Table 3: Representations in $USp(N_{f})$ and their dimensions  Reps. | Dims.
---|---
$a$ | $\dfrac{1}{2}(N_{f}-2)(N_{f}+1)$
$\mathrm{Ad}$ | $\dfrac{1}{2}N_{f}(N_{f}+1)$
$Y$ | $\dfrac{1}{12}(N_{f}-2)(N_{f}-1)N_{f}(N_{f}+3)$
$S$ | $\dfrac{1}{24}(N_{f}^{2}-1)N_{f}(N_{f}-6)$
$T$ | $\dfrac{1}{8}(N_{f}^{2}-1)(N_{f}-4)(N_{f}+2)$

We construct the projection operators for the product of broken generators as
with the previous subsections. We define the generator with the symplectic
metric as $(\tau^{a})_{ij}\equiv(T^{a})_{i}^{j^{\prime}}J_{jj^{\prime}}$,
which satisfies $(\tau^{a})^{T}=-\tau^{a}$. The Fierz identity for the broken
generator is given by

$\begin{split}(\tau^{a})_{ij}(\tau^{a})_{kl}=\frac{1}{4}\left(J_{ik}J_{jl}-J_{il}J_{jk}-\frac{2}{N_{f}}J_{ij}J_{kl}\right)\,.\end{split}$
(115)

The projection operators with the labels of the broken generators are defined
by

$\begin{split}P_{R}^{ab;cd}=4(\tau^{a})_{ij}(\tau^{b})_{kl}(P_{R})^{ijkl;i^{\prime}j^{\prime}k^{\prime}l^{\prime}}(\tau^{c})_{i^{\prime}j^{\prime}}(\tau^{d})_{k^{\prime}l^{\prime}}\,.\end{split}$
(116)

Now, we decompose the product $(\tau^{a})_{ij}(\tau^{b})_{kl}$ into its
irreducible representations as with previous subsection. We use the short-hand
notation $a_{ij}\equiv(\tau^{a})_{ij}$. We define the traceless part of the
product $a_{ij}b_{kl}$ as follows.

$\begin{split}&\tilde{a}_{ij}\tilde{b}_{kl}\equiv a_{ij}b_{kl}\\\
&\quad-\frac{1}{(N_{f}^{2}-4)(N_{f}+1)(N_{f}-4)}\left[2(3N_{f}+2)J_{ij}J_{kl}+(N_{f}^{2}+4)(J_{il}J_{jk}-J_{ik}J_{jl})\right]\mathrm{tr}(ab)\\\
&\quad+\frac{2}{(N_{f}-4)(N_{f}+2)}(J_{ij}J^{mn}a_{km}b_{nl}+J_{ij}J^{mn}a_{lm}b_{nk}+J_{kl}J^{mn}a_{im}b_{nj}+J_{kl}J^{mn}a_{jm}b_{ni})\\\
&\quad-\frac{4}{(N_{f}^{2}-4)(N_{f}-4)}(J_{il}J^{mn}a_{km}b_{nj}-J_{jl}J^{mn}a_{km}b_{ni}-J_{ik}J^{mn}a_{lm}b_{nj}+J_{jk}J^{mn}a_{lm}b_{ni})\\\
&\quad+\frac{N_{f}^{2}-2N_{f}-4}{(N_{f}^{2}-4)(N_{f}-4)}(J_{il}J^{mn}a_{jm}b_{nk}-J_{jl}J^{mn}a_{im}b_{nk}-J_{ik}J^{mn}a_{jm}b_{nl}+J_{jk}J^{mn}a_{im}b_{nl})\,.\end{split}$
(117)

Here, the invariant tensor is $J_{ij}$, but not $\delta_{ij}$, and the trace
is defined by $\mathrm{tr}(ab)\equiv J^{ik}J^{jl}a_{ij}b_{kl}$. This trace
coincides with the trace of $SU(N_{f})$ generators $T^{a}$ as follows:

$\begin{split}\mathrm{tr}(ab)&=J^{ik}J^{jl}(\tau^{a})_{ij}(\tau^{b})_{kl}=-J^{ik}J^{jl}J_{jj^{\prime}}J_{kk^{\prime}}(T^{a})_{i}^{j^{\prime}}(T^{b})_{l}^{k^{\prime}}\\\
&=\delta^{i}_{k^{\prime}}\delta^{l}_{j^{\prime}}(T^{a})_{i}^{j^{\prime}}(T^{b})_{l}^{k^{\prime}}=\mathrm{tr}(T^{a}T^{b})\,.\end{split}$
(118)

Using the traceless part, we can decompose $a_{ij}b_{kl}$ into their
irreducible parts. We introduce $(P_{R}^{ab})_{ijkl}$ as a projection of
$a_{ij}b_{kl}$ into its irreducible representation $R$,

$\begin{split}a_{ij}b_{kl}=\sum_{R}(P_{R}^{ab})_{ijkl}\,,\end{split}$ (119)

and each representation takes the form

$\displaystyle(P_{S}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{6}(\tilde{a}_{ij}\tilde{b}_{kl}-\tilde{a}_{ik}\tilde{b}_{jl}+\tilde{a}_{il}\tilde{b}_{jk}+\tilde{a}_{jk}\tilde{b}_{il}-\tilde{a}_{jl}\tilde{b}_{ik}+\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(120) $\displaystyle(P_{Y}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{6}(2\tilde{a}_{ij}\tilde{b}_{kl}+\tilde{a}_{ik}\tilde{b}_{jl}-\tilde{a}_{il}\tilde{b}_{jk}-\tilde{a}_{jk}\tilde{b}_{il}+\tilde{a}_{jl}\tilde{b}_{ik}+2\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(121) $\displaystyle(P_{T}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{2}(\tilde{a}_{ij}\tilde{b}_{kl}-\tilde{a}_{kl}\tilde{b}_{ij})\,,$
(122) $\displaystyle(P_{a}^{ab})_{ijkl}$
$\displaystyle=-\frac{2}{(N_{f}-4)(N_{f}+2)}\biggl{[}J_{ij}\\{a,b\\}_{kl}+J_{kl}\\{a,b\\}_{ij}$
$\displaystyle\qquad-\frac{N_{f}}{4}\left(J_{ik}\\{a,b\\}_{jl}+J_{jl}\\{a,b\\}_{ik}-J_{il}\\{a,b\\}_{jk}-J_{jk}\\{a,b\\}_{il}\right)$
$\displaystyle\qquad\left.+\frac{4}{N_{f}}J_{ij}J_{kl}\mathrm{tr}(ab)-(J_{il}J_{jk}-J_{ik}J_{jl})\mathrm{tr}(ab)\right]\,,$
(123) $\displaystyle(P_{\mathrm{Ad}}^{ab})_{ijkl}$
$\displaystyle=\frac{1}{2(N_{f}-2)}\left(J_{ik}[a,b]_{jl}-J_{jk}[a,b]_{il}-J_{il}[a,b]_{jk}+J_{jl}[a,b]_{ik}\right)\,,$
(124) $\displaystyle(P_{1}^{ab})_{ijkl}$
$\displaystyle=-\frac{1}{(N_{f}+1)(N_{f}-2)}\left[J_{il}J_{jk}-J_{ik}J_{jl}-\frac{2}{N_{f}}J_{ij}J_{kl}\right]\mathrm{tr}(ab)\,.$
(125)

We obtain the projection operators in terms of the pion indices by multiplying
$4c^{ij}d^{kl}$. Similarly to the previous subsection, the broken generators
are parts of the $SU(N_{f})$ generators, and hence we may write the projection
operators in terms of the group-theoretical constants for the $SU(N_{f})$
generators.

$\displaystyle P_{1}^{ab;cd}$
$\displaystyle=\frac{2}{(N_{f}-2)(N_{f}+1)}\delta^{ab}\delta^{cd}\,,$ (126)
$\displaystyle P_{\mathrm{Ad}}^{ab;cd}$
$\displaystyle=\frac{2}{N_{f}-2}f^{ab\bar{e}}f^{cd\bar{e}}\,,$ (127)
$\displaystyle P_{a}^{ab;cd}$
$\displaystyle=\frac{2N_{f}}{(N_{f}-4)(N_{f}+2)}d^{abe\bar{}}d^{cd\bar{e}}\,,$
(128) $\displaystyle P_{T}^{ab;cd}$
$\displaystyle=\frac{1}{2}(\delta^{ac}\delta^{bd}-\delta^{ad}\delta^{bc})-\frac{2}{N_{f}-2}f^{ab\bar{e}}f^{cd\bar{e}}\,,$
(129) $\displaystyle P_{Y}^{ab;cd}$
$\displaystyle=\frac{N_{f}+2}{3N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})-\frac{2(N_{f}+2)}{3N_{f}(N_{f}+1)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad+\frac{1}{3}(d^{ac\bar{e}}d^{bd\bar{e}}+d^{ad\bar{e}}d^{bc\bar{e}})-\frac{N_{f}+4}{3(N_{f}+2)}d^{ab\bar{e}}d^{cd\bar{e}}\,,$
(130) $\displaystyle P_{S}^{ab;cd}$
$\displaystyle=\frac{N_{f}-4}{6N_{f}}(\delta^{ac}\delta^{bd}+\delta^{ad}\delta^{bc})+\frac{2(N_{f}-4)}{3N_{f}(N_{f}-2)}\delta^{ab}\delta^{cd}$
$\displaystyle\qquad-\frac{1}{3}(d^{ac\bar{e}}d^{bd\bar{e}}+d^{ad\bar{e}}d^{bc\bar{e}})+\frac{N_{f}-8}{3(N_{f}-4)}d^{ab\bar{e}}d^{cd\bar{e}}\,.$
(131)

Here, the indices without bars denote the labels of broken generator, and the
indices with bars run over all $SU(N_{f})$ generators. These are consistent
with the projection operators given in Ref. [88].

Again, the trace of projection operators gives the dimension of the
representation,

$\begin{split}\sum_{a,b,c,d}P_{R}^{ab;cd}P_{R}^{cd;ab}=d_{R}\,,\end{split}$
(132)

where $d_{R}$ refers to the dimensions of the representation $R$ (see Table 3
for the dimensions). We confirm properties of the projection operators by
taking the trace of the projection operators:

$\begin{split}\sum_{R}\sum_{a,b}P_{R}^{ab;ab}=\frac{1}{4}(N_{f}-2)^{2}(N_{f}+1)^{2}\,,\qquad\sum_{a,b,c,d}P_{R}^{ab;cd}P_{R^{\prime}}^{ab;cd}=\delta_{RR^{\prime}}d_{R}\,.\end{split}$
(133)

Here, as we discussed in the previous section, we use the Fierz identity for
the broken generators [see Eq. 115] instead of Eq. 88.

###  C.4 Partial-wave amplitudes in 2-to-2 vertices

We summarize the tree-level partial-wave amplitudes $T_{\ell}^{R}$, which are
defined in Eq. 8, for each symmetry breaking. The symmetric representations
only have the $s$-wave amplitude ($T_{\ell=1}^{R}=0$ for a symmetric
representation $R$), while the antisymmetric representations only have the
$p$-wave amplitude ($T_{\ell=0}^{R}=0$ for an antisymmetric representation
$R$). For $G=SU(N_{f})\times SU(N_{f})$ and $H=SU(N_{f})$, the partial-wave
amplitudes are

$\begin{split}T_{\ell=0}^{Y}&=\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{f_{\pi}^{2}}\,,\qquad
T_{\ell=0}^{S}=-\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{f_{\pi}^{2}}\,,\\\
T_{\ell=0}^{\mathrm{Ad}_{S}}&=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(\frac{N_{f}}{2}s-\frac{4}{N_{f}}m_{\pi}^{2}\right)\,,\qquad
T_{\ell=0}^{1}=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(N_{f}s-\frac{2}{N_{f}}m_{\pi}^{2}\right)\,,\\\
T_{\ell=1}^{\mathrm{Ad}_{A}}&=\frac{1}{96\pi}\frac{2N_{f}}{f_{\pi}^{2}}\frac{s\beta^{2}}{4}\,,\qquad
T_{\ell=1}^{I}=0\quad(I=T\,,\overline{T})\,.\end{split}$ (134)

For $G=SU(N_{f})$ and $H=SO(N_{f})$, the partial-wave amplitudes are

$\begin{split}T_{\ell=0}^{Y}&=\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{2f_{\pi}^{2}}\,,\qquad
T_{\ell=0}^{S}=-\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{f_{\pi}^{2}}\,,\\\
T_{\ell=0}^{s}&=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(\frac{N_{f}}{4}s+\frac{N_{f}-4}{N_{f}}m_{\pi}^{2}\right)\,,\qquad
T_{\ell=0}^{1}=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(\frac{N_{f}}{2}s+\frac{N_{f}-2}{N_{f}}m_{\pi}^{2}\right)\,,\\\
T_{\ell=1}^{\mathrm{Ad}}&=\frac{1}{96\pi}\frac{N_{f}+2}{f_{\pi}^{2}}\frac{s\beta^{2}}{4}\,,\qquad
T_{\ell=1}^{T}=0\,.\end{split}$ (135)

For $G=SU(N_{f})$ and $H=USp(N_{f})$, the partial-wave amplitudes are

$\begin{split}T_{\ell=0}^{Y}&=-\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{2f_{\pi}^{2}}\,,\qquad
T_{\ell=0}^{S}=\frac{1}{32\pi}\frac{s-2m_{\pi}^{2}}{f_{\pi}^{2}}\,,\\\
T_{\ell=0}^{a}&=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(\frac{N_{f}}{4}s-\frac{N_{f}+4}{N_{f}}m_{\pi}^{2}\right)\,,\qquad
T_{\ell=0}^{1}=\frac{1}{32\pi}\frac{1}{f_{\pi}^{2}}\left(\frac{N_{f}}{2}s-\frac{N_{f}+2}{N_{f}}m_{\pi}^{2}\right)\,,\\\
T_{\ell=1}^{\mathrm{Ad}}&=\frac{1}{96\pi}\frac{N_{f}-2}{f_{\pi}^{2}}\frac{s\beta^{2}}{4}\,,\qquad
T_{\ell=1}^{T}=0\,.\end{split}$ (136)

We note that the projection operators for the higher-dimensional antisymmetric
representation do not appear in this vertex function.

###  C.5 Group Factors

The group-theoretical constants used in four-pion interactions [see Eq. 4] are
defined by

$\begin{split}r_{abcd}&\equiv\frac{2}{3}\mathrm{Tr}\left([T^{a},T^{c}][T^{b},T^{d}]+[T^{a},T^{d}][T^{b},T^{c}]\right)=-\frac{1}{3}(f_{ac\bar{e}}f_{bd\bar{e}}+f_{ad\bar{e}}f_{bc\bar{e}})\,,\\\
c_{abcd}&\equiv\frac{1}{3}\mathrm{Tr}\left(T^{\\{a}T^{b}T^{c}T^{d\\}}\right)\\\
&=\frac{1}{3}(d_{ab\bar{e}}d_{cd\bar{e}}+d_{ac\bar{e}}d_{bd\bar{e}}+d_{ad\bar{e}}d_{bc\bar{e}})+\frac{2}{3N_{f}}(\delta_{ab}\delta_{cd}+\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc})\,.\end{split}$
(137)

Here, the brace in the definition of $c_{abcd}$ denotes the sum of the
symmetric permutations. Again, the indices without bars run only over the
broken generators, while the indices with bars run over all $SU(N_{f})$
generators. The self-scattering cross section of dark pions depends on a
group-theoretical factor $a^{2}$.

$\begin{split}\sigma_{\mathrm{scatter}}=\frac{m_{\pi}^{2}}{32\pi
f_{\pi}^{4}}\frac{a^{2}}{N_{\pi}^{2}}\,.\end{split}$ (138)

Here, $a^{2}$ is given by the partial-wave amplitudes in the non-relativistic
limit as follows.

$\displaystyle SU(N_{f})$ $\displaystyle:\qquad
a^{2}=\frac{1}{2^{4}}\frac{(N_{f}^{2}-1)(3N_{f}^{4}-2N_{f}^{2}+6)}{2N_{f}^{2}}\,,$
(139) $\displaystyle SO(N_{f})$ $\displaystyle:\qquad
a^{2}=\frac{1}{2^{4}}\frac{(N_{f}-1)(N_{f}+2)(3N_{f}^{4}+7N_{f}^{3}-2N_{f}^{2}-12N_{f}+24)}{16N_{f}^{2}}\,,$
(140) $\displaystyle USp(N_{f})$ $\displaystyle:\qquad
a^{2}=\frac{1}{2^{4}}\frac{(N_{f}+1)(N_{f}-2)(3N_{f}^{4}-7N_{f}^{3}-2N_{f}^{2}+12N_{f}+24)}{16N_{f}^{2}}\,.$
(141)

These are consistent with the results in the original paper [51], which are
directly computed from the group-theoretical constants.

$\begin{split}a^{2}=\frac{1}{64}\sum_{abcd}(c_{abcd}^{2}+16r_{abcd}^{2})\,.\end{split}$
(142)

Last but not least, we show the coefficients $t_{R}^{2}$, which are used in
Eq. 62, for different symmetries. We compute $t_{R}^{2}$ by using the Fierz
identity for the generators instead of computing the group-theoretical
constants, $d$ and $f$, since all pion indices are contracted. These are given
by

$\displaystyle SU(N_{f})$ $\displaystyle:$ $\displaystyle t_{\mathrm{Ad}}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}(N_{f}^{2}-1)(N_{f}^{2}-4)}{48}\,,$
$\displaystyle t_{T,\overline{T}}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}(N_{f}^{2}-1)(N_{f}^{2}-4)}{96}\,,$
(143) $\displaystyle SO(N_{f})$ $\displaystyle:$ $\displaystyle
t_{\mathrm{Ad}}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}^{2}(N_{f}^{2}-1)(N_{f}-2)}{768}\,,$
$\displaystyle t_{T}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}(N_{f}+4)(N_{f}^{2}-1)(N_{f}-2)}{768}\,,$
(144) $\displaystyle USp(N_{f})$ $\displaystyle:$ $\displaystyle
t_{\mathrm{Ad}}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}^{2}(N_{f}^{2}-1)(N_{f}+2)}{768}\,,$
$\displaystyle t_{T}^{2}$
$\displaystyle=\frac{1}{2^{5}}\frac{N_{f}(N_{f}-4)(N_{f}^{2}-1)(N_{f}+2)}{768}\,,$
(145)

and others not listed here are zero. The sum of them is consistent with the
values given by [51].

$\displaystyle SU(N_{f})$ $\displaystyle:\qquad
t^{2}=\frac{1}{4^{5}}\frac{4N_{f}(N_{f}^{2}-1)(N_{f}^{2}-4)}{3}\,,$ (146)
$\displaystyle SO(N_{f})$ $\displaystyle:\qquad
t^{2}=\frac{1}{4^{5}}\frac{N_{f}(N_{f}^{2}-1)(N_{f}^{2}-4)}{12}\,,$ (147)
$\displaystyle USp(N_{f})$ $\displaystyle:\qquad
t^{2}=\frac{1}{4^{5}}\frac{N_{f}(N_{f}^{2}-1)(N_{f}^{2}-4)}{12}\,.$ (148)

## References

  * [1] Y. Hochberg, E. Kuflik, T. Volansky, and J. G. Wacker, “Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 113 (2014) 171301, arXiv:1402.5143 [hep-ph].
  * [2] D. N. Spergel and P. J. Steinhardt, “Observational evidence for selfinteracting cold dark matter,” Phys. Rev. Lett. 84 (2000) 3760–3763, arXiv:astro-ph/9909386.
  * [3] S. Tulin and H.-B. Yu, “Dark Matter Self-interactions and Small Scale Structure,” Phys. Rept. 730 (2018) 1–57, arXiv:1705.02358 [hep-ph].
  * [4] S. B. Gudnason, C. Kouvaris, and F. Sannino, “Dark Matter from new Technicolor Theories,” Phys. Rev. D74 (2006) 095008, arXiv:hep-ph/0608055 [hep-ph].
  * [5] D. D. Dietrich and F. Sannino, “Conformal window of SU(N) gauge theories with fermions in higher dimensional representations,” Phys. Rev. D75 (2007) 085018, arXiv:hep-ph/0611341 [hep-ph].
  * [6] M. Yu. Khlopov and C. Kouvaris, “Strong Interactive Massive Particles from a Strong Coupled Theory,” Phys. Rev. D77 (2008) 065002, arXiv:0710.2189 [astro-ph].
  * [7] M. Yu. Khlopov and C. Kouvaris, “Composite dark matter from a model with composite Higgs boson,” Phys. Rev. D78 (2008) 065040, arXiv:0806.1191 [astro-ph].
  * [8] R. Foadi, M. T. Frandsen, and F. Sannino, “Technicolor Dark Matter,” Phys. Rev. D80 (2009) 037702, arXiv:0812.3406 [hep-ph].
  * [9] J. Mardon, Y. Nomura, and J. Thaler, “Cosmic Signals from the Hidden Sector,” Phys. Rev. D80 (2009) 035013, arXiv:0905.3749 [hep-ph].
  * [10] G. D. Kribs, T. S. Roy, J. Terning, and K. M. Zurek, “Quirky Composite Dark Matter,” Phys. Rev. D81 (2010) 095001, arXiv:0909.2034 [hep-ph].
  * [11] R. Barbieri, S. Rychkov, and R. Torre, “Signals of composite electroweak-neutral Dark Matter: LHC/Direct Detection interplay,” Phys. Lett. B688 (2010) 212–215, arXiv:1001.3149 [hep-ph].
  * [12] M. Blennow, B. Dasgupta, E. Fernandez-Martinez, and N. Rius, “Aidnogenesis via Leptogenesis and Dark Sphalerons,” JHEP 03 (2011) 014, arXiv:1009.3159 [hep-ph].
  * [13] R. Lewis, C. Pica, and F. Sannino, “Light Asymmetric Dark Matter on the Lattice: SU(2) Technicolor with Two Fundamental Flavors,” Phys. Rev. D85 (2012) 014504, arXiv:1109.3513 [hep-ph].
  * [14] Lattice Strong Dynamics (LSD) Collaboration, T. Appelquist et al., “Lattice calculation of composite dark matter form factors,” Phys. Rev. D88 no. 1, (2013) 014502, arXiv:1301.1693 [hep-ph].
  * [15] A. Hietanen, R. Lewis, C. Pica, and F. Sannino, “Composite Goldstone Dark Matter: Experimental Predictions from the Lattice,” JHEP 12 (2014) 130, arXiv:1308.4130 [hep-ph].
  * [16] J. M. Cline, Z. Liu, G. Moore, and W. Xue, “Composite strongly interacting dark matter,” Phys. Rev. D90 no. 1, (2014) 015023, arXiv:1312.3325 [hep-ph].
  * [17] Lattice Strong Dynamics (LSD) Collaboration, T. Appelquist et al., “Composite bosonic baryon dark matter on the lattice: SU(4) baryon spectrum and the effective Higgs interaction,” Phys. Rev. D89 no. 9, (2014) 094508, arXiv:1402.6656 [hep-lat].
  * [18] A. Hietanen, R. Lewis, C. Pica, and F. Sannino, “Fundamental Composite Higgs Dynamics on the Lattice: SU(2) with Two Flavors,” JHEP 07 (2014) 116, arXiv:1404.2794 [hep-lat].
  * [19] G. Krnjaic and K. Sigurdson, “Big Bang Darkleosynthesis,” Phys. Lett. B751 (2015) 464–468, arXiv:1406.1171 [hep-ph].
  * [20] W. Detmold, M. McCullough, and A. Pochinsky, “Dark Nuclei I: Cosmology and Indirect Detection,” Phys. Rev. D90 no. 11, (2014) 115013, arXiv:1406.2276 [hep-ph].
  * [21] W. Detmold, M. McCullough, and A. Pochinsky, “Dark nuclei. II. Nuclear spectroscopy in two-color QCD,” Phys. Rev. D90 no. 11, (2014) 114506, arXiv:1406.4116 [hep-lat].
  * [22] M. Asano and R. Kitano, “Partially Composite Dark Matter,” JHEP 09 (2014) 171, arXiv:1406.6374 [hep-ph].
  * [23] J. Brod, J. Drobnak, A. L. Kagan, E. Stamou, and J. Zupan, “Stealth QCD-like strong interactions and the $t\bar{t}$ asymmetry,” Phys. Rev. D91 no. 9, (2015) 095009, arXiv:1407.8188 [hep-ph].
  * [24] O. Antipin, M. Redi, and A. Strumia, “Dynamical generation of the weak and Dark Matter scales from strong interactions,” JHEP 01 (2015) 157, arXiv:1410.1817 [hep-ph].
  * [25] E. Hardy, R. Lasenby, J. March-Russell, and S. M. West, “Big Bang Synthesis of Nuclear Dark Matter,” JHEP 06 (2015) 011, arXiv:1411.3739 [hep-ph].
  * [26] T. Appelquist et al., “Stealth Dark Matter: Dark scalar baryons through the Higgs portal,” Phys. Rev. D92 no. 7, (2015) 075030, arXiv:1503.04203 [hep-ph].
  * [27] T. Appelquist et al., “Detecting Stealth Dark Matter Directly through Electromagnetic Polarizability,” Phys. Rev. Lett. 115 no. 17, (2015) 171803, arXiv:1503.04205 [hep-ph].
  * [28] O. Antipin, M. Redi, A. Strumia, and E. Vigiani, “Accidental Composite Dark Matter,” JHEP 07 (2015) 039, arXiv:1503.08749 [hep-ph].
  * [29] E. Hardy, R. Lasenby, J. March-Russell, and S. M. West, “Signatures of Large Composite Dark Matter States,” JHEP 07 (2015) 133, arXiv:1504.05419 [hep-ph].
  * [30] R. T. Co, K. Harigaya, and Y. Nomura, “Chiral Dark Sector,” Phys. Rev. Lett. 118 no. 10, (2017) 101801, arXiv:1610.03848 [hep-ph].
  * [31] K. R. Dienes, F. Huang, S. Su, and B. Thomas, “Dynamical Dark Matter from Strongly-Coupled Dark Sectors,” Phys. Rev. D95 no. 4, (2017) 043526, arXiv:1610.04112 [hep-ph].
  * [32] H. Ishida, S. Matsuzaki, and Y. Yamaguchi, “Bosonic-Seesaw Portal Dark Matter,” PTEP 2017 no. 10, (2017) 103B01, arXiv:1610.07137 [hep-ph].
  * [33] S. J. Lonsdale, M. Schroor, and R. R. Volkas, “Asymmetric Dark Matter and the hadronic spectra of hidden QCD,” Phys. Rev. D96 no. 5, (2017) 055027, arXiv:1704.05213 [hep-ph].
  * [34] J. M. Berryman, A. de Gouvêa, K. J. Kelly, and Y. Zhang, “Dark Matter and Neutrino Mass from the Smallest Non-Abelian Chiral Dark Sector,” Phys. Rev. D96 no. 7, (2017) 075010, arXiv:1706.02722 [hep-ph].
  * [35] M. I. Gresham, H. K. Lou, and K. M. Zurek, “Nuclear Structure of Bound States of Asymmetric Dark Matter,” Phys. Rev. D96 no. 9, (2017) 096012, arXiv:1707.02313 [hep-ph].
  * [36] M. I. Gresham, H. K. Lou, and K. M. Zurek, “Early Universe synthesis of asymmetric dark matter nuggets,” Phys. Rev. D97 no. 3, (2018) 036003, arXiv:1707.02316 [hep-ph].
  * [37] A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, “Dark Matter as a weakly coupled Dark Baryon,” JHEP 10 (2017) 210, arXiv:1707.05380 [hep-ph].
  * [38] M. I. Gresham, H. K. Lou, and K. M. Zurek, “Astrophysical Signatures of Asymmetric Dark Matter Bound States,” Phys. Rev. D98 no. 9, (2018) 096001, arXiv:1805.04512 [hep-ph].
  * [39] M. Ibe, A. Kamada, S. Kobayashi, and W. Nakano, “Composite Asymmetric Dark Matter with a Dark Photon Portal,” JHEP 11 (2018) 203, arXiv:1805.06876 [hep-ph].
  * [40] E. Braaten, D. Kang, and R. Laha, “Production of dark-matter bound states in the early universe by three-body recombination,” JHEP 11 (2018) 084, arXiv:1806.00609 [hep-ph].
  * [41] A. Francis, R. J. Hudspith, R. Lewis, and S. Tulin, “Dark Matter from Strong Dynamics: The Minimal Theory of Dark Baryons,” arXiv:1809.09117 [hep-ph].
  * [42] Y. Bai, A. J. Long, and S. Lu, “Dark Quark Nuggets,” Phys. Rev. D 99 no. 5, (2019) 055047, arXiv:1810.04360 [hep-ph].
  * [43] X. Chu, C. Garcia-Cely, and H. Murayama, “Finite-size dark matter and its effect on small-scale structure,” Phys. Rev. Lett. 124 no. 4, (2020) 041101, arXiv:1901.00075 [hep-ph].
  * [44] E. Hall, T. Konstandin, R. McGehee, and H. Murayama, “Asymmetric Matters from a Dark First-Order Phase Transition,” arXiv:1911.12342 [hep-ph].
  * [45] Y.-D. Tsai, R. McGehee, and H. Murayama, “Resonant Self-Interacting Dark Matter from Dark QCD,” arXiv:2008.08608 [hep-ph].
  * [46] P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, “Accidentally Asymmetric Dark Matter,” Phys. Rev. Lett. 127 no. 21, (2021) 211101, arXiv:2103.09822 [hep-ph].
  * [47] M. Zhang, “Leptophilic composite asymmetric dark matter and its detection,” Phys. Rev. D 104 no. 5, (2021) 055008, arXiv:2104.06988 [hep-ph].
  * [48] S. Bottaro, M. Costa, and O. Popov, “Asymmetric accidental composite dark matter,” JHEP 11 (2021) 055, arXiv:2104.14244 [hep-ph].
  * [49] E. Hall, R. McGehee, H. Murayama, and B. Suter, “Asymmetric Dark Matter May Not Be Light,” arXiv:2107.03398 [hep-ph].
  * [50] G. D. Kribs and E. T. Neil, “Review of strongly-coupled composite dark matter models and lattice simulations,” Int. J. Mod. Phys. A31 no. 22, (2016) 1643004, arXiv:1604.04627 [hep-ph].
  * [51] Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky, and J. G. Wacker, “Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 115 no. 2, (2015) 021301, arXiv:1411.3727 [hep-ph].
  * [52] J. Wess and B. Zumino, “Consequences of anomalous Ward identities,” Phys. Lett. B 37 (1971) 95–97.
  * [53] E. Witten, “Global Aspects of Current Algebra,” Nucl. Phys. B 223 (1983) 422–432.
  * [54] M. Hansen, K. Langæble, and F. Sannino, “SIMP model at NNLO in chiral perturbation theory,” Phys. Rev. D 92 no. 7, (2015) 075036, arXiv:1507.01590 [hep-ph].
  * [55] U. Aydemir, M. M. Anber, and J. F. Donoghue, “Self-healing of unitarity in effective field theories and the onset of new physics,” Phys. Rev. D 86 (2012) 014025, arXiv:1203.5153 [hep-ph].
  * [56] M. E. Peskin, “The Alignment of the Vacuum in Theories of Technicolor,” Nucl. Phys. B 175 (1980) 197–233.
  * [57] J. Preskill, “Subgroup Alignment in Hypercolor Theories,” Nucl. Phys. B 177 (1981) 21–59.
  * [58] E. Witten, “Current Algebra, Baryons, and Quark Confinement,” Nucl. Phys. B 223 (1983) 433–444.
  * [59] D. A. Kosower, “SYMMETRY BREAKING PATTERNS IN PSEUDOREAL AND REAL GAUGE THEORIES,” Phys. Lett. B 144 (1984) 215–216.
  * [60] S. R. Coleman, J. Wess, and B. Zumino, “Structure of phenomenological Lagrangians. 1.,” Phys. Rev. 177 (1969) 2239–2247.
  * [61] C. G. Callan, Jr., S. R. Coleman, J. Wess, and B. Zumino, “Structure of phenomenological Lagrangians. 2.,” Phys. Rev. 177 (1969) 2247–2250.
  * [62] A. Katz, E. Salvioni, and B. Shakya, “Split SIMPs with Decays,” JHEP 10 (2020) 049, arXiv:2006.15148 [hep-ph].
  * [63] A. Kamada, H. Kim, and T. Sekiguchi, “Axionlike particle assisted strongly interacting massive particle,” Phys. Rev. D 96 no. 1, (2017) 016007, arXiv:1704.04505 [hep-ph].
  * [64] A. Manohar and H. Georgi, “Chiral Quarks and the Nonrelativistic Quark Model,” Nucl. Phys. B 234 (1984) 189–212.
  * [65] H. Georgi, “Generalized dimensional analysis,” Phys. Lett. B 298 (1993) 187–189, arXiv:hep-ph/9207278.
  * [66] G. ’t Hooft, “A Planar Diagram Theory for Strong Interactions,” Nucl. Phys. B 72 (1974) 461.
  * [67] G. ’t Hooft, “A Two-Dimensional Model for Mesons,” Nucl. Phys. B 75 (1974) 461–470.
  * [68] E. Witten, “Baryons in the 1/n Expansion,” Nucl. Phys. B 160 (1979) 57–115.
  * [69] E. Witten, “Current Algebra Theorems for the U(1) Goldstone Boson,” Nucl. Phys. B 156 (1979) 269–283.
  * [70] S. R. Coleman and E. Witten, “Chiral Symmetry Breakdown in Large N Chromodynamics,” Phys. Rev. Lett. 45 (1980) 100.
  * [71] E. Witten, “Large N Chiral Dynamics,” Annals Phys. 128 (1980) 363.
  * [72] R. S. Chivukula, M. J. Dugan, and M. Golden, “Electroweak corrections in technicolor reconsidered,” Phys. Lett. B 292 (1992) 435–441, arXiv:hep-ph/9207249.
  * [73] R. S. Chivukula, M. J. Dugan, and M. Golden, “Analyticity, crossing symmetry and the limits of chiral perturbation theory,” Phys. Rev. D 47 (1993) 2930–2939, arXiv:hep-ph/9206222.
  * [74] E. Kuflik, M. Perelstein, N. R.-L. Lorier, and Y.-D. Tsai, “Phenomenology of ELDER Dark Matter,” JHEP 08 (2017) 078, arXiv:1706.05381 [hep-ph].
  * [75] M. H. Namjoo, T. R. Slatyer, and C.-L. Wu, “Enhanced n-body annihilation of dark matter and its indirect signatures,” JHEP 03 (2019) 077, arXiv:1810.09455 [astro-ph.CO].
  * [76] D. Bhatia and S. Mukhopadhyay, “Unitarity limits on thermal dark matter in (non-)standard cosmologies,” JHEP 03 (2021) 133, arXiv:2010.09762 [hep-ph].
  * [77] E. Braaten and H. W. Hammer, “Universal Two-body Physics in Dark Matter near an S-wave Resonance,” Phys. Rev. D 88 (2013) 063511, arXiv:1303.4682 [hep-ph].
  * [78] X. Chu, C. Garcia-Cely, and H. Murayama, “A Practical and Consistent Parametrization of Dark Matter Self-Interactions,” JCAP 06 (2020) 043, arXiv:1908.06067 [hep-ph].
  * [79] S.-M. Choi, H. M. Lee, and M.-S. Seo, “Cosmic abundances of SIMP dark matter,” JHEP 04 (2017) 154, arXiv:1702.07860 [hep-ph].
  * [80] K. Saikawa and S. Shirai, “Precise WIMP Dark Matter Abundance and Standard Model Thermodynamics,” JCAP 08 (2020) 011, arXiv:2005.03544 [hep-ph].
  * [81] A. Dobado and J. R. Pelaez, “The Inverse amplitude method in chiral perturbation theory,” Phys. Rev. D 56 (1997) 3057–3073, arXiv:hep-ph/9604416.
  * [82] A. Soni and Y. Zhang, “Hidden SU(N) Glueball Dark Matter,” Phys. Rev. D 93 no. 11, (2016) 115025, arXiv:1602.00714 [hep-ph].
  * [83] L. Forestell, D. E. Morrissey, and K. Sigurdson, “Non-Abelian Dark Forces and the Relic Densities of Dark Glueballs,” Phys. Rev. D 95 no. 1, (2017) 015032, arXiv:1605.08048 [hep-ph].
  * [84] S.-M. Choi, Y. Hochberg, E. Kuflik, H. M. Lee, Y. Mambrini, H. Murayama, and M. Pierre, “Vector SIMP dark matter,” JHEP 10 (2017) 162, arXiv:1707.01434 [hep-ph].
  * [85] N. P. Mehta, S. T. Rittenhouse, J. P. D’Incao, J. von Stecher, and C. H. Greene, “General theoretical description of $n$-body recombination,” Phys. Rev. Lett. 103 (Oct, 2009) 153201. https://link.aps.org/doi/10.1103/PhysRevLett.103.153201.
  * [86] L. M. Delves, “Tertiary and general-order collisions,” Nuclear Physics 9 no. 3, (1958) 391–399.
  * [87] L. M. Delves, “Tertiary and general-order collisions (ii),” Nuclear Physics 20 (1960) 275–308.
  * [88] J. Bijnens and J. Lu, “Meson-meson Scattering in QCD-like Theories,” JHEP 03 (2011) 028, arXiv:1102.0172 [hep-ph].