Pion dark matter in a $\theta$ vacuum: a thermal relic with sharp velocity-dependent self-interactions

TL;DR

This paper proposes a dark matter model involving a non-zero topological angle and a dark photon portal, which explains the observed relic abundance and predicts velocity-dependent self-interactions relevant for astrophysical phenomena.

Contribution

It introduces a novel dark matter scenario with a sharp velocity-dependent self-interaction mechanism via a Breit-Wigner resonance, connecting QCD-like theories with astrophysical observations.

Findings

Accounts for observed dark matter relic abundance

Predicts velocity-dependent self-interactions in halos

Explains strong-lensing system observations

Abstract

As recently proposed, a non-vanishing topological angle may play a central role in QCD-like theories of dark matter (DM). In this work, we introduce a dark photon portal to the Standard Model in order to establish thermal equilibrium in the early Universe, and discuss the ensuing phenomenological constraints, including the stability of DM. The resulting dynamics accounts for the observed DM relic abundance and yields velocity-dependent DM self-interactions in astrophysical halos. Due to the sharp velocity dependence arising from a Breit-Wigner resonance, dedicated studies are required to assess the gravothermal evolution in detail, especially in the core-collapse regime. This is particularly timely in light of self-interacting DM interpretations of strong-lensing systems such as SDSS J0946+1006, which can be naturally explained within our framework.

Tables

Table 1. Table 1: Charge assignments and associated properties for benchmark models. The two-step thermalization between the DM and the SM bath is achieved by means of scatterings with the meson P ± P^{\pm} . See text for details.

Charge Assignment	P±	DM Lifetime
$Q_d^{\left(1\right)}=\text{diag}(+1,+1,-1)$	$K^{\pm }$	${\tau }_{\mathrm{DM}}$ in Eq. (42)
$Q_d^{\left(2\right)}=\text{diag}(+1,-1,-1)$	${\pi }^{\pm }$	${\tau }_{\mathrm{DM}}$ in Eq. (42) times ${\sin }^2{\theta }_{\pi \eta }$

Equations

$$\mathcal{L}=-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }+\frac{g^2\theta }{32{\pi }^2}{F}_{\mu \nu }{\tilde{F}}_{\mu \nu }+\bar{q}i/Dq-\left({\bar{q}}_{L}M{q}_R+\text{h.c.}\right),$$

$$q_L\to e^{-i\theta Q/2}{q}_L,$$

$$\theta \to \theta (1-\text{Tr}Q),$$

$$\text{Tr}Q=1.$$

$${\mathcal{L}}_{\mathrm{eff}}=\frac{f_{\pi }^2}{4}\text{Tr}[{\partial }_{\mu }{U}^{†}{\partial }^{\mu }U]+\frac{f_{\pi }^2}{2}{B}_0\text{Tr}[M_{\theta }^{†}U+U^{†}{M}_{\theta }],$$

$$\sin (\theta Q)=\alpha (\theta )M^{-1},$$

$$M_{\theta }=\cos (\theta Q)M+i\alpha (\theta )11,$$

$${\mathcal{L}}_{\mathrm{eff}}$$

$$\theta =\sum _i\arcsin \left(\frac{\alpha (\theta )}{m_i}\right),$$

$$\alpha (\theta )=m\sin \left(\theta /N_f\right),$$

$$\alpha (\theta )=\frac{\theta }{\mathrm{Tr}{M}^{-1}},$$

$$(M_{\text{meson}}^2{)}^{ab}=\frac{1}{2}{B}_0\alpha (\theta )\mathrm{Tr}\left(\cot (\theta Q)\{{\lambda }^a,{\lambda }^b\}\right),$$

$$(M_{\text{meson}}^2{)}^{ab}=\frac{1}{2}{B}_0\mathrm{Tr}\left(\sqrt{M^2-\alpha (\theta {)}^2}\{{\lambda }^a,{\lambda }^b\}\right).$$

$$m_i\to m_i^{\mathrm{eff}}\equiv \sqrt{m_i^2-\alpha (\theta {)}^2}.$$

$${\mathcal{L}}_{\text{eff, odd terms}}=-\frac{B_0\alpha (\theta )}{3f_{\pi }}\left(d_{abc}{\pi }_a{\pi }_b{\pi }_c-\frac{c_{abcde}}{10f_{\pi }^2}{\pi }_a{\pi }_b{\pi }_c{\pi }_d{\pi }_e\right)+\mathcal{O}(f_{\pi }^{-5}).$$

$${\mathcal{S}}_{\mathrm{WZW}}$$

$$\{\begin{array}{l}{m}_{{\pi }^{\pm }}^2=B_0(m_1+m_2),\\ m_{K^{\pm }}^2=B_0(m_1+m_3),\\ m_{K,{\bar{K}}^0}^2=B_0(m_2+m_3),\end{array}$$

$$\bigg{(}\begin{array}[]{cc}B_{0}(m_{1}+m_{2})&\frac{B_{0}}{\sqrt{3}}(m_{1}-m_{2})\\ \frac{B_{0}}{\sqrt{3}}(m_{1}-m_{2})&\frac{B_{0}}{3}(m_{1}+m_{2}+4m_{3})\end{array}\bigg{)}\,.$$

$$\bigg{(}\begin{array}[]{c}\pi^{0}\\ \eta\end{array}\bigg{)}=\bigg{(}\begin{array}[]{cc}\cos{\theta_{\eta\pi}}&\sin{\theta_{\eta\pi}}\\ -\sin{\theta_{\eta\pi}}&\cos{\theta_{\eta\pi}}\end{array}\bigg{)}\bigg{(}\begin{array}[]{c}\pi_{3}\\ \pi_{8}\end{array}\bigg{)}\,,$$

$$\tan (2{\theta }_{\eta \pi })=\frac{\sqrt{3}(m_2-m_1)}{2m_3-m_1-m_2}.$$

$$v_R=2\sqrt{\frac{m_{\eta }-2m_{{\pi }^0}}{m_{{\pi }^0}}},$$

$$\Gamma (\eta \to {\pi }^0{\pi }^0)=\frac{\alpha (\theta {)}^2{B}_0^2{\cos }^2(3{\theta }_{\eta \pi })}{24\pi f_{\pi }^2{m}_{\eta }}\sqrt{1-\frac{4m_{{\pi }^0}^2}{m_{\eta }^2}}.$$

$$Y=\frac{n}{s(T)},$$

$$f=\exp \left(\frac{\mu -E}{T}\right).$$

$$n=\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}f=\exp \left(\frac{\mu }{T}\right)n_{\mathrm{eq}}$$

$$\dot{n}+3Hn$$

$$\dot{n}+3Hn=sHx\frac{dY}{dx},$$

$$sHx\frac{dY}{dx}$$

$$\gamma \left(i\to f\right)=\int \left(\prod _{a\in i}\frac{d^3{\mathbf{p}}_a}{(2\pi {)}^{3}2E_a}{e}^{-\frac{E_a}{T}}\right)\left(\prod _{b\in f}\frac{d^3{\mathbf{p}}_b}{(2\pi {)}^{3}2E_b}\right){\left|{\mathcal{M}}_{i\to f}\right|}^2(2\pi {)}^4{\delta }^4\left(p_i-p_f\right).$$

$$\{\begin{array}{l}\gamma (1\to 23)=n_{1,\mathrm{eq}}\frac{K_1(m_1/T)}{K_2(m_1/T)}\Gamma (1\to 23),\\ \gamma (12\to 34)=n_{1,\mathrm{eq}}{n}_{2,\mathrm{eq}}⟨{\sigma }_{12\to 34}v⟩.\end{array}$$

Full text

1 Introduction

Understanding dark matter (DM) remains one of the most pressing open problems in modern science. DM in the form of dark pions Hochberg:2014kqa constitutes a compelling solution given its similarity with the QCD sector of the Standard Model (SM). In a recent work Garcia-Cely:2024ivo , it was pointed out that the topological theta angle of the dark sector can substantially modify the dynamics of DM pions. Concretely, Ref. Garcia-Cely:2024ivo , hereafter referred to as Paper I, demonstrated that a non-vanishing theta angle induces number-changing processes among the DM pions. As a result, their freeze-out can thermally give rise to the observed DM relic density, particularly in the presence of resonant effects. See also Ref. Kamada:2017tsq , which identified the number-changing processes associated with the $\theta$ angle, but reproduces the relic density with a mechanism independent of $\theta$, as well as Refs. Hochberg:2015vrg ; Kuflik:2015isi ; Bernal:2015bla ; Bernal:2015xba ; Bernal:2015ova ; Choi:2015bya ; Choi:2016hid ; Soni:2016gzf ; Kamada:2016ois ; Bernal:2017mqb ; Cline:2017tka ; Choi:2017mkk ; Kuflik:2017iqs ; Heikinheimo:2018esa ; Choi:2018iit ; Hochberg:2018rjs ; Bernal:2019uqr ; Choi:2019zeb ; Katz:2020ywn ; Smirnov:2020zwf ; Xing:2021pkb ; Braat:2023fhn ; Bernreuther:2023kcg ; Garani:2021zrr ; Dey:2016qgf ; Frumkin:2025iit for related work on MeV DM, and Refs. Redi:2016kip ; Draper:2018tmh ; Abe:2024mwa for the effect of the $\theta$ angle in QCD-like DM models in other contexts.

In this paper, we revisit the framework proposed in Ref. Garcia-Cely:2024ivo and address an aspect that remained unexplored, namely, the interactions between DM and ordinary particles of the SM. A thorough investigation is crucial for the following reasons: on the one hand, such interactions are essential for understanding the thermal processes leading to the relic abundance, which inherently assume that the DM sector thermalizes with the SM sector. On the other hand, with the portal interaction specified, the number-changing processes in the DM sector can lead to indirect detection signatures of DM, potentially including its decay. One must therefore ensure the DM meta-stability.

In contrast to the previous work Garcia-Cely:2024ivo , this analysis is necessarily model-dependent, as multiple realizations are possible Hochberg:2015vrg ; Hochberg:2018rjs ; Kamada:2017tsq ; Katz:2020ywn depending on the portal that communicates DM and the SM fields. To carry out our study, and for the sake of demonstrating that there exists at least one model that satisfies all current constraints, we focus on the benchmark model BM1 introduced in the Paper I. We extend this model by connecting it to the SM via a massive dark photon Holdom:1985ag , endowing the dark sector with an additional $U(1)_{d}$ gauge interaction.

We examine the Boltzmann equations relevant to determining the DM relic density and justify several assumptions made in the Paper I. In fact, we will show that, although correct, such assumptions were overly restrictive and, to a certain extent, unnecessary. This refinement has implications for the allowed theta-angle parameter space and may impact predictions for DM halo structures. For instance, we find a new family of solutions reproducing the observed relic density.

Moreover, the theta angle can also lead to modifications in the structure of small-scale DM halos today Spergel:1999mh ; Tulin:2017ara , while still agreeing with the predictions of collisionless dark matter (CDM) on larger scales Harvey:2015hha ; Bondarenko:2017rfu ; Harvey:2018uwf ; Sagunski:2020spe ; DES:2023bzs . Thus, the scenario naturally realizes the self-interacting dark matter (SIDM) paradigm Spergel:1999mh . When both effects are taken into account, one obtains a viable model of SIDM wherein the relic density is governed by the same physics responsible for altering DM halo structures. This result is highly non-trivial, as it is often believed that strongly coupled models of DM, while potentially predicting DM self-interactions, cannot simultaneously account for velocity-dependent scattering and therefore cannot capture the size-dependent effects in halo dynamics. Comprehensive reviews can be found in Refs. Tulin:2017ara ; Adhikari:2022sbh .

In particular, we now include a discussion of systems such as SDSS J0946+1006 Vegetti:2009cz , whose properties suggest an SIDM halo undergoing gravothermal collapse, see e.g Minor:2020hic , and which can be naturally explained within our framework through its sharp velocity dependence with enhanced cross sections at intermediate scales and suppressed interactions in clusters.

The structure of this paper is as follows. We begin by reviewing QCD-like theories of DM in Section 2. In Section 3, we discuss the Boltzmann equations in full generality, maintaining a model-independent perspective. The portal to the Standard Model is introduced in Section 4, where we examine thermalization, DM stability, and indirect detection constraints. In Section 5, we discuss the implications in DM halos. Finally, we present our conclusions in Section 6.

2 Pion DM in a $\theta$ vacuum

2.1 Pion DM for unitary gauge groups

We consider QCD-like theories with gauge group $SU(N_{c})$ and $N_{f}$ flavors of light Dirac quarks $q$ in the fundamental representation. Our results can be easily generalized to other gauge group choices such as $SO(N_{c})$ and $Sp(N_{c})$. Then, the Lagrangian of the theory is given by

[TABLE]

where $M=\text{diag}(m_{1},\ldots,m_{N_{f}})$ is the quark mass matrix with $m_{q}\neq 0$ and real for all flavors. Strong interactions confine at some energy scale $\Lambda$ and generate a fermion condensate, $\langle\bar{q}q\rangle\sim\Lambda^{3}$ Bai:2013xga , in analogy with standard QCD. If $m_{q}\ll\Lambda$, this Lagrangian exhibits an approximate $SU(N_{f})_{L}\times SU(N_{f})_{R}$ symmetry, corresponding to independent chiral transformations of the left- and right-handed fermion components in flavor space. In particular, an axial transformation takes the form

[TABLE]

where $Q$ is an arbitrary $N_{f}\times N_{f}$ hermitian matrix. Being anomalous under the $SU(N_{c})$ gauge interactions, such a transformation induces not only a change in the mass matrix but also a shift in the $\theta$ angle. Concretely,

[TABLE]

Hence, one can effectively move the $\theta$ parameter from the $F\tilde{F}$ term into the mass matrix by imposing

[TABLE]

In particular, one may take $Q$ diagonal such that $Q=\text{diag}(q_{1},..,q_{N_{f}})$.

It is believed that the fermion condensate spontaneously breaks the chiral symmetry to its diagonal subgroup, $SU(N_{f})_{L}\times SU(N_{f})_{R}\rightarrow SU(N_{f})_{V}$. This leads to the emergence of $N_{f}^{2}-1$ pseudo-Goldstone bosons, $\pi^{a}$, whose dynamics at sufficiently low energies is described in chiral perturbation theory Pich:1991fq ; Scherer:2002tk by

[TABLE]

We normalize the generators of $SU(N_{f})$ such that $\operatorname{tr}[\lambda^{a}\lambda^{b}]=2\delta^{ab}$, and parametrize the fermion condensate in terms of the dark meson constant, $f_{\pi}$, as $\langle\bar{q}q\rangle=-B_{0}f_{\pi}^{2}$. Note that it is expected that the confinement scale is of the order $\Lambda\sim 4\pi f_{\pi}/\sqrt{N_{c}}$ Kamada:2022zwb .

Let us emphasize a few points. The broken generators of the axial $SU(N_{f})_{A}$ subgroup, corresponding to the chiral transformations, $q_{L,R}\to\exp[\mp i\alpha^{a}\lambda^{a}]q_{L,R}$, are in a one-to-one correspondence with the light pseudo-scalar mesons $\pi_{a}$. Also, the $SU(N_{f})_{V}$ symmetry is an exact symmetry only in the limit of degenerate quarks; in general, it is explicitly broken by the differing values of the quark masses and is only approximately conserved.

Cancellation of tadpoles.

The vacuum is unstable unless the linear terms in $\pi^{a}$ are absent from $\mathcal{L}_{\rm eff}$. These terms cancel when $\text{Tr}[M_{\theta}\lambda^{a}-\lambda^{a}M_{\theta}^{\dagger}]=0$, or equivalently, $\text{Tr}[\sin(\theta Q)M\lambda^{a}]=0$ for every $a$. Since the matrices $\lambda^{a}$ generate the space of traceless $N_{f}\times N_{f}$ hermitian matrices, $\sin(\theta Q)M$ must be proportional to the identity. This happens if there exists a function $\alpha(\theta)$ such that

[TABLE]

In particular, one can cast $M_{\theta}$ in Eq. (3) as

[TABLE]

Likewise, the $M$-dependent part of the effective Lagrangian is

[TABLE]

The function $\alpha(\theta)$ can be found imposing the condition in Eq. (4). One implicit solution111There are more solutions corresponding to the different branches of the $\arcsin$. is

[TABLE]

Explicit solutions can be found in two important cases:

• •

For a quark degenerate mass spectrum, $m_{i}=m$, for which

[TABLE]

This regime was studied in Ref. Kamada:2017tsq . Since all meson masses are identical in this case, several effects are not captured; in particular, no resonant effects can arise.

• •

For $\theta\ll 1$, according to Eq. (6), $\alpha(\theta)\ll 1$ and hence

[TABLE]

This regime for $\theta$ was assumed in the Paper I; however, it was merely a simplifying assumption. As we shall see, allowing for arbitrary $\theta$ does not lead to significant differences.

Note that the $N_{f}+1$ parameters $\theta$ and $m_{i}$ can be exchanged for $\theta$, $\alpha$ and $q_{i}$ ($i=1,\ldots,N_{f}-1$) taking them as independent.

Meson mass spectrum.

According to Eq. (2.1), the meson masses follow from

[TABLE]

with $a,b=1,\dots,N_{f}$. In the particular cases for which $|q_{i}\theta|<\pi/2$,

Eq. (6) allows to write $\cos(q_{i}\theta)=\sqrt{1-(\alpha(\theta)/m_{i})^{2}}$. This leads to $\alpha(\theta)\cot(q_{i}\theta)=m_{i}\sqrt{1-(\alpha(\theta)/m_{i})^{2}}$, and hence

[TABLE]

Since the function $\alpha(\theta)$ is a number, in such cases, the sole effect of the $\theta$ angle on the meson masses with respect to a theory with $\theta=0$ is to replace the quark masses by

[TABLE]

In particular, assuming $m_{1}<m_{2}<\ldots<m_{N_{f}}$, we have $m^{\rm eff}_{1}<m^{\rm eff}_{2}<\ldots<m^{\rm eff}_{N_{f}}$.

Number-changing interactions.

While the first term in Eq. (2.1) describes meson masses and interactions involving an even number of mesons, the second term leads to interactions with an odd number of mesons. The leading two terms in the chiral expansion of the latter give

[TABLE]

With $d_{abc}=\frac{1}{4}{\rm Tr}(\{\lambda_{a},\lambda_{b}\}\lambda_{c})$, and $c_{abcde}=(\delta_{ab}d_{cde}+\delta_{cd}d_{abe})/N_{f}+\frac{1}{2}d_{abm}d_{cdn}d_{mne}$, which are non-zero only for $N_{f}\geq 3$. As already mentioned, in Paper I it was assumed that $\theta\ll 1$, which justifies using the expression for $\alpha(\theta)$ given in Eq. (11). However, formally speaking, the effective Lagrangian with the odd terms in $\pi^{a}$ remains the same for arbitrary $\theta$. Consequently, the results presented in that work concerning number-changing processes remain largely valid even for large values of $\theta$.

Finally, Eq. (15) must be compared against the Wess-Zumino-Witten (WZW) term Wess:1971yu ; Witten:1983tw , an interaction of topological origin that reads

[TABLE]

As Eq. (15), it gives rise to 3-to-2 meson annihilations.

2.2 Our benchmark model

The general properties outlined above were noted in Paper I and illustrated using two benchmark models. In this paper, we focus on BM1. Let us briefly summarize this model. It is based on a unitary gauge group with $N_{f}=3$ and quark mass matrix $M=\text{diag}(m_{1},m_{2},m_{3})$. The resulting meson spectrum resembles that of the SM, featuring mesons $\pi^{0}$, $\pi^{\pm}$, $K^{0}$, $\bar{K}^{0}$, $K^{\pm}$, and $\eta$.222Throughout this work, we will refer to the dark mesons with this notation; unless specified otherwise, SM mesons will not be referred to in this way. In particular, the superscript does not refer to any charge. The DM candidate is identified with the lightest meson, $\pi^{0}$, with the rest of the mesons eventually decaying or annihilating into DM.

Mass spectrum for $\theta=0$.

In this case,

[TABLE]

while $\pi_{3}$ and $\pi_{8}$ mix with their masses squared given by the eigenvalues of the matrix

[TABLE]

The physical states $\pi^{0}$ and $\eta$ are given by

[TABLE]

with a mixing angle

[TABLE]

Assuming $0<m_{1}\leq m_{2}\leq m_{3}$, the resulting mass spectrum exhibits a well-defined hierarchy: the $\eta$ meson is the heaviest state, while the $\pi^{0}$ is the lightest. The $\pi^{\pm}$ mesons are the second lightest states, with the mass splitting $m_{\pi^{\pm}}-m_{\pi^{0}}$ vanishing in the limit $m_{1}=m_{2}$. This is illustrated in the left panel of Fig. 1 (solid line) for which we assume $m_{\eta}=2m_{\pi^{0}}$ when $\theta=0$. Note that then $v_{R}|_{\theta=0}=0$ where

[TABLE]

which, as we shall see, characterizes the resonant scatterings of DM Chu:2018fzy . More precisely, the process $\pi^{0}\pi^{0}\to\pi^{0}\pi^{0}$ receives a resonant enhancement from the exchange of $\eta$ meson for non-relativistic velocities matching $v_{R}$.

As in Paper I, we remain agnostic of the origin of values of $v_{R}$ of the order of non-relativistic velocities in natural units, which arises from assuming $v_{R}|_{\theta=0}=0$ above. For examples of model-building in QCD-like models giving rise to this resonant behavior, see e.g., Refs. Tsai:2020vpi ; Csaki:2022xmu (see also Ref. Lee:2025lko ).

The case of a non-vanishing $\theta$ angle.

As explained above, for $|q_{i}\theta|<\pi/2$, the meson spectrum can be readily computed from the formulae of $\theta=0$ using Eq. (14) . Moreover, the four parameters $m_{1}$, $m_{2}$, $m_{3}$, and $\theta$ can be expressed as functions of the following four quantities: $\theta$, $\alpha$, the velocity $v_{R}|_{\theta=0}$, and the mass ratio $r_{12}=m_{1}/m_{2}$. Notably, since $\alpha$ is the only parameter out of the latter with mass dimension, any dimensionless combination of the meson masses may depend only on $\theta$, $r_{12}$, $v_{R}|_{\theta=0}$. Motivated by the resonant effect associated with $\eta$, throughout we will focus on $v_{R}|_{\theta=0}=0$. In Fig. 1 (left panel), we show the corresponding mass spectrum of the mesons taking into account the contribution of the $\theta$ angle.

Taking $v_{R}|_{\theta=0}=0$, all dimensionless combinations of masses will be functions of $\theta$ and $r_{12}$. Three dimensionless combinations of particular importance are the resonance velocity $v_{R}$, the mass splitting between $\pi^{\pm}$ and $\pi^{0}$ parametrized by $\delta\equiv m_{\pi^{\pm}}/m_{\pi^{0}}-1$, as well as the mixing angle $\theta_{\pi\eta}$. They are respectively shown in the Figs. 1 (right panel), 2 (left panel) and 2 (right panel). In the white regions of Fig. 2 the assumption $|q_{i}\theta|<\pi/2$ does not apply.

The role of $\eta$ for number-changing processes and resonant scattering.

The trilinear interaction in Eq. (15) –only present for a non-vanishing $\theta$ angle– induces the decay of $\eta$ meson into two pions as well as resonant elastic scatterings between pions, $\pi^{0}\pi^{0}\to\eta\to\pi^{0}\pi^{0}$. For these processes, the rate for $\eta\rightarrow\pi^{0}\,\pi^{0}$ enters in an important way, and is given by

[TABLE]

In addition to this decay and its inverse, the interactions in Eq. (15) induce $3\leftrightarrow 2$ self-annihilation processes. As such, the $\theta$ angle plays a central role in enabling efficient number-changing dynamics within BM1.

3 The relic density in the presence of the $\theta$ vacuum

In our benchmark scenario, as will be explained in section 5, observations of galaxy clusters impose an upper bound on the DM self-scattering cross section of approximately $0.5~\text{cm}^{2}/\text{g}$ Harvey:2015hha ; Bondarenko:2017rfu ; Harvey:2018uwf ; Sagunski:2020spe ; DES:2023bzs . In the region where chiral perturbation theory is valid, namely $m_{\pi^{0}}/f_{\pi}<4\pi/\sqrt{N_{c}}$, and where this bound is satisfied, the maximum strength allowed for the $3\to 2$ processes does not yield sufficient DM annihilation if $v_{R}\gg 1$, resulting in a relic abundance significantly larger than the observed DM density. This issue becomes even more severe in the limit $\theta\to 0$, where the WZW term is the sole source of number-changing interactions, as the resulting $3\leftrightarrow 2$ processes are velocity suppressed and the scattering in halos is velocity independent. Hochberg:2014kqa ; Kamada:2022zwb . In contrast, in the regime $v_{R}\lesssim 1$, the exchange of the $\eta$ meson resonantly enhances the 3-to-2 processes, reconciling the predicted relic density with the observed value. Concretely, the $\eta$ effectively acts as a catalyzer Chu:2024rrv , mediating number-changing interactions such as $\pi^{0}\,\pi^{0}\rightarrow\eta$ followed by $\eta\,\pi^{0}\rightarrow\pi^{0}\,\pi^{0}$, effectively inducing a resonant 3-to-2 process.

This section examines how via these number-changing processes a non-zero $\theta$ angle impacts the relic abundances of the mesons. Before discussing the corresponding Boltzmann equations, let us define notation. As usual, $H$ is the Hubble rate, which, in a radiation-dominated Universe, is given by $H\simeq\sqrt{4\pi^{3}g_{*}/45}\,T^{2}/M_{\rm Pl},$ with $g_{*}$ being the energy density degrees of freedom. Likewise, we find it convenient to normalize the number density of mesons to the entropy density of the Universe

[TABLE]

where $s(T)=2\pi^{2}g_{*s}\,T^{3}/45$, and $g_{*s}$ is the effective number of degrees of freedom in entropy.

3.1 The evolution with $T_{\rm dark}=T$

We start by assuming that the dark sector and the SM are in kinetic equilibrium at the same temperature. Later on, we will show that deviations from this assumption have a minor effect on the relic density and do not significantly alter the results for benchmark points of interest. Neglecting degeneracy effects, $T_{\rm dark}=T$ means that the dark mesons are distributed in phase space according to

[TABLE]

The corresponding densities are given by

[TABLE]

Hence, the number densities are determined by the chemical potentials, and vice versa, at a given SM temperature. For simplicity, let us first consider a given meson, which could be the DM or a heavier state. Quite generically, the Boltzmann equation determining its number density, when one particle is destroyed in an individual reaction $i\to f$, is given by

[TABLE]

where the second line accounts for the inverse reaction. In terms of the yield in Eq. (23)

[TABLE]

By substituting Eq. (24) into the Boltzmann equation, and using the principle of detailed balance –which relates the squared amplitudes of a reaction and its inverse– we obtain

[TABLE]

Here, the interaction rate density is

[TABLE]

This quantity is related to the ordinary decay rates and thermally averaged cross sections by

[TABLE]

For two identical particles in the initial state, an additional $1/2$ must be included. The analogue factor for identical particles in the final state is included in the definition of $\Gamma$ and $\langle\sigma v\rangle$. Adding several reactions and particles is straightforward. For simplicity, throughout we use the following notation

[TABLE]

We are now in a position to write down the Boltzmann equations that determine the DM density, including the effects of coannihilations. Concretely, the evolution of $\pi_{i}$ and $\eta$ populations is described by the following coupled Boltzmann equations

[TABLE]

We neglect the effect of kaons because, although they are lighter than the $\eta$, there are no resonant processes involving them that would impact the final DM abundance. The remaining contributions are limited to conversion processes, which can be safely neglected since $m_{K}\gtrsim 1.7\,m_{\pi^{0}}$. Furthermore, we also neglect the $3\leftrightarrow 2$ self-annihilations originated from both the WZW term and the 5-point interactions induced by $\theta$, as well as further number-changing processes, such as 4-to-2 annihilations. This stems from the fact that such processes are subdominant with respect to the resonant $3\to 2$ interactions.

The thermally-averaged cross sections entering in the rates are the following. For the semi-annihilations $\eta\,\pi_{a}\to\pi_{b}\,\pi_{c}$ ($a,b,c=0,1,2$)

[TABLE]

where the nonzero values for $\beta_{abc}$ are

[TABLE]

For the (co-)annihilations

[TABLE]

In Fig. 3, we present the results of our numerical analysis for a representative parameter choice: $(m_{\pi^{0}},f_{\pi})=(20~\text{MeV},34~\text{MeV})$. We numerically solve the Boltzmann equations (3.1)–(3.1), varying $\theta$ and $r_{12}$ within the ranges discussed above with the relic density fixed to the observed value. Several comments are in order.

First, we can recognize two distinct families of solutions: those where the DM abundance depends only on the value of $r_{12}$ and is independent of $\theta$ (on the vertical line), and those where the abundance is sensitive to both parameters (on the diagonal line).

• •

The former, studied in Paper I, are characterized by the fact that the processes $\eta\leftrightarrow\pi^{0}\pi^{0}$ are much faster than the Universe expansion when $\eta\pi^{0}\to\pi^{0}\pi^{0}$ freezes out, so that they keep chemical equilibrium among $\eta$ and $\pi^{0}$, which enforces $\mu_{\eta}=2\mu_{\pi}$ after freeze-out. See also Refs. Frumkin:2021zng ; Garcia-Cely:2024ivo . We illustrate this in the left panel of Fig. 4 for a representative choice of parameters. Concretely, we show the evolution of the abundances of $\pi^{0},\pi^{\pm},\eta$ (upper left panel) and their chemical potentials (lower left panel). Notice that $\mu_{\eta}=\mu_{\pi^{0}}=0$ when semi-annihilation and (inverse) decays are simultaneously active, while $\mu_{\eta}=2\mu_{\pi^{0}}$ after the decoupling of semi-annihilations. These relations between chemical potentials allow us to write a simplified Boltzmann equation for the combination $Y=Y_{\pi}+2Y_{\eta}$:

[TABLE]

where $Y_{\pi}=Y_{\pi^{0}}+2Y_{\pi^{\pm}}$ and $\langle\sigma_{\eta\pi}v\rangle$ is the proper combination of cross sections. Notice that, since the mass splitting between $\pi^{0}$ and $\pi^{\pm}$ is less than $10\%$, the latter contribute to the DM relic abundance, converting as $\pi^{+}\pi^{-}\to\pi^{0}\pi^{0}$ some time after freeze-out. Most importantly, this simplified equation is independent of $\theta$, in agreement with our results in Fig. 3. As discussed in Paper I, the simplified equation can be solved numerically, or alternatively, an approximate analytical solution of the form $Y_{\pi,\rm eq}(x_{\rm fo})$ can be obtained solving for $x_{\rm fo}$ from $H(x_{\rm fo})\simeq n_{\eta,\rm eq}(x_{\rm fo})\langle\sigma_{\eta\pi}v\rangle$. Both of these approaches are in excellent agreement with the numerical solution of the full set of Boltzmann equations, thus corroborating the results of Paper I.

• •

On the other hand, the points on the diagonal lines of Fig. 3 correspond to a new family of solutions, not discussed in Paper I. In such points, the inverse decay processes decouple before semi-annihilations. In the right panel of Fig. 4, we show the abundances (upper right panel) and the chemical potentials (lower right panel) of the different species for a representative choice of parameters. Once again, we observe that all the chemical potentials vanish when both semi-annihilations and inverse decays are in equilibrium. When the inverse decays become inefficient, the chemical equilibrium between $\eta$ and $\pi$ is kept by the semi-annihilation, leading to $\mu_{\eta}=\mu_{\pi^{0}}$. This lasts until semi-annihilation freezes out. After that, both processes involved are inactive and their rates are comparable in magnitude, such that the chemical equilibrium is not achieved. Eventually, the inverse decay becomes active again, restoring the chemical equilibrium and fixing the chemical potentials to $\mu_{\eta}=2\mu_{\pi}$. A similar behavior of the chemical potentials has been identified in other DM scenarios Garny:2017rxs ; DAgnolo:2017dbv ; Puetter:2022ucx .

In this case, a simplified Boltzmann equation, such as Eq. (38), fails to provide results in quantitative agreement with the full numerical solution. Nevertheless, we signal that a simplified equation involving only the (inverse) decay processes and enforcing $\mu_{\pi}=\mu_{\eta}$ provides results in qualitative agreement with the full numerical solution.

3.2 Departure from kinetic equilibrium

Before discussing the thermalization details—and thus the SM portal—let us emphasize that the results so far assume $T=T_{\text{dark}}.$ On the one hand, this requires strong interactions for the dark mesons so that the dark sector develops a thermal distribution. On the other hand, this also requires that the dark sector and the SM are kept in kinetic equilibrium, occurring if the elastic scatterings among the dark mesons and SM particles are sufficiently strong to allow an efficient energy transfer from one sector to the other.

If kinetic equilibrium between the two sectors is never realized, the dark sector develops its own temperature $T_{\rm dark}$ in view of elastic scattering among dark sector particles, which evolves independently of the SM temperature $T$. We consider the more predictive scenario in which the interactions between the SM and the dark sectors are strong enough to keep kinetic equilibrium among them initially, thus enforcing $T_{\rm dark}=T$ until the Universe cools down to some temperature, $T_{\text{dec}}$. At that time, known as kinetic decoupling, elastic scatterings of the dark mesons off SM particles become slower than the expansion of the Universe and the dark sector evolves with its own temperature. If number-changing processes within the dark sector are inefficient at the time of kinetic decoupling and afterwards, the temperature evolves as $T_{\rm dark}\sim 1/a$ if DM is relativistic or $T_{\rm dark}\sim 1/a^{2}$ if DM is non-relativistic Pappadopulo:2016pkp ; Katz:2020ywn . See Ref. Bernal:2015ova for a discussion in the context of SIDM. Notice that in this scenario, the evolution of the dark temperature is predicted in terms of that of the SM.

In the following, we denote the time of kinetic decoupling as $x_{\rm dec}=m_{\pi^{0}}/T_{\rm dec}$ and define the time of DM freeze-out as $x_{\rm fo}=m_{\pi^{0}}/T_{\rm fo}$. In this way, for $x<x_{\rm fo}$ the dark mesons have a vanishing chemical potential enforced by number-changing processes, while chemical decoupling occurs for $x>x_{\rm fo}$ with the dark mesons acquiring a chemical potential. For the family of solutions corresponding to the vertical line of Fig. 3, in which semi-annihilations decouple before inverse decays, $x_{\rm fo}$ corresponds to the usual definition of the time at which $\eta\pi^{0}\to\pi^{0}\pi^{0}$ becomes inefficient.

In the next section, after introducing an explicit SM portal and the interactions between the SM and the dark mesons, we will be interested in the region of the parameter space where kinetic decoupling occurs after DM freeze-out, $x_{\rm dec}>x_{\rm fo}$. Since $x_{\rm fo}\sim\mathcal{O}(20)$ and $\pi^{0}$ is the lightest dark particle, the dark sector is non-relativistic at the time of kinetic decoupling, so that for $x>x_{\rm dec}$ its temperature evolves as $T_{\rm dark}=T\,x_{\rm dec}/x\sim 1/a^{2}$. We will take this effect into account by solving the Boltzmann equations for the dark meson abundances in two steps:

• •

for $x<x_{\rm dec}$ we solve the set of coupled Boltzmann equations in Eq. (3.1)-(3.1);

• •

for $x>x_{\rm dec}$ we solve the same set of equations replacing the dark meson equilibrium distributions as $n_{a,{\rm eq}}(T)\to n_{a,{\rm eq}}(T_{\rm dark})$, so that $Y_{a,{\rm eq}}\to n_{a,{\rm eq}}(T_{\rm dark})/s(T)$, with $T_{\rm dark}=Tx_{\rm dec}/x$.

This typically leads to a reduction of the DM relic abundance of a factor of $\mathcal{O}(10-30)\%$.333In the opposite regime, i.e $x_{\rm dec}<x_{\rm fo}$, number-changing processes are still active after kinetic decoupling enforcing $\mu_{i}=0$ for all species and heating up the dark sector, whose temperature evolves logarithmically with the scale factor, $T_{\rm dark}\sim 1/\ln a$ Carlson:1992fn ; Pappadopulo:2016pkp .

4 Portal to the SM

4.1 Dark photon portal and DM decays

As mentioned above, we will assume a dark sector initially in thermal equilibrium with the SM thermal bath. This is a crucial assumption as, in the absence of an explicit portal, number changing processes within the dark sector would modify the evolution of the dark sector temperature and the corresponding DM population.

In this section, we substantiate this assumption by constructing explicit realizations of BM1 when it is coupled to the SM by means of a dark photon. This portal Holdom:1985ag has been studied in the context of models of strongly interacting massive particle (SIMP) DM where the relic abundance is obtained by 3-to-2 annihilations induced by the WZW term, see e.g. Refs. Hochberg:2015vrg and Katz:2020ywn . Although other portals can be considered –notably those induced by axion-like particles Hochberg:2018rjs ; Kamada:2017tsq , potentially related to the dark theta angle– we focus on dark photons to demonstrate the existence of one possible realization of QCD-like DM produced by the theta angle, as well as the internal consistency and phenomenological viability of the framework.

For simplicity, in this section we restrict our analysis to small values of the $\theta$ angle, namely $\theta\ll 1$. As discussed in the previous sections, we do not expect significant differences for larger values, $\theta\sim\mathcal{O}(1)$. We extend the model by introducing a $U(1)_{d}$ gauge symmetry, associated with a coupling constant $e_{d}$. The corresponding dark photon, with mass $m_{V}$ and kinetic mixing $\varepsilon$, is described by

[TABLE]

where $B_{\mu\nu}$ is associated with the SM hypercharge. We remain agnostic regarding the origin of the dark photon mass, which could arise either through the Stückelberg mechanism Stueckelberg:1938hvi or via interactions with a dark Higgs field Englert:1964et ; Higgs:1964pj .

On the one hand, as usual, the SM sector carries no dark charge and communicates with the DM only through the parameter $\varepsilon$, that induces a mixing between the dark photon and the SM photon. After diagonalization, the SM photon only couples to the SM electromagnetic current $J_{\mu}^{\rm SM}$ as $eA^{\mu}J_{\mu}^{\rm SM}$, while the physical dark photon couples to both the dark current $J_{\mu}^{\rm dark}$ and the SM current with coupling $V^{\mu}(e_{d}\,J_{\mu}^{\rm dark}+e\,\varepsilon\,c_{W}\,J_{\mu}^{\rm SM})$, where $c_{W}$ denotes the cosine of the weak mixing angle. Since we are primarily interested in the MeV scale, we ignore the small mixing between the dark photon and the $Z$ gauge boson, which is suppressed by $\varepsilon\,(m_{V}/m_{Z})^{2}$ Babu:1997st .

On the other hand, the dark quarks are charged under the $U(1)_{d}$ gauge group, with charges $Q_{d}={\rm diag}\,(\tilde{q}_{1},\tilde{q}_{2},\tilde{q}_{3})$. The chiral Lagrangian in Eq. (5) in presence of the dark photon is extended to  Scherer:2002tk ; Katz:2020ywn ; Braat:2023fhn

[TABLE]

Here, the covariant derivative is defined as $D_{\mu}U=\partial_{\mu}U-ie_{d}V_{\mu}\,[Q_{d},U]$. The last term of the first line is responsible for the mass corrections of the dark mesons charged under $U(1)_{d}$, with $c$ being a constant. Assuming a heavy dark photon, $m_{V}>\Lambda$, it scales as $c\sim(\Lambda/m_{V})^{2}\ln(m_{V}^{2}/\Lambda^{2})$ Katz:2020ywn .

In the presence of the portal, DM stability is no longer guaranteed, as the number-changing processes might induce its decay. Note that this does not hold if a conserved symmetry exists in the dark sector, but this is not the case for BM1, in which the flavor symmetry is broken by quark mass differences, $m_{1}\neq m_{2}\neq m_{3}$. Thus, the DM particle is meta-stable, and its interactions with SM particles determine whether it is sufficiently long-lived.

DM decay is induced by the operator $\pi^{0}\,V_{\mu\nu}\tilde{V}^{\mu\nu}$ and, in presence of a non-vanishing $\pi^{0}-\eta$ mixing, by $\eta\,V_{\mu\nu}\tilde{V}^{\mu\nu}$. For a generic choice of dark quark charges, both operators are generated via the axial anomaly: the approximate (spontaneously broken) $SU(3)_{A}$ global symmetry of the Lagrangian is anomalous under the dark $U(1)_{d}$ interactions.444More precisely, even in the limit in which the dark quarks are massless, a subset of the currents associated with the $SU(3)_{A}$ transformations are not conserved because of the quantum anomaly, namely $\partial^{\mu}J_{\mu,A}^{3,8}\propto e_{d}^{2}V_{\mu\nu}\tilde{V}^{\mu\nu}{\rm Tr}[\lambda^{3,8}Q_{d}^{2}]$ with $J_{\mu,A}^{3,8}=\bar{q}\,\gamma_{\mu}\gamma_{5}\lambda^{3,8}\,q$. Note that this closely resembles the anomalous decay of neutral mesons into two photons observed in the SM. The presence of the anomaly is encoded in the second line of Eq. (40), which provides interaction terms proportional to $(e_{d}^{2}/f_{\pi})\,\pi_{3,8}\,V_{\mu\nu}\tilde{V}^{\mu\nu}{\rm Tr}[\lambda^{3,8}Q_{d}^{2}]$. These, once expressed in terms of the mass eigenstates, induce the aforementioned operators with couplings suppressed by only one power of $f_{\pi}$.

However, it is possible to choose the charge matrix $Q_{d}$ in such a way that the anomaly of the axial currents vanishes, ${\rm Tr}[\lambda_{3,8}Q_{d}^{2}]=0$. In such a case, the corresponding currents are conserved in the massless limit, and the anomalous contribution to the operators discussed in the previous paragraph vanishes. This requirement can only be fulfilled by assigning the dark charges so that $\tilde{q}_{1}^{2}=\tilde{q}_{2}^{2}=\tilde{q}_{3}^{2}$, and hence $Q_{d}\propto(\pm 1,\pm 1,\pm 1)$. Without loss of generality, there are only three charge assignments (up to a global normalization). Two of them are listed in Table 1, while the third one, $Q_{d}^{(3)}\propto{\rm diag}(+1,-1,+1)$, is qualitatively similar to $Q_{d}^{(2)}$ and provides analogous results.

Even for a vanishing anomaly, the interaction terms $\pi^{0}V_{\mu\nu}\tilde{V}^{\mu\nu}$ and $\eta V_{\mu\nu}\tilde{V}^{\mu\nu}$ are generated by higher-order operators in the chiral Lagrangian expansion. More precisely, the relevant operators arise at $\mathcal{O}(p^{6})$, that is, next-to-next-to-leading order. Following the approach in Ref. Katz:2020ywn , as a representative operator we consider

[TABLE]

where $\alpha_{d}=e_{d}^{2}/4\pi$. This operator induces the $\pi^{0}V_{\mu\nu}\tilde{V}^{\mu\nu}$ and $\eta V_{\mu\nu}\tilde{V}^{\mu\nu}$ interactions with couplings suppressed by an extra factor of $(m_{\pi^{0}}/f_{\pi})^{2}$ compared to the anomalous contribution discussed before.555Notice that the presence of a non-vanishing $\theta$ angle may induce additional CP-violating operators in the chiral Lagrangian, potentially giving rise to $\theta\,\eta\,V_{\mu\nu}V^{\mu\nu}$ and $\theta\,\pi^{0}\,V_{\mu\nu}V^{\mu\nu}$ interactions. Their contribution to the DM decay rate is nevertheless suppressed for small values of $\theta$. One could expect $\mathcal{O}(1)$ corrections for large $\theta\sim\mathcal{O}(1)$. Thus, in the case of vanishing anomaly, the DM decay rate is significantly smaller close to the chiral limit, $m_{\pi^{0}}\ll f_{\pi}$, providing a much larger DM lifetime.

Finally, DM is absolutely stable in the limit in which the global $SU(2)_{V}$ symmetry involving the first two flavors of dark quarks is conserved, i.e. dark isospin. This requires both $r_{12}=0$ and $\tilde{q}_{1}=\tilde{q}_{2}$. For the cases in which $\tilde{q}_{1}=\tilde{q}_{2}$ but $r_{12}\neq 0$, the decay of $\pi^{0}$ is induced by dark isospin breaking, or, in other words, by the $\pi^{0}-\eta$ mixing, as given in Eq. (20). As shown in Fig. 2 (right panel), this is fixed by $r_{12}$ and $\theta$.

In the following, we assume $m_{V}>m_{\eta}$ so that the decays $\pi^{0},\eta\to VV$ are kinematically forbidden. As the dark mesons do not couple to the physical SM photon, the dark decays $\pi^{0},\eta\to\gamma\gamma$ are absent and DM decays to SM leptons. For $m_{\eta}\lesssim 200$ MeV, the relevant DM decay channel is $\pi^{0}\to e^{+}e^{-}e^{+}e^{-}$, corresponding to a tree level diagram with two dark photon propagators, see Fig. 5. This process involves the electron-dark photon coupling twice, so that both DM and $\eta$ decay rates are suppressed by $\varepsilon^{4}e_{d}^{4}$. An analogous decay mode arises for the $\eta$ meson. However, its branching ratio is extremely suppressed, as $\eta$ mostly decays to $2\pi^{0}$ via the $\theta$ angle. We will briefly comment on this below.

In the following, we will describe the possible charge assignments presented in Table 1 for which the axial anomaly vanishes and the DM is long-lived.

• •

Charge assignment 1. The dark mesons $\pi^{0}$, $\pi^{\pm}$ and $\eta$ are neutral under $U(1)_{d}$, whereas $K^{+}$ and $K^{0}$ carry a charge of $+2$, while $K^{-}$ and $\overline{K^{0}}$ carry a charge of $-2$. Being charged, the four kaons are stable. Furthermore, the resulting Lagrangian is invariant under the accidental global symmetry generated by $Q^{\prime}={\rm diag}\,(1,-1,0)$. Being the lightest particles charged under this symmetry, the $\pi^{\pm}$ are stable. In contrast, $\eta$ and $\pi^{0}$ are not protected by any symmetry, and therefore decay. DM decays are induced via the operator in Eq. (41), with a lifetime estimated up to $\mathcal{O}(1)$ factors as

[TABLE]

For reference, as shown in Fig. 2 (right panel), the range of values for the mixing angle is $0.23\gtrsim\sin\theta_{\pi\eta}\gtrsim 0.005$ for $\theta\ll 1$. The decay is strongly suppressed close to the chiral limit $m_{\pi^{0}}\ll f_{\pi}$, as expected. Furthermore, the DM is stable in the limit $\sin\theta_{\pi\eta}=0$ (equivalent to $r_{12}=0)$ in view of the isospin-preserving charge assignement, $\tilde{q}_{1}=\tilde{q}_{2}=1$. The corresponding constraints from DM indirect searches will be discussed in Section 4.3.

Moreover, the kaon masses receive a small correction due to the interaction with the dark photon (similarly to the SM). As long as $m_{V}>\Lambda$, this can be estimated as $m_{K}^{2}\to m_{K}^{2}+\Delta m_{K}^{2}$ with $\Delta m_{K}^{2}\sim 8\,c\,e_{d}^{2}\,f_{\pi}^{2}$. Recall, nevertheless, that the DM relic density is only marginally affected by this.

• •

Charge assignment 2. The charge of the dark mesons is analogous to the ordinary electric charge in the visible sector: the mesons $\pi^{0}$, $\eta$, $K^{0}$, and $\overline{K^{0}}$ are uncharged, whereas $\pi^{\pm}$ and $K^{\pm}$ carry charge $\pm 2$. While the $\pi^{\pm}$ states are stable as a result of this, $\pi^{0}$ decays with a lifetime estimated –up to factors of $\mathcal{O}(1)$– using Eq. (42) without the $(1/\sin^{2}\theta_{\pi\eta})$ factor. In this case, as $\tilde{q}_{1}\neq\tilde{q}_{2}$, the charge matrix explicitly breaks isospin, so that the operator in Eq. (41) directly induces the decay of $\pi^{0}$ regardless of its mixing with $\eta$. For the same parameters used in Eq. (42), the lifetime is shorter, namely, $\tau_{\rm DM}\sim 10^{27}$ s.

Along with this, the ${\pi^{\pm}}$ interaction with dark photon shift their mass squared: $m_{\pi^{\pm}}^{2}\to m_{\pi^{\pm}}^{2}+\Delta m_{\pi^{\pm}}^{2}$. We can express this in terms of the pion splitting defined in the previous sections as

[TABLE]

where \delta(r_{12})\big{|}_{e_{d}=0} is the correction in the absence of the portal, shown in Fig. 2 (left panel). The other term can be estimated for $m_{V}>\Lambda$ as Katz:2020ywn

[TABLE]

where we used $\Delta m_{\pi^{\pm}}^{2}\sim 8ce_{d}^{2}f_{\pi}^{2}$.

Notice that the additional source of mass splitting for the charged pions affects DM relic density when $\delta_{V}\gtrsim 0.01$. Indeed, the charged pions become sufficiently heavy that their number density is significantly Boltzmann-suppressed at freeze-out, $n_{\pi^{\pm}}\ll n_{\pi^{0}}$, thereby reducing the total DM abundance, $n_{\rm DM}=n_{\pi^{0}}+2n_{\pi^{\pm}}$. Such a reduction in the DM number density requires a larger value of $f_{\pi}$ to match the cosmological relic abundance for fixed DM mass and $r_{12}$ compared to the case of $\delta_{V}=0$. Furthermore, if $\delta\gtrsim 0.5$, the total cross section for $\eta\pi\to\pi\pi$ becomes smaller, as the $\eta\pi^{0}\to\pi^{+}\pi^{-}$ becomes kinematically closed.

We will take these effects into account in the following section while solving the Boltzmann equations for the dark meson abundances for a representative benchmark point.

4.2 Thermalization

In the rest of this section, we discuss the conditions for kinetic equilibrium between the SM and the dark sector. Kinetic equilibrium is established via elastic scatterings of the dark mesons with SM fermions mediated by the dark photon. Although this discussion shares some similarities with previous studies of other SIMP candidates (see e.g. Refs. Hochberg:2015vrg ; Katz:2020ywn ), our DM candidate is neutral under $U(1)_{d}$ and therefore does not interact directly with SM fermions. Thus, kinetic equilibrium is established in two steps:

• I.

The lightest dark mesons charged under $U(1)_{d}$ elastically scatter off SM leptons (mainly electrons and positrons) via the $t$-channel exchange of a dark photon. If these interactions are sufficiently strong, the charged mesons remain in kinetic equilibrium with the SM bath, thereby maintaining a common temperature. The lightest charged mesons are the kaons $K^{\pm}$ within the charge assignment 1 and the pions $\pi^{\pm}$ within charge assignment 2, see Table 1. The relevant elastic cross section is given by

[TABLE]

where $p_{f}$ is the momentum of the SM fermion $f=e^{\pm}$ and $Q_{K}=Q_{\pi}=2$ is the dark charge of the corresponding meson. They thermalize with the SM as long as Hochberg:2015vrg ; Katz:2020ywn

[TABLE]

where the scattering rate is given by

[TABLE]

• II.

The DM particle $\pi^{0}$ thermalizes with the lightest meson charged under $U(1)_{d}$, $P^{\pm}$, thus sharing their temperature, which is the same of the SM bath if the condition in Eq. (46) is satisfied. The relevant thermal cross sections are given by

[TABLE]

We numerically find that these expressions are of the same order of magnitude. Once established, the kinetic equilibrium is maintained as long as the elastic scatterings are more rapid than the Hubble expansion, that is,

[TABLE]

The main difference between the two charge assignments is that the charged pions are nearly degenerate in mass with the neutral pion, i.e., $m_{\pi^{\pm}}\sim m_{\pi^{0}}$, resulting in a significantly larger abundance of charged pions at freeze-out compared to that of kaons, $n_{\pi^{\pm}}\gg n_{K}$ at $T=T_{\rm fo}$. Consequently, the thermalization condition is more easily satisfied for the charge assignment 2. Thus, for simplicity in the numerical analysis that follows, we choose benchmark points within the charge assignment 2.

Finally, we compute the largest temperatures at which Eq. (46) or Eq. (50) no longer holds, which corresponds to the kinetic decoupling $T_{\rm dec}$ as defined in Section 3.2.For concreteness, we will focus on parameter space where kinetic decoupling occurs after DM freeze-out, that is $T_{\rm dec}<T_{\rm fo}$.

4.3 Viable parameter space

Taking charge assignment 2, Fig. 6 shows the region of the parameter space $(m_{V},\epsilon)$ where thermalization is achieved for a representative choice of the parameters $m_{\pi^{0}}$, $f_{\pi}$, $r_{12}$ and $\alpha_{d}$. In the purple region, denoted as early thermal decoupling, at least one of Eqs. (46) and (50) is violated before DM freeze-out and kinetic equilibrium is lost. The yellow and gray regions are excluded by indirect detection and collider searches to be discussed below. Here, the annihilation channel $\pi^{+}\pi^{-}\to e^{+}e^{-}$ is slower than the Universe expansion rate at $T=T_{\rm fo}$ and is irrelevant for the evolution of the DM density666The corresponding annihilation rate is given by

$\Gamma_{\rm ann}=\frac{\sum_{\pi}Q_{\pi}^{2}}{N^{2}_{\pi}}\langle\sigma v(\pi^{+}\,\pi^{-}\ \to e^{+}e^{-})\,n_{\pi}\rangle=\frac{2\,e^{2}e_{d}^{2}c_{W}^{2}m_{\pi}\varepsilon^{2}T}{\pi(m_{V}^{2}-4m_{\pi}^{2})^{2}}sY_{\rm DM}\,,$

(51)

and satisfies $\Gamma_{\rm ann}<H|_{T_{\rm fo}}$ for all the points displayed in Fig. 6.. Furthermore, as $m_{V}>2m_{\pi}$, the semi-annihilation channel $\pi\pi\to\pi V$ is kinematically closed.

The Boltzmann equations for the dark meson abundances are solved for a representative benchmark point, indicated by a star in Fig. 6, with the corresponding results shown in the left panel of Fig. 7. Although the behavior observed in the left panel of Fig. 7 is similar to what has been studied in the literature for standard SIMP scenarios Hochberg:2015vrg ; Katz:2020ywn ; Braat:2023fhn , our value of $m_{\pi^{0}}/f_{\pi}\sim 0.6$ lies well within the perturbative regime of chiral expansion. This opens up parameter space in agreement with indirect detection constraints. In the numerical solution of the Boltzmann equations, we take into account the kinetic decoupling discussed in the previous section as well as the correction to the charged pion masses, which gives $\delta_{V}\sim 0.05$. Such corrections have a significant impact on the DM relic density. In particular, the abundance of the charged pions (and therefore the total DM number density) at $x_{\rm fo}$ is reduced compared to the case without a portal. This effect, together with the kinetic decoupling, explains the fact that this benchmark point –which belongs to the vertical family of solutions discussed in Section 3– does not reproduce the DM relic abundance for the value of $r_{12}$ along the yellow vertical lines of Fig. 3, but for a slightly larger one. Nevertheless, this effect on the relic density is at most at the tens-of-percent level.

Constraints on DM decays.

Current bounds on MeV DM decaying to electrons require a lifetime $\tau_{\rm DM}\gtrsim 10^{26}\mathrm{\,s}$ Cirelli:2024ssz . More precisely, gas-rich dwarf galaxies with a low cooling rate are extremely sensitive to heat injections from DM decays. The strongest constraint in the mass range 1-15 MeV is set by observations of the Leo T dwarf galaxy Wadekar:2021qae . Furthermore, a DM candidate decaying to electrons would produce an imprint on charged cosmic rays. While the solar magnetic field prevents their local observation, interstellar cosmic rays can nonetheless be detected by Voyager Boudaud:2016mos , which sets the most stringent bound in the mass range 15-50 MeV. Finally, larger masses in the sub-GeV range can be probed by X-ray telescopes, which detect photons produced via Inverse Compton scattering of background radiation off electrons originating from DM decay, with the strongest constraints provided by XMM-Newton Koechler:2023ual .

Thus, for the choice of DM mass of 20 MeV adopted in Fig. 6, the strongest constraint is set by Voyager, $\tau_{\rm DM}\gtrsim 3\times 10^{26}\mathrm{\,s}$. This is shown in Fig. 6, as the yellow exclusion region. Notice that for the point indicated with the star, the DM lifetime is $\tau_{\rm DM}\simeq 5\times 10^{26}\mathrm{\,s}$. Concerning future prospects, HERA’s future measurements of the 21-cm power spectrum will probe DM lifetimes, $\tau_{\rm DM}\lesssim 10^{28}\mathrm{\,s}$ Sun:2023acy , potentially exploring most of the parameter space in Fig. 6. We show this as a dotted orange line in the figure.

It is important to mention that all these bounds are obtained assuming that the DM particle decays to an $e^{+}e^{-}$ pair. In our case, since the DM decays into $e^{+}e^{-}e^{+}e^{-}$, a correction comparable in magnitude to the lifetime itself is expected Slatyer:2016qyl . Nevertheless, given the large exponent of the DM mass in Eq. (42), this translates into a negligible correction to the DM mass and it is therefore inconsequential for Fig. 6.

Finally, the $\eta$ meson can be produced as a $s$-channel resonance from DM scatterings in galactic halos, as we discuss in more detail in Section 5. If its branching ratio to SM particles is large enough, severe constraints apply from DM indirect searches, see e.g. Ref. Chu:2018fzy . However, the decay rate of the $\eta$ meson is dominated by the $\theta$-induced process $\eta\to 2\pi^{0}$, while the branching ratio into electrons is extremely suppressed, yielding no relevant experimental constraint. For the benchmark point shown in Fig. 7, we obtain $\text{BR}(\eta\to\text{leptons})\sim\mathcal{O}(10^{-30})$. Larger values of the $\theta$ angle further reduce this branching ratio. Such minuscule sensitivities are far beyond the reach of any current experiment.

Bounds from dark photon searches.

The most stringent bounds on the parameter space under consideration arise from collider, fixed-target, and beam-dump experiments. The constraints change depending on the branching ratio between visible and invisible decays. For the choice of parameters of Fig. 6, we have $m_{V}>2m_{\eta}$ and $\alpha_{d}\gg\alpha_{\rm em}\epsilon$, hence the dark photon mainly decays to the dark mesons charged under $U(1)_{d}$, $V\to\pi^{+}\pi^{-},\,K^{+}K^{-}$ for charge assignment 2.777 For charge assignment 1, the dark photon decays as $V\to K^{+}K^{-},\,K^{0}/\overline{K^{0}}$, and the corresponding constraints remain unchanged. Thus, the relevant limits come from searches for invisible decays. The BaBar collaboration BaBar:2017tiz currently provides the most stringent bounds on invisibly decaying dark photons in the mass range of 500 MeV to 8 GeV. Future data from Belle II Belle-II:2018jsg are expected to significantly improve these bounds, reaching almost $\varepsilon\simeq 10^{-4}$. On the other hand, among the beam dump experiments, currently the most competitive constraint in the region relevant for thermalization below 500 MeV comes from NA64 NA64:2023wbi In the future, both NA64 NA64:2024nwj and LDMX LDMX:2018cma will explore even further this parameter space, see e.g. Ref. Antel:2023hkf .

5 Impact of the $\theta$ vacuum in DM halos

5.1 Reconciling the relic density with velocity-dependent scattering

Models in which DM production occurs via $3\to 2$ annihilations, commonly referred to as SIMP models, were initially proposed in the early 1990s Carlson:1992fn . These models experienced a resurgence of interest roughly a decade ago, as they naturally predict the observed DM relic abundance within regions of parameter space where the elastic self-interaction cross section per unit mass is of the order $\mathcal{O}(\text{cm}^{2}/\text{g})$. Scenarios with such interaction strengths are known as SIDM models Spergel:1999mh , which attracted considerable attention due to their ability to reduce the central densities of DM halos. This effect appears to alleviate the tension between predictions from collisionless CDM simulations Dubinski:1991bm ; Navarro:1995iw ; Navarro:1996gj and certain astrophysical observations Dave:2000ar ; Vogelsberger:2012ku ; Rocha:2012jg ; Peter:2012jh ; Elbert:2014bma ; Fry:2015rta .

However, it has become increasingly evident that SIDM models with constant elastic cross sections face significant tension with observational data from galaxy clusters, which impose upper bounds on the self-interaction cross section Harvey:2015hha ; Sagunski:2020spe ; Harvey:2018uwf ; Bondarenko:2017rfu ; DES:2023bzs :

[TABLE]

As typical DM velocities in galaxy clusters are larger compared to those in smaller halos, viable SIDM models now frequently incorporate velocity-dependent interactions Tulin:2017ara ; Adhikari:2022sbh . This poses a challenge for the original SIMP framework Hochberg:2014kqa , which, in its simplest form, does not accommodate velocity-dependent self-interactions. For example, in pion DM scenarios with $\theta=0$, the elastic scattering cross section remains constant at non-relativistic velocities. In contrast, the $\theta$-induced DM production mechanism considered in this work naturally yields velocity-dependent self-interaction cross sections that can evade cluster constraints while remaining sizable in smaller halos. These interactions are mediated by the $\eta$ meson acting as a resonance when $\theta\neq 0$, and are described by the non-relativistic Breit-Wigner formula Chu:2018fzy ; Chu:2019awd :

[TABLE]

where $v$ is the relative velocity of DM particles, $v_{R}$ is the resonant velocity defined in Eq. (21) and $\Gamma$ is the decay rate of the mediator. For BM1

[TABLE]

and $\Gamma$ is given in Eq. (22). In BM1, all DM today consists of $\pi^{0}$ due to the efficient conversion of $\pi^{\pm}$ (and heavier stable mesons) into $\pi^{0}$ in the early universe. Furthermore, we are assuming that the production of $\pi^{\pm}$ from $\pi^{0}$ in astrophysical halos is kinematically forbidden, which holds when $\delta\gg v^{2}\sim 10^{-5}$ today, excluding the gray region shown in the left panel of Fig. 2. As described in Section 4.1, in the presence of the dark photon portal, $\pi^{\pm}$ may become heavier due to radiative corrections. This makes its production in halos even harder, enlarging the viable parameter space.

Following the analysis of Ref. Chu:2018fzy , and for the benchmark point indicated by a star in Fig. 6, Fig. 7 demonstrates that the velocity dependence encoded in $\sigma(v)$ enables our QCD-like DM model to realize SIDM while remaining consistent with both cluster constraints and the observed relic density.

5.2 The SIDM parameter space

The ability of SIDM to explain cored density profiles arises from the predicted elastic self-scattering between DM particles in halos, which thermalizes their velocities and enables heat transport from the outer halo to the center, thereby flattening the central density profile—a mechanism absent in standard collisionless cold DM. This thermal evolution typically proceeds through two distinct phases. During the core formation phase, heat transfer from the hotter outer halo to the cooler central region leads to the development of an isothermal, constant-density core that replaces the steep central cusp predicted by collisionless models. Over longer timescales, continued energy transport causes the core to lose thermal support, transferring heat outward and initiating the core collapse phase, also known as gravothermal collapse. In this runaway process, the core contracts while its central density rises rapidly, analogous to the gravothermal catastrophe in globular clusters. The transition between these two phases depends on not only $\sigma(v)/m_{\rm DM}$ but also the halo mass and assembly history, and can be influenced by environmental effects such as tidal interactions. For a review with references to the original literature, we refer the reader to Ref. Adhikari:2022sbh .

Recent advances, largely based on non-resonant self-scattering cross sections, have identified the preferred values for $\sigma(v)/m_{\rm DM}$. In particular, successful models require strong velocity dependence, producing core formation in high-mass halos and gravothermal core collapse in low-mass halos Shah:2023qcw ; Ando:2024kpk , with outcomes influenced by assembly history and halo concentration Yang:2022mxl ; Nadler:2023nrd . Observational constraints span a broad range: as already mentioned, galaxy clusters with $v\sim 1000~\mathrm{km\,s^{-1}}$ imply $\sigma/m_{\rm DM}\lesssim 0.1$–$1~\mathrm{cm^{2}\,g^{-1}}$ Rocha:2012jg ; Harvey:2015hha ; Kaplinghat:2015aga ; Sagunski:2020spe ; Andrade:2020lqq ; galaxies with $v\sim 100~\mathrm{km\,s^{-1}}$ favor $\sigma/m_{\rm DM}\gtrsim 3~\mathrm{cm^{2}\,g^{-1}}$ Ren:2018jpt ; Nadler:2023nrd ; Roberts:2024uyw ; and sub-galactic systems with $v\sim 10~\mathrm{km\,s^{-1}}$ require $\sigma/m_{\rm DM}\gtrsim 5$–$10~\mathrm{cm^{2}\,g^{-1}}$ as well as core collapse in some halos Valli:2017ktb ; Read:2018pft ; Sameie:2019zfo ; Silverman:2022bhs . Furthermore, other analyses assuming the core formation phase Hayashi:2020syu ; Kamada:2023dse find that ultra-faint dwarf galaxies impose a stringent limit of $\sigma/m_{\rm DM}\lesssim 0.1~\mathrm{cm^{2}\,g^{-1}}$ at very low velocities, though uncertainties remain due to limited stellar kinematic data.

In addition to this, in the last few years there is growing interest in certain diffuse satellites and dense substructures suggesting $\sigma/m_{\rm DM}\sim 100~\mathrm{cm^{2}\,g^{-1}}$ or greater. This includes the lens system SDSS J0946+1006 Vegetti:2009cz , as it appears to host an unusually compact and over-concentrated dark subhalo, one of the rare dark halos smaller than a galaxy that has been directly detected by means of lensing Enzi:2024ygw . Its steep density profile sharply deviates from the predictions of collisionless DM, prompting investigations within the SIDM framework Minor:2020hic , particularly in the core-collapse phase. For instance, Ref. Tajalli:2025qjx has analyzed this system and found the observations could be explained by an SIDM halo undergoing gravothermal collapse, provided the self-interaction cross section is much larger than $\sigma/m_{\rm DM}\sim 100~\mathrm{cm^{2}\,g^{-1}}$. Another study has revisited the hypothesis of a core-collapsed SIDM halo, modeling its lensing signature via an isothermal Jeans approach. This showed that only a very massive core-collapsed halo (mass of order $10^{11}~M_{\odot}$) can reproduce the observations, but such a halo would be expected to host a luminous galaxy, which is not observed, posing a challenge to the SIDM interpretation Li:2025kpb . Overall, SDSS J0946+1006 remains a compelling but debated candidate for observing gravothermal collapse in SIDM.

Motivated by this, in Fig. 7 we indicate a lower bound Nadler:2023nrd ; Li:2025kpb of about $\sigma/m_{\rm DM}\simeq 150~\mathrm{cm^{2}\,g^{-1}}$ at intermediate velocities of order $60~\mathrm{km\,s^{-1}}$, consistent with our benchmark point. In fact, once $v_{R}\sim 100~\mathrm{km\,s^{-1}}$ is fixed, our QCD-like framework naturally predicts cross sections of this magnitude at $v\sim v_{R}$ as an unavoidable outcome of the model. Larger values of the $\theta$ angle further enhance the likelihood of producing these effects as they lead to larger cross sections. Interestingly, cross sections in this range have been proposed as a mechanism for seeding black holes, particularly in the context of SDSS J0946+1006, which may represent a striking phenomenological consequence of resonant SIDM in the early Universe Jiang:2025jtr ; Roberts:2024wup .

Taken at face value, these results suggest $\sigma/m_{\rm DM}\lesssim 0.1~\mathrm{cm^{2}\,g^{-1}}$ at both $v\gtrsim 1000~\mathrm{km\,s^{-1}}$ and $v\ll 30~\mathrm{km\,s^{-1}}$. This pattern implies a velocity-dependent cross section with a peak near $v\sim 100~\mathrm{km\,s^{-1}}$ and a steep drop toward lower velocities, consistent with models in which DM self-scattering proceeds via a resonant intermediate state by means of a Breit-Wigner cross section Kamada:2023dse ; Chu:2018fzy . While Fig. 7 is in qualitative agreement with the aforementioned effects, in this work we do not attempt a comprehensive global fit, including the relic density constraint, due to the ongoing investigations associated with resonant self-interacting DM. For instance, at present there is no $N$-body simulation for models exhibiting the velocity dependence of the Breit-Wigner form given in Eq. (53). Moreover, while a study based on a gravothermal fluid approximation has been conducted Kamada:2023dse , the core-collapse phase has been examined only for models with velocity-independent cross sections or for scenarios motivated by light-mediator interactions Tran:2024vxy ; Yang:2022hkm ; Yang:2022zkd ; Outmezguine:2022bhq , remaining insufficiently understood for Breit-Wigner cross sections. A detailed investigation will be presented elsewhere.

Thus, QCD-like theories with a non-zero $\theta$ parameter provide a novel mechanism for DM production through number-changing interactions, while naturally accommodating velocity-dependent self-interactions. These features allow such models to evade stringent cluster constraints and remain compatible with small-scale structure observations, potentially revitalizing the SIDM framework in light of recent astrophysical data.

6 Conclusions

Understanding the nature of DM remains one of the most compelling and unresolved challenges in contemporary physics. Within this context, QCD-like theories that feature dark pions as DM candidates offer a particularly attractive framework due to their structural similarity with the visible sector. A recent advance in this direction Garcia-Cely:2024ivo demonstrated that a non-zero topological $\theta$ angle can qualitatively alter the dynamics of dark pions, notably by inducing number-changing interactions that can drive thermal freeze-out consistent with the observed relic abundance. Moreover, the same mechanism naturally leads to resonant self-interactions that may influence the distribution of DM in galactic halos, thereby realizing the SIDM paradigm in a theoretically coherent manner.

In this work, we revisited the framework introduced in Ref. Garcia-Cely:2024ivo and addressed several open issues, particularly those related to the connection between the dark sector and the SM. By systematically analyzing the Boltzmann equations that govern the thermal history of DM, we have shown that some of the simplifying assumptions made in the original treatment, while valid, were unnecessarily restrictive.

To provide a concrete realization of the proposed scenario, we extended the benchmark model BM1 introduced in Ref. Garcia-Cely:2024ivo by introducing a dark photon mediator associated with an additional $U(1)_{d}$ gauge symmetry. This portal interaction not only ensures thermal contact between the dark and visible sectors but also has an impact on DM stability and indirect detection signatures. We demonstrated that the extended model remains consistent with current observational constraints, including those from indirect detection experiments and astrophysical probes.

Importantly, our analysis confirms that the strongly coupled dark sector can simultaneously account for the observed relic abundance and the small-scale structure anomalies that motivate SIDM models. The presence of resonant enhancements in the self-scattering cross section enables the required velocity dependence: enhanced interactions on galactic scales and suppressed interactions on cluster scales.

Looking forward, our results open several avenues for future research. From a theoretical standpoint, it would be valuable to explore alternative portal interactions and gauge group structures, particularly solving the Boltzmann equations for the DM distribution in phase space. On the observational front, upcoming data from indirect detection experiments, gravitational lensing studies, and precise measurements of dwarf galaxy dynamics may offer critical tests of the velocity-dependent SIDM framework outlined in this work. This is particularly timely in light of systems such as SDSS J0946+1006, which are claimed to host an SIDM halo undergoing gravothermal collapse. A thorough assessment of this requires a detailed study of the gravothermal evolution during the core-collapse phase of DM halos, in which particles self-scatter via a Breit–Wigner resonance.

In summary, our findings demonstrate that QCD-like DM theories with a non-zero $\theta$ angle furnish a viable and predictive framework that naturally unifies thermal freeze-out with small-scale structure phenomenology. By establishing a connection between relic abundance and halo dynamics through a common underlying mechanism, these models provide a compelling DM sector.

Acknowledgments

We are indebted to Ayuki Kamada for insightful and valuable discussions. GL is supported by the Generalitat Valenciana APOSTD/2023 Grant No. CIAPOS/2022/ 193. C.G.C. is supported by a Ramón y Cajal contract with Ref. RYC2020-029248-I, the Spanish National Grant PID2022-137268NA-C55 and Generalitat Valenciana through the grant CIPROM/22/69. L.M. is supported by the Generalitat Valenciana through the grant CIACIF/2023/91. O.Z. has been partially supported by Sostenibilidad-UdeA, the UdeA/CODI Grant 2022-52380, and Ministerio de Ciencias Grant CD 82315 CT ICETEX 2021-1080. O.Z. also acknowledges the support of the Fundación Carolina and the hospitality of the Instituto de Física Corpuscular UV-CSIC in Valencia, Spain.