# Even SIMP miracles are possible

Xiaoyong Chu<EMAIL_ADDRESS>Institute of High Energy Physics,
Austrian Academy of Sciences, Georg-Coch-Platz 2, 1010 Vienna, Austria  Marco
Nikolic<EMAIL_ADDRESS>Institute of High Energy Physics, Austrian
Academy of Sciences, Georg-Coch-Platz 2, 1010 Vienna, Austria  University of
Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Vienna, Austria  Josef
Pradler<EMAIL_ADDRESS>Institute of High Energy Physics, Austrian
Academy of Sciences, Georg-Coch-Platz 2, 1010 Vienna, Austria  University of
Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Vienna, Austria

###### Abstract

Strongly interacting massive particles $\pi$ have been advocated as prominent
dark matter candidates when they regulate their relic abundance through odd-
numbered $3\pi\to 2\pi$ annihilation. We show that successful freeze-out may
also be achieved through even-numbered interactions $XX\to\pi\pi$ once bound
states $X$ among the particles of the low-energy spectrum exist. In addition,
$X$-formation hosts the potential of also catalyzing odd-numbered $3\pi\to
2\pi$ annihilation processes, turning them into effective two-body processes
$\pi X\to\pi\pi$. Bound states are often a natural consequence of strongly
interacting theories. We calculate the dark matter freeze-out and comment on
the cosmic viability and possible extensions. Candidate theories can encompass
confining sectors without a mass gap, glueball dark matter, or $\phi^{3}$ and
$\phi^{4}$ theories with strong Yukawa or self-interactions.

††preprint: UWThPh 2024-5

#### Introduction.

The several decades-long efforts to detect dark matter (DM) non-
gravitationally have, to a significant degree, been fueled by a relic density
argument: the couplings of DM to the Standard Model (SM) that allow for
phenomenological exploration also successfully generate the DM abundance in
the early Universe through thermal two-body annihilation of DM into SM states.
As has been realized more recently, when the link between DM and SM becomes
too weak, DM may still regulate its abundance through $3\to 2$ DM-only
processes Hochberg  _et al._ (2014); also Dolgov (1980); Carlson  _et al._
(1992); de Laix  _et al._ (1995). This offers a pathway of DM as a thermal
relic even when being partially secluded from SM.

The odd-numbered $3\to 2$ reaction is naturally realized through the Wess-
Zumino-Witten (WZW) five-point interaction of a strongly interacting dark
sector Hochberg  _et al._ (2015). The low-energy DM states are the massive
pseudo-Goldstone bosons of the confining theory. Their even-numbered
interactions are described by the chiral Lagrangian, which, for suitable
topological structure, is amended by the WZW term. The possibility of a DM
number-depleting process that proceeds without participation of additional
degrees of freedom is hence very attractive and has led to a flurry of further
exploration, see Hochberg  _et al._ (2016); Kuflik  _et al._ (2016); Bernal
_et al._ (2015); Bernal and Chu (2016); Bernal  _et al._ (2016); Kuflik  _et
al._ (2016); Choi and Lee (2015, 2016); Soni and Zhang (2016); Kamada  _et
al._ (2016); Bernal  _et al._ (2017); Cline  _et al._ (2017); Choi  _et al._
(2017); Kuflik  _et al._ (2017); Heikinheimo  _et al._ (2018); Choi  _et al._
(2018); Hochberg  _et al._ (2018); Bernal  _et al._ (2020); Choi  _et al._
(2019); Katz  _et al._ (2020); Smirnov and Beacom (2020); Xing and Zhu (2021);
Braat and Postma (2023); Bernreuther  _et al._ (2023) among others.

Unconsidered within this paradigm is the possibility that strongly interacting
massive particles (SIMPs) may also form bound states among themselves. Calling
the SIMP DM candidate $\pi$ (flavor/internal indices suppressed) and their
bound state $X\equiv[\pi\pi]$ raises the question of the presence of $X$
impacting previous predictions and potentially changing the narrative of the
“SIMP mechanism.” In particular, bound states may lead to a catalysis of
reactions,

standard $3\to 2$ annihilation: $\displaystyle\ \pi+\pi+\pi\to\pi+\pi\,,$ (1a)
catalysed $3\to 2$ annihilation: $\displaystyle\ \pi+X\to\pi+\pi\,,$ (1b)
catalysed $4\to 2$ annihilation: $\displaystyle\ X+X\to\pi+\pi\,.$ (1c)

The last two are effective $2\to 2$ processes and compete with the free $3\to
2$ and $4\to 2$ counterpart reactions in depleting the overall DM mass
density. Moreover, whereas (1a) and (1b) are related through the same
underlying odd-numbered interaction, the final process $XX\to\pi\pi$ can be
entirely due to _even-numbered_ interactions, such as the four-point self-
interaction. This releases a requirement on the interaction structure of the
theory and opens the door to a SIMP mechanism without relying on anomaly-
mediated interactions.

Figure 1: Left to right: $X$-formation with $t-$channel enhancement;
$X$-assisted annihilation via two 4-point interactions; catalyzed $3\to 2$
annihilation.

Of course, the prospect of catalyzed reactions (1b) and (1c) taking place
requires $X$ to be part of the low-energy spectrum. Once the theory allows for
the existence of sufficiently long-lived $X$, their formation is guaranteed
through the radiationless exoergic process

guaranteed $X$ formation: $\displaystyle\quad\pi+\pi+\pi\to\pi+X\,.$ (2)

This reaction may also be mediated through even-numbered interactions, and in
its effective strength, it is not suppressed relative to the standard $3\to 2$
process. Hence, $X$ may form efficiently, and it shows that already for models
of SIMPs in isolation, the role of bound states calls to be studied; see Fig.
1 for illustration.

#### Exemplary SIMP model.

To allow for a paralleling exposition close to the original papers on the
SIMP-mechanism Hochberg  _et al._ (2014, 2015) we shall consider the low-
energy effective theory of $a=1,\dots,N_{\pi}$ massive pseudo-Goldstone bosons
$\pi_{a}$ as the DM candidates emerging from a confining dark non-Abelian
gauge group of $N_{f}$ fermion fields. The dynamics discussed is not exclusive
to this choice and further possibilities will be commented on.

In the construction of the chiral Lagrangian, $\pi_{a}$ are written as
fluctuations of the orientation of the chiral condensate $\Sigma_{0}$,
$\Sigma=e^{i\pi/f_{\pi}}\Sigma_{0}e^{i\pi^{T}/f_{\pi}}$, with
$\pi=\sum_{n}\pi_{n}T^{n}$ where $T^{n}$ are the $N_{\pi}$ broken generators
of the flavor group with normalization
$\operatorname{Tr}[T^{a}T^{b}]=\delta^{ab}/2$. The chiral Lagrangian is given
by,

$\displaystyle\mathcal{L}_{\chi}$
$\displaystyle=\frac{f_{\pi}^{2}}{4}\text{Tr}\left[\partial_{\mu}\Sigma\partial^{\mu}\Sigma^{\dagger}\right]-\frac{\mu^{3}}{2}\left(\operatorname{Tr}\left[M\Sigma\right]+h.c.\right).$
(3)

Expanding (3) in terms of $\Sigma$ yields the canonically normalized kinetic
terms, masses, and even-numbered interactions of $\pi$. Considering a flavor-
degenerate quark mass matrix $M$ with entries $m$, their universal mass is
given by $m_{\pi}^{2}=\pm{2\mu^{3}m}/{f_{\pi}^{2}}$; $f_{\pi}$ is the decay
constant, and the plus (minus) sign applies to $Sp(2N_{f})$ ($SU(N_{f})$ or
$SO(N_{f})$) residual flavor symmetry. Interactions are given by,

$\displaystyle\mathcal{L}^{\rm even}_{\rm int}$
$\displaystyle\supset-\frac{1}{3f_{\pi}^{2}}\operatorname{Tr}\left([\pi,\partial_{\mu}\pi][\pi,\partial^{\mu}\pi]\right)+\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\operatorname{Tr}\left[\pi^{4}\right]$
(4)

plus higher order terms $O(\pi^{6}/f_{\pi}^{6})$. Odd-numbered interactions in
form of a non-vanishing WZW term are only present for symmetry-breaking
pattern with coset spaces with non-trivial fifth homotopy groups Witten
(1983). The leading order WZW Lagrangian then reads,

$\displaystyle\mathcal{L}_{\rm int}^{\text{odd}}$
$\displaystyle=\frac{2N_{c}}{15\pi^{2}f_{\pi}^{5}}\epsilon^{\mu\nu\rho\sigma}\operatorname{Tr}\left[\pi\partial_{\mu}\pi\partial_{\nu}\pi\partial_{\rho}\pi\partial_{\sigma}\pi\right].$
(5)

In the picture of strongly interacting theories, $X$ would be a “meson
molecule” or “tetraquark” of mass $m_{X}=2m_{\pi}-E_{B}$ and $E_{B}>0$. For
the exposition of our ideas, we assume a shallow bound molecule with
$\kappa\equiv E_{B}/m_{\pi}\sim 0.1$ so that it can be treated as a non-
relativistic bound state Petraki  _et al._ (2015). This points to a theory
with $m/\Lambda\lesssim 1$ where $\Lambda\simeq 2\pi f_{\pi}$ Georgi (1993).
Theories with $m/\Lambda\ll 1$ such as for light quarks in QCD have a mass-gap
whereas for $m/\Lambda\gg 1$, the lowest lying states are expected to be
gluonia. The general ideas presented here also apply to deeper bound systems,
but their treatment requires advanced field theoretical tools, greatly
complicating matters. Even if we are far from the chiral limit, using (3)
allows for a most direct comparison with the original SIMP idea. In the
following, we shall take an $Sp(4)_{c}$ gauge theory ($N_{c}=4$) for
concreteness with $N_{f}=2$ fundamental Dirac fermions. After chiral symmetry
breaking, the vacuum alignment is $\Sigma_{0}=E$ where $E$ is the unitized
invariant matrix of the remaining $Sp(4)$ flavor symmetry group. A detailed
study of this choice is presented in Kulkarni  _et al._ (2022). Another
strongly-interacting option is taking $m/\Lambda\gg 1$ with glueball DM
Faraggi and Pospelov (2002); Soni and Zhang (2016) and to consider their bound
states Giacosa  _et al._ (2022). This is left for future work Chu. and Pradler
.

#### Bound state-assisted SIMP annihilation.

Before a detailed analysis, simple estimates may convince us that bound state
formation and $X$-assisted annihilation are both efficient and may even
supersede odd-numbered interactions. The parametric ratio of rates of $X$
formation to odd-numbered annihilation reads,

$\displaystyle\frac{\Gamma_{3\pi\to\pi X}}{\Gamma_{3\pi\to 2\pi}}$
$\displaystyle=\frac{\langle\sigma_{3\pi\to\pi
X}v^{2}\rangle}{\langle\sigma_{3\pi\to
2\pi}v^{2}\rangle}\approx\frac{|\psi(0)|^{2}f_{\pi}^{2}}{m_{\pi}^{5}}x_{\rm
f}^{2}.$ (6)

Here $\langle\sigma_{i}v^{2}\rangle$ are the thermally averaged collision
integrals of the respective processes. The dimensionful factor that relates
both “cross sections” is the square of the bound state wave function at the
origin $|\psi(0)|^{2}$; $x^{2}_{\rm f}\equiv(m_{\pi}/T_{\rm f})^{2}\sim 400$
is an enhancement factor that accounts for the different velocity scalings of
rates at non-relativistic freeze-out temperature $T_{\rm f}\sim m_{\pi}/20$.
Taking as an estimate $|\psi(0)|^{2}=1/a_{B}^{3}$ with the Bohr radius
$a_{B}=2/(\alpha_{s}^{\prime}m_{\pi})$ given in terms of the dark strong
coupling constant $\alpha_{s}^{\prime}=O(1)$ shows that the ratio in (6)
easily exceeds unity on account of $m_{\pi}/f_{\pi}=O(1)$ in strongly
interacting theories.111Alternatively, considering a spherical well
$V(r)=-V_{0}\theta(r-r_{0})$, setting $V_{0}\simeq 3.16m_{\pi}$ and
$r_{0}=m_{\pi}^{-1}$, results in a binding energy $E_{B}=0.1m_{\pi}$ and
yields $R(0)=\sqrt{4\pi}|\psi(0)|\simeq m_{\pi}^{3/2}$ Mahbubani  _et al._
(2021).

We may also convince ourselves that $X$-assisted annihilation competes with
its free counterpart. Naive dimensional analysis suggests that, for a pair of
$\pi$-particles, it is more likely for them to meet as constituents of a bound
state than as free particles,

$\displaystyle\frac{n_{X}|\psi(0)|^{2}}{n_{\pi}^{2}}\approx
2\sqrt{2}\pi^{3/2}x_{\rm f}^{3/2}e^{\kappa x_{\rm
f}}\frac{|\psi(0)|^{2}}{m_{\pi}^{3}}\,,$ (7)

where we used non-relativistic Maxwell-Boltzmann statistics for $n_{\pi}$ and
$n_{X}$. Note that $2\sqrt{2}\pi^{3/2}x_{\rm f}^{3/2}\approx 10^{3}$ for
$x_{\rm f}=20$. Thus, the ratio in (7) may easily exceed unity and suggests on
general grounds that ${XX\to\pi\pi}$ dominates over ${4\pi\to 2\pi}$, and that
${\pi X\to\pi\pi}$ dominates over ${3\pi\to 2\pi}$ when odd-numbered
interactions are present. This is what we mean by “catalysis.”

#### Cross sections.

We now calculate the relevant cross sections from (4) and (5) and first
consider the $X$-formation cross section $(\sigma_{3\pi\to\pi X}v^{2})$. In
the approximation that we are working in, the amplitude for a bound state
process is obtained from the free amplitude as follows. Calling
$\langle\vec{k}_{1},\vec{k}_{2},\vec{k}_{3}|\mathcal{T}\\{\dots\\}|\vec{p}_{1},\vec{p}_{2},\vec{p}_{3}\rangle$
the matrix-element of $3\pi\to 3\pi$ with respective incoming and outgoing
momenta $\vec{p}_{i}$ and $\vec{k}_{i}$, the amplitude for the bound state
formation $3\pi\to\pi X$ is obtained from

$\displaystyle\frac{\sqrt{2M_{X}}}{2m_{\pi}}\\!\\!\int\\!\\!\frac{d^{3}\vec{q}}{(2\pi)^{3}}\Big{\langle}\frac{\vec{K}}{2}-\vec{q},\frac{\vec{K}}{2}+\vec{q},\vec{k}_{3}\Big{|}\mathcal{T}\\{\dots\\}\Big{|}\vec{p}_{1},\vec{p}_{2},\vec{p}_{3}\Big{\rangle}\tilde{\psi}(\vec{q}).$
(8)

Here, $\vec{K}$ is the three-momentum of $X$; $\tilde{\psi}(\vec{q})$ is the
Fourier transform $\tilde{\psi}(\vec{q})=\int d^{3}\vec{x}\
\psi(\vec{x})e^{-i\vec{q}\cdot\vec{x}}$ of the wave function that is a
solution to the non-relativistic Schrödinger equation with confining potential
$V(r)$ where $r=|\vec{x}|$ is the separation of the constituent SIMPs.

On general grounds, many terms contribute to the integral in (8). They can be
classified by the power of the Cartesian components of $q_{i}$ that enter
through the matrix element, enforcing a selection rule on $\psi(\vec{x})$.
Decomposing the latter into angular and radial parts, $\psi(\vec{x})=Y_{\ell
m}(\hat{x})R(r)$, constant terms $(q_{i})^{0}$ yield for the integral
$\psi(\vec{x}=0)$ which is only non-vanishing for $S$-states $\ell=0$ where
$R(0)=\sqrt{4\pi}|\psi(0)|\neq 0$. This is our most important case, as it
concerns the ground state of $X$. Terms proportional to $q_{i}$ yield the
derivative of the radial wave function at the origin, $R^{\prime}(0)\equiv
dR/dr|_{r=0}$, so the lowest contributing angular momentum state is $P$-wave
with $\ell=1$. We will encounter this case for the WZW interaction below.
Finally, terms quadratic in $q_{i}$ yield nominally divergent integrals in
(8), probing the short distance behavior of the theory. In dimensional
regularization one may show that $\nabla^{2}R(0)=-m_{\pi}E_{B}R(0)$ holds
Brambilla  _et al._ (2003); Biondini and Shtabovenko (2022). Since we take
$E_{B}\ll m_{\pi}$, this becomes subleading in processes involving $S$-states,
and we are hence allowed to neglect such contributions.

For the bound state formation process $3\pi\to\pi X$, there are then various
diagrams to consider. During non-relativistic freeze-out, the dominant
processes are $t$\- and $u$-channel type diagrams of the sort depicted in Fig.
1 where two 4-point interactions are connected via a SIMP propagator. The
denominators of the propagators are enhanced by a matching kinematic condition
$k^{2}-m_{\pi}^{2}\propto-E_{B}m_{\pi}$. This renders other diagrams
irrelevant. The cross section is then given by

$\langle\sigma_{3\pi\to\pi
X}v^{2}\rangle\simeq\frac{\numprint{57041}}{\numprint{1310720}\sqrt{3}\pi^{2}}\frac{R_{S}^{2}(0)}{f_{\pi}^{8}}\left(\frac{m_{\pi}}{E_{B}}\right)^{3/2}\,.$
(9)

The prefactor depends on the gauge group and symmetry breaking pattern and is
averaged all possible incoming and summed over all outgoing flavor-
combinations so that $n_{\pi}^{3}\langle\sigma_{3\pi\to\pi X}v^{2}\rangle$
yields the total number change per time. In obtaining the result, we used a
non-relativistic expansion, assuming that the typical incoming kinetic energy
satisfies $m_{\pi}\langle v^{2}\rangle/2\lesssim E_{B}$. This is equivalent to
demanding $T_{\rm f}\lesssim E_{B}$. The thermal average over a Maxwell-
Boltzmann ensemble was taken in the final step. When comparing with (6) we
observe an additional enhancement by a factor of $({m_{\pi}}/{E_{B}})^{3/2}$
in the ratio of rates relative to the WZW-mediated $3\pi\to 2\pi$
annihilation. Detailed calculations of all cross sections are provided in the
supplementary material to this letter.

Similarly, we may proceed with the calculation of the annihilation cross
section for $XX\to\pi\pi$. A general computation would be a formidable
challenge, but we may again profit from the imposed selection rules, focusing
on $S$-wave initial bound states. There are six $t$\- and $u$-channel
diagrams, which become related in the limit that we neglect the internal
motion of the constituents. The final result reads

$\displaystyle\langle\sigma_{XX\to\pi\pi}v\rangle\simeq\frac{\numprint{2529757}}{\numprint{424673280}\sqrt{3}\,\pi^{3}}\frac{R_{S}^{4}(0)}{f_{\pi}^{8}}\,,$
(10)

where we have again summed over all flavor combinations. Finally, when
considering odd-numbered $3\to 2$ interactions enabled by (5), we obtain the
cross section for the related $\pi X\to\pi\pi$ process as

$\displaystyle\langle\sigma_{\pi
X\to\pi\pi}v\rangle\simeq\frac{\sqrt{5}N_{c}^{2}m_{\pi}^{3}}{512\pi^{6}f_{\pi}^{10}x}R_{P}^{\prime
2}(0)\,.$ (11)

Importantly, the process requires $X$ to be in a $P$-wave state $X_{P}$ with
$\ell=1$.

Figure 2: Evolution of the DM abundance, $Y_{\pi}+2Y_{X}\simeq Y_{\pi}$ shown
by the thick solid lines for $m_{\pi}=0.2\,$GeV, $m_{\pi}/f_{\pi}=5$,
$\kappa_{S}=4\kappa_{P}=0.1$, $R(0)=0.3m_{\pi}^{3/2}$, and
$dR(0)/dr=0.03m_{\pi}^{5/2}$. Thin solid and dotted lines show the equilibrium
and detailed balance abundances, respectively. Top: even-numbered case where
$XX$ annihilation freezes out at $x_{1}\simeq 20$ (left vertical line) and
bound-state formation freezes out at $x_{2}\simeq 115$ (right vertical line);
corresponding rates normalized to $H$ are shown as labeled. Middle: odd-
numbered WZW interactions included. $X_{P}$ annihilation freezes out
comparatively later (left vertical line), maintaining longer chemical
equilibrium. Free $3\to 2$ freeze-out is also shown for comparison. Bottom:
associated bound state abundances with matching colors. For the choice of
parameters, the observed DM abundance is reached for the odd case; using
instead $R(0)=1.4m_{\pi}^{3/2}$ yields the correct abundance for the even
case.

#### Abundance evolution.

We are now in a position to solve the evolution equations for the two
populations, free $\pi$ and $X$. Their respective total comoving number
densities, normalized to the total entropy density $s$, are given by
$Y_{\pi,X}=n_{\pi,X}/s$, where $s$ is the total entropy density of the
Universe and a sum over all flavors is implicit. We assume that kinetic
equilibrium with SM is maintained; we develop on this in the next section. At
high temperatures (small $x$), $\pi$ and $X$ follow their equilibrium
distributions, $Y_{\pi,X}=Y_{\pi,X}^{\rm eq}$, due to fast number-changing
processes. Their chemical decoupling happens at $x_{1}\sim 20$ when
$n_{X}^{2}\langle\sigma_{XX\to\pi\pi}v\rangle/n_{\pi}\simeq H(x_{1})$.
Subsequently, considering only even-numbered interactions, assuming dominance
of $XX\to\pi\pi$ over the free $4\pi\to 2\pi$ counterpart and neglecting the
inverse process, a particularly simple form of the Boltzmann equation is found
for the combination $Y_{\pi}+2Y_{X}$,222The approximate evolution of the
overall mass density in the dark sector is $\rho_{\rm
DM}\simeq(Y_{\pi}+2Y_{X})m_{\pi}s$ up to corrections $O(E_{B}/m_{\pi})$.

$\displaystyle\frac{d(Y_{\pi}+2Y_{X})}{dx}=-\frac{2\langle\sigma_{XX\to\pi\pi}v\rangle
Y_{X}^{2}s(x)}{xH(x)}\quad(x>x_{1})\,.$ (12)

The immediate evolution that ensues for $x>x_{1}$ is non-trivial because bound
state formation $3\pi\leftrightarrow\pi X$ is still operative. This maintains
a detailed balance between the $X$ and $\pi$ populations,

$Y_{X}=\frac{Y_{\pi}^{2}Y_{X}^{\rm eq}}{(Y_{\pi}^{\rm
eq})^{2}}=Y_{\pi}^{2}\frac{N_{X}}{N_{\pi}^{2}}(2\pi x)^{3/2}e^{\kappa
x}\left(\frac{m_{X}}{m_{\pi}}\right)^{3/2}\frac{s(x)}{m_{\pi}^{3}}\,,$ (13)

where $N_{\pi}=5$ and $N_{X}=(N_{\pi}+1)N_{\pi}/2=15$ are the possible flavor
combinations. Together, (13) and (12) determine the evolution of $Y_{\pi}$ and
$Y_{X}$ for as long as their detailed balance holds until bound state
formation freezes out at $x_{2}>x_{1}$ when $\Gamma_{3\pi\to\pi X}\equiv
n_{\pi}^{2}\langle\sigma_{3\pi\to\pi X}v^{2}\rangle=H(x_{2})$.

Using (13) in (12) with $Y_{X}\ll Y_{\pi}$ and neglecting $dY_{X}/dx$, the
Boltzmann equation becomes one for $dY_{\pi}/dx$ that can be integrated. To
leading order in $x_{1}/x_{2}$ we obtain the following solution,

$\displaystyle
Y_{\pi}^{-3}(x_{2})\simeq\frac{256\sqrt{2}\pi^{8}g_{*}^{5/2}m_{\pi}M_{P}\langle\sigma_{XX\to\pi\pi}v\rangle}{6075\sqrt{5}x_{2}^{4}}\frac{N_{X}^{2}}{N_{\pi}^{4}}$
$\displaystyle\times\left[8(\kappa x_{2})^{4}\operatorname{Ei}(2\kappa
x_{2})-e^{2\kappa x_{2}}(3+2\kappa
x_{2}+2\kappa^{2}x_{2}^{2}+4\kappa^{3}x_{2}^{3})\right],$ (14)

where $\operatorname{Ei}$ is the exponential integral function and $M_{P}$ is
the reduced Planck mass and $g_{*}$ are the effective degrees of freedom at
$x_{2}$. This approximation works for $\kappa x_{1}\gtrsim 1$ which is
congruent with assuming $E_{B}>T$ at chemical decoupling. In writing the
solution, we have also taken $Y^{-3}_{\pi}(x_{1})\ll Y^{-3}_{\pi}(x_{2})$. To
within a factor of two, we may then put the solution in suggestive form,

$\displaystyle\Omega^{\rm even}_{\pi}\sim
0.2\left(\frac{10^{3}}{\kappa^{4}e^{\kappa x_{2}}}\frac{{\rm
bn/GeV}}{\langle\sigma_{XX\to\pi\pi}v\rangle/m_{\pi}}\frac{m_{\pi}}{\text{GeV}}\right)^{1/3}.$
(15)

This is a central result. First, note that the relic density depends on
$x_{2}$, i.e., the moment of freeze-out of bound state formation, and not
$x_{1}$. Second, we observe a strong dependence on $\kappa$, and with $x_{2}$
typically between 50 to 100, $\kappa=0.1$ suggests sub-GeV DM with cm2/gram
self-interactions—in the same ballpark as the odd-numbered SIMP case. Finally,
when compared to ordinary $2\to 2$ freeze-out, the inverse dependence on the
annihilation cross section is softened by the cubic root.

The top panel of Fig. 2 shows the full numerical solution for $Y_{\pi}$ and
$Y_{X}$ as a function of $x$ for the parameters given in the caption. The
evolution contains three steps as discussed above. First, chemical decoupling
happens at $x_{1}\simeq 20$ when $XX\leftrightarrow\pi\pi$ freezes out. Both
$X$ and $\pi$ develop a chemical potential and keep a detailed balance via
$3\pi\leftrightarrow\pi X$ until $x_{2}\simeq 115$. For $x>x_{2}$ free $\pi$
have reached their relic density value while $Y_{X}$ keeps decreasing further,
as the rate of $XX\to\pi\pi$ with respect to $n_{X}$ is still larger than the
Hubble rate, $n_{X}\langle\sigma_{XX\to\pi\pi}v\rangle>H(x_{2})$. The
evolution of $Y_{X}$ is shown in the bottom panel.

#### Odd-numbered case.

We now turn our attention to the scenario when we are additionally afforded
odd-numbered interactions. As calculated in (11), the efficiency of $\pi
X_{P}\to\pi\pi$ entirely hinges on the availability $X_{P}$ which must be
present in the low energy spectrum. Since the path to collisional excitation
is open, one may consider the detailed balancing relation
$n_{X_{P}}/n_{X_{S}}=3\exp[-|\kappa-\kappa_{P}|x]$ as an estimate for the
number density $n_{X_{P}}$ of excited states; $\kappa_{P}=E_{P}/m_{\pi}$ with
$E_{P}$ being the $P$-wave binding energy. As we shall see now, the impact of
bound states $X$ can also be substantial. Even with the bottleneck of $P$-wave
states for WZW interactions, $\pi X_{P}\to\pi\pi$ supersede the free $3\pi\to
2\pi$ scenario, and is generally stronger than $XX\rightarrow\pi\pi$.

For the odd-numbered case, the chemical decoupling happens when
$n_{X_{P}}\langle\sigma_{\pi X_{P}\rightarrow\pi\pi}v\rangle\simeq H(x_{1})$.
For $x\geq x_{1}$, the right hand side of (12) is replaced by
$-s\langle\sigma_{\pi X_{P}\rightarrow\pi\pi}v\rangle Y_{\pi}Y_{X_{P}}/(xH)$.
Since $Y_{X_{P}}/Y_{X_{S}}$, and thus the rate of $\pi X_{P}\to\pi\pi$,
decrease exponentially for $x>x_{1}$, $Y_{\pi}$ already freezes out at
$x_{1}$. If $3\pi\to\pi X$ decouples later, $x_{2}>x_{1}$, (13) allows us to
estimate the yield from $n_{X_{P}}^{\rm eq}\langle\sigma_{\pi
X_{P}\rightarrow\pi\pi}v\rangle/H(x_{1})\simeq(Y_{\pi}^{\rm
eq})^{2}/Y^{2}_{\pi}(x_{1})$. This gives

$\Omega^{\rm odd}_{\pi}\simeq
0.2\,\left(\frac{x_{1}}{20}\right)^{5/4}\left(\frac{e^{-\kappa_{P}x_{1}}\,10^{-3}\,{\rm
bn/GeV}}{\langle\sigma_{\pi
X_{P}\rightarrow\pi\pi}v\rangle/m_{\pi}}\right)^{1/2}\,.$ (16)

The middle panel of Fig. 2 shows the full numerical freeze-out solution. In
comparison with the even-numered case, the stronger $\pi X_{P}\to\pi\pi$
reaction maintains longer chemical equilibrium. At the same time, with
$p$-wave states diminishing more rapidly with temperature, there is no
distinct intermediate phase, and freeze-out happens in one step (unless
considering large values for $dR(0)/dr$). Also shown is the standard SIMP
scenario through free $3\to 2$ reactions, and we observe an order of magnitude
smaller freeze-out yield for the chosen set of parameters when $X$ is
considered (catalysis). This softens the notorious tension between maximum
permissible elastic scattering cross section and relic density requirement in
SIMP models Hansen  _et al._ (2015), and for the shown case, both requirements
are indeed satisfied. The leading order elastic scattering cross section is
$\sigma_{\rm el}^{\rm LO}/m_{\pi}=0.2$ bn/GeV, receiving higher order
corrections in the chiral expansion Hansen  _et al._ (2015) as well as from
the bound state in the spectrum. For the latter, we may estimate a $s$-wave
scattering length of the order $1/(\kappa_{S}m_{\pi}^{2})^{1/2}$ Braaten and
Hammer (2013), leading to resonant-induced $\sigma_{\rm el}/m_{\pi}\sim
1/(\kappa_{S}m_{\pi}^{3})$. For $m_{\pi}=0.2\,$GeV and $\kappa_{S}\sim 0.1$
adopted here, it suggests an elastic scattering in the same ballpark and below
bn/GeV.

Within the exemplary scenario of pseudo-Goldstone bosons making $X$, we may
also comment on the influence of additional low-lying states, such as
$\rho$-mesons of mass $m_{\rho}\lesssim 2m_{\pi}$. Additional annihililation
channels become available, such as $3\pi\to\pi\rho^{*}\to\pi\pi$ Choi  _et
al._ (2018) or $3\pi\to\pi\rho$ Bernreuther  _et al._ (2023). For the
parameter region of interest and unless one considers a finely-tuned resonance
region $m_{\rho}\simeq 2m_{\pi}$, we find that these processes are generally
subleading to the $X$-mediated ones and we are allowed to neglect them.

#### Coupling to SM and longevity of $X$.

As is pertinent to all SIMP scenarios that freeze out through self-depletion,
kinetic equilibrium with radiation must be maintained to achieve a cold DM
scenario. An elastic scattering process $\pi\,{\rm SM}_{i}\to\pi\,{\rm
SM}_{i}$ with rate $\Gamma_{\pi\,{\rm SM}}=\langle\sigma_{\pi\,{\rm
SM}}c\rangle n_{i}>H$, that brings SM and dark sector into kinetic equilibrium
during freeze-out, generally also enables $\pi\pi\to{\rm SM}_{i}\overline{\rm
SM}_{i}$ annihilation. The SIMP mechanism then requires $\Gamma_{\rm
ann}=n_{\pi}\langle\sigma_{\rm ann}v\rangle<H$, where $\sigma_{\rm ann}$ is
the cross section for $\pi\pi\to{\rm SM}_{i}\overline{\rm SM}_{i}$. On the
account of $n_{i}/n_{\pi}\gg 1$, where $n_{i}$ is the number density of a
relativistic SM species, both conditions are generically satisfied Hochberg
_et al._ (2014). In the current context, interactions of $\pi$ with SM may
also destabilize $X$ through $X=[\pi\pi]\to{\rm SM}_{i}\overline{\rm SM}_{i}$
and we must ensure that for the decay rate $\Gamma_{X}<H$ holds until after
freeze-out.

Assuming $\sigma_{\rm ann}v\simeq const.$ and noting that $|\psi(0)|^{2}v$ has
units of particle flux, we may estimate the induced decay width of $X$ as
$\Gamma_{X}\sim|\psi(0)|^{2}(\sigma_{\rm ann}v)$ where $(\sigma_{\rm ann}v)$
is the $\pi\pi\to{\rm SM}_{i}\overline{\rm SM}_{i}$ annihilation cross
section. The longevity requirement $\Gamma_{X}/H<1$ thereby translates into an
upper limit on the annihilation cross section,

$\displaystyle\sigma_{\rm ann}v\lesssim 10^{-3}{\rm pb}\
x^{-2}\left(\frac{m_{\pi}}{100\
\mathrm{MeV}}\right)^{2}\frac{\mathrm{MeV}^{3}}{|\psi(0)|^{2}}\,.$ (17)

In the simplest cases, such as contact interactions through a heavy mediator,
elastic and annihilation cross sections are additionally related and in the
same ballpark, $\sigma_{\pi\,{\rm SM}}c\sim\sigma_{\rm ann}v$. We may then use
the bound in (17) to estimate the implied ceiling on the elastic scattering
rate,

$\displaystyle 1\lesssim\frac{\Gamma_{\pi\,{\rm
SM}}}{H}\lesssim\frac{10^{6}}{x^{3}}\left(\frac{m_{\pi}}{100\
\mathrm{MeV}}\right)^{3}\frac{\mathrm{MeV}^{3}}{|\psi(0)|^{2}}\,.$ (18)

This can be easily satisfied at freeze-out for $|\psi(0)|<m_{\pi}^{3/2}$. We
hence conclude that it is possible to retain kinetic equilibrium while
maintaining sufficient longevity of $X$ and paired with sub-Hubble two-body
annihilation. Therefore, the model-building requirements for coupling the dark
sector to the SM are not escalated compared to the standard SIMP mechanism,
and one may use the options already entertained in the original work Hochberg
_et al._ (2014).

Finally, additional $X$ formation and breakup reactions may open when
introducing couplings to SM. It is important to note that the detailed
balancing condition (13) between $X$ and $\pi$—being a Saha equation—remains
unaltered. If the new processes dominate over $3\pi\leftrightarrow\pi X$, (13)
retains its validity longer, $x_{2}$ will be larger, and the relic density
smaller. Hence, if anything, the introduction of SM-interactions harbors the
potential to make $X$-assisted freeze-out even more efficient, adding a level
of richness, without jeopardizing the overall picture.

#### Conclusions.

Bound-state-assisted self-depletion offers a novel approach to DM relic
density generation. It supports an even-numbered relic SIMP mechanism and
enhances odd-numbered counterparts. Both are realized in strongly interacting
theories. A broader study of the many aspects mentioned in this work, as well
as the exploration of other particle-physics realizations giving rise to $X$,
such as glueballs or strong Yukawa forces, will be the subject of upcoming
work Chu. and Pradler .

#### Acknowledgements.

JP thanks M. Pospelov for a useful discussion on QCD-like theories. This work
was supported by the Research Network Quantum Aspects of Spacetime (TURIS) and
by the FWF Austrian Science Fund research teams grant STRONG-DM (FG1).
Funded/Co-funded by the European Union (ERC, NLO-DM, 101044443).

## References

  * Hochberg  _et al._ (2014) Y. Hochberg, E. Kuflik, T. Volansky, and J. G. Wacker, Phys. Rev. Lett. 113, 171301 (2014), arXiv:1402.5143 [hep-ph] .
  * Dolgov (1980) A. D. Dolgov, Yad. Fiz. 31, 1522 (1980).
  * Carlson  _et al._ (1992) E. D. Carlson, M. E. Machacek, and L. J. Hall, Astrophys. J. 398, 43 (1992).
  * de Laix  _et al._ (1995) A. A. de Laix, R. J. Scherrer, and R. K. Schaefer, Astrophys. J. 452, 495 (1995), arXiv:astro-ph/9502087 .
  * Hochberg  _et al._ (2015) Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky, and J. G. Wacker, Phys. Rev. Lett. 115, 021301 (2015), arXiv:1411.3727 [hep-ph] .
  * Hochberg  _et al._ (2016) Y. Hochberg, E. Kuflik, and H. Murayama, JHEP 05, 090 (2016), arXiv:1512.07917 [hep-ph] .
  * Kuflik  _et al._ (2016) E. Kuflik, M. Perelstein, N. R.-L. Lorier, and Y.-D. Tsai, Phys. Rev. Lett. 116, 221302 (2016), arXiv:1512.04545 [hep-ph] .
  * Bernal  _et al._ (2015) N. Bernal, C. Garcia-Cely, and R. Rosenfeld, JCAP 04, 012 (2015), arXiv:1501.01973 [hep-ph] .
  * Bernal and Chu (2016) N. Bernal and X. Chu, JCAP 01, 006 (2016), arXiv:1510.08527 [hep-ph] .
  * Bernal  _et al._ (2016) N. Bernal, X. Chu, C. Garcia-Cely, T. Hambye, and B. Zaldivar, JCAP 03, 018 (2016), arXiv:1510.08063 [hep-ph] .
  * Choi and Lee (2015) S.-M. Choi and H. M. Lee, JHEP 09, 063 (2015), arXiv:1505.00960 [hep-ph] .
  * Choi and Lee (2016) S.-M. Choi and H. M. Lee, Phys. Lett. B 758, 47 (2016), arXiv:1601.03566 [hep-ph] .
  * Soni and Zhang (2016) A. Soni and Y. Zhang, Phys. Rev. D 93, 115025 (2016), arXiv:1602.00714 [hep-ph] .
  * Kamada  _et al._ (2016) A. Kamada, M. Yamada, T. T. Yanagida, and K. Yonekura, Phys. Rev. D 94, 055035 (2016), arXiv:1606.01628 [hep-ph] .
  * Bernal  _et al._ (2017) N. Bernal, X. Chu, and J. Pradler, Phys. Rev. D 95, 115023 (2017), arXiv:1702.04906 [hep-ph] .
  * Cline  _et al._ (2017) J. M. Cline, H. Liu, T. Slatyer, and W. Xue, Phys. Rev. D 96, 083521 (2017), arXiv:1702.07716 [hep-ph] .
  * Choi  _et al._ (2017) S.-M. Choi, H. M. Lee, and M.-S. Seo, JHEP 04, 154 (2017), arXiv:1702.07860 [hep-ph] .
  * Kuflik  _et al._ (2017) E. Kuflik, M. Perelstein, N. R.-L. Lorier, and Y.-D. Tsai, JHEP 08, 078 (2017), arXiv:1706.05381 [hep-ph] .
  * Heikinheimo  _et al._ (2018) M. Heikinheimo, K. Tuominen, and K. Langæble, Phys. Rev. D 97, 095040 (2018), arXiv:1803.07518 [hep-ph] .
  * Choi  _et al._ (2018) S.-M. Choi, H. M. Lee, P. Ko, and A. Natale, Phys. Rev. D 98, 015034 (2018), arXiv:1801.07726 [hep-ph] .
  * Hochberg  _et al._ (2018) Y. Hochberg, E. Kuflik, R. Mcgehee, H. Murayama, and K. Schutz, Phys. Rev. D 98, 115031 (2018), arXiv:1806.10139 [hep-ph] .
  * Bernal  _et al._ (2020) N. Bernal, X. Chu, S. Kulkarni, and J. Pradler, Phys. Rev. D 101, 055044 (2020), arXiv:1912.06681 [hep-ph] .
  * Choi  _et al._ (2019) S.-M. Choi, H. M. Lee, Y. Mambrini, and M. Pierre, JHEP 07, 049 (2019), arXiv:1904.04109 [hep-ph] .
  * Katz  _et al._ (2020) A. Katz, E. Salvioni, and B. Shakya, JHEP 10, 049 (2020), arXiv:2006.15148 [hep-ph] .
  * Smirnov and Beacom (2020) J. Smirnov and J. F. Beacom, Phys. Rev. Lett. 125, 131301 (2020), arXiv:2002.04038 [hep-ph] .
  * Xing and Zhu (2021) C.-Y. Xing and S.-H. Zhu, Phys. Rev. Lett. 127, 061101 (2021), arXiv:2102.02447 [hep-ph] .
  * Braat and Postma (2023) P. Braat and M. Postma, JHEP 03, 216 (2023), arXiv:2301.04513 [hep-ph] .
  * Bernreuther  _et al._ (2023) E. Bernreuther, N. Hemme, F. Kahlhoefer, and S. Kulkarni, (2023), arXiv:2311.17157 [hep-ph] .
  * Witten (1983) E. Witten, Nucl. Phys. B 223, 422 (1983).
  * Petraki  _et al._ (2015) K. Petraki, M. Postma, and M. Wiechers, JHEP 06, 128 (2015), arXiv:1505.00109 [hep-ph] .
  * Georgi (1993) H. Georgi, Phys. Lett. B 298, 187 (1993), arXiv:hep-ph/9207278 .
  * Kulkarni  _et al._ (2022) S. Kulkarni, A. Maas, S. Mee, M. Nikolic, J. Pradler, and F. Zierler, (2022), arXiv:2202.05191 [hep-ph] .
  * Faraggi and Pospelov (2002) A. E. Faraggi and M. Pospelov, Astropart. Phys. 16, 451 (2002), arXiv:hep-ph/0008223 .
  * Giacosa  _et al._ (2022) F. Giacosa, A. Pilloni, and E. Trotti, Eur. Phys. J. C 82, 487 (2022), arXiv:2110.05582 [hep-ph] .
  * (35) X. Chu. and J. Pradler, In preparation.
  * Mahbubani  _et al._ (2021) R. Mahbubani, M. Redi, and A. Tesi, JCAP 02, 039 (2021), arXiv:2007.07231 [hep-ph] .
  * Brambilla  _et al._ (2003) N. Brambilla, D. Eiras, A. Pineda, J. Soto, and A. Vairo, Phys. Rev. D 67, 034018 (2003), arXiv:hep-ph/0208019 .
  * Biondini and Shtabovenko (2022) S. Biondini and V. Shtabovenko, JHEP 03, 172 (2022), arXiv:2112.10145 [hep-ph] .
  * Hansen  _et al._ (2015) M. Hansen, K. Langæble, and F. Sannino, Phys. Rev. D 92, 075036 (2015), arXiv:1507.01590 [hep-ph] .
  * Braaten and Hammer (2013) E. Braaten and H. W. Hammer, Phys. Rev. D 88, 063511 (2013), arXiv:1303.4682 [hep-ph] .
  * Kamada  _et al._ (2023) A. Kamada, S. Kobayashi, and T. Kuwahara, JHEP 02, 217 (2023), arXiv:2210.01393 [hep-ph] .

##  Appendix A Boltzmann equations

Here we collect the relevant Boltzmann equations for $\pi$ and $X$. The right-
hand side thereby provides the definition of the thermally averaged cross
sections $\langle\sigma_{i}v\rangle$ and collision terms
$\langle\sigma_{i}v^{2}\rangle$ that are computed in the subsequent section.
Our yield variables are defined as flavor-summed quantities
$Y_{\pi}\equiv\sum_{a}{n_{\pi_{a}}}/{s}$, $Y_{X}=\sum_{a\leq
b}{n_{X_{ab}}}/{s}$ where the sum is over flavors $a,b,\dots=1,\dots N_{\pi}$.
The evolution of $\pi$ and $X$ is then governed by

$\displaystyle\frac{dY_{\pi}}{dx}$
$\displaystyle=\frac{1}{3H(x)}\frac{ds}{dx}\left[-2\langle\sigma_{XX\to
2\pi}v\rangle\left(Y_{X}^{2}-(\frac{Y_{\pi}}{Y^{\rm eq}_{\pi}}Y_{X}^{\rm
eq})^{2}\right)+2s\langle\sigma_{3\pi\to\pi
X}v^{2}\rangle\left(Y_{\pi}^{3}-Y_{\pi}(Y_{\pi}^{\rm
eq})^{2}\frac{Y_{X}}{Y_{X}^{\rm eq}}\right)\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad-\langle\sigma_{\pi X\to
2\pi}v\rangle\left(Y_{\pi}Y_{X}-\frac{Y_{\pi}^{2}}{Y_{\pi}^{\rm eq}}Y_{X}^{\rm
eq}\right)+s\langle\sigma_{3\pi\to
2\pi}v^{2}\rangle\left(Y_{\pi}^{3}-Y_{\pi}^{2}Y_{\pi}^{\rm eq}\right)\right],$
(19) $\displaystyle\frac{dY_{X}}{dx}$
$\displaystyle=\frac{1}{3H(x)}\frac{ds}{dx}\left[2\langle\sigma_{XX\to
2\pi}v\rangle\left(Y_{X}^{2}-(\frac{Y_{\pi}}{Y^{\rm eq}_{\pi}}Y_{X}^{\rm
eq})^{2}\right)-s\langle\sigma_{3\pi\to\pi
X}v^{2}\rangle\left(Y_{\pi}^{3}-Y_{\pi}(Y_{\pi}^{\rm
eq})^{2}\frac{Y_{X}}{Y_{X}^{\rm eq}}\right)\right.$
$\displaystyle\left.\qquad\qquad\qquad\qquad+\langle\sigma_{\pi X\rightarrow
2\pi}v\rangle\left(Y_{\pi}Y_{X}-\frac{Y_{\pi}^{2}}{Y_{\pi}^{\rm eq}}Y_{X}^{\rm
eq}\right)\right]\ ,$ (20)

where the first (second) lines of each equation is due to by even (odd)
numbered interactions. In our principal mass region of interest
$m_{\pi}=O(100\,\mathrm{MeV})$ freeze-out happens late enough that to
reasonable approximation, we may neglect the variation of the relativistic
degrees of freedom (including DM) and use $ds/dx=-3s/x$ when seeking
analytical solutions; in our numerical treatment, we account for the full
evolution. The dark matter mass density is then obtained from the mass-
weighted sum of their freeze-out yields,

$\displaystyle\Omega_{\rm dm}=\frac{s_{0}}{\rho_{\rm
c}}(m_{X}Y_{X,0}+m_{\pi}Y_{\pi,0})\,,$ (21)

where the “0” superscript denotes the quantities today. In practice,
$Y_{X,0}\ll Y_{\pi,0}$ or $Y_{X,0}$ is vanishing altogether because $X$ has
decayed.

##  Appendix B Cross sections and collision integrals

###  B.1 Interaction terms and Feynman Rules

The leading-order chiral Langrangian (3) of the main text, expanded to six-
point interaction terms and dropping mass and canonically normalized kinetic
terms reads,

$\displaystyle\mathcal{L}_{\rm int}^{\rm even}$
$\displaystyle=-\frac{2}{3f_{\pi}^{2}}\operatorname{Tr}\left[\pi^{2}\partial^{\mu}\pi\partial_{\mu}\pi-\pi\partial^{\mu}\pi\pi\partial_{\mu}\pi\right]+\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\operatorname{Tr}\left[\pi^{4}\right]$
$\displaystyle\quad+\frac{16}{45f_{\pi}^{4}}\operatorname{Tr}\left[\pi^{2}\left(\pi^{2}\partial^{\mu}\pi\partial_{\mu}\pi-\pi\partial^{\mu}\pi\pi\partial_{\mu}\pi\right)\right]-\frac{2m_{\pi}^{2}}{45f_{\pi}^{4}}\operatorname{Tr}\left[\pi^{6}\right]+\dots$
(22)
$\displaystyle=-\frac{1}{12f_{\pi}^{2}}\sum_{kn}\left(\pi_{k}\pi_{k}\partial_{\mu}\pi_{n}\partial^{\mu}\pi_{n}-\pi_{k}\partial_{\mu}\pi_{k}\pi_{n}\partial^{\mu}\pi_{n}\right)+\frac{m_{\pi}^{2}}{48f_{\pi}^{2}}\left(\sum_{k}\pi_{k}^{2}\right)^{2}$
$\displaystyle\quad+\frac{1}{180f_{\pi}^{4}}\sum_{jkn}\pi_{j}^{2}\left(\pi_{k}\pi_{k}\partial_{\mu}\pi_{n}\partial^{\mu}\pi_{n}-\pi_{k}\partial_{\mu}\pi_{k}\pi_{n}\partial^{\mu}\pi_{n}\right)-\frac{m_{\pi}^{2}}{2880f_{\pi}^{4}}\left(\sum_{k}\pi_{k}^{2}\right)^{3}+\dots\,.$
(23)

The sums run over $i,j,\dots=1,\dots,N_{\pi}$. For $\pi\in SU(4)/Sp(4)$, i.e.
specializing to $N_{f}=2$, there are $N_{\pi}=5$ pseudo-Nambu-Goldstone
fields.

A technical aspect of computing the relevant cross sections is the
combinatorics of fields only distinguished by their flavor index. We may
obtain the relevant Feynman rules (up to a global factor $i$) for the four-
and six-point interactions directly by contracting the respective operators
with the multiparticle state $|p_{1,a_{1}}p_{2,a_{2}}\dots
p_{m,a_{m}}\dots\rangle$, which is to be understood as an ordered list of
momenta: $m$ labels the particle by its position in the multiparticle state
and $a_{m}$ denotes its flavor. Each momentum is defined with its direction
pointing towards the interaction vertex. We start with the simplest term, the
quartic interaction with constant coupling $m_{\pi}^{2}/3f_{\pi}^{2}$,

$\displaystyle\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\langle
0|\operatorname{Tr}\left[\pi^{4}\right]|p_{1,a_{1}}p_{2,a_{2}}p_{3,a_{3}}p_{4,a_{4}}\rangle$
$\displaystyle=\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\sum_{ijkl}\operatorname{Tr}\left[T^{i}T^{j}T^{k}T^{l}\right]\langle
0|\pi_{i}\pi_{j}\pi_{k}\pi_{l}|p_{1,a_{1}}p_{2,a_{2}}p_{3,a_{3}}p_{4,a_{4}}\rangle$
$\displaystyle=\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\sum_{ijkl}\operatorname{Tr}\left[T^{i}T^{j}T^{k}T^{l}\right]\,\sum_{\sigma}\delta_{i,\sigma(a_{1})}\delta_{j,\sigma(a_{2})}\delta_{k,\sigma(a_{3})}\delta_{l,\sigma(a_{4})}\,,$
$\displaystyle=\frac{m_{\pi}^{2}}{3f_{\pi}^{2}}\sum_{\sigma}\operatorname{Tr}\left[T^{\sigma(a_{1})}T^{\sigma(a_{2})}T^{\sigma(a_{3})}T^{\sigma(a_{4})}\right]\,,$
(24)

where $\sigma=\sigma(a_{1})\,\sigma(a_{2})\,\sigma(a_{3})\,\sigma(a_{4})$ sums
up all $4!$ permutations of $a_{1}\,a_{2}\,a_{3}\,a_{4}$. For the kinetic
quartic interaction, derivatives $\partial^{\mu}\pi_{\sigma(a_{m})}$ are
replaced by $ip^{\mu}_{m,\sigma(a_{m})}\pi_{\sigma(a_{m})}$ and we write
$p_{m,\sigma(a_{m})}\equiv p_{m}$ below for brevity. We obtain,

$\displaystyle-\frac{2}{3f_{\pi}^{2}}\langle
0|\operatorname{Tr}\left[\pi^{2}\partial^{\mu}\pi\partial_{\mu}\pi-\pi\partial^{\mu}\pi\pi\partial_{\mu}\pi\right]|p_{1,a_{1}}...\rangle$
$\displaystyle=-\frac{2}{3f_{\pi}^{2}}\sum_{ijkl}\operatorname{Tr}\left[T^{i}T^{j}T^{k}T^{l}\right]i^{2}(p_{k}\cdot
p_{l}-p_{j}\cdot p_{l})\langle
0|\pi_{i}\pi_{j}\pi_{k}\pi_{l}|p_{1,a_{1}}...\rangle$
$\displaystyle=\frac{2}{3f_{\pi}^{2}}\sum_{\sigma}\operatorname{Tr}\left[T^{\sigma(a_{1})}T^{\sigma(a_{2})}T^{\sigma(a_{3})}T^{\sigma(a_{4})}\right]\,(p_{3}\cdot
p_{4}-p_{2}\cdot p_{4})\,,$ (25)

For the concrete model of $SU(4)/Sp(4)$ the trace of generators of the coset
space evaluates to

$\operatorname{Tr}\left[T^{a_{1}}T^{a_{2}}T^{a_{3}}T^{a_{4}}\right]={1\over
16}\left(\delta_{a_{1},a_{2}}\delta_{a_{3},a_{4}}+\delta_{a_{1},a_{4}}\delta_{a_{2},a_{3}}-\delta_{a_{1},a_{3}}\delta_{a_{2},a_{4}}\right)\,.$
(26)

This allows us to write the four-point vertex explicitly as

$\displaystyle\Gamma_{\left\\{p_{1},a_{1}\right\\},\left\\{p_{2},a_{2}\right\\},\left\\{p_{3},a_{3}\right\\},\left\\{p_{4},a_{4}\right\\}}=\frac{1}{6f_{\pi}^{2}}\big{\\{}$
$\displaystyle\delta_{a_{1},a_{2}}\delta_{a_{3},a_{4}}\left(2p_{1}\cdot
p_{2}+2p_{3}\cdot p_{4}-p_{1}\cdot p_{3}-p_{1}\cdot p_{4}-p_{2}\cdot
p_{3}-p_{2}\cdot p_{4}+m_{\pi}^{2}\right)$ $\displaystyle+$
$\displaystyle\delta_{a_{1},a_{3}}\delta_{a_{2},a_{4}}\left(2p_{1}\cdot
p_{3}+2p_{2}\cdot p_{4}-p_{1}\cdot p_{2}-p_{1}\cdot p_{4}-p_{2}\cdot
p_{3}-p_{2}\cdot p_{4}+m_{\pi}^{2}\right)$ $\displaystyle+$
$\displaystyle\delta_{a_{1},a_{4}}\delta_{a_{2},a_{3}}\left(2p_{1}\cdot
p_{4}+2p_{2}\cdot p_{3}-p_{1}\cdot p_{2}-p_{1}\cdot p_{3}-p_{2}\cdot
p_{4}-p_{2}\cdot p_{4}+m_{\pi}^{2}\right)\big{\\}}\,,$ (27)

where the kinetic interaction induces the terms involving momentum scalar
products, while the last terms proportional to $m_{\pi}^{2}$ are from ${\rm
Tr}[\pi^{4}]$.

Regarding Feynman rules of six-point interactions, we similarly obtain

$\sum_{\sigma}\left[-\frac{2m_{\pi}^{2}}{45f_{\pi}^{4}}-\frac{16}{45f_{\pi}^{4}}(p_{5}\cdot
p_{6}-p_{4}\cdot
p_{6})\right]\operatorname{Tr}\left[T^{\sigma(a_{1})}T^{\sigma(a_{2})}T^{\sigma(a_{3})}T^{\sigma(a_{4})}T^{\sigma(a_{5})}T^{\sigma(a_{6})}\right]\,,$
(28)

where
$\sigma=\sigma(a_{1})\,\sigma(a_{2})\,\sigma(a_{3})\,\sigma(a_{4})\,\sigma(a_{5})\,\sigma(a_{6})$
sums up the $6!=720$ permutations of
$a_{1}\,a_{2}\,a_{3}\,a_{4}\,a_{5}\,a_{6}$. The prefactor is insensitive to
switching the order of $a_{1}$ and $a_{2}$, i.e., for $SU(4)/Sp(4)$ we can use
the equality

$\sum_{\sigma}\operatorname{Tr}\left[{T^{a_{1}}T^{a_{2}}+T^{a_{2}}T^{a_{2}}\over
2}T^{a_{3}}T^{a_{4}}T^{a_{5}}T^{a_{6}}\right]=\sum_{\sigma}{\delta_{a_{1},a_{2}}\over
8}\operatorname{Tr}\left[T^{a_{3}}T^{a_{4}}T^{a_{5}}T^{a_{6}}\right]$ (29)

to reduce the number of generators in the trace. Then, replacing $a_{i}$
$(i=1,..,6)$ by $(a,b,c,d,f,g)$ results in

$\displaystyle\Gamma_{\rm 6-point}|_{\rm
contact}=-\frac{m_{\pi}^{2}}{60f_{\pi}^{4}}$
$\displaystyle(\delta_{a,g}\delta_{b,f}\delta_{c,d}+\delta_{a,f}\delta_{b,g}\delta_{c,d}+\delta_{a,g}\delta_{b,d}\delta_{c,f}+\delta_{a,d}\delta_{b,g}\delta_{c,f}+\delta_{a,f}\delta_{b,d}\delta_{c,g}+\delta_{a,d}\delta_{b,f}\delta_{c,g}$
$\displaystyle+\delta_{a,g}\delta_{b,c}\delta_{d,f}+\delta_{a,c}\delta_{b,g}\delta_{d,f}+\delta_{a,b}\delta_{c,g}\delta_{d,f}+\delta_{a,f}\delta_{b,c}\delta_{d,g}+\delta_{a,c}\delta_{b,f}\delta_{d,g}+\delta_{a,b}\delta_{c,f}\delta_{d,g}$
$\displaystyle+\delta_{a,b}\delta_{c,d}\delta_{f,g}+\delta_{a,d}\delta_{b,c}\delta_{f,g}+\delta_{a,c}\delta_{b,d}\delta_{f,g})\,$
(30)

for the contact interactions originating from the mass term in the chiral
Lagrangian. For the six-point interaction originating from the kinetic term,
we further simplify the expression by tailoring it to the process
$X_{ab}+X_{cd}\to\pi_{f}+\pi_{g}$, taking the non-relativistic limit for
initial states (and thus $p_{f}\cdot p_{g}=2m_{X}^{2}-m_{\pi}^{2}$ for final
states)

$\displaystyle\Gamma_{\rm 6-point}|_{\rm
kinetic}=\frac{13m_{X}^{2}-4m_{\pi}^{2}}{180f_{\pi}^{4}}$
$\displaystyle(\delta_{a,g}\delta_{b,f}\delta_{c,d}+\delta_{a,f}\delta_{b,g}\delta_{c,d}+\delta_{a,g}\delta_{b,d}\delta_{c,f}+\delta_{a,d}\delta_{b,g}\delta_{c,f}+\delta_{a,f}\delta_{b,d}\delta_{c,g}+\delta_{a,d}\delta_{b,f}\delta_{c,g}$
$\displaystyle+\delta_{a,g}\delta_{b,c}\delta_{d,f}+\delta_{a,c}\delta_{b,g}\delta_{d,f}+\delta_{a,b}\delta_{c,g}\delta_{d,f}+\delta_{a,f}\delta_{b,c}\delta_{d,g}+\delta_{a,c}\delta_{b,f}\delta_{d,g}+\delta_{a,b}\delta_{c,f}\delta_{d,g}$
$\displaystyle-4\delta_{a,b}\delta_{c,d}\delta_{f,g}-4\delta_{a,d}\delta_{b,c}\delta_{f,g}-4\delta_{a,c}\delta_{b,d}\delta_{f,g})\,.$
(31)

Regarding the odd-numbered case, when the WZW interaction is present, there
may exist additional mass-depletion processes, such as $\pi\pi\pi\to\pi\pi$
and $\pi X\to\pi\pi$. Concretely, the WZW interaction is given by

$\displaystyle\mathcal{L}_{\rm int}^{\rm odd}=$
$\displaystyle\frac{2N_{c}}{15\pi^{2}f_{\pi}^{5}}\epsilon_{\mu\nu\rho\sigma}\operatorname{Tr}\left[\pi\partial^{\mu}\pi\partial^{\nu}\pi\partial^{\rho}\pi\partial^{\sigma}\pi\right]+\dots$
(32) $\displaystyle=$
$\displaystyle\frac{N_{c}}{10\sqrt{2}\pi^{2}f_{\pi}^{5}}\epsilon_{\mu\nu\rho\sigma}\left(\pi_{1}\partial^{\mu}\pi_{2}\partial^{\nu}\pi_{3}\partial^{\rho}\pi_{4}\partial^{\sigma}\pi_{5}+\pi_{2}\partial^{\mu}\pi_{3}\partial^{\nu}\pi_{4}\partial^{\rho}\pi_{5}\partial^{\sigma}\pi_{1}+\pi_{3}\partial^{\mu}\pi_{4}\partial^{\nu}\pi_{5}\partial^{\rho}\pi_{1}\partial^{\sigma}\pi_{2}\right.$
$\displaystyle~{}\left.\qquad\qquad\qquad\qquad+\pi_{4}\partial^{\mu}\pi_{5}\partial^{\nu}\pi_{1}\partial^{\rho}\pi_{2}\partial^{\sigma}\pi_{3}+\pi_{5}\partial^{\mu}\pi_{1}\partial^{\nu}\pi_{2}\partial^{\rho}\pi_{3}\partial^{\sigma}\pi_{4}\right)+\dots\,,$
(33)

where the ellipses stand for higher order odd-numbered interactions. The
Lagrangian above leads to non-vanishing mass-depletion processes if all
involved $\pi$’s have different flavors. In this case, its Feynman rule is
given by

$\displaystyle\Gamma_{\left\\{p_{1},a\right\\},\left\\{p_{2},b\right\\},\left\\{p_{3},c\right\\},\left\\{p_{4},d\right\\},\left\\{p_{5},e\right\\}}=\frac{N_{c}}{2\sqrt{2}\pi^{2}f_{\pi}^{5}}\epsilon_{\mu\nu\rho\sigma}p_{2}^{\mu}p_{3}^{\nu}p_{4}^{\rho}p_{5}^{\sigma}\quad(a\neq
b\neq c\neq d\neq e)\,.$ (34)

Again, all momenta are defined with their direction pointing towards the
interaction vertex.

###  B.2 Bound state formation $3\pi\to\pi X$

Figure 3: Feynman diagrams contributing to $3\pi\to\pi X$. The $t$\- and
$u$-channel type diagrams (left and middle) are enhanced in the non-
relativistic regime, relative to $s$-channel (top left) and final 6-point
contact interaction (bottom left) diagrams.

With all Feynman rules at hand, we now compute the collision term for bound
state formation $3\pi\to\pi X$. The process is dominated by the first and
second $t$\- and $u$-channel diagrams of both rows of Fig. 3. For the first
two, the amplitudes read,

$\displaystyle\mathcal{M}_{t}=$
$\displaystyle\frac{1}{k_{t}^{2}-m_{\pi}^{2}}\sum_{d}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\sqrt{2M_{X}}}{2m_{\pi}}\widetilde{\psi}^{*}(\vec{q})\,\Gamma_{\left\\{a,p_{1}\right\\},\left\\{b,p_{2}\right\\},\left\\{e,-k_{1}\right\\},\\{d,-k_{t}\\}}\Gamma_{\\{d,k_{t}\\},\left\\{c,p_{3}\right\\},\left\\{f,-\frac{K}{2}-q\right\\},\left\\{g,-\frac{K}{2}+q\right\\}}\,,$
(35) $\displaystyle\mathcal{M}_{u}=$
$\displaystyle\frac{1}{k_{u}^{2}-m_{\pi}^{2}}\sum_{d}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\sqrt{2M_{X}}}{2m_{\pi}}\widetilde{\psi}^{*}(\vec{q})\,\Gamma_{\left\\{a,p_{1}\right\\},\left\\{b,p_{2}\right\\},\left\\{f,-\frac{K}{2}-q\right\\},\\{d,-k_{u}\\}}\Gamma_{\\{d,k_{u}\\},\left\\{c,p_{3}\right\\},\left\\{e,-k_{1}\right\\},\left\\{g,-\frac{K}{2}+q\right\\}}\,,$
(36)

where $k_{t}=p_{1}+p_{2}-k_{1}$ and $k_{u}=p_{1}+p_{2}-\frac{K}{2}-q$ are the
four-momenta of the $t$ and $u$-channel propagators, respectively. In the non-
relativistic limit and using 3-body kinematics in the center of momentum (CM)
frame, it can be shown that $k_{t}^{2}-m_{\pi}^{2}\simeq-8E_{B}m_{\pi}/3$ and
$k_{u}^{2}-m_{\pi}^{2}\simeq-E_{B}m_{\pi}/3$, respectively, leading to near
on-shell enhancements of both diagrams. In all simplifications we assume that
$E_{B}>T$ so that $m_{\pi}v^{2}/2E_{B}$ remains a small parameter where $v$ is
a typical DM velocity, $v\sim\sqrt{T/m_{\pi}}$. In addition, the wave-function
of each bound state, assumed to be flavor-blind here, is normalized as

$\int\frac{d^{3}q}{(2\pi)^{3}}\>\widetilde{\psi}^{*}(\vec{q})\widetilde{\psi}(\vec{q})=\int
d^{3}r\ \psi^{*}(\vec{r})\psi(\vec{r})=1\,,\text{~{}~{}together
with~{}~{}}\widetilde{\psi}({\vec{q}})\equiv\int d^{3}r\
\psi({\vec{r}})e^{-i{\vec{q}\cdot\vec{r}}}\,.$ (37)

Thus, for constant terms of $\Gamma$, one can simplify the amplitude with the
equality

$\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\sqrt{2M_{X}}}{2m_{\pi}}\widetilde{\psi}^{*}(\vec{q})=\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\sqrt{2M_{X}}}{2m_{\pi}}\int
d^{3}r\psi^{*}({\vec{r}})e^{i{\vec{q}\cdot\vec{r}}}=\frac{\sqrt{2M_{X}}}{2m_{\pi}}\int
d^{3}r\psi^{*}({\vec{r}})\delta^{3}({\vec{r}})=\frac{\sqrt{2M_{X}}}{2m_{\pi}}\psi^{*}(0)\,.$
(38)

The $3\pi\rightarrow\pi X$ cross section of negative mass dimension five may
then be written in the CM frame as

$\displaystyle\sigma_{3\pi\rightarrow\pi X}v^{2}\,$
$\displaystyle\equiv\frac{1}{2E_{p_{1}}\,2E_{p_{2}}\,2E_{p_{3}}}\int\frac{d^{3}k_{1}}{(2\pi)^{3}2E_{k_{1}}}\frac{d^{3}K}{(2\pi)^{3}2E_{K}}(2\pi)^{4}\delta^{4}\left(p_{1}+p_{2}+p_{3}-k_{1}-K\right)|\mathcal{M}_{3\pi\rightarrow\pi
X}|^{2}$ (39)
$\displaystyle\simeq\frac{1}{128\pi^{2}E_{p_{1}}E_{p_{2}}E_{p_{3}}}\frac{|\vec{k}_{1}|}{(E_{p_{1}}+E_{p_{2}}+E_{p_{3}})}\int
d\Omega\,|\mathcal{M}_{3\pi\rightarrow\pi X}|^{2}\,,$ (40)

where $\mathcal{M}_{3\pi\rightarrow\pi X}$ is obtained from the sum of the
four $t$\- and $u$-channel-type diagrams of Fig. 3 (the specific flavor
combinations will be discussed momentarily); $d\Omega$ is the differential
final state solid angle. The overall collision term or “annihilation cross
section” is defined by the physical number-changing rate of free $\pi$ as it
enters our Boltzmann equation in the previous section, i.e., in the absence of
Hubble expansion

${dn_{\pi}}=\sum_{k}{dn_{\pi_{k}}}=2n_{\pi}^{3}\langle\sigma_{3\pi\rightarrow\pi
X}v^{2}\rangle{dt}\,.$ (41)

By adding up all distinct processes that contribute to $dn_{\pi}/dt$, we
obtain

$\displaystyle n_{\pi}^{3}\langle\sigma_{3\pi\rightarrow\pi
X}v^{2}\rangle=\sum_{a>b>c}n_{\pi_{a}}n_{\pi_{b}}n_{\pi_{c}}\langle\sigma_{\pi_{a}\pi_{b}\pi_{c}\rightarrow\pi
X}v^{2}\rangle+\sum_{a\neq b}{n_{\pi_{a}}^{2}\over
2!}n_{\pi_{b}}\langle\sigma_{\pi_{a}\pi_{a}\pi_{b}\rightarrow\pi
X}v^{2}\rangle+\sum_{a}{n_{\pi_{a}}^{3}\over
3!}\langle\sigma_{\pi_{a}\pi_{a}\pi_{a}\rightarrow\pi X}v^{2}\rangle\,,$ (42)

where the prefactors $2!$ and $3!$ avoid double-counting of identical initial
states. As in the thermal bath
$n_{\pi_{a}}=n_{\pi_{b}}=n_{\pi_{c}}=n_{\pi}/N_{\pi}$ holds, our thermally-
averaged cross section can be written as

$\displaystyle\langle\sigma_{3\pi\rightarrow\pi X}v^{2}\rangle={1\over
N_{\pi}^{3}}{1\over 2!}\sum_{e,f,g}\left[\sum_{a\neq b\neq c}{1\over
3!}\langle\sigma_{\pi_{a}\pi_{b}\pi_{c}\rightarrow\pi_{e}X_{fg}}v^{2}\rangle+\sum_{a\neq
b}{1\over
2!}\langle\sigma_{\pi_{a}\pi_{a}\pi_{b}\rightarrow\pi_{e}X_{fg}}v^{2}\rangle+\sum_{a}{1\over
3!}\langle\sigma_{\pi_{a}\pi_{a}\pi_{a}\rightarrow\pi_{e}X_{fg}}v^{2}\rangle\,\right]\,,$
(43)

where $N_{\pi}=5$, and $\frac{1}{2!}\sum_{e,f,g}$ sums up the final states,
taking into account that $X_{fg}=X_{gf}$. For the specific flavor
combinations, the following matrix elements enter

$\displaystyle\mathcal{M}_{\pi_{a}\pi_{b}\pi_{c}\rightarrow\pi_{e}X_{fg}}$
$\displaystyle=(\mathcal{M}_{t}+\mathcal{M}_{u})+({\left\\{a,p_{1}\right\\}\leftrightarrow\left\\{c,p_{3}\right\\}})+({\left\\{b,p_{2}\right\\}\leftrightarrow\left\\{c,p_{3}\right\\}})\,,$
(44)
$\displaystyle\mathcal{M}_{\pi_{a}\pi_{a}\pi_{c}\rightarrow\pi_{e}X_{fg}}$
$\displaystyle=(\mathcal{M}_{t}+\mathcal{M}_{u})|_{b\to
a}+({\left\\{a,p_{1}\right\\}\leftrightarrow\left\\{c,p_{3}\right\\}})|_{b\to
a}\,,$ (45)
$\displaystyle\mathcal{M}_{\pi_{a}\pi_{a}\pi_{a}\rightarrow\pi_{e}X_{fg}}$
$\displaystyle=(\mathcal{M}_{t}+\mathcal{M}_{u})|_{b\to a,\,c\to a}\,,$ (46)

where $\mathcal{M}_{t,u}$ have been defined above, and
$({\left\\{...\right\\}\leftrightarrow\left\\{...\right\\}})$ stands for
switching the corresponding two particles in the first term on the right hand
side (RHS) of each corresponding equation. Putting everything together, we
obtain the effective cross section for the bound state formation rate as

$\langle\sigma_{3\pi\rightarrow\pi
X}v^{2}\rangle\simeq\frac{\numprint{57041}|\psi(0)|^{2}}{\numprint{327680}\sqrt{3}\pi
f_{\pi}^{8}}\left({m_{\pi}\over E_{B}}\right)^{3/2}+O\left(\sqrt{m_{\pi}\over
E_{B}}\right)=\frac{\numprint{57041}R(0)^{2}}{\numprint{1310720}\sqrt{3}\pi^{2}f_{\pi}^{8}}\left({m_{\pi}\over
E_{B}}\right)^{3/2}+O\left(\sqrt{m_{\pi}\over E_{B}}\right)\,.$ (47)

This equation is applicable in the non-relativistic limit, where the ground-
state wavefunction of the two-particle bound state is
$\psi(r)=R(r)/\sqrt{4\pi}$. The numerical prefactor in the last equality is
approximately $3\times 10^{-3}$. Note that we have neglected the leftmost
diagrams of Fig. 3, which are not enhanced by inverse powers of binding
energy.

Figure 4: Feynman diagrams contributing to $3\pi\to\pi X$. Some of the
diagrams are identical when $a=b$ or $c=d$.

###  B.3 Annihilation $XX\to\pi\pi$

Now, we turn to the mass-depleting process $XX\to\pi\pi$ mediated by even-
numbered interactions. Here, we take into account all diagrams shown in Fig.
4, among which the amplitude of the first one is given by

$\mathcal{M}_{abf,cdg}=\frac{1}{k^{2}-m_{\pi}^{2}}\sum_{e}\int\frac{d^{3}\vec{q}_{1}d^{3}\vec{q}_{2}}{(2\pi)^{6}}\frac{2M_{X}\widetilde{\psi}(\vec{q}_{1})\widetilde{\psi}(\vec{q}_{2})}{(2m_{\pi})^{2}}\,\Gamma_{\left\\{a,{K_{1}\over
2}+q_{1}\right\\},\left\\{b,{K_{1}\over
2}-q_{1}\right\\},\left\\{f,-k_{1}\right\\},\\{e,-k\\}}\Gamma_{\left\\{c,{K_{2}\over
2}+q_{2}\right\\},\left\\{d,{K_{2}\over
2}-q_{2}\right\\},\left\\{e,k\right\\},\left\\{g,-k_{2}\right\\}}\,,$ (48)

where the subscripts $abf$ and $cdg$ denote the incoming and outgoing states
that are connected by the same vertex. The 4-momentum of the intermediate
state is given by $k=(K_{1}/2+q_{1})+(K_{1}/2-q_{1})-k_{1}=K_{1}-k_{1}$,
pointing from the $f$\- to the $g$-associated vertex. The amplitudes of the
following five diagrams can be obtained by switching the momenta and indices
in both $\Gamma$’s and $k$, labelled as $\mathcal{M}_{acf,bdg}$,
$\mathcal{M}_{bcf,adg}$, $\mathcal{M}_{adf,bcg}$, $\mathcal{M}_{cdf,abg}$,
$\mathcal{M}_{bdf,acg}$, respectively. The amplitude of the last one,
$\mathcal{M}_{abcdfg}$, is the sum of $\Gamma_{\rm 6-point}|_{\rm contact}$
and $\Gamma_{\rm 6-point}|_{\rm kinetic}$ that have been given above.

Again we define the thermally-averaged cross section such that it yields the
physical mass-depletion rate, i.e.,
${dn_{X}}/{dt}=2n_{X}^{2}\langle\sigma_{2X\rightarrow 2\pi}v\rangle$ in the
absence of Hubble expansion. Here, $n_{\pi}$ and $n_{X}$ sum up all flavors of
free and bound states. Denoting by $N_{X}=N_{\pi}(N_{\pi}+1)/2=15$ the number
of flavor-blind bound states we have

$\displaystyle\langle\sigma_{2X\rightarrow 2\pi}v\rangle=$
$\displaystyle{1\over
N_{X}^{2}}\sum_{f>g}\left[\sum_{a=c>b>d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle+\sum_{b>a=c>d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle+\sum_{b>d>a=c}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle\right]$
$\displaystyle+{1\over
N_{X}^{2}}\sum_{f>g}\left[\sum_{a=b>c>d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle+\sum_{c>a=b>d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle+\sum_{c>d>a=b}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle\right]$
$\displaystyle+{1\over
N_{X}^{2}}\sum_{f>g}\left[\sum_{a=b=c>d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle+\sum_{d>a=b=c}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{g}}v\rangle\right]$
$\displaystyle+{1\over
N_{X}^{2}}\sum_{f}\left[\sum_{a=b>c=d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{f}}v\rangle+{1\over
2!}\sum_{a=b=c=d}\langle\sigma_{X_{ab}X_{cd}\rightarrow\pi_{f}\pi_{f}}v\rangle\right]\,,$
(49)

where there are no $f=g$ ($f>g$) terms in the first three lines (last line),
and $a,\,b,\,c,\,d$ cannot be all different from each other, because all even-
numbered interaction are flavor-conserving in this concrete model. In fact,
the final states are uniquely fixed by the flavors of initial states once
$f\geq g$ is chosen. In the last line, a factor of $1/2!$ avoids double-
counting of the identical initial states. Meanwhile, the scattering angle
ranges only from $0$ to $\pi/2$ for identical final states. Finally, note that
$\sum_{a=c>b>d}+\sum_{b>a=c>d}+\sum_{b>d>a=c}$ in the first line can be
unified to $\sum_{a=c\neq b\neq d}/2$. However, it then becomes less obvious
that there is no double-counting issue with the latter expressions. A similar
simplification applies to the second and third lines. For instance, the third
line corresponds to $\sum_{a=b=c\neq d}$.

As shown in Fig. 4, some diagrams are the same when $a=b$, as the effect of
switching two identical particles in a bound state is absorbed by the
corresponding wavefunction. Moreover, when $a=b$ and $c=d$, only the first,
second, fifth, and last diagrams are distinct. Taking this into account, we
express the general amplitude for Eq. (B.3) as

$\mathcal{M}_{X_{ab}X_{cd}\to\pi_{e}\pi_{f}}=\mathcal{M}_{abf,cdg}+\mathcal{M}_{acf,bdg}+(1-\delta_{a,b})\mathcal{M}_{bcf,adg}+(1-\delta_{c,d})\mathcal{M}_{adf,bcg}+\mathcal{M}_{cdf,abg}+(1-\delta_{a,b})\mathcal{M}_{bdf,acg}+\mathcal{M}_{abcdfg}\,.$

With this amplitude, the effective cross section for the mass-depletion rate
can be calcualted to yield

$\displaystyle\langle\sigma_{2X\rightarrow 2\pi}v\rangle\simeq$
$\displaystyle\frac{|\psi(0)|^{4}}{f_{\pi}^{8}}\sqrt{1-{m_{\pi}^{2}\over
m_{X}^{2}}}\left(\frac{\numprint{118790261}\,m_{X}^{8}-\numprint{104839888}\,m_{\pi}^{2}m_{X}^{6}+\numprint{36685056}\,m_{\pi}^{4}m_{X}^{4}-\numprint{1573120}\,m_{\pi}^{6}m_{X}^{2}+\numprint{4444800}\,m_{\pi}^{8}}{\numprint{23887872000}\,\pi\,m_{\pi}^{4}m_{X}^{4}}\right)$
$\displaystyle\simeq$
$\displaystyle\frac{\numprint{2529757}}{\numprint{26542080}\sqrt{3}\pi}\frac{|\psi(0)|^{4}}{f_{\pi}^{8}}=\frac{\numprint{2529757}}{\numprint{424673280}\sqrt{3}\pi^{3}}\frac{R(0)^{4}}{f_{\pi}^{8}}\,$
(50)

in the non-relativistic limit, where we have taken $m_{X}\simeq 2m_{\pi}$ and
$\psi(0)=R(0)/\sqrt{4\pi}$ for the last two equalities. The numerical
prefactor in the last equality is $10^{-4}$.

###  B.4 Standard $3\pi\to 2\pi$ and catalyzed $\pi X\to\pi\pi$ process

The final set of cross sections used in this work are the ones induced by the
WZW term. Both, the free $3\pi\to 2\pi$ and its bound-state version $\pi
X\to\pi\pi$, are velocity-dependent in the non-relativistic freeze-out regime.
We therefore have to include momentum-dependent terms in the amplitudes. The
case involving only free states, $3\pi\to 2\pi$, is simple, as its amplitude
is directly given by the vertex factor,

$\mathcal{M}_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}=\Gamma_{\left\\{p_{1},a\right\\},\left\\{p_{2},b\right\\},\left\\{p_{3},c\right\\},\left\\{-k_{1},d\right\\},\left\\{-k_{2},e\right\\}}\,,$
(51)

where $p_{i}$ and $k_{i}$ are the momenta of incoming and outgoing states,
respectively. Taking the definition of a $3\to 2$ cross section above, Eq.
(39), one obtains in the CM frame

$\displaystyle\sigma_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}v^{2}$
$\displaystyle=\frac{1}{128\pi^{2}E_{p_{1}}E_{p_{2}}E_{p_{3}}}\frac{|\vec{k}_{1}|}{\sqrt{s}}\int
d\Omega\,|\mathcal{M}_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}|^{2}$
$\displaystyle=\frac{N_{c}^{2}}{3072\pi^{5}f_{\pi}^{10}}\frac{s(s-4m_{\pi}^{2})^{3/2}\left[s_{12}s_{13}s_{23}-m_{\pi}^{2}(s-m_{\pi}^{2})^{2}\right]}{(s+m_{\pi}^{2}-s_{12})(s+m_{\pi}^{2}-s_{13})(s+m_{\pi}^{2}-s_{23})}\,,$
(52)

where $s=(p_{1}+p_{2}+p_{3})^{2}$, $s_{ij}=(p_{i}+p_{j})^{2}$. Note that
$a,b,c,d,e$ are all different. The factor
$\left[s_{12}s_{13}s_{23}-m_{\pi}^{2}(s-m_{\pi}^{2})^{2}\right]$ vanishes in
the limit $s_{ij}\to 4m_{\pi}^{2}$ and $s\to 9m_{\pi}^{2}$. Therefore, one
requires an expansion of the Mandelstam variables $s$ and $s_{ij}$ to higher
order in velocities.

Without including relativistic effects, we calculate the thermally-averaged
cross section as follows:

$\left\langle\sigma_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}v^{2}\right\rangle=\frac{\int
d^{3}\vec{p}_{1}d^{3}\vec{p}_{2}d^{3}\vec{p}_{3}\,f^{\rm
eq}(\vec{p}_{1})f^{\rm eq}(\vec{p}_{2})f^{\rm
eq}(\vec{p}_{3})\times(\sigma_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}v^{2})}{\int
d^{3}\vec{p}_{1}d^{3}\vec{p}_{2}d^{3}\vec{p}_{3}\,f^{\rm
eq}(\vec{p}_{1})f^{\rm eq}(\vec{p}_{2})f^{\rm eq}(\vec{p}_{3})}\,,$ (53)

where we adopt non-relativistic Maxwell-Boltzmann equilibrium statistics for
the three initial states distribution functions $f^{\rm eq}$. To simplify this
integral, the key is that for each set of
($\vec{p}_{1},\,\vec{p}_{2},\,\vec{p}_{3}$) we have the freedom to choose a
concrete frame where we fix the $y$-$z$ plane with
$\vec{p}_{1},\,\vec{p}_{2}$, and set $\vec{p}_{1}$ along the z-direction, so
that

$\displaystyle p_{1}=\left(E_{p_{1}},\,0,\,0,\,|\vec{p}_{1}|\right)\,,$ (54)
$\displaystyle
p_{2}=\left(E_{p_{2}},\,0,\,|\vec{p}_{2}|\sin\theta_{2},\,|\vec{p}_{2}|\cos\theta_{2}\right)\,,$
(55) $\displaystyle
p_{3}=\left(E_{p_{3}},\,|\vec{p}_{3}|\sin\theta_{3}\cos\phi_{3},\,|\vec{p}_{3}|\sin\theta_{3}\sin\phi_{3},\,|\vec{p}_{3}|\cos\theta_{3}\right)\,,$
(56)

where $E_{p_{i}}=\sqrt{p_{i}^{2}+m_{\pi}^{2}}$, $\theta_{2}$ is the angle
between ($\vec{p}_{1},\,\vec{p}_{2}$), independent of the frame after
neglecting second-order relativistic corrections, and the direction of
$\vec{p}_{3}$ is randomly distributed. That is, the distribution functions of
all six variables, ($|\vec{p}_{1}|$, $|\vec{p}_{2}|$, $|\vec{p}_{3}|$,
$\theta_{2}$, $\theta_{3}$, $\phi_{3}$), are fixed independently by the three-
dimensional Maxwell-Boltzmann distribution functions.

Those six variables can then describe the squared amplitude, in which the
Lorentz-invariant Mandelstam variables become

$\displaystyle
s_{12}=2m_{\pi}^{2}+2E_{p_{1}}E_{p_{2}}-2\sqrt{E_{p_{1}}^{2}-m_{\pi}^{2}}\sqrt{E_{p_{2}}^{2}-m_{\pi}^{2}}\cos\theta_{2}\,,$
(57) $\displaystyle
s_{13}=2m_{\pi}^{2}+2E_{p_{1}}E_{p_{3}}-2\sqrt{E_{p_{1}}^{2}-m_{\pi}^{2}}\sqrt{E_{p_{3}}^{2}-m_{\pi}^{2}}\cos\theta_{3}\,,$
(58) $\displaystyle
s_{23}=2m_{\pi}^{2}+2E_{p_{2}}E_{p_{3}}-2\sqrt{E_{p_{2}}^{2}-m_{\pi}^{2}}\sqrt{E_{p_{3}}^{2}-m_{\pi}^{2}}\left(\cos\theta_{2}\cos\theta_{3}+\sin\theta_{2}\sin\theta_{3}\sin\phi_{3}\right)\,,$
(59) $\displaystyle s=s_{12}+s_{13}+s_{23}-3m_{\pi}^{2}\,.$ (60)

As a result, the momentum-dependent part of the integral in Eq. (53) now
simplifies to

$\int
d^{3}\vec{p}_{1}d^{3}\vec{p}_{2}d^{3}\vec{p}_{3}=\int_{0}^{\infty}4\pi|\vec{p}_{1}|^{2}\,d|\vec{p}_{1}|\int_{0}^{\infty}2\pi|\vec{p}_{2}|^{2}\,d|\vec{p}_{2}|\int^{1}_{-1}d\cos\theta_{2}\int_{0}^{\infty}|\vec{p}_{3}|^{2}d|\vec{p}_{3}|\int^{1}_{-1}d\cos\theta_{3}\int^{2\pi}_{0}d\phi_{3}\,,$
(61)

which can be solved analytically. We obtain

$\left\langle\sigma_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}v^{2}\right\rangle=\frac{5\sqrt{5}N_{c}^{2}m_{\pi}^{3}T^{2}}{1024\pi^{5}f_{\pi}^{10}}+O\left({T^{4}\over
m_{\pi}^{4}}\right)\,,$ (62)

where $T$ is the temperature of the dark sector. Since this cross section is
only non-vanishing for $a\neq b\neq c$, the additional prefactor after
averaging over all flavors of initial states is the same as the first term on
the RHS of Eq. (43), resulting in

$\langle\sigma_{3\pi\rightarrow 2\pi}v^{2}\rangle={1\over N_{\pi}^{3}}{1\over
2!\,3!}\sum_{a\neq b\neq c\neq d\neq
e}\langle\sigma_{\pi_{a}\pi_{b}\pi_{c}\to\pi_{d}\pi_{e}}v^{2}\rangle={1\over
N_{\pi}^{3}}{5!\over
2!\,3!}\frac{5\sqrt{5}N_{c}^{2}m_{\pi}^{3}T^{2}}{1024\pi^{5}f_{\pi}^{10}}\,,$
(63)

where the new factor $5!$ is the number of non-equal ($a,b,c,d,e$)
combinations. The numerical factor in the last expression is approximately
$5\times 10^{-5}$. This result has an additional factor of $1/3$ with respect
to the thermally averaged cross section in Hochberg  _et al._ (2015) after
taking into account the different normalizations
$f_{\pi}=f_{\pi}^{\text{Ref.\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Hochberg:2014kqa}{\@@citephrase{(}}{\@@citephrase{)}}}}}/2$
and $T^{n}=T^{n}_{\text{Ref.\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Hochberg:2014kqa}{\@@citephrase{(}}{\@@citephrase{)}}}}}/2$,
but is in agreement with Kamada  _et al._ (2023).

Similarly, the mass-depletion process involving the bound state can be
calculated from the squared amplitude as

$\displaystyle|\mathcal{M}_{\pi_{a}X_{bc}\to\pi_{d}\pi_{e}}|^{2}=\bigg{|}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{\sqrt{2M_{X}}\widetilde{\psi}^{*}(\vec{q})}{2m_{\pi}}\Gamma_{\left\\{p,a\right\\},\left\\{P/2+q,b\right\\},\left\\{P/2-q,c\right\\},\left\\{-k_{1},d\right\\},\left\\{-k_{2},e\right\\}}\bigg{|}^{2}$
$\displaystyle\quad\quad\quad\simeq\frac{N_{c}^{2}(s-4m_{\pi}^{2})\left[m_{X}^{4}+(s-4m_{\pi}^{2})^{2}-2m_{X}^{2}(s+4m_{\pi}^{2})^{2}\right]\cos\theta_{k}^{2}}{32\pi^{4}f_{\pi}^{10}}\left({2m_{X}\over
4m_{\pi}^{2}}\,\bigg{|}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\widetilde{\psi}^{*}(\vec{q}){q_{z}}\bigg{|}^{2}\right)\,,$
(64)

in the CM frame, where the $z$ direction is defined to be perpendicular to the
plane spanned by the two vectors $\vec{p}$ and $\vec{k}_{1}$; the angle
between the two latter 3-momenta is the scattering angle $\theta_{k}$. We have
set $a\neq b\neq c\neq d\neq e$, as the cross section vanishes for all other
combinations. The leading contribution contains $q_{z}$, requiring to solve

$\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\widetilde{\psi}^{*}(\vec{q})q_{z}=\int{d^{3}r}\psi^{*}(\vec{r})\int\frac{d^{3}q}{(2\pi)^{3}}\,e^{i\vec{q}\cdot\vec{r}}{q_{z}}=\int{d^{3}r}\psi^{*}(\vec{r})\left(-i{\partial\over\partial
z}\int\frac{d^{3}q}{(2\pi)^{3}}\,e^{i\vec{q}\cdot\vec{r}}\right)=-i{\partial\over\partial
z}\psi^{*}(\vec{r})|_{\vec{r}\to 0}\,.$ (65)

Here, the first non-vanishing contribution comes from the $p$-wave bound state
with wavefunctions in the form of $\psi_{1m}(\vec{r})=Y_{1m}R(r)$, where
$m=0,\pm 1$. After averaging over the value of $m$ and using $R(r)/r|_{r\to
0}=R^{\prime}(0)$, we obtain

${1\over 3}\sum_{m=0,\pm
1}\left|\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\widetilde{\psi}_{1m}^{*}(\vec{q}){q_{z}}\right|^{2}={1\over
4\pi}R^{\prime}(0)^{2}\,.$ (66)

In addition, we have also calculated this quantity without fixing a predefined
$z$-direction, and reach

$\displaystyle\frac{1}{2\ell+1}\sum_{m}(\nabla_{i}\psi^{*}_{\ell
m})(\nabla_{j}\psi_{\ell
m})\stackrel{{\scriptstyle\ell=1}}{{=}}\frac{1}{4\pi}\begin{pmatrix}R^{\prime}(0)^{2}&0&0\\\
0&R^{\prime}(0)^{2}&0\\\ 0&0&R^{\prime}(0)^{2}\\\
\end{pmatrix}_{\text{Cartesian tensor}}\,,$ (67)

where $i,j=x,y,z$, and it agrees with Eq. (66). With this at hand, we may now
express the $2\to 2$ cross section as

$\displaystyle\sigma_{\pi_{a}X_{bc}\to\pi_{d}\pi_{e}}v=$
$\displaystyle\frac{|\vec{k}_{1}|}{64\pi^{2}E_{p}E_{P}(E_{p}+E_{P})}\int
d\Omega|\mathcal{M}_{\pi_{a}X_{bc}\to\pi_{d}\pi_{e}}|^{2}$
$\displaystyle\simeq$
$\displaystyle\frac{N_{c}^{2}R^{\prime}(0)^{2}m_{X}s^{2}}{6144\pi^{6}f_{\pi}^{10}m_{\pi}^{2}}\left[\frac{m_{X}^{4}+(s-m_{\pi}^{2})^{2}-2m_{X}^{2}(s+m_{\pi}^{2})}{s^{2}-(m_{X}^{2}-m_{\pi}^{2})^{2}}\right]\left(1-\frac{4m_{\pi}^{2}}{s}\right)^{\frac{3}{2}}\,.$
(68)

In the non-relativistic limit, we take the expansion
$s=(m_{\pi}+m_{X})^{2}+m_{\pi}m_{X}v^{2}+\mathcal{O}(v^{4})$ in terms of the
initial relative velocity $v$, leading to

$\displaystyle\sigma_{\pi_{a}X_{bc}\to\pi_{d}\pi_{e}}v=$
$\displaystyle\frac{N_{c}^{2}R^{\prime}(0)^{2}(m_{X}^{2}+2m_{X}m_{\pi}-3m_{\pi}^{2})^{\frac{3}{2}}m_{X}^{2}}{6144\pi^{6}f_{\pi}^{10}m_{\pi}(m_{X}+m_{\pi})}v^{2}\,.$
(69)

Finally, the thermally-averaged cross section obtained from the requirement
${dn_{X}}/{dt}=n_{\pi}n_{X}\langle\sigma_{\pi X\rightarrow\pi\pi}v\rangle$
(neglecting Hubble expansion) is given by

$\displaystyle\langle\sigma_{\pi X\rightarrow 2\pi}v\rangle=$
$\displaystyle{1\over N_{\pi}N_{X}}{5!\over
2!\,2!}\langle\sigma_{\pi_{a}X_{bc}\to\pi_{d}\pi_{e}}v\rangle$
$\displaystyle=$ $\displaystyle{1\over N_{\pi}N_{X}}{5!\over
2!\,2!}\frac{N_{c}^{2}R^{\prime}(0)^{2}(m_{X}^{2}+2m_{X}m_{\pi}-3m_{\pi}^{2})^{\frac{3}{2}}m_{X}^{2}}{6144\pi^{6}f_{\pi}^{10}m_{\pi}(m_{X}+m_{\pi})}\langle
v^{2}\rangle$ $\displaystyle\simeq$
$\displaystyle\frac{N_{c}^{2}R^{\prime}(0)^{2}(m_{X}^{2}+2m_{X}m_{\pi}-3m_{\pi}^{2})^{\frac{3}{2}}m_{X}}{5120\pi^{6}f_{\pi}^{10}m_{\pi}^{2}}T\,$
$\displaystyle=$
$\displaystyle\frac{\sqrt{5}N_{c}^{2}R^{\prime}(0)^{2}m_{\pi}^{2}}{512\pi^{6}f_{\pi}^{10}}T\,,$
(70)

where in the first equality we include the following symmetry factors: $5!$
for the non-equal flavor combinations of the process, $1/2!$ accounting for
$X_{bc}=X_{cb}$ and another $1/2!$ for the final states. For the next-to-last
equality we use $\langle v^{2}\rangle=3T(m_{X}^{-1}+m_{\pi}^{-1})$ and in the
last equality we make the approximation $m_{X}\simeq 2m_{\pi}$. Setting
$N_{c}=4$ yields for the numerical prefactor $7\times 10^{-5}$.

##  Appendix C Additional details on abundance evolution

In the main text we adopt a benchmark model fixing $m_{\pi}=5f_{\pi}=0.2$ GeV,
as well as $R(0)=0.3m_{\pi}^{3/2}$, $dR(0)/dr=0.03m_{\pi}^{5/2}$ to
demonstrate our numerical solutions to the abundance evolution as a function
of $x=m_{\pi}/T$. Here we provide additional quantitative information and
study the dependence on various parameters.

###  C.1 Even-numbered case: analytical and numerical solutions

The analytical solutions provided in the main text, Eq. (14) and (16), to
even- and odd-numbered cases, indicate the parametric dependencies of model
parameters. In Fig. 5, we test numerically how the evolution of the total DM
abundance changes with varying parameters and compare with our analytical
solutions. Concretely, the parameters adopted for our benchmark model in the
main text yield

$\displaystyle\sigma^{\rm LO}_{\rm el}/m_{\pi}$
$\displaystyle\simeq{0.2}\,\text{bn}/\text{GeV}\simeq{0.1}\,\text{cm}^{2}/\text{g}\,,$
(71) $\displaystyle\langle\sigma_{2X\to 2\pi}v^{2}\rangle/m_{\pi}$
$\displaystyle\simeq{6.5}\,\text{bn}/\text{GeV}\,,~{}~{}\langle\sigma_{\pi
X_{p}\to 2\pi}v\rangle/m_{\pi}\simeq{0.03x^{-1}}\,\text{bn}/\text{GeV}\,.$
(72)

For the even-numbered case, where the final chemical decoupling of bound state
formation, $3\pi\to\pi X$, happens at $x_{2}\simeq 115$, our analytical
solution, Eq. (14), yields $Y(x_{2})\simeq 0.35Y|_{\rm num}$, while its
approximation, Eq. (15), yields $Y(x_{2})\simeq 0.28Y|_{\rm num}$; here
$Y|_{\rm num}$ denotes the final yield of the full numerical solution. We see
that the simple analytical expression works to within a factor of 2-3. More
parameter choices for the even-numbered case are shown in the left panel of
Fig. 5, where we change $R(0)/m_{\pi}^{3/2}=0.35$ to $0.7$ and $1.4$, of which
the last yields approximately the observed DM relic abundance. For the last
parameter choice, one obtains $x_{2}\simeq 111$, so Eq. (14) yields
$Y(x_{2})\simeq 0.39Y|_{\rm num}$, while the approximation yields
$Y(x_{2})\simeq 0.27Y|_{\rm num}$. Interestingly, the relative differences
between the analytical and numerical results are similar in both sets of
parameters.

To study the parameter-dependence quantitatively, increasing $R(0)$ by a
factor of $2$ enhances the relevant process $\langle\sigma_{2X\to
2\pi}v^{2}\rangle$ by a factor of $2^{4}=16$. Our analytical expression
suggests that this change reduces the final abundance by $16^{1/3}\simeq 2.5$
if the freeze-out point $x_{2}$ remains. Indeed, our full numerical solution
shows an decrease of the final abundance by about a factor of $2.2$ in the
left panel of Fig. 5. Finally, we have also numerically calculated how the
relic abundance changes by doubling the binding energy to $\kappa=0.2$.
Although not demonstrated explicitly with figures here, our numerical results
show that this reduces $x_{2}$ by roughly a factor of $2$. The analytical
solution then suggests a suppression factor of $2^{-4/3}\simeq 0.40$ of the
final DM abundance. Indeed, the full numerical results yield a very similar
factor $0.42$. In summary, our analytical expression of the even-number case
estimates the final abundance within a factor of $2-3$ and shows excellent
scaling behavior with respect to a variation of model parameters.

Figure 5:  Evolution of DM abundance, $Y_{\pi}+2Y_{X}\simeq Y_{\pi}$ with
$m_{\pi}=5f_{\pi}=0.2$ GeV and binding energy $\kappa_{S}=4\kappa_{P}=0.1$.
_Left_ : even-numbered cases varying $R(0)/m_{\pi}^{3/2}=0.35,\,0.70,\,1.4$
(from top to bottom). _Middle_ : odd-numbered cases varying
$R(0)/m_{\pi}^{3/2}=0.03,\,0.12,\,0.48$ (from top to bottom) and fixing
$R^{\prime}(0)/m_{\pi}^{3/2}=0.03$. _Right_ : odd-numbered cases fixing
$R(0)/m_{\pi}^{3/2}=0.3$ and varying
$R^{\prime}(0)/m_{\pi}^{3/2}=0.01,\,0,04,\,0.16$ (from top to bottom). The
gray curve labeled “SIMP” shows the conventional WZW case without bound state
formation and the lower dashed gray line corresponds to the observed DM
abundance. Vertical lines give the corresponding values of $x_{1}$ (decoupling
point of mass depletion processes) and $x_{2}$ (decoupling point of
$3\pi\leftrightarrow\pi X$) with matching colors. In all panels $x_{1}\leq
x_{2}$ holds.

###  C.2 Odd-numbered case: analytical and numerical solutions

Now we turn to the odd-numbered case. Here, the DM abundance essentially
already freezes out at $x_{1}$, since $\pi X_{p}\to\pi\pi$ decouples with
exponential temperature dependence. The relic abundance can be estimated from
$n_{X_{p}}\langle\sigma_{\pi X\rightarrow 2\pi}v\rangle\simeq H(x_{1})$. When
$3\pi\leftrightarrow\pi X$ is strong enough at $x_{1}$, there exists the
detailed balance condition $Y_{X}={Y_{\pi}^{2}Y_{X}^{\rm eq}}/{(Y_{\pi}^{\rm
eq})^{2}}$. Using this equality at $x_{1}$, we obtain

$Y_{\pi}(x_{1})\simeq\left(\frac{91125N_{\pi}^{4}e^{-2\kappa_{P}x_{1}}x_{1}^{5}}{2048\pi^{9}g_{*}^{3}N_{X_{p}}^{2}M_{P}^{2}m_{\pi}^{2}\langle\sigma_{\pi
X_{P}\rightarrow 2\pi}v\rangle^{2}}\right)^{1/4}\,,$ (73)

which, upon evaluation, yields Eq. (16) of the main text. Comparing this
analytical solution to the numerical results is straightforward, as we show
explicitly the vertical line where $n_{X_{p}}\langle\sigma_{\pi
X_{P}\rightarrow 2\pi}v\rangle\simeq H(x_{1})$ is reached in the middle and
right panels of Fig. 5. We find that $Y_{\pi}(x_{1})$ is already within a
factor of two of its asymptotic value at $x\gg x_{1}$. Another way to estimate
the DM relic abundance in the odd-numbered case is to integrate the both sides
of the approximate differential equation

$\frac{dY_{\pi}}{dx}\simeq-\frac{s}{Hx}\langle\sigma_{\pi X_{P}\rightarrow
2\pi}v\rangle\frac{Y^{\rm eq}_{X}}{(Y^{\rm eq}_{\pi})^{2}}Y_{\pi}^{3}\,,$ (74)

which actually leads to an analytical expression very similar to Eq. (73) for
$\kappa_{P}x_{1}\lesssim 1$.

Note that there also exists another distinct possibility for very large values
of $R^{\prime}(0)$, where the detailed balance at $x=x_{1}$ is rather
established between the formation of $X$, via $3\pi\to\pi X$, and the
depletion of $X$, via $\pi X_{p}\to 2\pi$. In this case, we, in turn, have

$Y_{X_{(P)}}^{\rm eq}\left(\frac{Y_{\pi}}{Y_{\pi}^{\rm
eq}}\right)^{2}>Y_{X_{(P)}}>Y_{X_{(P)}}^{\rm
eq}\left(\frac{Y_{\pi}}{Y_{\pi}^{\rm eq}}\right)\,.$ (75)

This happens for the parameter choice of
$R(0)/m_{\pi}^{3/2}=R^{\prime}(0)/m_{\pi}^{5/2}=0.03$ (red dotted curve) in
the middle panel of Fig. 5, leading to $x_{1}=x_{2}$. Nevertheless, in reality
we expect $R(0)/m_{\pi}^{3/2}\gg R^{\prime}(0)/m_{\pi}^{5/2}$ for bound states
to typically hold. We leave a quantitative study of this possibility for
future work.

###  C.3 Prospective DM mass range

Finally, we explore the prospective higher mass DM range. The simple model
considered here is very restrictive as it only has a few principal parameters
such as $m_{\pi}$ and $m_{\pi}/f_{\pi}$. Larger masses $m_{\pi}$ require
larger interactions, i.e., larger ratios of $m_{\pi}/f_{\pi}$ to achieve the
observed relic abundance. One enters a regime where next-to-leading order and
unitarizing corrections become important. Here we simply explore the general
trends based on our leading-order formulation. In Fig. 6 we show examples for
$m_{\pi}=0.5\,\mathrm{GeV}$ (left panel), $m_{\pi}=1\,\mathrm{GeV}$ (middle
panel) and $m_{\pi}=2\,\mathrm{GeV}$ (right panel). Note that in all examples,
the perturbative limit, $m_{\pi}/f_{\pi}<4\pi$, and the condition for the
longevity of bound states, $|\Psi(0)|=R(0)/\sqrt{4\pi}<m_{\pi}^{3/2}$, are
satisfied.

Figure 6:  Evolution of DM abundance, $Y_{\pi}+2Y_{X}\simeq Y_{\pi}$ for
various DM masses with parameters chosen to achieve the observed relic
abundance. _Left_ : $m_{\pi}=8f_{\pi}=0.5$ GeV, for even case
$R(0)/m_{\pi}^{3/2}=1.1$, and for odd case
$(R(0)/m_{\pi}^{3/2},R^{\prime}(0)/m_{\pi}^{5/2})=(0.2,\,0.01)$. _Middle_ :
$m_{\pi}=10f_{\pi}=1$ GeV, even case $R(0)/m_{\pi}^{3/2}=1.2$, and odd case
$(R(0)/m_{\pi}^{3/2},R^{\prime}(0)/m_{\pi}^{5/2})=(0.2,\,0.01)$. _Middle_ :
$m_{\pi}=10f_{\pi}=2$ GeV, even case $R(0)/m_{\pi}^{3/2}=2.6$, and odd case
$(R(0)/m_{\pi}^{3/2},R^{\prime}(0)/m_{\pi}^{5/2})=(0.2,\,0.03)$. We fix the
binding energies to $\kappa_{S}=4\kappa_{P}=0.1$. The leading-order DM self-
scattering cross section per mass $\sigma_{\rm el}^{\rm LO}/m_{\pi}=0.08$,
$0.03$, $0.003$ bn/GeV, respectively, from left to right. Other curves are the
same as in Fig. 5.