SIMPly add a dark photon

Pions of a dark sector gauge group can be strongly interacting massive particle (SIMP) dark matter, produced by the freeze-out of $3 \to 2$ interactions, with naturally large self-interactions. We study if adding a dark photon to the set-up can do it all: i) maintain thermalization with the visible sector, ii) resonantly enhance the $3\to2$ interactions, thus allowing for a perturbative pion description, and iii) provide a velocity dependent self-interaction that can affect small scale structure formation. For $N_f=3$ this minimal setup is marginally excluded, as the required kinetic mixing is too small to maintain thermal equilibrium with the SM. Adding an extra dark quark opens up parameter space, and -- perhaps somewhat surprisingly -- we find that all bounds can be satisfied for dark pion masses $m_\pi \sim 250-600\,$MeV. Dropping the small scale structure requirement iii), a viable setup is reproduced for dark charges of $\alpha_d = 0.01-1$ and a dark pion mass $m_\pi \geq 30$ MeV. Late time annihilations are non-negligible making the SIMP dark pion a bit WIMPy.


Introduction
Despite ongoing experimental and theoretical efforts, the nature of DM remains elusive. An attractive possibility is that DM is a thermal relic, and its abundance is determined by freeze-out from the thermal plasma. Although most attention has been on weakly interacting particles (WIMPs), their parameter space is increasingly constrained [1][2][3], and other explanations have come to prominence. One such alternative is strongly interactive massive particles (SIMPs) [4][5][6]. In the SIMP scenario, dark matter freeze-out occurs in a dark sector via 3 → 2 interactions, which are typically stronger than in the WIMP scenario. As a result, SIMPs have a lower mass (typically MeV-GeV scale), to which direct detection experiments are less sensitive [7]. To avoid overheating the dark sector during freeze-out, a portal coupling maintains thermal equilibrium with the visible Standard Model sector. Cosmological observations question the collisionless dark matter (CDM) paradigm. For example, observations of dark matter halo density profiles do not match the expected NFWprofile [8,9] of CDM [10][11][12][13][14][15]. This discrepancey is known as the cusp vs. core problem, see e.g. [16] for a recent review. Although the inclusion of baryonic effects in the numerical simulations may resolve the discrepancy [17][18][19][20][21], it is also possible that the resolution lies in the dark matter properties. Indeed, dark matter with strong self-interactions naturally alleviate the problem by transferring heat from the inner to the outer parts of the halo, thus smoothening the density profile [22,23]. The required interactions are scale dependent -a factor 10 difference between galaxies and galaxy clusters -pointing to a velocity-dependent self-interaction [23].
The archetypical SIMPs are the pseudo Nambu-Goldstone bosons -the dark pionsof a condensed dark Yang-Mills theory, with the Wess-Zumino-Witten term providing the five-point interactions [24][25][26]. The dark pions have naturally large self-interactions, and may address the small scale problems of CDM as well [27], except that the interactions are not velocity dependent. Moreover, satisfying both the relic density and the self-interaction constraints (or more conservatively, the upper bound on the self-interactions from the Bullet cluster observations), is only possible for non-perturbatively large pion couplings, invalidating the chiral perturbation theory approach [4]. Both of these issues can be overcome with extra vector bosons in the model [28], and in this paper we will consider adding a massive dark photon. For finetuned dark photon mass, almost twice the dark pion mass, the photon mediated self-interactions can be on resonance, thus giving rise to a velocity-dependent effect [29]. Moreover, the WZW-interactions may likewise be resonantly enhanced, increasing the freeze-out interactions [30,31].
In addition, the dark photon can maintain equilibrium with the SM through kinetic mixing [32,33] with the SM photon [5,6]. In this case dark pion annihilations to SM final states should be included in the relic density calculations as well. The annihilation rate grows at late times, as the dark pions lose kinetic energy and the annihilation cross section gets more and more resonantly enhanced. As a result, even after freeze-out of the (resonantly enhanced) WZW interactions, the annihilations can still be important and affect the final relic density. Late time annihilations in photons and electrons are bounded by nucleosynthesis and cosmic microwave background observations.
In this paper we will study if one particle -the dark photon -can do it all, solve the small scale structure problems of CDM, enhance the freeze-out interactions such that the relic density is obtained for non-perturbative pion couplings, and maintain thermal equilibrium with the SM -this is dubbed the resonant self-interacting dark matter (RSIDM) scenario. We also include collider bounds on dark photon kinetic mixing, as well as cosmological constraints from nucleosynthesis and the cosmic microwave background. We find that RSIDM can affect small scale structure formation while maintaining thermal equilibrium with the SM for the non-minimal setup with four dark quarks. The minimal setup with three quarks is only marginally excluded by the thermalisation requirement. In both setups more precise calculations may be needed to make definite statements. We will also study how parameter space opens up if the requirement that dark matter affects small scale structure formation is dropped.
This paper is structured as follows. Section 2 introduces the dark sector, with the dark pions and dark photon. This is followed by a discussion of the dark matter self-interactions and Bullet cluster bound in section 3; analytical estimates for the freeze-out temperature and final relic density, including both WZW interactions and annihilations, in section 4; and the constraints on the kinetic mixing parameter in section 5. In section 6 we then discuss the parameter space for which the the relic density can be obtained in a perturbative set-up, the self-interactions can be resonantly enhanced to address the small scale structure problems of CDM, and the dark photon can keep the dark and visible sector in thermal equilibrium during freeze-out. In addition to the analytical estimates we will also provide numerical results. We end with concluding remarks in section 7. For completeness, we have also added the computation of the various (thermally averaged) cross sections in the appendix. Our results for the WZW and pion self-interactions agree with the literature; new is the photon mediated resonant contributions to the various cross section.

Lagrangian
Strongly interacting massive particles (SIMPs) freeze out via 3 → 2 dark matter interactions [7,34]. The large required number changing interactions can naturally be obtained in a dark sector with a non-abelian symmetry, with dark pions playing the role of the DM [4]. In this paper we study the phenomenology of this set-up if we add a dark photon [5,6]. The dark photon can provide a portal between the dark and SM sectors, and -for tuned masses -can resonantly enhance the dark pion interactions.
We thus consider a dark sector with an SU (N c ) × U (1) gauge symmetry, with N f dark quarks in the fundamental representation of the gauge group. The quark mass matrix is assumed diagonal M = m q 1, allowing for a Wess-Zumino-Witten (WZW) term in the action for N c ≥ 3 colors [24][25][26]. A dimension-four kinetic mixing operator connects the dark U (1) group with the SM hypercharge [32,33], and provides a portal between the dark and visible sectors.
At a scale Λ the non-abelian gauge group condenses, and the approximate flavor symmetry of the light left-and right-handed quarks is broken down to the diagonal subgroup. The (N 2 f − 1) dark pions are the pseudo-Goldstone bosons of the this symmetry breaking. They can be naturally lighter than the scale Λ, which sets the mass of the baryons in the theory.
Depending on their couplings to other sectors, including the SM, the dark pions can be stable on cosmological timescales and thus are a good dark matter candidate.
At energies below the condensation scale the effective action is The first two terms are the leading order operators of the chiral effective Lagrangian describing the dark pion dynamics. Here U = e 2iπ/fπ , π = π a T a with T a generators of SU (N f ), and f π ∼ √ N c Λ/(4π) the pion decay constant. The covariant derivative is D µ U = ∂ µ U + ig d [Q, U ]V µ , with V µ the dark photon field and g d the dark U(1) gauge coupling. We choose the charge matrix [6,35] Q = Diag(1, −1, 1, −1, ...), (2.2) with N f entries. With this charge assignment Tr(Q 2 T a ) = 0 and the mixed anomaly vanishes, avoiding decay of the neutral dark pion into to two (dark) photons [36,37]. The next two terms are the dark photon kinetic term, where we defined the gauge field strength V µν = ∂ µ V ν −∂ ν V µ , and the Stückelberg mass term for the dark photon. To avoid pion decay the dark photon mass should exceed twice the pion mass; in the limit that the photon mass is close to that threshold, the photon-pion interactions can be resonantly enhanced. We The last term in eq. (2.1) between the square brackes couples the dark photon to the SM vector current J µ SM = q ff γ µ f , and q f the electric charge of the SM fermion f . This term arises as a consequence of kinetic mixing between the dark U (1) group and SM hypercharge; 1 Here mπ is the mass of the charged -with unit charge -dark pions, which interact with the dark photon.
The mass of the charged dark pions receive loop corrections, and is slightly larger than the mass of the neutral pions.
after redefining the fields to make the kinetic terms canonical, and diagonalizing the mass matrix, the result is the coupling in eq. (2.1) [38]. Here we used that 1 and have dropped the O( 2 ) terms, and we have neglected the coupling to the Z-boson, valid if the dark photon mass is small compared to the electroweak scale.
Finally, the last term is the WZW action, which is present if the 5th homotopy group of the coset space π 5 (G/H) is non-trivial; this is the case for N f ≥ 3 mass degenerate flavors.
Expanding the action in pion fields, the Lagrangian is the sum of the Lagrangian for chiral perturbation theory (χPT), the dark photon Lagrangian, and the terms from the WZW action: L = L χPT + L V + L WZW . The relevant dark pion and and dark photon interactions are: The first two terms in the chiral lagrangian are the kinetic and mass term for the pion fields, with mass m 2 π = 2ζf π m q . Chiral perturbation theory is perturbative for The 3rd term gives the 4pnt pion interactions, and the last term the pion-dark photon coupling (V 2π interaction). The first two terms in the dark photon Lagrangian are the kinetic and mass term for the dark photon field, and the last term the coupling of the dark photon to the SM fermions from kinetic mixing. Finally, the WZW lagrangian contains the 5pnt pion interaction, and additional dark photon-pion couplings with an odd number of pions (V 3π and 2V π interactions). We have only included the most relevant, lowest dimensional operators. Higher dimensional operators scale with the inverse cutoff scale 1/f π , which is suppressed compared to the pion mass scale for small ξ. For odd N f these terms can make the neutral pion unstable [39], which can be avoided by introducing a small mass splitting between the quarks [40]. We will instead assume that for odd N f the Wilson coefficients of the higher order contributions are sufficiently small for the dark pion to be stable on cosmological timescales, and present results for the minimal setup with N f = 3 as well as for N f = 4.

Self-interactions
The dark pion can scatter via a 4-point contact interaction appearing in the χPT Lagrangian, and via the exchange of a dark photon. The cross sections for these contributions are calculated in the non-relativistic limit in appendix B. The photon exchange contribution is subdominant, unless enhanced by an s-channel resonance which can appear for fine-tuned dark photon masses eq. (2.3). The cross section is can then be approximated by a sum of the velocity independent contact interaction and velocity dependent resonance contribution σ SI ≈ σ 4pnt SI + σ res SI . For SIMP dark matter the self-interactions are naturally large, and the Bullet cluster observations [41][42][43] put a strong constraint on the cross section. In the resonant selfinteracting dark matter (RSIDM) scenario [29,44], with a judicial choice of parameters, the self-interactions can affect structure formation on small scales in a velocity dependent way, and thus may address the putative problems of collisionless dark matter [22,[45][46][47].

Bullet cluster bound
The s-wave part of the dark pion self-interaction cross section is eq. (B.12) 2 × 10 5 cm 2 /g. The cross section is not very sensitive to the number of quark flavors N f . The Bullet cluster observation puts an absolute upper bound on the self-interaction cross section, given by a int ≈ 1. In the RSIDM scenario, where the resonant interactions from dark photon exchange become important at small scales, the data is fit by a smaller s-wave contribution and a int ≈ 0.11 (and the r.h.s. of eq. (3.1) becomes an equality) [29]. The bound on the cross section can be translated in a condition on the dark pion mass m π ≥ 14.9 MeV

Resonant self-interactions
The velocity dependent contribution to the self-interactions from nearly on-shell dark photon exchange can be written in Breit-Wigner form eq. (B.15) with S = 3S f /N 2 π the ratio of multiplicities of the resonance dark photon (3 polarizations) and DM particles (N π = N 2 f − 1) times a symmetry factor S f = 2 that takes into account that there are two identical particles in the final state of both the self-interaction process and in dark photon decay. E(v) = m π v 2 /4 is the kinetic energy in terms of the relative velocity v, and E(v R ) = m π δm the resonant kinetic energy. Further, Γ d (v) is the running decay rate of the dark photon into dark sector pions, and Γ = Γ d + Γ v the total running decay rate, which includes the decay into the visible SM sector fermions and pions via kinetic mixing. The velocity dependence of the decay widths can be parameterized as with α d = g 2 d /(4π) and α = e 2 /(4π) the dark sector and SM fine-structure constants. The resonance is peaked for v = v R , with v R = 3.6 × 10 −4 c from small scale structure data [29]. This determines δm The height of the peak is fit by m π = 4000 MeVS 1/3 (B d γ d ) 1/3 , with B d = Γ d /Γ the branching ratio for decay into the dark sector. This determines the dark photon gauge coupling: The resonant enhancement of the cross section is large, and small α d is required to avoid too large self-interactions. For RSIDM both the mass and the dark gauge coupling are given in terms of ξ, which together with the kinetic mixing parameter are the only free parameters left. The combined constraints eq. (3.2) with a int = 0.11 and eqs. (3.5) and (3.6) give

Relic density
The dark pions can freeze out via 3 → 2 number changing SIMP interactions and via annihilation into SM fermions. In the SIMP scenario thermal equilibrium with the SM sector should be maintained through freeze-out, to avoid entropy production and heating up of the dark sector. We will assume this is the case in this section, and return to the question whether the dark photon can be responsible for this in the next section. The SIMP interactions get a contribution form the 5-point coupling in the WZW Lagrangian, and from diagrams with dark photon exchange allowed by the (3π)V coupling in the WZW Lagrangian; with our choice of dark charges eq. (2.2) the π(2V ) coupling vanishes.
The thermally averaged cross section is calculated in appendix E, and is given in eqs. (E.10) and (E.27). Introducing the 'time' variable it can be written in the form The first term comes from the 5pnt pointlike interaction. The 2nd term is from dark photon exchange, which is dominated by an s-channel resonance if the dark photon is nearly on shell δm 1, and we have we used the narrow width approximation to evaluate the thermally averaged cross section. We have calculated the latter term only for N f = 3 and N f = 4 flavors. The effective couplings are Dark pions can annihilate into SM electrons (and depending on their mass, into heavier charged SM particles) via the kinetic mixing portal. Annihilation is also dominated by the s-channel resonance and the thermally averaged cross section is eq. (D.9) σv ann = α ann with as before B d = Γ d /Γ the branching ratio for decay into dark sector states. We note that both the resonant part of the WZW interactions and the resonant annihilations are independent of the mass splitting δm, except for the exponential factor, which determines below which temperature x 1/δm these interactions are 'turned off'. The Boltzmann equation for the dark pions readṡ n + 3Hn = − σv ann (n 2 − n 2 eq ) − σv 2 3→2 (n 3 − n 2 n eq ), (4.5) which in terms of the number density fraction Y = n/s becomes Here with s(m π ) = (2π 2 g * s m 3 π )/45 and H(m π ) = (π √ g * m 2 π )/(3 √ 10m pl ) the entropy density and Hubble constant at T = m π .
The relic dark matter density matches observations [48] for with Y ∞ = n/s the asymptotic number density fraction, and s 0 the entropy density today.

SIMP freeze-out
Consider first the case that freeze out of the dark pion is determined by the WZW interactions, either by the 5pnt contact interaction or the dark photon mediated contribution. This means annihilation are subdominant at freeze out σv ann σv 2 3→2 n eq at x = x f . We will estimate the bound this condition gives on in the next section. Note, however, that the annihilation cross section grows at late times, and even if negligible at freeze out, annihilation may still affect the final relic density significantly.
To describe freeze out we thus set the annihilation contribution to zero in the Boltzmann eq. (4.6). At late times the equilibrium distributions can be dropped, which allows to solve for the asymptotic distribution where we introduced the notation Y f for the asymptotic number density after freeze-out of the SIMP reactions. The freeze-out temperature can be estimated from n 2 π σv 2 3→2 H which gives where we used that the non-relativistic number density is n = N π m 3 π (2πx) −3/2 e −x . In the limit that the 5pnt interaction respectively the resonant contribution dominates the cross section we can estimate the freeze-out temperature using that x n e 2x = c gives x ≈ ln √ c − n 2 ln(ln √ c). The SIMP interactions are negligible after freeze out, but annihilations may still have an effect. To estimate this we solve the Boltzmann equation with the boundary condition where we assumed that (δm x f ) 1. The annihilation rate increases at late time, and is most efficient just before the exponential cutoff at x = 1/δm kicks in; this is how the δmdependence appears in the estimate for the relic density. It follows that annihilations are that is, only for very small kinetic mixing.

Annihilation scenario
In the opposite limit that 3 → 2 interactions are always subdominant, σv ann ≥ σv 2 3→2 n eq at freeze-out, the relic density is set by annihilation reactions only. We can still use eq. (4.11) for the relic density, but now Y f (x f ) is the number density as annihilations freeze out. The freeze-out temperature can be estimated from n π σv ann H which gives We estimate the freeze out density where we used eq. (4. 13). If x f in eq. (4.14) is smaller than that for SIMP reactions eq. (4.10), it follows freeze-out is dominated by annihilations. An earlier freeze out means a larger density Y f . Hence we can write the asymptotic number density as with the freeze out density for 3 → 2 reactions and annihilation given in eqs. (4.9) and (4.10), and eqs. (4.13) and (4.14) respectively.

Kinetic mixing
Kinetic mixing between the dark and SM photons provides a portal between the dark and visible sector. In the SIMP scenario, in which the relic density is determined by the freezeout of 3 → 2 dark pion number changing interactions, both sectors need to be in thermal equilibrium during freeze-out to avoid heating up the dark sector. The SIMP scenario further requires that dark pion annihilation into SM particles is subdominant during freeze-outalthough, as we have seen in the previous subsection, annihilation still may affect the relic density at late times. In this subsection we determine the constraints on the mixing parameter that these two requirements give. The relevant cross sections are computed in appendix D. We also quickly review the relevant cosmological and collider bounds on the kinetic mixing parameter.

Thermal equilibrium between the dark and visible sector
The dark and visible sector can be kept in thermal equilibrium via pion scattering with SM electrons and positrons. The non-relativistic cross section for this process is eq. (D.16) This can be straightforwardly generalized to include muon and SM pion scattering as well, 2 see below eq. (D. 16), in case freeze-out occurs at temperatures exceeding the muon and pion mass. Here p ≈ E e is the incoming electron momentum, and C 4 = 4 (8) for N f = 3 (4) the same color factor as appearing in the dark photon decay rate eq. (3.4). The scattering rate can be estimated as Γ scat ≈ n e E 2 e (σv) scat /E 2 e [5], and demanding that it exceeds the Hubble rate Γ scat > H at the time of freeze out gives the bound If the dark photon is to maintain thermal equilibrium with the SM, annihilations cannot be neglected for the relic density calculation if (comparing eq. 4.12 and eq. 5.2) The last equality applies to RSIDM, for which annihilations thus always play a role, depleting the dark matter abundance at late times.

Annihilations subdominant during freeze-out
Annihilations are subdominant at freeze-out if σv ann < σv 2 3→2 n π at x = x f eq. (4.10), which translates to α 3→2 (1 + α res x f . This gives an upper bound on the kinetic mixing parameter where we set κ 3→2 , B d to unity.

CMB bound and other bounds on kinetic mixing
BBN and CMB observations place bounds on the energy injected in the photon fluid from late time dark pion annihilations to SM final states. These constraints can be stringent as the resonant contribution to the thermally averaged cross section grows as σv ∝ x 3/2 e −δmx . At late times x > δm −1 , when the (average) dark matter velocity falls below the resonance velocity, resonant annihilations into SM final states are exponentially suppressed. Due to the momentum dependent coupling of the dark pions to the dark photon the remaining nonresonant contribution is p-wave suppressed, and evades all CMB bounds. For self-interacting resonant DM the mass splitting δm ∼ 3×10 −8 is small eq. (3.5), and the exponential suppression of the cross section only kicks in shortly before or after the CMB is formed, depending on the mass of the dark pions. The thermally averaged cross section thus peaks at this time. CMB observations are more stringent for s-wave annihilations (as opposed to more stringent BBN bounds for p-wave annihilations) [49]. Since our thermally averaged cross sections scales inversely proportional to the velocity, CMB bounds are stronger than bounds from BBN, and we thus only consider the former. The CMB bound is given by and f (z) = 0.01 − 1 is a function that quantifies the efficiency of energy injection in the CMB [48]. Recasting this to a bound on the mixing parameter , the constraint is given by (setting which vanishes above dark pion masses of 150 MeV. This bound was derived for s-wave scattering. In our model, the annihilation cross section increases for later times until the exponential cut-off kicks in, so the bound is expected to be stronger at late times. At the same time, the energy injection into the CMB is maximized at z ∼ 600, where f (z) ≈ 1 [50,51]. Although applying the CMB bound naively for our case at the different choices of z affects the exact constraint on the mixing parameter somewhat, this has no consequences on parameter space of the SIMP scenarios discussed in the next section, as this is dominated by the constraints from thermalisation and beam dump experiments. For simplicity, then, we imposed the CMB bound at z ∼ 600, and that is what is shown in our figures in section 6. For larger δm only the BBN bound applies. Following [49], energy injection during BBN is most efficient in the range 1/T ∼ 10 2 − 10 3 MeV −1 . Requiring the annihilations to be suppressed at this time (δmx 1), the BBN bound can evaded for mass splittings δm 10 −3 .
Finally, there are bounds from dark photon searches at beam dump or fixed target experiments at electron or proton colliders. In these experiments large number of dark photons can be produced from Bremsstrahlung or secondary meson decays. The experiments typically search for highly displaced vertices in the detector [52,53]. For our region of interest, 10 −4 10 −8 and 10 MeV m γ 1 GeV we consider bounds from NuCal [54,55], CHARM [56], and E137 [57]. Given the p-wave nature of our scattering cross section, scattering at late times is heavily suppressed. Bounds on millicharged particles from direct detection experiments like XENON are therefore too weak to constrain the kinetic mixing parameters and are therefore not considered.

SIMP scenarios
In the 'standard' SIMP scenario the dark pion relic density results from the freeze out of the 5pnt 3 → 2 WZW processes. We reproduce this set-up by turning off the dark photon interactions α d = 0. The Bullet cluster observation puts an upper bound on the pion selfinteractions, and consequently on the pion mass eq. (3.2). The correct relic density is obtained for ξ ∼ 4π for N f = N c = 3, uncomfortably close to the perturbativity bound ξ 4π [4]. Increasing the number of colors improves the situation slightly, but a large number of colors is needed to be within the perturbative regime.
The problem with only the 3 → 2 interactions, and no resonance enhancement, is that dark matter is over produced. The dark pions contribute a fraction to the DM density eq. (4.8) The freeze-out temperature only depends logarithmicly on the model parameters, and ranges from x f 20 = 0.68 − 1.1 for ξ = 1 − 10. The relic density is reproduced for R = 1, which requires large ξ 4π. This is illustrated in fig. 1, which shows the numerical solution to the Boltzmann equation without (red) and with (blue) adding the 3 → 2 resonance, where α d is chosen from eq. (3.7).
The dashed lines correspond to the estimate of eq. (4.9) with α res = 0 (α res = 0) for the red (blue) curve. The shaded areas are excluded by the perturbativity cutoff on ξ and the Bullet cluster bound on the self-interaction.
Without the dark photon resonance, the required value of ξ is at or above the perturbativity cutoff. Including the resonance, but not considering annihilations, the situation improves significantly as the increased 3 → 2 interactions reduce the dark matter density. The observed relic density is obtained for a larger dark pion mass for a fixed ξ. For such a large pion mass, however, the kinetic mixing parameter is too small to maintain kinetic equilibrium with the SM eq. (5.3) and either an additional portal interaction is required or annihilations should be considered.
In the following subsections we discuss two SIMP scenarios that give the correct relic density, satisfy self-interaction constraints, and the kinetic mixing portal interaction maintains thermal equilibrium with the SM during freeze-out. First, we focus on the possibility that the dark pions are RSIDM. The 3 → 2 freeze-out interactions can be resonantly enhanced for large enough α d /ξ 4 10 −6 . Given the small mass difference δm in eq. (3.5), annihilations become important at late times, and reduce the relic density further. Second, we consider the more classical SIMP scenario, and drop the requirement that the self-interactions can affect small scale structure formation. The self-interactions should still satisfy the upper bound from the bullet cluster. Both 3 → 2 interactions and annihilations may be resonantly enhanced, by tuning the dark photon mass, but now have more freedom in the resonant condition δm and the dark gauge coupling α d to satisfy all constraints.

Self-interacting resonant SIMP DM
Consider first the RSIDM scenario that the relic density is produced via the SIMP mechanism, i.e. via the freeze out of 3 → 2 reactions, and that resonant DM self-interactions can address the small scale structure problems. The requirements on the self-interactions eq. (3.7) fixes the parameters m π , α d , δm in terms of ξ, which in turn is determined from the correct relic density eqs. (4.8) to (4.10).
For large enough kinetic mixing the relic density will be reduced by (late-time) annihilations, which reduces the required ξ value. Assuming annihilations are subdominant during freeze-out and B d ≈ 1, then Y ∞ is given by eq. (4.11) and R now becomes where R res ≡ 1 + 30x  Figure 2 shows the value of as a function of the dark pion mass m π for which R = 1 from numerically solving the Boltzmann equations for N f = 3 (orange) and N f = 4 (red), as well as the analytical estimate eq. (6.3) above (dashed curve). In addition, the bounds from the CMB, colliders, and from the requirement of thermal equilibrium between the dark and visible sector during SIMP freeze out eq. (5.2) are shown.
The observed relic density can be obtained with perturbative couplings, e.g. ξ = 1 for kinetic mixing = 3 × 10 −8 . Note, however, that there is a maximum value ≤ 2.9 × 10 −7 (6.4) to get the observed relic density, at a dark pion mass m π ∼ 600 − 650 MeV. For larger dark pion masses, the 3 → 2 interactions alone underproduce DM given the imposed relations on the self-interaction eq. All constraints can be satisfied and the RSIDM scenario can be viable for N f = 4 and dark pion masses in the range m π 250 − 600 MeV, while for N f = 3 this is only marginally possible around m π ∼ 500 MeV. In our analysis we have analytically estimated the bound on the required kinetic mixing with the SM. To make more precise statements requires inclusion of the bath effects, which also allows for (partial) heating of the dark bath, in our numerical calculations. This is left for future research.

SIMP DM
We now consider the SIMP scenario, in which the relic density is determined by 3 → 2 interactions and possibly additional annihilations, but self-interactions are too weak to affect small scale structure formation. As we have seen eq. (6.1), with just the 5pnt WZW interactions and given the Bullet cluster bound, too much DM is produced for perturbative couplings ξ. The DM density can be reduced by a resonant enhancement of the WZW interactions and by annihilations. No longer constrained by the small scale structure data, the value of δm can now be larger. This immediately avoids CMB and BBN constraints as the thermally averaged annihilation cross section ∝ e −δm x eq. (4.4) is exponentially suppressed in these eras. Moreover, for larger δm dark photons will predominantly decay to dark pions, thus evading dark photon searches at beam dump experiments. For concreteness we will fix the mass splitting to δm = 10 −3 throughout this subsection. This choice avoids the cosmological constraints, while it can still give rise to interesting phenomenology at late times.
In the parameter space region where the dark photon maintains thermal equilibrium with the SM sector during freeze-out of the WZW-interactions, the kinetic mixing parameter is always large enough that late time -after freeze-out -annihilations cannot be neglected. For general m π , ξ, δm and α d the required mixing parameter that reproduces the correct relic density is with Y f from eq. (4.9). Parameter space opens up significantly compared to the RSIDM scenario, as all other parameters are free except for the Bullet cluster bound on the dark matter self-interactions. Figure 3 shows the numerical results for the kinetic mixing parameter as a function of the dark pion mass that reproduces the observed relic density for N f = 3 (left) and N f = 4 (right). The red, orange and purple curves correspond to different values of ξ = 2, 5, 10. In all plots the mass splitting is set to δm = 10 −3 . The numerical results are in excellent agreement with the analytical estimate eq. (6.5). Along the curves three different regions can be identified; i) a part where the self-interactions are larger than allowed by the Bullet cluster constraint, which is excluded (dotted); ii) a part where the 3 → 2 interactions dominate freeze-out, and annihilations are only important at later times times (solid); and iii) annihilations are the dominant freeze-out process (dashed).
The blue region in the plots is excluded by the thermalisation requirement eq. (5.2). For larger α d a smaller kinetic mixing parameter is required to maintain thermal equilibrium with the SM. The results are very similar for N f = 3 and N f = 4, although the kinetic mixing (bullet cluster) constraints are slightly weaker (stronger) for N f = 4.
The top plots shows the result for α d = 0.01. For such small gauge couplings, the resonance enhancement of the WZW interaction is negligible. DM is over produced unless annihilations are important. In fact, we see that for the parameter space allowed by the Bullet cluster constraints, the annihilations are actually so large that they always dominate freeze out. Hence, in this scenario the dark pions are WIMP rather than SIMP dark matter.
In the middle plots the gauge coupling is increased α d = 0.1, but still resonance effects on the WZW interactions are small except for large ξ. Indeed, for ξ ∼ 10 the dark pion can be SIMP DM. Although late time annihilations reduce the relic density somewhat -this is what generates the slope of the solid curve as a function of mixing parameter -the effect is not strong enough to allow for much smaller ξ than in the 'standard' scenario.
Finally, the bottom plots are for sizeable couplings α d , and the WZW interactions are significantly enhanced for smaller ξ = 1 − 5 as well, allowing SIMP dark matter in a perturbative set-up. The slope of the solid part of the curves shows the impact of late time annihilations as a function of kinetic mixing. The curves asymptote to a constant value for small mixing and annihilations are negligible at all times, thus providing a lower bound on the dark pion mass for a given ξ.

Conclusion
We have studied a dark sector containing dark pions and a dark photon. The dark pions are stable and can be SIMP dark matter, that is produced by freeze-out of 3 → 2 WZWinteractions. The dark photon mixes kinetically with the SM sector, and can maintain thermal equilibrium during freeze out. For a fine-tuned dark photon mass m V ≈ 2m π the WZW are resonantly enhanced, which opens up the possibility that the observed relic density is produced in the perturbative regime ξ = m π /f π 4π of the effective chiral Lagrangian. In addition, the pion self-interactions are resonantly enhanced and become velocity dependent, which opens up the possibility that it can address the small scale structure problems of collisionless dark matter -this scenario is dubbed resonant self-interaction dark matter RSIDM.
We found that the RSIDM scenario is possible for dark pion masses in the range m π ∼ 250 − 600 MeV for N f = 4 dark quark flavours. For N f = 3 flavours the RSIDM scenario is (marginally) excluded for all dark pions masses considered, because of the highly efficient dark pion annihilations. In this case the value of the mixing parameter required to reproduce the relic density is too small to maintain kinetic equilibrium with the SM. For both setups more precise calculations are required to assess the viability of the model in this region. In particular, one could allow for (partial) heating of the dark bath to study the effect on the dark bath. This is left for future research.
If we give up the demand that the self-interactions have an effect on small scale structure formation, and consequently are only constraint by an upper bound from observations of the Bullet cluster, then parameter space opens up and smaller ξ-values become possible for sufficiently dark gauge couplings α d ∼ 1.

A Cross sections
We list here the definitions of the (thermally averaged) cross sections. This appendix also serves to set the notation.

A.1 Scattering/annihilation cross section
The cross section for scattering with two particles in both the initial and final state, labeled by α and β respectively, is . We use P µ for 4-momenta, and p i for 3-momenta. The second expression is valid in the center of mass frame (CM), with p in = |p α | (p out = |p β |) the absolute value of the three-momentum of either incoming (outgoing) particle. S β = N ! for N identical particles in the final state, to avoid overcounting in the phase space integral. |M| 2 is the amplitude averaged over initial and summed over final state particles. The flux factor can be written as the center of mass energy squared.

A.2 Thermally averaged cross section
The thermally averaged cross sections can be defined in terms of the scattering rates appearing in the Boltzmann equation for the DM particle [34] n + 3Hn = − α,β ∆ αβ (γ α→β −γ β→α ) = − σv ann (n 2 − n 2 eq ) − σv 2 3→2 (n 3 − n 2 n eq ), (A.2) with n = n DM the number density of dark matter. We have included both DM annihilation and 3 → 2 interactions. ∆ αβ = (N dm α − N dm β ) is the difference between the number of DM particles in the initial (N dm α ) and final state (N dm β ). The collision terms arê with f α , N α the distribution functions and degrees of freedom of the initial states, and S α (S β ) = N ! for N identical particles in the initial (final) state. Assuming kinematic equilibrium for the DM and chemical equilibrium for all other particles gives the relationŝ with γ α→β ≡γ α→β (f eq α ) and f eq α = e −Eα/T the Maxwell-Boltzmann distribution. We can then express the thermally averaged cross sections appearing in the Boltzmann equation eq. (A.2) in terms of the collision rates as follows For annihilations the momentum integrations in γ ann can be partially done, and the final expression is given in terms of one remaining integral over the center of mass energy [58,59]. The thermally averaged cross section is with the Møller velocity related to the flux factor as F = (E 1 E 2 )v møl , and the factor ∆ ann /S α = 1 is set to unity. The equilibrium number density is defined as with the last expression valid in the non-relativistic limit. For m 1 = m 2 ≡ m we can rewrite this in dimensionless variables

A.3 S-channel resonance and narrow width approximation
If the DM interactions are mediated by a massive meditor particle -in our case, the dark photon -there will be an s-channel resonance for momenta that the mediator is nearly on shell s ≈ m 2 V . To incorporate this effect, we include the decay rate in the dark photon propagator where we used that Γ 2 m 2 V . Γ = Γ d + Γ v is the total decay width of the dark photon, which is the sum of the decay rate into pions and decay rate into SM fermions, i.e. into the dark and visible sector. In the resonance limit the most enhanced terms in the cross section will be ∝ D µν (s) 2 , which can be evaluated in the narrow width approximation

B Pion self-interactions
In this appendix we calculate the pion self-interaction cross section σ SI = σ(ππ → ππ), which has contributions from 4pnt self-interactions and from dark photon exchange.

B.1 Amplitude
Dark photon mediated self-interaction The 4pnt pion interaction follows from eq. (2.4) There are 4! possible orderings of the pions a, b, c.d. Consider first the amplitude for the {acbd}-term plus the cyclic permutations: where we took P a , P b as incoming momenta, and P c , P d as outgoing (∂π a → −iP a , ∂π c → iP c ), and on the 2nd line we used the Mandelstam variables The results are similar for the other possible permutations, and the total amplitude is [27] M 4pnt Dark photon mediated self-interaction The pion-dark photon interactions that follow from the covariant derivatives in the chiral Lagangian eq. (2.4) can be written in the form The (V 2π)-vertex interaction and dark photon propagator (in Lorentz gauge) are then with both momenta P a , P b incoming. The amplitude for ππ → ππ scattering via dark photon exchange is with color factor Resonant scattering arises for P 2 ≈ m 2 V , and we include the decay width in the propagator eq. (A.9) to describe this.

B.2 Cross section
The matrix element squared summed over final and averaged over initial states is |M| 2 SI = 1 N 2 π abcd |M abcd | 2 with N π = (N 2 f − 1) the number of pions. The amplitude is the sum of the 4pnt interaction eq. (B.4) and the photon exchange contribution eq. (B.7). The latter is subdominant, except for momenta near the s-channel resonance s ≈ m 2 V ; we can then neglect interference terms and the non-resonant contributions, and approximate the amplitude the s-channel resonance contribution from the photon exchange diagram.
4pnt self-interaction Starting with the 4-pnt contribution, the amplitude squared can be calculated using the Feyncalc Mathematica program [60][61][62] The cross section eq. (A.1) for pion scattering is with S f = 2 for two identical particles in the final state, and p = |p| the incoming momentum of either pion in the CM frame. The last expression is valid in the non-relativistic limit p 2 m 2 π in which the velocity-independent s-wave contribution dominates, and s ≈ 4m 2 π . The result agrees with [28], but differs a factor of 8 with [4] (after rescaling f them π = 2f π ). For N f = 2 it reproduces the results of [27] (except for the sign in front of the 5 6 p 4 term) but not for other N f .
Dark photon mediated self-interaction The diagram for photon exchange is negligble except near the s-channel resonance, where the amplitude can be approximated by the schannel diagrams, and ) with S f = 2 amd Γ(p) the running photon decay width, computed in appendix C. The color factor eq. (B.8) summed over flavors is C 2 4 = |C abcd | 2 = 4 2 (8 2 ) for N f = 3 (4). The cross section can be rewritten in the more familiar Breit-Wigner form. To do so, expand the momentum as p = 1 with v the relative velocity. The resonance velocity follows from eq. (C.4), which gives v R ≈ 2 √ δm for dark photon masses eq. (2.3). We further define the velocity dependent and resonant kinetic energies [29] (B.14) We can then rewrite the resonant cross section in Breit-Wigner form where we used the non-relativistic approximation s = m 2 V + O(v 2 ). The numerator is written in terms of the decay width into dark pions Γ d using the explicit expression eq. (C.3).
Ref. [29] list a numerical factor S = N V /N 2 π for the ratio of mulitplicities of the resonance dark photon and DM particles; we find here an additional S f = 2 symmetry factor.

B.3 Thermally averaged cross section
The thermally averaged self-interaction cross section is [29] σv SI = σ 4pnt with the velocity distribution f (v, v 0 ) = (4v 2 )/( √ πv 3 0 )e −v 2 /v 2 0 taken as a Maxwell-Boltzmann distribution truncated at the halo escape velocity v max , and v 0 a constant implicitly defined via v 2v 0 / √ π. Here we used that the self-interaction cross section is the sum of the (dominantly) s-wave 4pnt interaction and s-channel resonance contribution eqs. (B.12) and (B.15). The resonance contribution can be calculated in the narrow width approximation eq. (B.14): with B d = Γ d /Γ the branching ratio for decay into dark sector pions. The decay rate into dark pions (p-wave) eq. (C.3) and SM fermions (s-wave) eq. (C.7) can be parameterized as for decay into dark and visible sectors, and m V ≈ 2m π . This factors out the explicit velocity dependence of the (Γ d B d )-factor in the thermally averaged cross section.

C Dark photon decay rate
The dark photon can decay into dark pions, and in SM fermions and pions via kinetic mixing with the SM photon.
The dark photon decay rate in dark sector particles is Γ d = Γ(V → ππ). The decay into dark pions is mediated by the vertex interaction eq. (B.6), with amplitude with P a , P b the 4-momenta of the outgoing pions. The amplitude squared (summed over final states, and averaged over intial photon polarizations states) in the CM frame is with p out = |p a | = |p b | the final state 3-momentum, and color factor C 2 2 = ab |C ab | 2 = C 4 . The decay rate becomes with α d = g 2 d /(4π) and S f = 2 for two identical particles in the final state. The resonance momentum (taking the intial state photon at rest) is where we used the parameterization eq. (2.3) in the last step, and assumed the dark photon mass is close to resonance δm 2 m 2 π . Evaluating the decay rate at resonance p out = p R gives C.2 Dark photon decay rate in SM particles.
The dark photon can decay into SM electrons via kinetic mixing and Γ s = Γ(V →f f ) with f the electron; for large enough DM mass the decay into muons and charged SM pions also becomes kinematically accessible. The decay rate into SM pions will be of the form eq. (C.3) with α d → 2 α and instead of C 4 → C QCD 4 = 1/2 the appropiate color factor for QCD. It is velocity suppressed compared to decay into fermions, which is s-wave, and we neglect it in the following.
The amplitude for dark photon decay into a SM fermion is with P 1 , P 2 the 4-momenta of the outgoing fermions, the kinetic mixing paramer appearing in the Lagrangian eq. (2.4), and e the electric charge of the electron. The amplitude squared (summed over final states, and averaged over photon polarizations) in the CM frame and decay rate are with α = e 2 /(4π), symmetry factor S f = 1, and s = m 2 V . Decay into the darks sector pions dominates and Γ d Γ v or where we set v → v R .

D Dark pion annilation and scattering
Dark pions can annihilate into and scatter off SM fermions. These processes are important for the final relic density and for keeping the dark sector in thermal equilibrium with the SM sector.

D.1 Dark pion annihiliation into SM fermions
Consider the annihiliation of dark pions into Standard Model electrons π a (P 1 )π b (P 2 ) → f (K 1 )f (K 2 ). For large enough dark pion mass, also the decay channel into muons (m π > 105 MeV) and SM charged pions (m π > 139.6 MeV) opens up. The amplitude for annihilation in a SM Dirac fermion pairf f is mediated by the vertex eq. (B.6) and the coupling to the SM fermion current via kinetic mixing.
Averaging over initial states and summing over the spins of the final state fermions gives with color factor C 4 = C 2 2 ≡ ab |Tr([T a , T b ]Q| 2 . Here we denoted the CM 3-momentum of the incoming pions with p and that of the outgoing fermions with k, and θ the scattering angle. The cross section eq. (A.1) becomes In the last line we neglected the SM fermion mass, and rewrote the cross section in terms of the dimensionless variables = s/(4m 2 π ). The amplitude for annihilation in a pair of SM pions π ± SM with is Averaging over intitial states and summing over the spins of the final state pions -c, d running over the generators corresponding to π ± SM -gives with color factor C QCD 4 = 1/2. The cross section can then be expressed as where we neglected the electron mass. Substituting eq. (D.3) in the thermally averaged cross section eq. (A.8) for annihilation into SM electrons is with ∆ ann /S α = 1. We are interested in the thermally averaged cross section during freezeout, in the limit x = m π /T 1. This will be dominated by the resonance contribution which we can calculate in the narrow width approximation eq. (A.10). Including the decay width in the propagator eq. (A.9) and definingΓ = Γ/(2m π ),m V = m V /(2m π ) the thermally averaged cross section eq. (D.7) becomes where in the last step we used that K 2 (x) = K 1 (x) = π/(2x)e −x + ... for x 1, and expanded (m 2 V − 1) ≈ δm. Now expand the dark photon mass in small δm defined in eq. (2.3), and use the explicit expression for the decay width eq. (C.5) to get with B d = Γ d /Γ the branching ratio for decay into the dark sector. We can then write the full thermally averaged annihilation cross section as with g ann incorporating the degrees of freedom the dark pion can annihilate in. Below the muon threshold this is only electrons and g ann = 1. Including muons and SM pions with Θ the Heaviside step function.

D.2 Pion-electron scattering
Consider scattering of dark pions with SM particles, which for light DM is dominated by electron scattering via the reaction π a (P 1 )f (K 1 ) → π a (P 2 )f (K f ).
Averaging over intitial and summing over final states gives The spinor sum now becomes (...) = 4 2(P 1 + P 2 ).K 1 (P 1 + P 2 ).K 2 − (P 1 + P 2 ) 2 (K 1 .K 2 − m 2 f ) ≈ 16m 2 π p 2 (1 + cos θ) + O(p 2 ) (D.14) with p the CM 3-momenta of the incoming particles and k of the outgoing particles, and θ the scattering angle between the ingoing and outgoing pion. In the last step we set the fermion mass to zero and took the non-relativistic limit. The non-relativistic cross section becomes [5] 15) with A scat = 2πC 4 α D α/N π and m V ≈ 2m π . The total scattering cross section is written as , where for simplicity we have neglected the muon and SM pion masses.

E Cross section for 3 → 2 dark pion interactions
The 3 → 2 dark pion amplitude has a contribution from the 5pnt pion interaction and from dark photon exchange. The necessary vertices with an odd number of pions arise from the WZW Lagrangian eq. (2.4).

E.1 Pion 5pnt interaction from the WZW-term
The 5pnt pion interaction from the WZW term can be written as [4] L WZW ⊃ A 5 f 5 π µνρσ Tr [π∂ µ π∂ ν π∂ ρ π∂ σ π] , with N c the number of colors of the dark QCD-like gauge group. There are 5! different contributions to the amplitude M abc→de . We can group them by their momentum dependence with f a coefficients that each are the sum of 4! terms, and all momenta are taken as incoming.
The color-coefficient f e is defined as That is, f e is the sum of all traces of five generators with T e fixed in the first position, and all possible permutations of the other generators (of the {a, b, c, d} indices); the sign of each term is determined by the number of permutations away from the {a, b, c, d, e}-sequence. Likewise, all other f i coefficients can be defined. The matrix element squared averaged over initial and summed over final states is (the calculation is done with the Mathematica package FeynCalc) with N π = N 2 f − 1 the number of pions and N f the number of flavors. F is a complicated momentum-dependent function, which we will evaluate for non-relativistic incoming momenta. We parameterize the momenta in the CM frame (we set P d → −P d and P e → −P e to make them outgoing momenta with positive energy as the zeroth component of the 4-vector): P a = (E 1 , p 1 , 0, 0), P b = (E 2 , p 2 cos θ, p 2 sin θ cos ϕ, p 2 sin θ sin ϕ), P c = (E 3 , −p 1 − p 2 cos θ, −p 2 sin θ cos ϕ, −p 2 sin θ sin ϕ), , p 4 cosθ, p 4 sinθ cosφ, p 4 sinθ sinφ), The momentum p 1 is aligned with the z-axis. Ω 2 and Ω 4 are then the solid angle of p 2 and p 4 respectively. p 3 and p 5 are fixed in center of mass frame, and E 4 , E 5 are fixed by energy conservation. The energies are E i = p 2 i + m 2 π for i = 1, 2, E 2 3 = (p 2 1 + 2p 1 p 2 cos θ + p 2 2 ) + m 2 π , and p 2 4 + m 2 π = 1/2(E 1 + E 2 + E 3 ). We can thus express E i , p 4 in terms of p 1 , p 2 . To get the leading order term in the limit of non-relativistic momentum for the incoming particles set p i = p i for i = 1, 2 and expand in small . The first non-zero term arises at 4th order: F = 4 F 4 (p 1 , p 2 , θ, ϕ,θ,φ) + O( 6 ) with The transition amplitude eqs. (A.3) and (A.4) then becomes 1 N π f eq 1 dp 2 p 4 2 N π f eq 2 (E. 8) with S α = 3! and S β = 2! to account for identical particles in initial and final states. The integral dΩdΩF 4 = 750π 2 m 4 π p 2 1 p 2 2 . To get the final expression, we further used that in the non-relativistic limit E i ≈ m π for i = 1, 2, 3 and E i ≈ 3 2 m π for i = 4, 5, which yields The integrations over the phase space densities give in the non-rel limit N π d 3 p 3 (2π) 3 f eq 3 = n eq π , N π dp 1 p 4 1 f eq 1 = 6π 2 m π T n eq π (E.9) The thermally averaged cross section eq. (A.5) is then with x = m π /T and ξ = m π /f π . This result matches the result in Ref [4] (taking into account the different definitions f them π = 2f π ).

E.2 Dark photon interactions from WZW term
In the presence of a dark photon, there are additional diagrams with photon exchange contributing to the 3 → 2 cross section. The WZW term also contains photon interactions with an odd number of pions eq. (2.4). The A(3π) and (2A)π interactions are with T Qabc = Tr QT a T b T c and T Qa = Tr Q 2 T a . Including the 5pnt pion interaction discussed in the previous subsection and the A(2π)-interaciton from the chiral Lagrangian, the relevant photon-pion couplings vertices are and all momenta are taken as incoming. For the A abc vertex we used that there are 3! different terms, corresponding to abc + cycl.. We have symmetrized the A a vertex in the two photon legs. Q 2 = 1 For our choice of charge matrix eq. (2.2), and thus T Qa = Tr(Q 2 T a ) = Tr(T a ) = 0, and the (2A)π interaction vanishes A µρ a = 0. The propagator is where in the last step we introduced the dimensionless propagator∆ = ∆/(4m 2 π ).

E.2.1 Amplitude
The photon interactions give rise to 6 additional diagrams contributing to the 3 → 2 interactions [28]. For our charge matrix eq. (2.2) the A µρ a = 0 vertex vanishes, and only the first three diagrams contribute. We will calculate the amplitude for N c = N f = 3. The full amplitude is iM = A abcde + a 1 −iA µ abc A deµ ∆ (de) + a 2 P abc −iA µ ab A cdeµ ∆ (ab) + a 3 P abc;de −iA µ ce A abdµ ∆ (ce) (E. 15) with ∆ (ij) = ∆(P i + P j ) the propagator of momenta P i and P j . P abc means all cyclic permutations of (abc) -hence there are three diagrams contributing to a 2 -, and P abc;de cyclic permutations of (abc) times cyclic permutations of (de) -hence there are six diagrams contributing to a 3 . The a i are included as note keeping devices and can be set to unity at any time during the calculation.The amplitude in terms of theĀ-vertices becomes Diagram a 2 has an s-channel resonance for m V ≈ 2m π , as the propagators∆ ij with i, j = a, b, c go nearly on shell. Diagram a 1 can become resonant for m V ≥ 3m π , and the propagator ∆ de can be put on shell. This case was analysed in [63]. We will here not consider it any further, and instead focus on lighter dark photon masses.