Strongly interacting massive particles through the axion portal

Dark matter could be a thermal relic comprised of strongly interacting massive particles (SIMPs), where 3 → 2 interactions set the relic abundance. Such interactions generically arise in theories of chiral symmetry breaking via the Wess-Zumino-Witten term. In this work, we show that an axionlike particle can successfully maintain kinetic equilibrium between the dark matter and the visible sector, allowing the requisite entropy transfer that is crucial for SIMPs to be a cold dark matter candidate. Constraints on this scenario arise from beam dump and collider experiments, from the cosmic microwave background, and from supernovae. We find a viable parameter space when the axionlike particle is close in mass to the SIMP dark matter, with strong-scale masses of order a few hundred MeV. Many planned experiments are set to probe the parameter space in the near future.


I. INTRODUCTION
Dark matter (DM) comprises the majority of the matter budget of the Universe, but its microphysical properties and origin remain unknown. One possibility is that DM is a thermal relic from the early Universe. The most wellstudied thermal scenario is that DM is comprised of weakly interacting massive particles (WIMPs). The number density of WIMPs is set by 2 → 2 annihilations of the DM into standard model (SM) particles, and the observed DM relic abundance is achieved when both the DM mass and coupling to SM particles are near the scales relevant for electroweak processes.
An alternative thermal setup was proposed in Ref. [1] where 3 → 2 DM self-interactions set its abundance. In this scenario, the observed relic density indicates that the DM mass and self-coupling should be near the strong scale. This mechanism of strongly interacting massive particles (SIMPs) was shown to be generic in strongly coupled theories of chiral symmetry breaking, where the pions play the role of DM [2]. The 3 → 2 interactions are then sourced by the well known Wess-Zumino-Witten (WZW) action [3][4][5]. This provides a simple and calculable realization of the SIMP mechanism, although by no means the only one [6][7][8][9][10][11][12][13][14][15][16].
In addition to providing a novel thermal mechanism for explaining the dark matter abundance, SIMPs also offer a possible explanation for issues related to small-scale structure formation. In particular, observed dark matter subhalos tend to be less dense than in simulations (see Ref. [17] for a recent review). While many of these issues may be resolved with better understanding of astrophysical processes (for instance, in Ref. [18]), it is also possible to mitigate these issues if the dark matter can self-scatter (see Ref. [19] for a recent review). The strong self-annihilations of SIMP dark matter imply that their self-scatterings are also large, such that they naturally address these small-scale puzzles [1,2].
The 3 → 2 DM annihilations would raise the temperature of the residual DM due to conservation of comoving entropy. Therefore, the DM must be in thermal equilibrium with a heat sink, such as the SM bath, until after freeze-out [1]. Otherwise, the 3 → 2 DM annihilations would cause the steady depletion of DM particles and heating of the remaining DM, a scenario referred to as cannibalization. While cannibalization was originally proposed to provide a class of DM models intermediate between hot and cold DM, such models are not observationally-viable [20,21]. Obtaining the observed DM abundance inevitably leads to an unacceptable washout of small-scale structure.
To allow for adequate thermalization between the DM and the SM, Refs. [8,9] explored the kinetically mixed hidden photon portal. Here, we explore the possibility of a pseudoscalar portal using axionlike particles to accomplish the entropy transfer to photons. For brevity, we refer to axionlike particles simply as "axions" throughout the paper. We note that Ref. [22] also considered an axion portal, but focused on the regime where semiannihilations set the relic abundance. In contrast, we focus on the SIMP regime where 3 → 2 annihilations determine the relic density. For concreteness, we will use the SIMPlest pion realization of the DM based on an Spð2N c Þ gauge theory with four doublet Weyl fermions following Ref. [2]. Spð2N c Þ gauge groups with a larger number of flavors or SUðN c Þ and SOðN c Þ gauge groups allow for semiannihilations which can control the relic abundance, although there may still be parameter space where 3 → 2 annihilations determine the dark matter density.
This article is organized as follows. In Sec. II, we describe the framework and identify the interactions responsible for setting the correct DM relic abundance and cooling the DM via the axion portal. In order to cool the DM effectively, the axion must be in thermal equilibrium with both the DM and the SM. In Sec. III, we illustrate the theoretical and empirical requirements of axion-pion thermal equilibrium, while in Sec. IV we do the same for axion-SM thermal equilibrium. Concluding remarks and discussions follow in Sec. V.

II. THE FRAMEWORK
Our starting point is an Spð2N c Þ gauge theory with 2N f Weyl fermions that couple to an axionlike field a as where m a is the axion mass, m Q is the quark mass matrix, q i are the confining quarks and J is the Spð2N f Þ group invariant. 1 Upon dynamical chiral symmetry breaking, the ground state is expected to be given by Any transformation by the flavor symmetry V ∈ SUðN f Þ would also be a ground state, and in general Switching the description to the chiral Lagrangian, a spacetime-dependent flavor rotation gives the low-energy excitations, where π ≡ π b T b , T b are the Spð2N f Þ generators and f π is the pion decay constant. We use the normalization TrT b T c ¼ 2δ bc for the generators. In terms of the pion fields, The theory has an SUð2N f Þ=Spð2N f Þ flavor structure, where the residual Spð2N f Þ is exact due to the quark masses' proportionality to J. For N f ≥ 2, the fifth homotopy group of the coset space is nonvanishing and the WZW term exists [3][4][5], Generally, both aTrðπ 3 Þ and a 2 Trðπ 2 Þ terms can appear in the interaction Lagrangian. However, the former introduces semiannihilations of pions into a pion and an axion which might contribute to determining the relic abundance of the dark matter. Here, we are interested in exploring the role of an axion mediator in the SIMP mechanism of 3 → 2 self-annihilations of pions. Consequently, we focus on an Spð2N c Þ gauge theory with N f ¼ 2 fermions, where the flavor symmetry is SUð4Þ=Spð4Þ and N π ¼ ðN f − 1Þ× ð2N f þ 1Þ ¼ 5 pions emerge. In this theory, the semiannihilation process is absent since Trðπ 3 Þ ¼ 0, and pure 3 → 2 annihilations of pions via the WZW term are guaranteed to control the relic abundance of DM. For more flavors, or for other gauge groups, 3 → 2 annihilation may still control the relic abundance, though in a smaller region of parameter space. To leading order in pion fields, the WZW term for our choice of gauge group takes the form The excess kinetic energy generated in the dark sector from 3 → 2 annihilations needs to be transferred out, which can be obtained through kinetic coupling of the pions to the axions and the axions to the SM bath. Since the semiannihilation term is absent for our flavor group of choice, the interaction Lagrangian between pions and axions is If the axion coupling to the pions arises in a similar manner to what occurs in QCD, as in Eq. (1), the mass term for the hidden quarks q in the Spð2N c Þ gauge theory gives rise to an axion potential: where m π is the pion mass. Using the normalization of Trπ 2 ¼ 2π b π c δ bc , we identify the Feynman rule for the aaπ b π c vertex in Eq. (8) as Meanwhile, the interaction Lagrangian between the axions and SM photons is As long as the kinetic equilibrium between the pions and the SM is maintained through the interactions of Eqs. (9) and (11), the preferred mass for the dark matter is m π ≈ 300 MeV [2] with m π ∼ 2πf π to set the observed relic abundance. We find below that viable ALP masses are around the same scale, 10 MeV ≲ m a ≲ 1 GeV. Couplings that satisfy m π ∼ 2πf π correspond to the strongly interacting regime of the theory, where self-interactions are important on astrophysical scales. In this regime, Oð1Þ corrections to perturbative results are expected, and therefore should be thought of as a proxy for the scales involved. Phenomenologically interesting pion masses lie at the edge of perturbativity, where higher order corrections and vector meson effects can impact the range of observationallyviable pion masses [16,23].

A. Theoretical requirements
The interaction Lagrangian in Eq. (9) leads to annihilations of pions into axions and to elastic scattering of axions off of pions. The SIMP mechanism requires the former process to be suppressed at the time of 3 → 2 freezeout while the latter is active.
The requirement in order for the 3 → 2 pion selfannihilations to dominate the 2 → 2 annihilations of pions into axions at the time of freeze-out is where hσvi ann is the thermally averaged cross section for the annihilation process ππ → aa. The Hubble parameter at freeze-out is given by where T F is the freeze-out temperature of 3 → 2 interactions (typically T F ∼ m π =20 in the SIMP setup), M Pl is the Planck mass, and g Ã;F is the effective number of relativistic degrees of freedom at the time of freeze-out.
We have verified numerically that this requirement on the annihilation rate does maintain the correct relic abundance as set by the 3 → 2 SIMP mechanism. The thermally averaged annihilation cross section that appears in Eq. (12) can be readily calculated. For a trivial matrix element, M, of a process 12 → 34 [as is relevant for the Lagrangian of Eq. (9)] in which all states obey Maxwell-Boltzmann statistics, the thermally averaged cross section entering the Boltzmann equations [24] is expressed in terms of where g i counts degrees of freedom for particle i, λðx; y; zÞ ≡ ð1 − ðy þ zÞ 2 =x 2 Þð1 − ðy − zÞ 2 =x 2 Þ; ð16Þ and K 1 is the first order modified Bessel function of the second kind. The amplitude jMj 2 which appears is averaged over all degrees of freedom. For pion-axion scattering and pion annihilation to axions, the relevant amplitude is therefore since the trace requires that the pions be the same. The thermally averaged cross section for the annihilation process ππ → aa is where n eq i denotes the number density of particle i in equilibrium, where K 2 is the second order modified Bessel function of the second kind. In Eq. (18), the phase-space factor of 1=2 for identical initial particles is cancelled because the number density changes by two particles per annihilation.
For m a ≪ m π , Eq. (18) simplifies to In addition to the suppression of the 2 → 2 annihilations, the SIMP mechanism requires that the rate of energy transfer in the scattering process aπ → aπ is fast enough to successfully cool the DM. We require that thermal decoupling occurs after freeze-out, T D < T F , so that the energy transfer is efficient for the entirety of the freeze-out process. In the limit of small m a (using Bose-Einstein statistics), we follow the analytic derivation of Ref. [25] which obtains the thermal decoupling temperature T D of the πa → πa elastic scattering process The details of this derivation, including a normalization factor of 4 −1=4 Γð3=4Þ, can be found in Ref. [25]. We have checked numerically that this requirement on thermalization between the axions and pions does keep the DM cool.
In the regime where T F < m a ≲ m π , we generalize the approach laid out in Ref. [26] for particles scattering in the limit of low momentum transfer. We estimate that in the viable parameter space, this approximation is valid to within 10% even though the pion mass is larger than the axion mass by only order unity factors. Working to second order in the momentum transfer and using Maxwell-Boltzmann statistics for the axion, we find that decoupling occurs when We find that the low-and high-mass axion decoupling temperatures match onto each other when numerical differences between Bose-Einstein and Maxwell-Boltzmann statistics are taken into account in the intermediate regime. The details of the derivation can be found in Appendix B.
In addition to the above requirements, the decay constant f aπ must be greater than the cutoff scale of chiral symmetry breaking. Otherwise, the description in Eq. (1) breaks down. We require that f aπ ≳ 2πf π , where f π is determined for a given m π from the solution to the Boltzmann equation. Since the Spð2N c Þ gauge theory with N f ¼ 2 we discuss here points to the strongly interacting regime where m π ∼ 2πf π , we require that f aπ ≳ m π . In practice, however, suppressing 2 → 2 annihilations at freeze-out is always a stronger requirement.
An additional preference, though not a requirement, comes from considering how chiral symmetry breaking contributes to the axion mass in Eq. (9), The natural range for the axion mass-squared is therefore where Δm 2 a ≲ m 2 a , such that no fine tuning is present against an unspecified negative contribution, possibly from another confining gauge theory with θ ≈ π.
Satisfying the above requirements on f aπ as a function of m a for a variety of dark matter masses m π yields the viable SIMP regions depicted in Fig. 1. We take g Ã;F ¼ 10.75 at freeze-out since for the DM masses we consider, freeze-out happens below the temperature of muon-antimuon annihilation. We learn that viable SIMP-axion thermalization is achieved over a broad range of axion masses and couplings f aπ .
We note that elastically decoupling relic (ELDER) dark matter [25,27] is obtained along the thermalization curve in FIG. 1. The parameter space for axions coupling to pions. The shaded regions correspond to the regions where the SIMP mechanism is theoretically viable for a given dark matter mass m π . The boundaries of this region are set by requiring that the 3 → 2 rather than the 2 → 2 process sets the relic abundance (labeled as "WIMP regime") and that the pions can transfer their excess heat to the axions and hence the SM sector (labeled as "Thermalization"). Note that the parameter space along the "Thermalization" curve corresponds to the scenario where dark matter is an elastically decoupled relic (ELDER). We also indicate the natural mass range where the axion mass is at least as large as its contribution from chiral symmetry breaking (labeled as "Naturalness"). Also shown are the empirical upper limits on pion annihilation from energy injection into the CMB (thick solid lines labeled as "CMB limits on annihilation"). Fig. 1. For ELDER DM, the kinetic decoupling between the DM and SM baths occurs before 3 → 2 pion selfannihilations freeze out. This causes the relic abundance to be exponentially sensitive to this elastic scattering while being relatively insensitive to the strength of the 3 → 2 pion process. On the thermalization curve in Fig. 1, the elastic scattering of pions off of axions dominates over the 3 → 2 pion self-annihilation process in setting the relic abundance.

B. Empirical requirements
Having established the theoretical requirements on the axion-SIMP parameter space, we now move to the observational constraints coming from the cosmic microwave background (CMB). In standard cold dark matter cosmology, the intergalactic medium (IGM) is almost entirely neutral after recombination and CMB photons free stream. If some fraction of the DM annihilates to SM particles and partially ionizes the IGM, this will cause some CMB photons to re-scatter which modifies CMB anisotropies in a characteristic way. For the scenario we consider in this work, the process of interest is ππ → aa → 4γ. In the parameter space where there is sufficient thermalization between axions and the SM (see Sec. IV), the decay of the intermediate-state axions happens immediately. Thus, the use of the narrow-width approximation is appropriate and the cross section for this process is set by the cross section for the annihilation process ππ → aa. We use limits derived in Ref. [28], which are not very sensitive to whether there are two final-state photons or four [29]. The resulting upper limits are shown in Fig. 1 as a set of thick, solid lines corresponding to the different depicted pion masses. We thus find a viable SIMP-axion parameter space below the CMB curve and above the thermalization curve.

IV. AXIONS AND PHOTONS IN EQUILIBRIUM
A. Theoretical requirements Fig. 1 presents the viable region where the SIMP and axion maintain thermal equilibrium. If the axion and SM also maintain thermal equilibrium via the axion-photon coupling in Eq. (11), then the pions will share a temperature with the SM. For most of the axion masses we consider, decays and inverse decays into SM photons are the most efficient processes for kinetic equilibrium with the SM at freeze-out. The rate for these decays at rest is For the axions to thermalize with the SM during freeze-out, two conditions must be satisfied. First, axion decays and inverse decays must be fast enough to thermalize the axions with the SM, This is the strongest condition in the regime where the axion is relatively light and abundant. In the regime where the axion mass is comparable to or larger than that of the pion, a second condition on their decay becomes stronger than just pure thermalization between the axions and SM: This second condition requires that the axions decay quickly enough to transfer entropy that has accumulated from the pion sector to the SM. This matters more for higher axion masses m a ≳ m π since the axion number density is lower than the pion number density, so that each axion decay must transfer several pions' entropy. The detailed derivations of these conditions can be found in Appendix C. Decays and inverse decays come into equilibrium at late times. A priori, this could suggest the need for Eqs. (25) and (26) to hold prior to freeze-out in order to sufficiently transfer entropy from the annihilating pions into the SM particles. However, we verified numerically that this is not the case: the SIMP relic abundance is unaffected if decays and inverse decays into SM particles only come into equilibrium close to the time of freeze-out.
For axions with masses below the freeze-out temperature, the scattering process ae → γe is more efficient than the decays considered above. This arises because the rate of scattering is enhanced relative to decays like ∼ðT=m a Þ 3 for axions that are relativistic at the time of freeze-out, while the rate of decays is suppressed due to the boost factor of ∼ðm a =TÞ. We find that the parameter space with m a < T F is tightly constrained for the pion masses we consider and therefore do not include the scattering process ae → γe in our analytics or numerics, since it is expected to be subdominant for axion masses with m a ≳ T F . Note that by including decays as the only channel to transfer entropy, we are being conservative since adding the ae → γe channel would only lower the required coupling strength between axions and the SM.
In Fig. 2, we depict the requirement on f aγ such that decays and inverse decays sufficiently transfer entropy between the sectors. Each solid curve corresponds to the lower bound on f −1 aγ to maintain thermal contact for a fixed pion mass. We use the full Boltzmann equations and full energy transfer rates. A crossover between two regimes occurs at m a ∼ m π , where the lower axion number density starts to matter and Eq. (26) becomes a stronger condition than Eq. (25). In the regime where m a ≳ m π , kinetic equilibrium is maintained by axions in the exponential tail of the distribution, which causes the precipitous increase in f −1 aγ . As is evident, kinetic equilibrium between the axion and the SM through decays and inverse decays is possible over a range of axion masses. The conditions outlined above amount to requiring that the axions and SM have the same temperature at freeze-out; however, there is another possibility which we outline here. . Solid curves correspond to the lower bound on the decay rate for thermalization between the axions and SM for various pion masses. Vertical, dashed lines correspond to the largest axion mass allowed by CMB constraints for a given pion mass (see Fig. 1). Below the thermalization curves, the SIMP mechanism may still be viable down to the black solid line if the axion decays out of equilibrium and reheats above the neutrino decoupling temperature. Solid, filled regions correspond to existing constraints from supernova 1987a [30,31], LEP and CDF [30,32,33], BABAR [30], beam dump experiments SLAC 137, SLAC 141, CHARM, and NuCal [30,[34][35][36][37][38], and kaon decay experiments E949, NA62, NA48=2, and KTeV [39][40][41][42]. Regions enclosed by dotted, black lines correspond to projected reach by SHiP [37,43], NA62 [37], Belle II 3γ [30,42] (noting that the projected Belle II constraint from γ þ invisible falls between NuCal and NA62 [30]), BABAR [42], SeaQuest [44] and FASER [45]. See the main text for more details.
with the SM, the pion-axion sector will be slightly hotter than the SM and the relic abundance will be slighter larger for the same value of 3 → 2 rate. Therefore if the axions and SM are not thermalized at freeze-out, the value of f π must be increased slightly to give the right relic abundance. For sufficiently large m a , the universe undergoes a brief matter-dominated phase where the axions dominate the energy density of the universe. When the axions decay, they reheat the SM components, the universe becomes radiationdominated again, and the pion abundance is diluted. This happens when We require that the reheat temperature be larger than the temperature of neutrino decoupling: if the reheat temperature is lower, then the photons get preferentially heated and the effective number of relativistic neutrinos (N eff ) becomes smaller than allowed by observations of the CMB [46]. We take the neutrino decoupling temperature to be ∼3 MeV [47], at which point g Ã ¼ 10.75. Having such a high reheat temperature also enforces that the decay products do not affect Big Bang Nucleosynthesis (BBN). Therefore, the region in Fig. 2 between the solid curves and straight line labeled "Reheat before neutrino decoupling" may be viable with a slightly different 3 → 2 cross section than in the standard SIMP scenario. More dedicated numerical studies are necessary in this case and will be explored in future work.

B. Empirical requirements
Having established a theoretically-viable parameter space, we must check whether it is allowed by current experiments and observations. Constraints arise from early universe cosmology, astrophysical bodies, and terrestrial experiments.

Light Degrees of Freedom
If light axions are in thermal equilibrium with the SM bath, a bound on their mass arises from their effect on the temperature ratio T ν =T γ after neutrinos have decoupled. This difference alters the effective number of neutrino species contributing to the radiation density, N eff , which can be measured in the CMB by comparing the photon diffusion scale to the sound horizon scale [46,48]. Such constraints are relevant for light particles in equilibrium with the photon or electron plasma beneath the temperature of neutrino decoupling unless the particle couples to neutrinos as well. When applicable, this bound is stronger than the BBN bound of ∼MeV, which comes from the fact that changes to N eff change the expansion history and hence modify the abundance of the light elements. Because of this bound, only values of m a > 2.6 MeV are shown Fig. 2.

Supernova 1987a
The direct coupling of the axion to photons can lead to excess emission from supernovae (SN) via the Primakoff scattering process [49]. When the coupling between axions and the SM is sufficiently strong, the scattering of eγ → ea produces axions in the stellar medium which leads to excess cooling if the axions escape the SN. However, if the coupling is too strong, then trapping occurs via the inverse ea → eγ process along with axion decays, in which case the axion does not carry any energy out of the star. SN cooling primarily proceeds through neutrino emission; due to the observed neutrino signal from SN 1987A, any new SN cooling process must carry away less energy than the neutrinos, ∼3 × 10 53 ergs. The region of parameter space excluded by the excess cooling of SN 1987A [30] is shown in Fig. 2. For photon couplings that are too weak to produce significant energy loss in the supernova, there are still constraints from escaping axions decaying into an observable burst of photons [31], which we also show in Fig. 2.

Terrestrial
The couplings between axions and SM particles are constrained by terrestrial experiments. However, these constraints often come with assumptions about how the axions interact with the SM. We classify constraints on the axion-photon coupling based on different assumptions about its fundamental origin, namely that the photon coupling arises from: (1) solely coupling to Uð1Þ Y ; (2) solely coupling to SUð2Þ W ; (3) equal couplings to Uð1Þ Y and SUð2Þ W , in which case, the aZγ coupling vanishes. Measurements from the LEP collider and CDF constrain the decay Z → γγðγÞ [30,32,33] as shown in Fig. 2. BABAR also constrains the decay Z → γ þ invisible [30]. In the third case above, in which the aZγ coupling vanishes due to equal couplings to Uð1Þ Y and SUð2Þ W , both of these Z decay constraints are alleviated. In their place, there is a LEP bound on e þ e − → γγ [33] and a BABAR bound on e þ e − → γ þ inv [30]. Constraints from electron beam dump experiments SLAC 137, SLAC 141 [30,34,35,38], and proton beam dump experiments CHARM and NuCal [36,37] apply for axions coupled to photons regardless of how the coupling arises. Constraints from K L → π 0 a and K AE → π AE a with a → γγ assume the axion couples to SUð2Þ W . These kaon results were obtained in Ref. [42] from analyses of fixedtarget kaon rare decay experiments by the E949 [39], NA62 and NA48=2 [40], and KTeV [41] collaborations.
In addition to existing constraints, we show the projected reach of several future experiments and analyses on the photon coupling to an axion, indicated by the dashed black curves in Fig. 2. We include the projected reach of SHiP [37,43], NA62 [37], BABAR [42], Belle II [30,42], SeaQuest [44] and FASER [45]. In principle there could be a constraint from a process involving the aZZ coupling for all three scenarios, though we expect it would be weaker than the constraints we present and at this time we are not aware of any existing or projected limits from such a process.

V. DISCUSSION
In this paper, we have considered the pion realization of SIMP dark matter in strongly coupled gauge theories, and have shown that it can be realized with axions as the thermalization portal between the dark matter and SM. Throughout this work, we have required that all three sectors-the SIMPs, axions, and SM-share the same temperature as the 3 → 2 annihilations freeze out. This requirement sets a target range of masses and couplings for this mechanism to be theoretically viable.
In examining the couplings between the SIMPs and axions, we have required that the coupling is strong enough to thermalize the two sectors via 2 → 2 scattering. At the same time, we require that the coupling not be strong enough for 2 → 2 annihilations to overwhelm the 3 → 2 process that is the hallmark of the SIMP mechanism. Combined, these requirements lead to a well-defined range of couplings between the pion dark matter and the axions such that the SIMP mechanism can work. Constraints on annihilation coming from the CMB narrow the allowed range of couplings, though a broad parameter space remains. It is possible that a future CMB spectral distortions experiment can probe this parameter space further, though exploring this possibility is beyond the scope of this work.
Considering the couplings between the axions and the SM, we focused on the coupling to photons. For a given pion mass, there is a range of axion masses allowed for maintaining SIMP-axion equilibrium. Within this axion mass range, the main requirement for axion-SM thermalization is that the axions decay quickly enough to successfully transfer the entropy from the pions to the SM, which can easily be achieved. The relevant couplings to photons can be probed in a multitude of ways. The range of axion masses considered here are at an energy scale that is relevant for supernovae, which constrains part of the parameter space. Additionally, terrestrial beam dump and collider experiments have probed complementary parameter space. We find that the SIMP mechanism can be realized in a broad region of parameter space that is not excluded by current constraints. Several upcoming experiments are forecast to probe much of the viable parameter space that is currently allowed, providing an excellent handle for testing the framework.
There are several possible ways to extend the parameter space for axion-mediated SIMPs. Some of these possibilities are already excluded by existing limits. For instance, heavy axions which mediate the entropy transfer through off-shell interactions (both through the interactions we consider here and through the CP-violating interaction L ⊃ f CPV 2 aπ b π b ) are excluded by the LEP constraint shown in Fig. 2. Another heavily-constrained scenario is that the axions couple to all fermions through a universal Yukawa coupling, which is almost entirely ruled out by SLAC 137, CHARM, kaon decays, B decays, supernova 1987a, BABAR and the muon anomalous magnetic moment [50][51][52]. A more promising possibility is that axions couple only to electrons or to charged leptons; in this case, there are known limits from SLAC 137 and the muon anomalous magnetic moment [50,52]. In addition, the axion-electron parameter space should be constrained by limits from supernova 1987a and from loop-induced couplings to photons. However, such constraints have not been explored in the literature and will be the subject of future work [53]. Another possibility is that the axions have a long enough lifetime that they decay out of equilibrium, dumping entropy to the SM and diluting the SIMPs-this too will be explored in future work [53].
where the left-hand side is the relativistic Liouville operator in a Friedmann-Robertson-Walker spacetime and C½f X is the collision term. In the regime where 2 → 2 annihilations are negligible, the relevant terms which appear in the collision term are the 3 → 2 interactions which set the relic abundance, the DM self-interactions which give it a thermal distribution, the elastic interactions which transfer energy between the pions and axions, and decays (and inverse decays) of axions to SM particles. We neglect axions converting to photons (and vice versa) via t-channel scattering off electrons, which is less efficient than decays (and inverse decays) at thermalizing the axions with the SM. In the parameter space of interest (with the exception of axions in the out-of-equilibrium decay scenario), all particles will be interacting sufficiently frequently so that they have thermal distributions, where the þsign (−sign) is for Fermi-Dirac (Bose-Einstein) statistics and μ is the chemical potential. At temperatures below the mass of a given particle, the effects of quantum statistics become negligible and the Maxwell-Boltzmann (MB) distribution is recovered.
To solve the Boltzmann equations, it is most useful to look at two moments of the phase space distribution, which correspond to the number density and energy density where g X is the number of degrees of freedom and đ 3 p ≡ d 3 p=ð2πÞ 3 . The Boltzmann equations for axionmediated SIMPs are ∂n π ∂t þ 3Hn π ¼ −hσ 3→2 v 2 iðn 3 π − n 2 π n eq;T π π Þ; ðA5Þ ∂ρ π ∂t þ 3Hðρ π þ p π Þ ¼ −hσ el vδEin π n a ; ðA6Þ ∂ρ a ∂t þ 3Hðρ a þ p a Þ ¼ hσ el vΔEin π n a − m a Γðn a − n eq;T SM a Þ; where p X is the pressure densities of species X (which is related to the energy density through the equation of state w X ), and n eq;T X is the thermal equilibrium density for species X at temperature T. Additionally, hΓi T ¼ Γhm a =E a i T is the thermally averaged decay rate of the axion at temperature T, hσ 3→2 v 2 i is the thermally averaged 3 → 2 cross section of the pions for this choice of gauge group, labeled i ¼ 1…5, and n π n a hσ el vδEi is the energy transfer rate between the pions and axions (with the initial and final states labeled as 1 and 2), n π n a hσ el vΔEi For MB statistics, the equilibrium values for the number, energy, and pressure densities for a particle of mass m and temperature T with g degrees of freedom are and the thermally averaged boost factor is m E

APPENDIX B: PION-AXION KINETIC DECOUPLING
We require that the pions and axions are in kinetic (thermal) equilibrium during the entire time over which the 3 → 2 process is active. This can be recast as a requirement that the thermal decoupling temperature between the two sectors is lower than the temperature at which the 3 → 2 process freezes out. In this range of temperatures, the pions are guaranteed to be nonrelativistic since T F ∼ m π =20. Therefore, the energy transfer rate can be rewritten as The form of the collision term in the integrand will depend on whether the axions are relativistic or not at around the time of freeze-out, as detailed below.

Relativistic axions
In the regime where the axions are still relativistic at the time of freeze-out, there are well-known methods for computing the collision term [26]. For pion-axion scattering in this regime, the collision term takes the form Integrating over the pion phase space in Eq. (B1) yields n π n a hσ el vΔEi ¼ πg 2 a m 5 π 120f 4 aπ T a m π 4 ðT π − T a Þ: ðB3Þ While the 3 → 2 is actively depleting the number density, the pions are nonrelativistic and follow MB statistics, which means that their energy density Boltzmann equation can then be expressed as The first term on the right-hand side comes from the 3 → 2 and causes the pion temperature to increase, while the second term comes from pion-axion elastic scattering and pushes T π → T a . The second term cannot keep up with the first as the temperature drops and the pions and axions decouple. The temperature of decoupling is

Nonrelativistic axions
The standard result for the collision term derived in Ref. [26] applies only when the axion is relativistic. However, as long as the momentum transfer is still small, we can still apply the same methods as Ref. [26] in deriving the collision term. We will be interested in the regime where the axion mass is still smaller than the pion mass (which kinematically enforces that the momentum transferred in a single collision is relatively small) but where the axion is sufficiently heavier than the freeze-out temperature T F ∼ m π =20 such that it becomes Boltzmann suppressed. Most generally, the collision term for 2 → 2 scattering of pions and axions can be written as where J is the relevant combination of phase space factors.
In the regime of interest, everything is MB distributed at thermal decoupling so The collision term can be written as an expansion in the momentum transfer, In this expansion, C 0 vanishes simply because if the momentum transfer is zero and the number of a species does not change, the collision term is identically zero. For a contact interaction which has no angular dependence (as in the scenario we consider here, jMj 2 ∼ const), C 1 also vanishes because the angular integral contains the integrand ð ⃗ p a 2 − ⃗ p a 1 Þ · ⃗ p π ≡ ⃗ q · ⃗p π 1 , which is odd over the angular domain. Therefore, the leading-order term is C 2 . The momentum transfer scales like Δp a ∼ ðm π p a − m a p π Þ= ðm a þ m π Þ, and plugging in thermal values for the typical momentum indicates that when truncated at OððΔp a =p π Þ 2 Þ, the expansion in is accurate at the ∼10% level when the axion mass is m a ≲ m π =3.
Following Ref. [26] and plugging in the matrix element of Eq. (17), the leading order piece of the collision term is then where J 0 ≡ e −E π 2 =T π e −E a 2 =T a . This feeds into the calculation of the energy transfer rate of Eq. (B1) under the assumption that the pions are nonrelativistic, By analogy to Eq. (B4), the term in the Boltzmann equations that equalizes temperatures between the two sectors is and decoupling happens when this is order unity. The requirement that thermal decoupling happens after pion freeze-out can be recast as a requirement on f aπ Equations (B13) and (B5) do not match exactly due to a difference in numerical prefactors for BE vs MB statistics. When that relative factor π 4 =90 is taken into account, then the two match exactly.

APPENDIX C: AXION-SM THERMALIZATION
In order for the pions to maintain thermal equilibrium with the SM, two conditions must be satisfied: the decays and inverse decays of the axions need to thermalize the axions with the SM, and the axions need to lose kinetic energy via decays faster than they gain energy from the pion 3 → 2 heating. We have verified numerically that as long as these conditions are satisfied at the freeze-out temperature of the pion, the relic abundance of DM is unaffected and the pions constitute cold DM.

Axions in thermal equilibrium with photons
To understand the requirement that axions maintain thermal contact with the SM, we can ignore the pions and consider only the relevant Boltzmann equations for the axions: ∂n a ∂t þ 3Hn a ¼ −Γ a m a ðhE −1 a i T a n a − hE −1 a i T SM n eq;T SM a Þ ðC1Þ ∂ρ a ∂t þ 3Hðρ a þ P a Þ ¼ −m a Γ a ðn a − n eq;T SM a Þ ðC2Þ where the average axion energy is hE a i ¼ ρ a =n a . To make the notation less cumbersome for the remainder of this Appendix, the label for equilibrium distributions will denote chemical equilibrium and kinetic equilibrium between the axions and SM, i.e., T a ¼ T SM ≡ T. With this notation, these equations can be re-expressed as −T ∂n a ∂T þ 3n a ¼ − m a Γ a H n a hE −1 a i − hE −1 a i eq n eq a n a ≡ − m a Γ a H n a c n ðC3Þ − T ∂hE a in a ∂T þ 3hE a in a ð1 þ w a Þ ¼ − m a Γ a H n a 1 − n eq a n a ≡ − m a Γ a H n a c ρ ; ðC4Þ where w a is the equation of state of the axion, which is a function of time and axion temperature w a ¼ w a ðT a Þ. In order to eliminate ∂n a =∂T, Eqs. (C3) and (C4) can be combined to give a differential equation for hE a i: In the first term, the expansion is driving the change in the average energy, while in the second term, the decay is driving the average energy. One can check that the first term matches the expectations for a decoupled particle. Meanwhile, for further examination of the second term, we define the variable α ¼ αðT a Þ such that The value of α changes monotonically α ∈ ½1; π 6 = ð360ζð3Þ 2 Þ as T a goes from 0 to ∞. Using this definition, we find The second term vanishes when the particle is in equilibrium, i.e., n a ¼ n eq a with T a ¼ T. If the particle is driven out of equilibrium (for instance by the expansion), this term will push it back into equilibrium. In order to overcome the expansion, Γ a needs to be large enough so that the second term is larger than the first.
First we consider the case that m a ≪ T such that w a ¼ 1=3 and α ¼ π 6 =ð360ζð3Þ 2 Þ. Assuming the axion is near equilibrium and expanding Eq. (C7) around T a ¼ T to leading order gives If the two temperatures differ, then the second term will drive the system back into equilibrium if it is comparable to the first, m a Γ a TH ≳ hE a i ð2α − 3ÞT ≃ 4 m a ≪ T: Below, we find the strongest requirement on kinetic equilibrium when m a ≳ m π . For intermediate masses T ≲ m a ≲ m π , we analytically continue the condition in Eq. (C9) and require m a Γ a TH e −m a =2T ≳ 4 m a ≲ m π : We have checked numerically that this requirement ensures thermal equilibrium between the axions and SM for the entire region m a ≲ m π . This is the condition listed in Eq. (25).

Energy transfer through decays faster than from cannibalization
The second condition requires that the kinetic energy transferred to the axions from pion 3 → 2 can be compensated by axion decays. This condition will only be important when the axion is heavier than the pion and has a smaller number density. The rate of kinetic energy density loss through decays for nonrelativistic axions is Γ a n a T 2 =m a . Meanwhile, for the axions to sufficiently cool the pions, we require that hσvΔEin a n π ≳ m π _ n π ∼ Hm 2 π n π =T at freeze-out when the pions are still barely in chemical equilibrium. Therefore, the requirement is Γ a T 2 n a m a ≳ Hm 2 π n π T ðC11Þ at freeze-out. We have numerically checked that this condition keeps the pion, axion, and SM at the same temperature in the regime m a ≳ m π .