Split SIMPs with Decays

We discuss a minimal realization of the strongly interacting massive particle (SIMP) framework. The model includes a dark copy of QCD with three colors and three light flavors. A massive dark photon, kinetically mixed with the Standard Model hypercharge, maintains kinetic equilibrium between the dark and visible sectors. One of the dark mesons is necessarily unstable but long-lived, with potential impact on CMB observables. We show that an approximate"isospin"symmetry acting on the down-type quarks is an essential ingredient of the model. This symmetry stabilizes the dark matter and allows to split sufficiently the masses of the other states to suppress strongly their relic abundances. We discuss for the first time the SIMP cosmology with sizable mass splittings between all meson multiplets. We demonstrate that the SIMP mechanism remains efficient in setting the dark matter relic density, while CMB constraints on unstable relics can be robustly avoided. We also consider the phenomenological consequences of isospin breaking, including dark matter decay. Cosmological, astrophysical, and terrestrial probes are combined into a global picture of the parameter space. In addition, we outline an ultraviolet completion in the context of neutral naturalness, where confinement at the GeV scale is generic. We emphasize the general applicability of several novel features of the SIMP mechanism that we discuss here.


Introduction
In recent years, the strongly interacting massive particle (SIMP) framework [1] has emerged as an attractive possibility for thermal dark matter, alternative to the traditional weakly interacting massive particle (WIMP) paradigm. The SIMP relic density is set by the freezeout of 3 → 2 self-annihilations, whose parametrics naturally point toward masses comparable to the strong scale, roughly between 10 MeV and 1 GeV, and strong coupling. If kinetic equilibrium between the dark matter and the Standard Model (SM) bath is maintained until the 3 → 2 processes freeze out at T ≈ m DM /20, the dark matter remains sufficiently cold to avoid conflict with structure formation bounds [1], which otherwise exclude [2] a completely secluded 3 → 2 freezeout [3]. The SIMP mechanism finds its most attractive realizations in the context of confining gauge theories with chiral symmetry breaking [4], where the pseudo Nambu-Goldstone bosons (pNGBs) play the role of dark matter, and 3 → 2 annihilations are mediated by the Wess-Zumino-Witten (WZW) action [5,6]. The pNGBs are also naturally characterized by strong 2 → 2 self-scattering [1,4], which may help to address possible shortcomings of collision-less cold dark matter (CDM) on small scales (as reviewed for example in Ref. [7]). While such issues may eventually be resolved without modifying the CDM paradigm, dark matter that self-interacts with cross sections in the range σ/m DM ∼ 0.1 -10 cm 2 /g is a very interesting possibility in this respect; see Ref. [8] for an extensive review.
Existing studies of pNGBs as SIMP dark matter generally consider scenarios where all the pNGB mesons are (approximately) mass-degenerate, and the (dominant) component of dark matter is stable. The aim of this work is to study a minimal realization of the SIMP mechanism where variations of such properties can be explored. We consider a dark copy of the SM QCD, consisting of an SU (N c ) gauge theory with N c = 3 and N f = 3 light hidden quark flavors, which is the smallest N f that admits a WZW action [6], necessary for the realization of the 3 → 2 processes central to the SIMP mechanism. Hidden electromagnetism is gauged by a massive dark photon kinetically mixed with the SM hypercharge, providing a viable mediation mechanism that maintains kinetic equilibrium between the hidden and SM sectors [9,10].
It was pointed out in Ref. [11] that when N f is odd, one of the pNGBs is necessarily unstable. We show that, for typical parameters, the unstable meson η decays to SM particles with lifetime comparable to the timescale of recombination, potentially leading to strong constraints from cosmic microwave background (CMB) anisotropies [12,13]. We show that an approximate SU (2) U global symmetry acting on the down and strange quarks, which we refer to as "isospin," is a crucial component in this setup. This symmetry plays two key roles. First, it stabilizes the lightest multiplet, a triplet of dark mesons with mass in the 100 -300 MeV range. Second, it allows for separation between the masses of the up quark and the degenerate down-type quarks, which in turn raises the η mass relative to the dark matter mass, suppressing the η abundance to a level allowed by CMB measurements. In this regime, the masses of all SU (2) U multiplets are separated by similar, sizable amounts, raising another important question: namely, whether the 3 → 2 freezeout remains effective even in this scenario of larger mass splittings. Our analysis provides a positive answer, opening up new parameter space for the SIMP mechanism. These results can be easily generalized to other models with odd N f .
In this context, we address another question that has remained surprisingly understudied in the literature: the absolute stability of SIMP dark matter. The neutral pion, which is one of the components of the dark matter triplet, decays through its mixing with η induced by small breaking of isospin. We analyze quantitatively the sensitivity of current dark matter indirect detection searches to the order parameter of isospin breaking, finding that it should not exceed O(10 −5 ), and present the projected reach at future experiments. We combine these astrophysical and cosmological probes with laboratory tests of the dark photon mediator, painting a global picture of the parameter space. We emphasize that many of our results have broader applicability, beyond the minimal model adopted here.
A brief discussion of the abundance and decays of unstable SIMP mesons was presented in Ref. [14], albeit in a setup with N f = 4. One of the main novelties of our work is a quantitative analysis of the SIMP cosmological history, focusing on larger mass splittings than previously considered in the literature and highlighting the key role of CMB anisotropy bounds on the unstable mesons. We work in pure chiral perturbation theory for the pNGBs, neglecting resonances such as the vector mesons, whose role has been extensively discussed in Ref. [11] (see also Ref. [15]).
As is characteristic of the SIMP framework, we find large dark matter self-scattering cross sections [4]. Our minimal choice of N c = 3 leads to σ/m DM ∼ few cm 2 /g, in some tension with bounds from the Bullet cluster [16] and halo shapes [17,18]. However, given the evolving status of the small-scale CDM puzzles, we believe that it would be premature to discard the minimal and theoretically appealing setup analyzed here. Furthermore, our results can serve as a useful basis for building models that feature smaller self-interaction cross sections.
Theoretical motivation for the scenario discussed here comes from neutral naturalness theories, such as the Twin Higgs [19], which address the little hierarchy problem by introducing top partner particles that are not charged under SM color. Neutral naturalness models typically single out N c = 3 in the hidden sector, as this allows the top partner to cancel the top quark loop correction to the Higgs mass. In addition, as argued in Ref. [20], two-loop naturalness considerations suggest that hidden color should confine at a scale similar to that of the SM QCD. Motivated by these arguments, we outline an embedding of the low-energy theory into a neutral naturalness model, along the lines of the vector-like Twin Higgs [21]. In contrast with the "Twin SIMPs" setup of Ref. [14], we do not introduce a full mirror copy of the SM, but instead propose a minimal construction for SIMP dark matter where the light dark quarks and the top partner(s) are more directly linked.
The remainder of our paper is structured as follows. In Section 2 we introduce the effective theory for the hidden mesons, discussing in detail its symmetries, mass spectrum, and leading interactions. We outline ultraviolet completions in the framework of neutral naturalness in Section 3, which can be skipped by readers who are only interested in dark matter phenomenology. The lifetimes of the unstable mesons are calculated in Section 4. In Section 5 we discuss in detail the cosmological history, while signatures and constraints are presented in Section 6. Finally, we conclude in Section 7 with a brief summary and outlook. Appendix A provides complete Boltzmann equations for our setup.

Effective theory for the hidden mesons
In this section we take a phenomenological approach and discuss the effective field theory (EFT) of hidden QCD, remaining agnostic about specific ultraviolet (UV) completions. Possible UV completions in the framework of neutral naturalness are addressed in Section 3.
We assume an SU (N c ) hidden color gauge group with N f = 3 light hidden quark flavors. All our numerical results assume N c = 3, although we occasionally comment on the effect of changing the number of colors. We also introduce a massive dark photon, kinetically mixed with the SM hypercharge, which keeps the hidden and SM sectors in kinetic equilibrium 1 via elastic scattering at least until the 3 → 2 processes among the dark matter (DM) particles freeze out [9,10], as required for SIMP DM [1]. 2 The Lagrangian for the dark photon is After diagonalization of the gauge kinetic and the mass terms, the physical dark photon A couples to the SM electromagnetic (EM) current at O(ε). For fermions, this interaction reads εec w Q f A µf γ µ f . We do not specify the origin of the A mass, which could arise from the Stückelberg mechanism or from the coupling to a dark Higgs field. If the three quarks u, d, s are light compared to the confinement scale, the pattern of low-energy spontaneous symmetry breaking is SU where λ a are the Gell-Mann matrices, satisfying Tr(λ a λ b ) = 2δ ab . Our normalization is such that f SM π ≈ 92.4 MeV, and the cutoff of the EFT is which is also our definition of the dark strong coupling scale. The terms of the chiral Lagrangian most relevant to our discussion are (see e.g. Ref. [25]) where Σ → LΣR † , and likewise for the quark mass spurion M . The interactions with external vector fields are described by spurions with formal transformation properties Q L → L Q L L † and Q R → R Q R R † , but in Eq. (2.4) we have already set Q L = Q R = Q as appropriate to describe the coupling to hidden EM. The covariant derivative of Σ reads 1 An axion-like particle [22,23] or vector mesons [11] have also been studied as mediators for pNGB SIMP dark matter. 2 If the elastic scattering decouples before the 3 → 2 processes, elastically-decoupling relic DM, or ELDER, can be realized [24].
. The second and third lines of L display pieces of the WZW action that arise due to the presence of the gauge fields. In particular, the second line is relevant to the calculation of π → A * A * decays, whereas the third line is responsible for the semi-annihilations ππ → πA . The fourth line shows the piece of the WZW action that controls 3 → 2 scattering among the mesons; it is nonzero only if all five participating mesons are different.
The SM values of the electric charges are More generally, assuming the underlying microscopic theory to be vector-like (so that the dark baryon number B is anomaly-free), we have the freedom to gauge a linear combination of Q and B. In particular, gauging results in the vanishing of the axial-vector-vector (AVV) anomalies as a consequence of Tr(Q 2 λ a ) = 0, ensuring that even singlet mesons do not decay through anomalous diagrams. In this paper we consider both possibilities, Q = Q and Q = Q , for the charges. We will require that the semi-annihilations ππ → πA be sub-leading to 3 → 2 processes during DM freezeout (see Section 5.5), which imposes a lower bound m A 2m π on the mass of the dark photon [10].

Mass spectrum
Setting the mass spurion to its physical value M = diag (m u , m d , m s ) gives the following pNGB masses, 7) where ∆m 2 em = 2cê 2 f 2 π is the electromagnetic correction, with c being a constant. The π 3 and π 8 mix as where B is an a-priori unknown parameter of order the strong scale Λ, defined in Eq. (2.3).
In the SM, the π 0 mass gives B SM = 2.7 GeV (B has dimension of mass, so the natural dimensionless parameter is b = B/Λ with b SM = 2.3) and the π ± -π 0 mass difference gives c SM = 0.8. In the following we simply take b ∼ O(1), whereas the EM correction c has a different scaling compared to the SM, since we assume the dark photon is heavy, Ref. [26]), implying that this correction is small throughout our parameter space.
For phenomenological reasons that will become clear momentarily, we focus on the scenario where the up quark is moderately heavier than the down-type quarks, which are approximately degenerate. Assuming m u + m d > 2m s , the above matrix is diagonalized as The rotation matrix is defined as R(ϕ) ≡ c ϕ s ϕ −s ϕ c ϕ , employing short-hand notations for sine and cosine that we use throughout the paper (in particular, s w and c w refer to the weak mixing angle). When the mass splittings are neglected, we denote the common octet mass simply as m π .
As we have anticipated, we focus on a scenario where the down-type quarks are nearly degenerate due to an approximate SU (2) symmetry. Hence, we parametrize the dark quark masses as (m u , m d , m s ) = (m + ∆m, m, m + dm), (2.11) where ∆m > 0 is a sizable splitting of O(0.1 -1)m, whereas dm can have either sign but is very small, |dm|/∆m 1. The meson masses in this limit are where in the ordering of the charged mesons we have assumed dm > 0 for concreteness. The mixing angle between the π 0 and η reads where we have defined δ ≡ dm ∆m . (2.14) This parameter measures the strength of the SU (2) U breaking compared to the chiral SU (3) breaking; however, since in most of our discussion we consider ∆m/m ∼ O(1), we can take δ as effectively measuring the strength of isospin breaking. We also define the relative splittings ∆ π,η ≡ m π + ,η − m π 0 m π 0 , (2.15) which will be used frequently in later discussions. In the isospin-symmetric limit δ = 0, the mesons transform as 3 0 , 2 ±1 and 1 0 under the SU (2) U × U (1) Q symmetry, which is exactly preserved by the quark masses and by rep. meson quark content mass stable?
no; decays via AVV anomaly or higher-order operators (Section 4.1) stable in isospin-symmetric limit; decays via ∝ δ mixing with η (Section 4.2) Table 1. Summary of the properties of the hidden mesons, ordered by decreasing mass (assuming δ > 0). The first column lists the representation under the SU (2) U × U (1) Q global symmetry. The quark contents correspond to the isospin-symmetric limit. Some of the masses are approximate; complete expressions are given in Section 2.1.
electromagnetism. 3 In particular, the SU (2) U triplet, which constitutes the DM, is formed by (π 0 , K 0 , K 0 ), the doublets are (π + , K + ) and (π − , K − ), and the singlet is η. In the SM context the SU (2) U is often called U-spin, justifying its name, but for simplicity we refer to it as "isospin." We assume that U (1) Q is a good low-energy global symmetry, despite a heavy A , implying, in particular, that all the charged pions and kaons are (almost) massdegenerate and stable. This is guaranteed if m A arises from a Stückelberg mechanism, whereas it could be a good approximation in certain realizations of a dark Higgs mechanism.
In Table 1 we present the explicit field contents of the dark mesons in the isospinsymmetric limit, as well as an overview of other salient properties, some of which will be analyzed later. With the exception of the η, for δ = 0 the mesons cannot decay at any order in the chiral Lagrangian. As we discuss in Section 4, the η does decay even in the isospinsymmetric limit, while the π 0 decays via small isospin-breaking effects. A stable SU (2) triplet of SIMP DM mesons composed of hidden d, s quarks was previously considered in Ref. [14], albeit in a theory also containing light and degenerate u, c quarks, leading to an extended pattern of chiral symmetry breaking.
Note that choosing ∆m < 0 in Eq. (2.11) leads to a spectrum where the singlet is the lightest meson and therefore comes to dominate the hidden sector abundance after freezeout. As the singlet is unstable even in the isospin-symmetric limit, with a lifetime much longer than one second but not arbitrarily long due to the requirement of thermalization between the hidden and SM sectors, this possibility is not viable.
The SIMP mechanism typically requires rather large values of m π /f π , and in this work we consider m π 0 /f π ∼ 8 -10 (for reference, in the SM m K /f π ≈ 5.4). As the stronglycoupled regime is approached, higher-order corrections in chiral perturbation theory be-come important [28]. Furthermore, when the pNGBs are heavy, additional resonances can play an important role in the dynamics. In particular, Refs. [11,15] explicitly introduced vector mesons in the effective theory, showing that this opens up additional parameter space for hidden meson DM. In general, the effects of the vector mesons can be neglected as long as their masses satisfy m V 2m π . Here, we assume this is the case and focus on a minimal framework, neglecting resonances.

Ultraviolet completions in neutral naturalness
We now outline possible UV completions of the chiral Lagrangian presented in the previous section in the context of neutral naturalness. The reader who is only interested in the phenomenological aspects of our work, or perhaps favors a different class of completions, may choose to proceed directly to Section 4 without loss of continuity.
The minimum requirements for a neutral natural theory that realizes the meson spectrum discussed in Section 2 are: • two light and degenerate down-type quarks; • one moderately heavier (but still lighter than the confinement scale) up-type quark; • one top partner that cancels the quadratic UV sensitivity of the Higgs mass induced by the top Yukawa coupling.
The simplest construction with these characteristics can be built along the lines of the vector-like Twin Higgs [21], where = iσ 2 . The vector-like masses for the fermion fields are assumed to be small perturbations to the Yukawas and not written explicitly. The fields have the following charges under the twin SU (2) L × U (1) Y symmetry, 2) where we have assumed the same hypercharge assignments as in the SM. The fields Q, u, d transform as triplets of the twin SU (3) c , and Q c , u c , d c transform as anti-triplets. We assume that the global (approximate) SU (4) is non-linearly realized and the radial mode is heavier than the cutoff. Therefore, in unitary gauge the twin Higgs doublet is parametrized as Ifŷ t = y t is enforced by a Z 2 symmetry andŷ t ,ŷ b ,ŷ b y t , the quadratic correction to the Higgs mass from the SM top loop is canceled by the twin top with mass (f /v)m t . In addition, gauging the twin SU (2) L with Z 2 -symmetric couplingĝ = g ensures that the leading gauge corrections to the Higgs mass cancel as well. This is a simple vector-like Twin Higgs scenario with light twin quarks, leading to chiral symmetry breaking in the hidden sector (see Ref. [29] for similar ideas).
Naively, takingŷ t ŷ b =ŷ b in Eq. (3.1) seems to provide exactly the light quark spectrum we desire. However, the fact that y t =ŷ t ŷ t leads at one loop level to different renormalization of the two down-sector Yukawas via diagrams involving the twin W . The natural value of the isospin-breaking parameter is therefore one loop factor, significantly exceeding δ 10 −5 , which is necessary to render the π 0 sufficiently long-lived to be a viable DM component. Thus, in this model an acceptably small δ can be achieved only at the price of fine tuning.
This problem can be fixed by requiring thatŷ t =ŷ t , thus considering instead the Lagrangian In this case the little hierarchy problem can be solved ifŷ t = y t / √ 2, with two heavy, degenerate top partners of mass (f /v)m t / √ 2 now canceling the SM top loop. We envisage that a UV completion of Eq. (3.4) may be constructed with Orbifold Higgs methods [30,31], although we do not attempt to do so here. We have introduced an additional vector-like fermion u , u c in order to have a light up-type quark in the spectrum. Our phenomenological study shows that only a moderate coincidence of scales, M u ŷ b f within an O(1) factor, is necessary for a viable SIMP scenario. Equation (3.4) provides a technically natural setup that preserves the isospin SU (2) U . Additional masses and interactions can then provide small breaking of isospin. 4 We can estimate the confinement scale for hidden QCD, Λ QCD , 5 by requiring [20] that the visible and hidden color gauge couplings are approximately equal at the scale Λ UV 4πf , where the theory needs to be extended. Including, in addition to the three light quark flavors, the two degenerate top partners with mass (f /v) m t / √ 2 , and allowing for |(ĝ s − g s )/g s | < 0.2 at Λ UV = 5 TeV, we obtain 0.12 < Λ QCD /GeV < 4.7, where we took f = 750 GeV for illustration. 6 The mass scale required for viable 3 → 2 freezeout of the dark mesons falls toward the lower end of this range.
While the Higgs portal interactions decouple early, at temperatures around a few GeV, the kinetic mixing between the hidden hypercharge gauge fieldB and its SM counterpart 4 A priori, another possibility to obtain the desired spectrum is to assume that the Yukawa couplings are negligible in the down sector and the leading contribution to the down-type quark masses comes from MQ = M d = M . In this case we have a single top partner withŷt = yt, the required coincidence of scales for the up quark mass readsŷ t f M , and u , u c are not needed. However, diagrams involving the twin EW gauge bosons introduce isospin breaking with one-loop size, requiring fine tuning to achieve a phenomenologically viable δ 10 −5 .
5 Note that ΛQCD is distinct from Λ defined in Eq. (2.3). 6 We used 2-loop running and αs(mZ ) = 0.118 as input. For reference, the same running procedure applied to the SM leads to ΛQCD, SM ≈ 370 MeV.
can maintain kinetic equilibrium between the two sectors until the DM freezes out. The gauge Lagrangian reads The diagonalization of the kinetic and mass Lagrangian for the four neutral gauge fieldŝ W 3 ,B, W 3 , B was performed in Ref. [32]. For mB < m Z the mass of the physical dark photon is m A ≈ĉ w mB , whereĉ w is the cosine of the twin weak mixing angle, and the A coupling to the SM fermions f is εec wĉw Q f . We do not specify whether mB comes from a Stückelberg mechanism (see Ref. [33] for a review) or from a dark Higgs. In either case, attention must be paid to avoid introducing an additional naturalness problem related to the A mass, which needs to satisfy Λ m A < m Z for viable SIMP phenomenology. An embedding of SIMP DM in the Twin Higgs framework was presented in Ref. [14], where a complete mirror spectrum was introduced and an exact SU (2) flavor symmetry relating the first two generations was imposed. The more minimal proposal sketched in Eq. (3.4) aims at a more direct connection between the heavy degrees of freedom essential for Higgs naturalness and the light quark spectrum dictating DM phenomenology. We leave a detailed study of this completion to future work.

Dark meson decays
In this section we calculate the lifetimes of the hidden mesons, which provide key inputs to the analysis of cosmological and astrophysical constraints discussed below in Section 6. First we focus on the singlet η, the fastest-decaying meson since it is not protected by any symmetry. We then consider decays of the DM components: we discuss π 0 decays induced by small isospin breaking, and show that the neutral kaons are accidentally stable even when the isospin is broken by the quark masses, due to a residual U (1) symmetry.

η decay
The η, being a singlet under SU (2) U × U (1) Q , is unstable even when this is an exact symmetry. The diagrams mediating its decay are shown in Fig. 1. The strength of the ηA * A * interaction, and therefore the η lifetime, depend on the choice of hidden electric charges.
For the standard choice of Q, the coupling is dominated by the AVV anomaly (see the second line of Eq. (2.4)). For the tree-level 7 decay to four electrically-charged SM fermions, since m A > m η , the amplitude can be matched to the effective operator Here we have taken two distinct fermion-antifermion pairs to avoid subtleties with identical particles; we will adjust the symmetry factors as appropriate when discussing the relevant case η → 4e. We have computed the decay width corresponding to the operator in Eq. (4.1) with the help of FeynRules [34] and MadGraph5 [35], assuming massless fermions. The result is Focusing on the decay to 4e, we find the lifetime where we have included an extra factor of 2 to approximately correct for the identical particles. For the helicity-suppressed, one-loop decay to two fermions we estimate where in the second line we have focused on the decay to µµ, which is kinematically open in part of the parameter space we consider. Note that when m η 2m µ , the eeµµ channel is also open.
With the Q charges, the AVV anomalies vanish but the ηA * A * interaction is still generated by higher-order, O(p 6 ) operators in the chiral Lagrangian [11]. Operators with one insertion of the quark mass matrix include for example , the operator in Eq. (4.5) gives a larger contribution. By matching to Eq. (4.1), we find where we have taken B = 4πf π and N c = 3. In addition, there are O(p 6 ) operators without mass insertions but containing additional derivatives, such as This operator leads to K , the same result as in Eq. (4.7) up to a sign. In the remainder of this paper we assume that the decays of η are mediated by d 1 , but it should be kept in mind that this coefficient actually represents a combination of several coefficients of O(p 6 ) operators.
The result in Eq. (4.7) illustrates that the decay through higher-order operators is strongly suppressed close to the chiral limit. However, in practice we consider mesons that are only moderately lighter than 4πf π , and the resulting lifetime is only mildly longer than in the anomalous case, For the two-body decay to µµ, our estimate is Note that the above amplitudes do not mediate the decay of η to SM π 0 + − , as axialvector couplings are not involved. The results for the decay to ff in Eqs. (4.4) and (4.10) only represent rough estimates and should therefore be taken with some caution.
In the anomalous case we find that the amplitude for π 0 → 4f is given by Eq. (4.1) with the replacement 2 (4.11) For the one-loop helicity-suppressed decay to µµ we similarly find (4.12) In the scenario without AVV anomalies we obtain iTr(Q M Σ † ) + h.c. ⊃ −mδπ 0 /f π , and from Eq. taking B = 4πf π . Our estimate for the decay to µµ is While detailed constraints on DM decay are discussed in Section 6.2, we anticipate that the π 0 lifetime must be longer than about 10 25 seconds to be phenomenologically viable.
For the neutral kaon K 0 , we note that the explicit isospin breaking caused by dm preserves a residual U (1) symmetry under which the quarks have charges Q r = diag (0, 1, −1), since Q r and the quark mass matrix M commute. Under this symmetry the complex mesons have charges Q r (π + ) = −1, Q r (K + ) = +1, and Q r (K 0 ) = +2 . It follows that K 0 , being the lightest particle charged under U (1) Qr , is accidentally stable.

Cosmological history
In this section we discuss various aspects of the cosmological history of our setup. We assume that the hidden and SM sectors are in equilibrium in the early Universe (the necessary conditions will be detailed below) and study the evolution and freezeout of the dark meson abundances. Where relevant, for convenience we use the labels L, M, H to refer to the light (π 0 , K 0 , K 0 ), middle (π ± , K ± ), and heavy (η) meson states, respectively, and m L,M,H to indicate their masses.

Freezeout and relic abundances
We begin by discussing the thermal freezeout and relic abundances of the dark mesons. As is well known, the 3 → 2 strong interaction processes mediated by the WZW term can remain efficient in keeping dark mesons in chemical equilibrium down to x ≡ m π 0 /T ≈ 20, thereby setting the relic density of the lightest multiplet to the observed value Y DM ≈ 4.1 × 10 −10 (GeV/m DM ), where Y ≡ n/s is the number density normalized to the total entropy. We will see that 2 → 2 scattering processes among dark mesons remain active down to much lower temperatures, strongly suppressing the abundances of the heavier multiplets, in particular of the unstable η. The complete Boltzmann equations that we solve to obtain the evolution of the meson abundances in the early Universe are provided in Appendix A.
The evolution of the meson abundances for these parameters is shown in the left panel of Fig. 2. Early on, for m π 0 /T < 20, number-changing 3 → 2 processes and dark meson-SM scatterings keep all dark mesons in chemical and kinetic equilibrium with the SM bath, so that their abundances follow the Boltzmann-suppressed thermal distributions with Y η < Y π + < Y π 0 . The 5-point interactions enable every meson to scatter with two states belonging to the lightest multiplet, which are the most abundant of the dark mesons; hence all 3 → 2 processes freeze out roughly when (n eq π 0 ) 2 1 6 As is characteristic of the SIMP mechanism, this freezeout occurs roughly at x ≈ 20, after which the abundances of π 0 , K 0 , K 0 remain approximately constant. The abundances of the heavier mesons continue to deplete even after the 3 → 2 processes freeze out, thanks to 2 → 2 annihilation processes such as HH → M M, M M → LL, M M → HL, which remain active. 9 These interactions maintain the densities of the heavier mesons on the "shifted" equilibrium curves Y eq i (Y π 0 /Y eq π 0 ) Y π 0 e −∆ i x [14]. Note that these annihilations provide negligible corrections to the frozen-out abundances of the L states, as Y η , Y π + Y π 0 at this stage. Assuming kinetic coupling between the hidden and SM sector, our numerical results show that the 2 → 2 annihilations begin to freeze out around x ≈ 400 (see solid curves in the left panel of Fig. 2). To estimate the freezeout abundance Y π + , we can apply the instantaneous freezeout approximation to the M M → LL process as where the thermally-averaged 2 → 2 cross section is defined in Eq. (A.7). For our parameters, this gives x fo + ≈ 420 and Y fo π + ≈ 4.5 × 10 −19 , in good agreement with the numerical results.
The evolution of the η abundance is more complex. After 3 → 2 freezeout, Y η decreases due to the HH → M M, HH → LL processes, as well as the thermally driven HL → M M . Around x ≈ 400, all of these processes become inefficient at depleting η; as a consequence, Y η departs from the shifted equilibrium curve and undergoes a short period of increase due to injections from the M M → HL processes, which are still active. The increased η density eventually leads to a re-coupling of HH → M M , and the final ratio of the η and π + abundances is determined by detailed balance between these two processes, For the spectrum under consideration, this ratio is ≈ 0.89. The key to understanding this behavior is to observe that both M M → HL and HH → M M are kinematically allowed at T = 0 ; hence the freezeout abundance is driven not by the familiar Boltzmann suppression, but by detailed balance between different processes. Furthermore, note that since Y η Y π + at x fo + (when the M states freeze out), such interplay only gives a mild correction to Y fo π + . The above estimates assume that the SM and dark sectors remain in kinetic equilibrium throughout, via the scattering of π ± , K ± on electrons. However, depending on the value ofŷ ≡ εê(m π 0 /m A ) 2 , the two sectors may decouple at some point in the evolution. After decoupling, the hidden sector redshifts like matter and thus its temperature decreases faster, T D ∝ a −2 (with a the scale factor), compared to the SM bath temperature, which decreases as T ∝ a −1 , giving a ratio T D /T x dec /x for x > x dec . To take this effect into account, Eq. (5.3) can be modified as Assuming x dec = 25, i.e. the two sectors decouple immediately after the light states L freeze out, results in x fo + ≈ 110 (corresponding to (x fo + ) 2 /x dec ≈ 450) and Y fo + ≈ 1.1 × 10 −19 , a factor ∼ 4 smaller than in the coupled scenario. These estimates are borne out in the numerical analysis, see the dotted curves in the left panel of Fig. 2. Earlier decoupling of the two sectors therefore results in further O(1) suppression of the M and H abundances.
In the right panel of Fig. 2, we consider a scenario where M M → HL scattering is kinematically closed at T = 0. This can occur if higher-order effects raise the η mass above (2m π + − m π 0 ). Corrections to the masses of the pNGBs arise from O(p 4 ) operators with two insertions of the quark mass matrix M [36], among which only affects the η mass [11]. The NDA size of the coefficient is c 7 ∼ O(1). In the SM, its sign is negative due to η -η mixing. In our framework, however, since we require m π /f π ∼ 8 -10, which is much larger than in the SM, we cannot a priori exclude the scenario where η would be lighter than η, resulting in c 7 > 0 and thus raising the η mass. In the right panel of Fig. 2 we illustrate this scenario, assuming the input parameters B = 265 MeV, m = 76 MeV, ∆m = 9 MeV, and c 7 = +2. The η abundance continues to be efficiently depleted by the HL → M M processes even at very low temperatures, resulting in an enormous suppression. While this alternative possibility deserves to be kept in mind, in the light of the above discussion it appears somewhat less likely on theoretical grounds. Therefore, we neglect it in the rest of the discussion, focusing solely on the leading-order meson spectrum.

5.2
Consequences of π 0 -K 0 mass splitting In the previous subsection, the small isospin breaking introduced by dm could be safely neglected. At later times, however, this effect can play a role in determining the relic abundances of the DM components through K 0 K 0 ↔ π 0 π 0 processes. According to Eq. (2.12), isospin breaking generates a small splitting between the π 0 and K 0 masses, For δ 10 −5 , necessary to avoid observational bounds on DM decays (see Section 6.2), this splitting is O(10 −4 ) eV. The splitting enables K 0 K 0 → π 0 π 0 annihilation to slowly convert part of the kaon population into π 0 's. We expect this process to freeze out when where we have made the assumption that this freezeout occurs after kinetic decoupling of the two sectors. The final K 0 abundance is Y fo K 0 = Y π 0 e −∆ K x 2 fo /x dec . For δ = 10 −5 and x dec = 20 -200, we find that freezeout occurs at T fo D ∼ 2 -6 × 10 −4 eV (corresponding to T fo ∼ 50 -20 eV), to be compared with the mass splitting of 2.5 × 10 −4 eV, and the kaon abundance is moderately depleted to Y fo K 0 /Y π 0 ∼ 0.3 -0.6. For smaller δ 10 −6 , the freezeout temperature in the dark sector is much larger than the mass splitting, and the kaon density is not appreciably suppressed.
This mild depletion of the kaon abundance can, however, be reversed in the late Universe, when kaons are regenerated via π 0 π 0 → K 0 K 0 up-scatterings in DM halos. The regenerated fractional density can be roughly estimated as [37] where τ int is the "integration time" over which scatterings occur. After relating the forward and backward reactions and assuming ∆ K v 2 , the right-hand side of Eq. (5.9) becomes where we have taken as reference the DM density and velocity dispersion relevant for galaxy clusters. Therefore, kaons and pions are expected to become equally distributed within ∼ 2 Gyr. While this is only a rough estimate, it suggests that it is appropriate to assume that DM is composed equally of π 0 , K 0 , and K 0 for late-Universe phenomena such as cluster mergers. Similar considerations apply to galaxies, where the typical velocity dispersion is v ∼ 250 km/s, and possibly even to dwarf galaxies, where the average velocity can be as low as 30 km/s but the DM density is larger.
The masses of π ± and K ± are also split by isospin breaking, but at O(δ), as can be read in Eq. (2.12). However, since the densities of the charged mesons are extremely suppressed, this effect does not have appreciable phenomenological consequences, and we neglect it.

Larger mass splittings
We now comment on scenarios where the mass splitting between the lightest (L) and nextto-lightest (M ) multiplets becomes appreciable. We will show that as long as 3m L > 2m M , the SIMP mechanism proceeds largely as in the degenerate case. We first provide simple analytical arguments supporting this conclusion, followed by detailed numerical results.
Recall that the L freezeout abundance is determined from detailed balance between the 3 → 2 process LLL → M M and its inverse 2 → 3 process M M → LLL. For comparable masses m L ≈ m M , the inverse process needs to make up energy equivalent to the mass of a single particle, incurring a corresponding Boltzmann suppression e −m L /T from the thermal tail. In other words, the rate of the inverse process is proportional to L e −m L /T , whereas the rate of the forward process is proportional to Y 3 L , hence detailed balance between the two gives the familiar Boltzmann-suppressed abundance Y L ∼ e −m L /T . When the mass splitting becomes larger than the bath temperature T but LLL → M M remains open without kinematic suppression, i.e. 3m L > 2m M , the inverse process M M → LLL needs to make up a smaller amount of energy, 3m L − 2m M , and the incurred thermal suppression e −(3m L −2m M )/T appears much weaker. However, note that the M abundance itself is further suppressed in this case, Y M ∼ e −m M /T , hence the rate of the inverse process is proportional to Y 2 M e −(3m L −2m M )/T = e −3m L /T , which is the same as in the degenerate scenario. Freezeout therefore largely proceeds as in the degenerate case, resulting in the same freezeout temperature T 3→2 fo , and the SIMP mechanism still produces the correct relic density.
We verify these findings numerically with the following mass spectrum, obtained by taking B = 4πf π = 222 MeV, m = 51 MeV, and ∆m = 45 MeV. The relative splittings with respect to the lightest multiplet are ∆ π = 0.20 and ∆ η 4∆ π /3 = 0.26, much larger than T 3→2 fo /m π 0 ≈ 0.05. We emphasize that the quark masses in this case roughly satisfy m u ∼ 2m d,s , namely, there is an O(1) mass splitting between the up and the down-type quarks. The results are illustrated in Fig. 3. For simplicity, in our numerical analysis we continue to employ the expression for σv 2 0 computed in the degenerate limit; the corrections to this approximation are expected to be subleading compared to the Boltzmann suppressions from the meson mass splittings. The constraints from η decay will be discussed later in Section 6.1, see Fig. 4. Incidentally, we notice that for large mass splittings the scenario where the M M → HL process is kinematically closed at T = 0 may become more plausible as the O(1) breaking of SU (3) by ∆m ∼ m lifts the η mass closer to the strong scale Λ. This in turn makes it more likely that m η > m η , resulting in c 7 > 0 in Eq. (5.6). In this case, the evolution of the relic abundances would be qualitatively similar to the one in the right panel of Fig. 2.
On the other hand, the outlook changes rapidly when 2m M > 3m L . The leading processes depleting the DM abundance in this case are LLL → M M with thermal suppression and LLM → LM . The (forward) rate for the former process is proportional to Y 3 L e −(2m M −3m L )/T , and for the latter process goes as Y 2 L in the degenerate case). It is then clear that the rates for these processes drop below the Hubble rate much earlier; depending on whether LLL → M M or LLM → LM is dominant, the new freezeout temperature can be estimated as valid for ∆ π > 0.5. In this regime, we therefore expect freezeout to occur much earlier, leading to a DM relic abundance several orders of magnitude larger than the desired value, signaling a breakdown of the SIMP mechanism as a viable method to produce thermal DM.

Thermalization of hidden and SM sectors
Viable SIMP freezeout requires the dark photon A to maintain kinetic equilibrium between the hidden and SM sectors down to T 3→2 fo m π /20. For m π GeV, this corresponds to T 3→2 fo 50 MeV, and it is a reasonable approximation to focus on dark meson scattering with electrons and neutrinos, as muons and pions are already somewhat nonrelativistic [9,10]. The cross section for, e.g., π + f → π + f scattering, where f is a Dirac SM fermion, is mediated by the exchange of neutral vector bosons V A,B = A , Z in the t-channel. Neglecting m f , we find where p f is the fermion three-momentum and we have averaged over the f spin states.
The dominant contribution comes from the exchange of A , which couples only weakly to neutrinos. Hence the total scattering rate is (5.14) where the ratio π Q 2 π /N π = 1/2 accounts for the fact that half of the N π = 8 mesons have unit charge under the dark EM and the other half are neutral (this is true regardless of whether we take Q or Q charges for the quarks). The relevant thermal average is , and g e − = g e + = 2 . The hidden sector is therefore thermalized down to freeze-out as long as 5ζ (5) 4 at T 3→2 fo [10]. This condition can be rewritten as a lower bound on ε, where the mass splittings among the mesons have been neglected. We note that the elastic scattering of the DM triplet on hidden-charged mesons is extremely efficient at T 3→2 fo , i.e. σv LM →LM n eq π + H. This remains true even if ∆ π is increased and n eq π + n π 0 at x 3→2 fo , maintaining kinetic equilibrium between the DM and the SM sector.

Suppressed semi-annihilation to dark photon and annihilation to SM
In order for 3 → 2 annihilations to drive DM freezeout, the semi-annihilation processes ππ → πA (mediated by the AAAV part of the WZW action that appears in the third line of Eq. (2.4)) should be sufficiently suppressed. This cannot be achieved by arbitrarily decreasingα, asαε 2 is bounded from below by the thermalization requirement in Eq. (5.16). Thus we restrict our analysis to m A 2m π , which, as discussed in Ref. [10], ensures the kinematic suppression of semi-annihilations, including the effect of the thermal tail, as Y π Y A Y 3 π with Y ∼ e −m/T at freezeout. In addition, we also require that dark meson annihilation to SM leptons is out of equilibrium at T 3→2 fo . The cross section for, e.g., π + π − → ff mediated by the exchange of neutral vectors ( V A,B = A , Z) in the s-channel is (5.17) where we have assumed that f is color neutral and neglected its mass. The exchange of A dominates, hence the thermally-averaged total annihilation rate is where s = 2π 2 g * s T 3 /45 and Y DM ≈ 4.1 × 10 −10 (GeV/m π ). Notice that we have neglected the annihilation to SM charged pions. Performing the average π Q 2 π /N 2 π = 1/8 and requiring Γ annihilation H at T 3→2 fo , we obtain the region above the dotted gray curve in the upper left portion of the (m A , ε) plane in Fig. 5. Here we have ignored the possible mass splittings among the mesons, which would make the charged mesons more Boltzmann suppressed at x 3→2 fo ≈ 20, mildly shifting the curve upwards. Finally, we comment that the decays of A do not affect any of the discussions in this section, as A freezes out with a very suppressed abundance and has a short lifetime. The A population decays dominantly to hidden charged mesons. Even in the narrow sliver of parameter space 2m π 0 < m A < 2m π + where this channel is forbidden by kinematics, A decays to SM particles with Γ A ∼ αc 2 w ε 2 m A , corresponding to a lifetime τ A ∼ 10 −9 s (10 −6 /ε) 2 (0.1 GeV/m A ), which is orders of magnitude too rapid to affect Big Bang nucleosynthesis (BBN).

Constraints and signatures
Our framework admits a wide variety of signals on several fronts, spanning early Universe cosmology, DM indirect detection, DM self-interactions, and dark photon searches. We now discuss these constraints and signatures in turn before summarizing them in Fig. 5.

η decays
As discussed in Section 4.1, the η decays even when the isospin symmetry is exact. If m η < 2m µ the dominant decay channel is η → 4e, whereas η → µµ dominates for heavier masses. The lifetime spans a large range of values in the allowed parameter space, but is generally much longer than one second, i.e. η decays after BBN, and shorter than the age of the Universe, so that η cannot be a significant component of the present DM.
We derive constraints on the decay of η in the early Universe based on the results presented in Ref. [13] (see also Ref. [12]), which provided limits on the energy fraction of a decaying particle as a function of its lifetime. In Fig. 4 we present bounds assuming m π 0 = 150 MeV and m π 0 /f π = 8.5, as in Fig. 3, for which only the η → 4e decay channel is relevant. We consider both anomalous and non-anomalous decay, with two different choices of the d 1 coefficient in Eq. in the (∆ η ,ŷ) plane, whose shapes reflect Fig. 5 of Ref. [13]. The parameter ∆ η controls the mass of η and therefore its freezeout abundance before it decays (see Fig. 3), as well as its lifetime to a moderate extent. The combinationŷ ≡ εê(m π 0 /m A ) 2 determines both the η lifetime and the kinetic decoupling; the latter has a mild impact on the η freezeout abundance. 10 For τ 10 12 s, constraints come from BBN observables and CMB spectral distortions, whereas CMB anisotropies provide the leading sensitivity for longer lifetimes. The bounds are strongest when the lifetime matches the timescale of recombination, τ ∼ 10 13 s, corresponding for example toŷ ∼ 10 −6 for decay mediated by the anomaly. Larger values of ∆ η lead to larger η -π 0 mass splitting and therefore to a greater suppression of the η freezeout density, relaxing such constraints. As seen in Fig. 4, the bounds completely disappear for ∆ η 0.2, corresponding to n η /n π 0 10 −11 , for which the lifetime becomes unconstrained [12,13]. For m η > 2m µ the exclusion from η decays has a different parametric dependence (see the middle panel of Fig. 5 and Section 6.6).

Dark matter decays
Next, we discuss constraints and prospects for DM decay. We mostly focus on the region m π 0 < 2m µ , where π 0 dominantly decays to 4e. As we learned in the previous sections, DM is composed of the K 0 and K 0 , which are stable, and the π 0 , which decays through its mixing with η, proportional to the small isospin-breaking parameter δ. For simplicity, in this subsection we quote lifetime bounds and projected sensitivities for a decaying species that constitutes all of DM. Our numerical results take into account that π 0 only forms a fraction r π 0 of the DM density, and as a consequence the limits are weaker by a factor 1/r π 0 . As we saw in Section 5.2, the π 0 -K 0 mass splitting can cause a depletion of the neutral kaons before the time of recombination. For realistic values of δ 10 −5 , however, the effect is mild, and r π 0 1/3 is appropriate for evaluating CMB constraints. Later, after structures form up-scattering equilibrates the three DM subcomponents, so that r π 0 = 1/3 at the present epoch.
The leading CMB sensitivity on DM decays comes from the anisotropies of the angular power spectra [12,13]. For O(100) MeV DM decaying to e + e − , the bound is τ DM 2 × 10 25 s [12], which we adopt here. A precise constraint would need to take into account the difference between our four-body decay and the two-body e + e − injection assumed in Refs. [12,13], likely resulting in an O(1) correction.
There also exist indirect detection bounds on DM decay in the present epoch. The O(100) MeV masses relevant to our setup fall in the so-called "MeV gap" of gamma astronomy, where the current constraints from gamma ray measurements are relatively weak. The leading bounds come from the COMPTEL experiment. Reference [39] used COMPTEL data to derive very conservative limits, without any background subtraction, of τ DM 10 25 s for DM decay into e + e − with the emission of final state radiation. 11 Stronger bounds of τ DM few × 10 26 s were derived in Ref. [41] from the measurement of electrons and positrons in the interstellar medium by Voyager 1. These e ± constraints come with the usual caveats about uncertainties in the modeling of cosmic ray propagation in the galaxy.
In the future, improved gamma ray limits are expected from the proposed AMEGO observatory [42]. To estimate the AMEGO sensitivity to π 0 decays we use the results of Ref. [43], which performed a detailed study of the future gamma ray reach on annihilating DM in the MeV-GeV range, using eASTROGAM [44] as benchmark. Taking into account that for m DM = O(100) MeV the ultimate AMEGO reach will be approximately a factor 2 weaker than for eASTROGAM [42,44], we apply the improvement factor for annihilating DM derived in Ref. [43] to our decaying DM, and estimate that the AMEGO sensitivity will be ∼ 50 times stronger than the current COMPTEL limit from Ref. [39], reaching τ DM 5 × 10 26 s. Thus AMEGO will become competitive with, and even surpass, the limits from Voyager 1 reported in Ref. [41], providing a strongly complementary exploration of the parameter space. Note that since the signal rate from DM decay scales as ρ DM , compared to ρ 2 DM for annihilation, uncertainties on the distribution of DM have smaller impact on indirect detection limits for decaying DM compared to annihilating DM.
We emphasize that the above indirect detection estimates are based on existing results for DM decaying (and annihilating) to e + e − , whereas the most relevant decay channel here is 4e, for which dedicated studies are not available. While a detailed treatment of this channel would be interesting, we expect that it would only induce O(1) corrections to our estimates. For example, Ref. [43] found the gamma ray constraint on a related process, DM DM → φ φ, φ → e + e − , to be about 4 times weaker than for DM DM → e + e − ; the limit on our direct four-body process would fall somewhere in between.
When π 0 is heavier than the dimuon threshold, as in the bottom panel of Fig. 5, for indirect detection bounds we make use of the results for the µµ channel presented in Refs. [41,43]. For the CMB constraints, as suggested in Ref. [13], we assume that, since neutrinos carry away an average fraction x ≈ 2/3 of the energy, the bound on τ DM→µµ is roughly given by (1 − x) times the bound on τ DM→ee [12] evaluated at (1 − x)m DM .

Dark matter self-interactions
One of the salient features of the SIMP paradigm is a large DM self-interaction cross section. The leading contributions to ππ → ππ scattering arise from the kinetic and the mass terms in the chiral Lagrangian, (6.1) where r abcd = 2(f ace f bde + f bce f ade ) and c abcd = 4 3N f (δ ab δ cd + δ ac δ bd + δ ad δ bc ) + 2 3 (d abe d cde + d ace d bde + d ade d cbe ) , with f abc the SU (N f ) structure constants and d abc its fully symmetric symbols. 12 We also defined m 2 π = 2Bm and ignored the small effect of isospin breaking. As discussed in Section 5.2, the very small mass splitting between π 0 and neutral kaon states induced by isospin breaking can cause an O(1) depletion of the latter population, but subsequent up-scattering of π 0 in halos re-equilibrates the densities of the three DM subcomponents. Therefore we assume the current DM abundance to be an equal admixture of π 0 , K 0 , and K 0 , for which the self-scattering cross section is [4] where, defining R abcd ≡ r abcd + r abdc + r bacd + r badc , we find and N π = N 2 f −1. 13 For our DM triplet, the cross section is given by Eq. (6.2) with N f = 2. The rate of self-scatterings is proportional to σ average /m DM , for which we obtain For the three benchmark scenarios presented in Fig. 5, this quantity ranges from ∼ 2 to ∼ 6 cm 2 /g, exhibiting some tension with the O(1) cm 2 /g constraints from the Bullet cluster [16] and halo shapes [17,18]. In view of the ongoing debate concerning the CDM small-scale puzzles and their possible DM explanations (see Ref. [8] for a review), we leave a conclusive statement about the viability of this cross section to the future. Incidentally, we note that the cross section in Eq. (6.4) also approximately applies to the setup of Ref. [14], where the DM is an SU (2) triplet in a 4-flavor hidden QCD theory.
If the neutral kaon and pion were split by an amount ∆ K v 2 0 , where v 0 ∼ 0.003 for the DM velocity on cluster scales, then K 0 K 0 → π 0 π 0 scattering would very efficiently deplete the kaons in the early Universe, and the current DM population would be dominated by π 0 's. These only undergo elastic π 0 π 0 → π 0 π 0 scattering, as up-scattering to heavier mesons is kinematically forbidden. Since r aaaa = 0, the cross section would be reduced to easily in agreement with constraints for our parameters. This observation was already made in Ref. [4], and the cross section in Eq. (6.5) was assumed in Ref. [11]. However, from Eq. (5.7) we find that the above scenario requires δ 10 −2 , corresponding to π 0 lifetimes that are orders of magnitude shorter than the experimental bounds. Thus, our results illustrate quantitatively that it is nontrivial to realize a setup where a single real pion constitutes all of SIMP DM and has a sufficiently long lifetime, motivating further investigation into this question.
Finally, we remind the reader that the DM self-scattering cross section decreases as the number of hidden colors is increased [4]: for fixed m π , N f and larger N c the correct relic density is obtained at larger f π , leading to a suppression σ ∝ N −4/5 c . Furthermore, a larger f π allows for heavier mesons within the perturbativity bound m π < 4πf π , which affords a further reduction of the self-scattering cross section. However, the large -N c scaling of the vector meson masses, m V ∼ 4πf π /N 1/2 c , indicates that when N c is increased the effects of these resonances become increasingly important [11], making pure chiral perturbation theory inapplicable. In addition, heavier mesons generically have shorter lifetimes, requiring a higher degree of symmetry to avoid conflict with experimental constraints.

Dark photon searches at laboratory experiments
We now discuss existing and future constraints on the A , which couples dominantly to dark quarks and therefore decays dominantly to invisible final states. These bounds and projections are collected in Fig. 5.
The strongest constraint comes from the BaBar search for e + e − → γA , A → invisible, based on 53 fb −1 of data [45] (see Refs. [46,47] for pioneering reinterpretations of preliminary BaBar results). The energy of the photon in the center-of-mass frame is given by where √ s corresponds to the Υ(2S) and Υ(3S) resonances, hence the experimental requirement E * γ 1.8 GeV implies sensitivity to m A 8 GeV. In addition, we show the expected reach attainable in the near future by Belle II with 20 fb −1 of data [48].
For small m A GeV, fixed-target experiments provide additional sensitivity. NA64 has performed a search for invisibly-decaying dark photons, exploiting bremsstrahlung production from a 100 GeV electron beam that scatters in an active dump and searching for the missing energy signature [49]. This result is based on 2.84×10 11 electrons on target, expected to increase up to 5 × 10 12 in the future [50]. Looking further ahead, the proposed LDMX experiment [51] will be able to extend the reach for light invisible dark photons by searching for the missing momentum signature. We consider two different projections for the long-term LDMX reach: an extended run beyond Phase I as discussed in Ref. [51] (see Fig. 79 therein), and the "ultimate" sensitivity attainable with a 16 GeV electron beam at the proposed eSPS facility at CERN [50]. For completeness, we also mention that the proposed KLEVER experiment could provide interesting sensitivity to lighter dark photons in the m A ∼ 0.1-0.3 GeV window [50] by searching for K L → π 0 A decays.
For m A 8 GeV the leading constraint comes from electroweak precision tests (EWPT), where the dark photon induces a shift in the Z mass and corrections to its couplings to SM fermions. We show these constraints as reported in Ref. [52]. In combination with the thermalization bound, Eq. (5.16), EWPT rule out the region m A 50 GeV. For m A m Z the main effect is the correction to the Z mass, m 2 Z m 2 Z 0 (1 + s 2 w ε 2 ) [53]. We also show the expected improvement after the completion of the LHC program, assuming in particular an 8 (440) MeV precision on m W (m t ) and a reduction of the uncertainty on ∆α (5) had (m Z ) by a factor 2, with the refined measurement of m W being largely responsible for the increased sensitivity of the fit [52].
In addition, we display the bound from a monophoton search at DELPHI, originally analyzed as a constraint on DM coupling to electrons in Ref. [54] and later recast to the dark photon case [55]. Monojet searches at hadron colliders also probe the invisiblydecaying A . In the region m A m Z the strongest bounds are still those derived from CDF data in Ref. [56]. ATLAS [57] and CMS [58] searches (based on 36 fb −1 of 13 TeV data) give comparable but weaker constraints due to much stricter selection criteria enforced in particular by trigger requirements, which imply a loss of sensitivity to the soft signal considered here [56]. In the light of this, it appears challenging for the LHC monojet searches to improve on EWPT limits, even in the high-luminosity phase. In the CMS monophoton search [59] the transverse momentum requirements are only moderately softer compared to the monojet channel, so we do not expect competitive sensitivity in our scenario. We remark that when the A decays dominantly to the dark sector, i.e. for α α ε 2 , the constraints from invisible final states are insensitive toα. Finally, searches in the + − final states (with = e or µ) at lepton and hadron colliders [60][61][62] rival EWPT for the best current sensitivity in the m A 8 GeV region, despite the strongly-suppressed branching ratio to SM particles of our dark photon. BaBar searched for e + e − → γA , A → ee, µµ in the mass range 0.02 GeV < m A < 10.2 GeV, using 514 fb −1 of data [60], whereas LHCb has recently searched for qq → A → µµ in the range 2m µ < m A < 70 GeV [61]. In our scenario these analyses are only competitive for m A 8 GeV, as the BaBar monophoton constraint is far stronger for smaller masses. In addition, CMS has recently performed a search for the dimuon signal in the 11.5 GeV < m A < 200 GeV region; for m A < 45 GeV, as relevant here, the search employs data scouting to enhance its sensitivity [62].

Dark matter direct detection and annihilation to SM
In our setup the (π 0 , K 0 , K 0 ) triplet is neutral under dark EM, whereas the charged doublet states (π ± , K ± ) have exponentially suppressed relic abundances, as required to avoid cosmological constraints from η decay. As a consequence, the signal of DM scattering on electrons via dark photon exchange [10] is beyond the foreseeable experimental reach.
For the same reasons, DM annihilation to SM leptons is extremely suppressed, avoiding any constraints from the CMB and indirect detection. Note that the annihilation cross section, Eq. (5.17), is p-wave suppressed, implying that such constraints are very weak even for degenerate mesons.

Parameter space overview
We present in Fig. 5 the current and projected constraints in the (m A , ε) parameter space based on the analysis in the previous subsections. Each of the three panels in this figure corresponds to a different representative benchmark for the meson masses and couplings.
Top panel: We consider the benchmark corresponding to Fig. 3, with Q charges for hidden EM (i.e., the AVV anomaly is present),α = 1/(4π), and isospin-breaking parameter δ = 10 −6 . The region shaded in gray is not viable for SIMP DM, as the thermalization condition Eq. (5.16) is not met ifŷ = εê(m π 0 /m A ) 2 is too small. In the upper left region, above the dotted gray curve labeled ππ → SM, direct annihilation to SM can be active at 3 → 2 freezeout, as discussed after Eq. (5.18), signaling a transition to a different regime. We recall that this curve was computed assuming degenerate mesons and would shift mildly upwards if mass splittings are included. Since m η < 2m µ the η decays to 4e, but the chosen mass splitting parameter ∆ π = 0.2 is large enough that the signal is below the Planck sensitivity on CMB anisotropies, see Fig. 4. The constraints on DM decays, i.e. π 0 → 4e with lifetime given by Eq. (4.11), are depicted as lines of constant y, corresponding to ε ∝ m 2 A . For laboratory tests of the dark photon, we show regions excluded by monophoton searches at BaBar [45] and DELPHI [55] (shaded in blue), as well as the projected near-term reach of Belle II [48] (dashed blue curve). We also display the EWPT bound (orange-shaded region) together with the improvement expected by the end of HL-LHC [52] (dashed orange), and the best current monojet constraint, based on CDF data [56] (red). For fixed-target experiments, we show the region already excluded by NA64 [49] (shaded in purple), the projected reach with 5 × 10 12 electrons on target [50] (dashed purple), and the projected sensitivity of LDMX in two possible scenarios (dotdashed purple curves), namely an extended run beyond Phase I [51] and a version with 16 GeV electron beam [50]. Notice that while the collider bounds depend weakly on m A (within the kinematically accessible range), the bremsstrahlung A production relevant to fixed-target experiments has a strong dependence on m A , resulting in the rather steep shape of the NA64 and LDMX curves. Finally, we show the current constraints derived from A → at BaBar [60] (dark yellow), as well as from A → µµ at LHCb [61] (dark red) and CMS [62] (light green). Although σ × BR ∝ ε 4 , as opposed to ε 2 for invisible final states, the strong sensitivity of the dilepton searches compensates for the suppression. For the sake of readability we do not show the BaBar and LHCb limits for m A < 8 GeV, where they are much weaker than the BaBar monophoton bound.
Middle panel: We adopt the benchmark corresponding to the left panel of Fig. 2, with the additional assumption of Q hidden electric charges so that the AVV anomaly is absent. The decays of η and π 0 are then mediated by the operator in Eq. (4.5), where we set d 1 = 0.1 andα = α. The isospin-breaking parameter is fixed to δ = 10 −5 . In contrast with the top panel, here the meson mass splittings are smaller (∆ π = 0.05), and the η density is not sufficiently suppressed to completely avoid CMB constraints (shaded in yellow). Since the η mass is (just) above the 2m µ threshold, the CMB exclusion is approximately delimited by ε ∝ m A lines, as can be read off the η → µµ lifetime expression in Eq. (4.10). On the other hand, the π 0 decays to 4e with lifetime given by Eq. (4.13). Compared to the top panel, here we consider a smaller fine-structure constant for hidden EM. This shifts the cosmological and astrophysical constraints, resulting for example in a stronger lower bound on ε from thermalization. Laboratory searches for A → invisible are essentially unaffected by this reduction inα, whereas the A → limits become stronger, benefiting from the increased branching ratio to SM particles. In particular, the BaBar and LHCb dilepton searches currently provide the leading sensitivity in the range 8 GeV m A 20 GeV, slightly outperforming EWPT.
Bottom panel: We consider a third benchmark with larger masses, m π 0 ,K 0 ,K 0 = 300 MeV, m π ± ,K ± = 330 MeV, m η = 339 MeV, (6.6) obtained by taking B = 4πf π = 343 MeV, m = 131 MeV, and ∆m = 55 MeV. The relative splittings with respect to the lightest multiplet are ∆ π = 0.10 and ∆ η 4∆ π /3 = 0.13. In addition, we assume Q hidden EM charges with d 1 = 0.5 andα = α. In comparison with the middle panel, larger mass splittings and a shorter η lifetime imply that there are no CMB constraints on η decays. The isospin-breaking parameter is set to δ = 10 −8 , significantly smaller than in the top and middle panels, because here π 0 decays dominantly to µµ with shorter lifetime, see Eq. (4.14). As a result, the bounds on DM decays have a different slope compared to the previous two panels, scaling like ε ∝ m A rather than ε ∝ m 2 A . For the chosen value of δ, current Voyager 1 constraints on decaying DM rule out the entire region m A 4 GeV. We emphasize that the values of the isospin-breaking parameter δ assumed in the three panels of Fig. 5 are simply meant to be illustrative. In all cases, changing δ → δ rescales the exclusion and projection lines for π 0 decay by ε → ε = ε (δ/ δ ) 1/2 at fixed m A .

Summary and outlook
In this paper we have studied a minimal realization of SIMP dark matter: an SU (3) hidden color gauge theory with N f = 3 light hidden quark flavors, the smallest N f that admits a Wess-Zumino-Witten action mediating 3 → 2 self-annihilation processes. An approximate isospin SU (2) U global symmetry among the down and strange quarks plays a crucial role, stabilizing the lightest hidden mesons, which form a triplet of dark matter particles. The dark photon, massive and kinetically mixed with the SM hypercharge, maintains kinetic equilibrium between the hidden and visible sectors. We summarize our novel results as follows: • We performed a detailed study of the evolution and fate of the singlet meson η, which is necessarily unstable in our setup, even if the axial-vector-vector anomaly is absent in the hidden sector. We found that η undergoes a peculiar freezeout process, driven by detailed balance between different 2 → 2 scattering processes. Its abundance and lifetime are subject to strong constraints from CMB anisotropy measurements, but these are avoided if the meson mass splittings are larger than approximately 20%.
• We studied the 3 → 2 freezeout process with mass splittings larger than the freezeout temperature, showing that the SIMP mechanism can produce the observed dark matter abundance for splittings as large as 50% of the lightest meson mass. This opens up new regions of parameter space where the CMB constraints on η decays are robustly evaded.
• We considered the possibility of decaying SIMP dark matter, as a consequence of the explicit breaking of the isospin symmetry. We analyzed indirect detection constraints in this scenario, quantifying the currently allowed amount of symmetry breaking δ 10 −5 , and discussed future prospects.
We conclude by emphasizing a few possible directions for future work. The viable parameter space, presented in Fig. 5, can be divided into two separate regions according to the mass of the invisibly-decaying dark photon. In the light region, 0.1 m A /GeV 5, existing or planned laboratory experiments such as Belle II and LDMX are set to test large swaths of parameter space in the future. The heavy region, 8 m A /GeV 50, appears more difficult to probe, motivating in particular further analysis of the LHC sensitivity through missing energy signatures.
As we have shown, dark matter decays can also provide powerful tests of the parameter space. Their relative importance can be put on firmer ground in the context of concrete ultraviolet completions, where a preferred size of δ may be predicted. As a step in this direction, we have outlined an embedding of our setup in the neutral naturalness framework, where a hidden SU (3) gauge theory with GeV-scale confinement is rather generic. We have only sketched the main guidelines, pointing out that ingredients beyond the minimal models are required, whereas the construction of a full completion and its detailed study are left for future work.
Finally, our dark matter SU (2) U triplet has a rather large self-scattering cross section, mediated by the kinetic term of the nonlinear sigma model. As is well-known, this cross section can be reduced in scenarios where the dark matter is composed of a single real meson species, whose self-interactions are then mediated by smaller explicit symmetry breaking effects. We have found, however, that this is not straightforward to realize in practice: a sufficient mass splitting of the triplet components requires values of δ that are in stark conflict with bounds on dark matter decay. We believe this aspect deserves further attention.

A Boltzmann equations
The 5-pion interaction in the last line of Eq. (2.4) can be written in the form [4] N c 240π 2 f 5 π µνρσ a < b < c < d < e T abcde π a ∂ µ π b ∂ ν π c ∂ ρ π d ∂ σ π e , (A.1) where T abcde = 60 Tr(λ a λ b λ c λ d λ e ), and explicit calculation gives  [4]. We now present the full Boltzmann equations for the three multiplets, written in terms of the yields per degree of freedom Y i (i = η, π + , π 0 ). Using x = m π 0 /T as the time variable, we have where we have defined Y + ≡ Y π + and Y 0 ≡ Y π 0 , as well as the quantities λ ≡ 4 √ 10 π 3 675 M Pl m 4 π 0 g 2 * s √ g * , κ ≡ 2 √ 10 π 15 g * s M Pl m π 0 √ g * , (A. 6) with M Pl being the reduced Planck mass, while σv 2 0 was defined in Eq. (5.2). We have assumed that all 3 → 2 annihilations are kinematically allowed at T = 0. The thermallyaveraged cross section for 2 → 2 scattering is defined as (A.7) We have taken 2m π + > m η + m π 0 , as verified for our leading-order spectrum. The equilibrium yields are defined as Y eq i (x) = Y eq (xm i /m π 0 ), where we employ the non-relativistic approximation Y eq (z) = 45 4π 4 g g * s z 2 K 2 (z) , g = 1 . (A.8) In addition, g * and g * s are also x-dependent. The above assumes that the hidden and SM sectors are in kinetic equilibrium at temperature T . If kinetic decoupling occurs at x dec (see e.g. the dotted curves in the left panel of Fig. 2 Fig. 3, which correspond to x dec = 25