Cosmology and Accelerator Tests of Strongly Interacting Dark Matter

A natural possibility for dark matter is that it is composed of the stable pions of a QCD-like hidden sector. Existing literature largely assumes that pion self-interactions alone control the early universe cosmology. We point out that processes involving vector mesons typically dominate the physics of dark matter freeze-out and significantly widen the viable mass range for these models. The vector mesons also give rise to striking signals at accelerators. For example, in most of the cosmologically favored parameter space, the vector mesons are naturally long-lived and produce Standard Model particles in their decays. Electron and proton beam fixed-target experiments such as HPS, SeaQuest, and LDMX can exploit these signals to explore much of the viable parameter space. We also comment on dark matter decay inherent in a large class of previously considered models and explain how to ensure dark matter stability.


I. INTRODUCTION
Despite the wealth of gravitational evidence for dark matter (DM), its nature remains a fundamental puzzle in particle physics. If DM possesses non-negligible interactions with the inflaton or Standard Model (SM) bath, a large thermal population is necessarily generated in the early universe. This abundance must be depleted to match the observed DM energy density. In the Weakly Interacting Massive Particle (WIMP) paradigm, this occurs through annihilations to SM particles. In general, however, once the temperature drops below the DM mass, any process that reduces the DM number density can potentially lead to a viable freeze-out scenario.
A simple and generic possibility is that DM is the lightest state of a strongly coupled hidden sector (HS) that is uncharged under SM forces. In this case, the dynamics of the HS opens up qualitatively new mechanisms for depleting the DM abundance. The best-known possibility is where annihilations of three DM particles into two DM particles deplete the thermal abundance [1], as realized in models of Strongly Interacting Massive Particles (SIMPs) [2,3]. Nonetheless, to be cosmologically viable, the HS must couple to the SM in order to maintain kinetic equilibrium, though the required coupling is significantly weaker than would be needed to achieve the DM abundance through pair-annihilation.
We focus on a class of strongly interacting hidden sectors with QCD-like particle content. Dark matter is composed of the lightest states in this sector, π D -the pseudo-Nambu-Goldstone bosons of a spontaneously broken chiral symmetry, similar to the pions of QCD [3][4][5][6][7][8][9]. Kinetic equilibrium between the dark sector and the SM can be maintained through interactions via a light mediator, such as a dark photon, A , the gauge boson of a broken U (1) D symmetry that kinetically mixes with SM hypercharge. The strongly interacting sector necessarily includes heavier excitations, such as vector mesons, V D , the analogues of SM ρ mesons, which are expected to mix with the dark photon.
These vector mesons, V D , have important consequences for both the cosmology and phenomenology of SIMPs. Past studies have largely ignored the impact of vector mesons on SIMP cosmology. This is justified in the chiral limit, where V D 's are parametrically heavier than DM pions. However, SIMP models rely on significant chiral symmetry breaking to achieve sufficiently large 3π D → 2π D cross-sections, so that one expects m V D ∼ m π D . In this case, V D → π D π D decays are kinematically forbidden, and the semi-annihilation [10] reaction, π D π D → π D V D followed by V D → SM, is often the primary process responsible for depleting the DM abundance. Accounting for this process consistently dramatically expands the viable range for the DM mass and the coupling to the SM. We find that pion masses ranging from ∼ 10 MeV to well above the GeV-scale are allowed, in contrast to the narrow mass window around ∼ 100 MeV when only 3 → 2 reactions are considered [2,3]. This process also allows the coupling between the HS and SM to be orders of magnitude smaller than previously thought.
A visibly decaying vector meson also opens up significant new avenues for accelerator-based dark sector searches. Some of these have been discussed in Refs. [11][12][13]. We focus here on the visible signals associated with dark photon decays to the strongly-coupled sector. A striking signature of a strongly interacting HS is the decay of a dark photon A → π D V D with subsequent resonant decay V D → + − (where is a SM lepton), as shown in the bottom-left diagram of Fig. 1. For natural parameter values, there is an observably displaced V D -decay vertex. This signal differs from the direct decay of a dark photon (or direct production of the V D ) in the energy spectrum of the decay products (see Fig. 6). More dramatically, the production rate for V D 's is parametrically larger than that of an A for a given lifetime. . Similar processes can also occur for non-singlet vector mesons which undergo a three-body decay (through an off-shell A ) into a hidden sector pion and a pair of Standard Model leptons (bottom-right). The inset shows a schematic hidden sector mass spectrum, with m A /2 mV D ∼ mπ D , which enables these decays.
A related process is the decay of the dark photon into vector mesons whose quantum numbers do not permit mixing with the dark photon, as shown in the bottomright diagram of Fig. 1. These vector mesons decay to π D + − final states with even longer lifetimes. These distinctive signatures can be searched for at beam dump and fixed-target experiments. Such searches are complementary to the minimal signals of HS DM, e.g., nuclear/electron recoils and invisible dark photon decays, the latter of which is shown in the top diagram of Fig. 1. Data from the E137 beam dump experiment is already able to probe interesting regions of parameter space, especially for ∼ 100 meter decay lengths. Complementary viable regions will be tested in the near future at the currently running Heavy Photon Search (HPS) experiment, an upgrade of the SeaQuest experiment, and at the proposed Light Dark Matter eXperiment (LDMX). Our main results are summarized in Fig. 5, where we show existing constraints as well as sensitivity of HPS, SeaQuest, and LDMX to cosmologically-motivated models that have not been tested otherwise. Similar signals are also observable above the muon threshold at the Bfactories BaBar and Belle-II and at the Large Hadron Collider (LHC). This paper is organized as follows. In Sec. II, we describe a benchmark model of a strongly interacting HS that we use throughout this work. We also show that HS vector mesons are long-lived for well-motivated parameter values and, therefore, can give rise to displaced vertex signals at fixed-target and collider experiments. In Sec. III, we discuss the cosmological importance of these vector mesons and clarify the issue of pion stability. We then demonstrate in Sec. IV that existing and future fixed-target, collider, and direct detection experiments are sensitive to cosmologically-motivated parameter space. We also briefly comment on various astrophysical and cosmological probes. Finally, we summarize our conclusions in Sec. V. Details of the model, cross-sections and decay rates, and Boltzmann equations are provided in Appendices A-C.

II. A STRONGLY INTERACTING SECTOR
We consider a strongly interacting HS described by a confining SU (N c ) gauge theory with N c = 3 colors, analogous to SM QCD. We also introduce N f light flavors of Dirac fermions in the fundamental representation. We are interested in the relative importance of 3π D → 2π D and π D π D → π D V D in dictating the DM abundance. We choose N f = 3, as this is the minimum number of flavors that is required to allow either process. In this section, we briefly outline the basics of the model, while a more detailed discussion is provided in Appendix A. Hereafter, we denote the HS pions and vector mesons as π and V , respectively (a subscript "D" is implied). For π and V , the superscripts, 0 and ±, denote charges under U (1) D , while for , they denote charges under U (1) em . The global chiral symmetry, SU (N f ) L × SU (N f ) R , is spontaneously broken by the hidden quark condensate to the diagonal subgroup, SU (N f ) V , during confinement. Thus, at low energies this is a theory of N 2 f − 1 pions, π, which constitute the DM of the universe. The low-energy pion selfinteractions are described by chiral perturbation theory; the strength of these interactions is characterized by the HS pion decay constant, f π . 1 The odd intrinsic-parity Wess-Zumino-Witten term of the chiral Lagrangian enables 3 → 2 pion annihilations [14,15]. If such processes are the dominant pion number-changing mechanism during thermal freeze-out, the HS pion abundance matches the observed DM energy density for m π /f π few [3]. On the other hand, observations of merging galaxy clusters lead to upper limits on the DM self-scattering crosssection, at the level of σ scatter /m π few×cm 2 /g [16][17][18]. As shown in previous studies, these considerations suggest that m π ∼ few × 100 MeV and m π /f π ∼ 4π [3].
Vector mesons, V , are expected to appear in the spectrum at a scale close to the cutoff of the theory, m V ∼ 4π f π / √ N c [19][20][21]. Hence, for large values of m π /f π , the DM pions are parametrically close in mass to the lightest spin-1 resonances of the HS. In particular, for the 3 → 2 cosmologically-motivated values, m π /f π 4, the vector mesons are expected to lie near the two pion threshold, m V ∼ 2 m π . If m V < 2 m π , the decays V → ππ are kinematically forbidden and instead V must decay to SM states. As a result, vector mesons play an important role in the cosmological depletion of DM pions.
Considerations of structure formation in the early universe demand that kinetic equilibrium between the HS and the SM bath must be maintained during freezeout [2]. A light mediator that couples weakly to the SM can allow for efficient heat transport between the two sectors. The dark photon, A , is an ideal candidate for this task since it interacts with the SM through a renormalizable operator; thus, if its interactions are not too feeble, it is able to maintain kinetic equilibrium between the two sectors even at low freeze-out temperatures relevant for light SIMPs (see also Ref. [22] for a recent investigation of HS-SM equilibration). The A is a gauge boson of a broken U (1) D symmetry that couples to the SM through kinetic mixing with the hypercharge gauge boson, B, where A µν and B µν are the U (1) D and U (1) Y field strengths, respectively, and θ W is the Weinberg angle [23,24]. The kinetic mixing parameter, , can be generated by loops of heavy particles charged under both U (1) D and U (1) Y . Its natural size is therefore 10 −2 ( 10 −4 ) if it is generated at one (two) loop(s) [24,25]. For m A 100 GeV, the A dominantly mixes with the SM photon. After electroweak symmetry breaking, SM fermions acquire a millicharge under U (1) D , given by e Q f , where e is the electromagnetic coupling and Q f is the electric charge of the SM fermion in units of e. Dark photon masses below the electroweak scale naturally arise in, e.g., supersymmetric theories in which the HS is sequestered from the source of supersymmetry 1 Note that our convention for fπ differs from that of Refs. [3,12] by a factor of two.
breaking [26]. In addition, m A ∼ × O(100) GeV is independently motivated in theoretical constructions where electroweak symmetry breaking triggers U (1) D breaking at the GeV-scale [27]. With A at the GeV-scale, different mass hierarchies for the other HS states are possible. From a bottom-up perspective, m A m π is required to prevent observable distortions of the cosmic microwave background (CMB) through late-time DM annihilations of the form ππ → A A followed by A → 2 [28]. Hence, we will focus on m π ∼ m V m A ∼ GeV. Interactions between the A and HS mesons are incorporated by embedding U (1) D in the unbroken SU (N f ) V subgroup of the global chiral symmetry. This embedding is specified by the HS quark U (1) D charge assignments, defined by the charge matrix, Q. A generic choice of Q enables decays of neutral HS pions into SM leptons via tree-level (π → A * A * → 4 ) and loop-induced (π → 2 ) processes. The leading contribution to these processes in chiral perturbation theory is controlled by the chiral anomaly, which is proportional to Tr(T π Q 2 ), where T π is the SU (N f ) generator corresponding to π. If the neutral pions are unstable during DM freeze-out at t ∼ H(T fo ) −1 , charged pions can scatter into these states, efficiently depleting the entire DM abundance to a negligible value [12]. Therefore, the viability of this DM scenario requires the neutral pions to be long-lived on timescales compared to freeze-out. The leading contribution to neutral pion decay vanishes if Q 2 ∝ 1 since Tr T π = 0. Following Refs. [12,29], we take It has been claimed that this charge assignment is sufficient to make all DM pions stable [12,29]. However, the anomalous term is only the leading contribution to the decay; higher dimensional operators in chiral perturbation theory also contribute to this process [30]. 2 Indeed, even in the absence of the electromagnetic chiral anomaly, the SM π 0 SM would be unstable but with a much longer lifetime [31]. Therefore, for general quark mass assignments, only the charged HS pions are stable since they are the lightest states protected by the global U (1) D symmetry. In Sec. III, we argue that these states can be viable DM candidates even if the neutral pions are unstable.
Hidden sector mesons can be created in terrestrial experiments through their interactions with the A . Gauge symmetry determines the interactions of A with the charged mesons, while the neutral vector mesons mix with the A through interactions of the form [12,32,33] where V µν is the vector meson field strength, e D is the U (1) D gauge coupling, and g is the V ππ interaction 2 We thank Markus Luty and Jack Gunion for emphasizing this point to us.
strength. We relate g to the vector meson mass and pion decay constant using the KSRF relation 3 [35,36], Anomalous decays, such as A → πV , are described by the gauged Wess-Zumino-Witten interactions [37]. An on-shell A can be radiated in an electron-or proton-nucleus collision as depicted in Fig. 1. It can also be produced through SM meson decays and Drell-Yan at proton beam fixed-target experiments. The observable signals are determined by its decays. For m A > 2 m π and 1, the A decays almost exclusively into HS mesons, where α D ≡ e 2 D /4π. Since the choice of Q in Eq. (2) is structurally similar to that of the SM, we adopt the standard meson naming scheme (see also Appendix A). 4 We present the branching ratios of the A to various HS final states in the upper panel of Fig. 2, as a function of m π /f π for m A /m π = 3, m V /m π = 1.8, and 1. For m π /f π 3, the anomalous decays, A → V π, dominate the total A width, while for m π /f π 3, the A has a significant invisible branching fraction into pairs of charged pions (and charged vector mesons if m A > 2 m V ).
The vector mesons are unstable. If V mixes with the A , i.e., Tr(Q T V ) = 0 (see Eq. (3)), and m V < 2 m π , V decays directly to the SM. The decay rate into a pair of SM leptons is of the form The corresponding decay length can easily be macroscopic since Γ is suppressed by the small parameters and α D /g 2 . The spin-1 U (1) D -neutral states ρ and φ decay to the SM in this manner, while ω − A mixing vanishes at leading order. As a result, the ω and charged vector mesons decay via three-body processes through an offshell A , e.g., ω → π A * → π + − and V ± → π ± A * → π ± + − . In the limit that m m π m V m A , these decay rates take the parametric form 3 In the Hidden Local Symmetry framework of vector meson-pion interactions, the KSRF relation follows from the Vector Meson Dominance limit [34], where the A couples to π only through mixing with an intermediate V .  and are further suppressed by three-body phase space factors. The corresponding decay lengths are typically much larger than those of the two-body decays in Eq. (6). 5 The decay lengths for two-and three-body decays are shown in the lower panel of Fig. 2 for α D = 10 −2 , = 10 −3 , and m π /f π = 3 as a function of the vector meson mass. We note that, even in the two-body case, V has a macroscopic decay length for m V GeV. Full expressions for the decays considered above are given in Appendix B.
As shown in the top panel of Fig. 2, HS vector mesons are produced with an O(1) branching fraction in A decays for m π /f π few. These vector mesons are parametrically long-lived compared to the A if m V < 2 m π . Thus, if such vectors are produced, their decays to the SM are naturally displaced, giving rise to unique experimental signatures. Higher-order corrections to meson masses from chiral symmetry breaking effects can give rise to various spectra where only a subset of the vector meson masses lie below 2 m π . For instance, if all vector mesons are lighter than 2 m π , then two-and three-body visible decays are present in the theory. The enhanced decay length from three-body processes favors high-current, long-baseline beam dump experiments for moderate val-ues of . 6 A qualitatively different scenario arises if the ω and all charged vector mesons are heavier than 2 m π . Then, only two-body decays of the remaining neutral vector mesons can give rise to visible signatures, with a decay length better suited for short-and medium-baseline experiments. We show in Sec. IV and in Fig. 5 that both possibilities can be tested with existing data, future runs of HPS and SeaQuest, and the proposed LDMX experiment. Alternatively, when m V > 2 m π , all vector mesons decay invisibly into pairs of DM pions. In this case, promising signals include DM-electron/nucleon scattering and missing energy signatures at accelerators [13].

III. COSMOLOGY
In this section, we discuss in more detail two aspects of SIMP cosmology mentioned in Secs. I and II. We previously noted that the singlet pions are unstable even if Q 2 ∝ 1 for general quark masses, M q . If these are the lightest of the DM pions, then their decays can efficiently deplete the entire DM energy density. Even if they are long-lived compared to the timescale of freeze-out, but with a lifetime greater than a second, late decays may be in conflict with the successful predictions of Big Bang nucleosynthesis or measurements of the CMB.
However, these pitfalls can be avoided. For example, if in addition to Q 2 ∝ 1, we enforce Tr Q = 0 and M q ∝ 1, then the singlet pions can be made absolutely stable through an enhanced symmetry. These requirements allow for the construction of a G-parity, where C is the U (1) D charge conjugation operator, Z A 2 sends A → −A , and U q ∈ SU (N f ) V is a unitary isospin transformation of the HS quarks such that U † q QU q = −Q T . Note that such a U q exists only if the numbers of positively-and negatively-charged HS quarks are equal since a unitary transformation cannot alter the determinant or trace of Q. The singlet pion (defined in Eq. (4.8) of Ref. [12]) is odd under G. Hence, the above choices of Q and M q forbid the construction of an operator relevant for π → A ( * ) A ( * ) → 4 or π → 2 . This can easily be seen from the fact that any candidate amplitude must be built out of traces of at least two powers of Q, the generator corresponding to the pion, T π , and possibly M q .
If M q ∝ 1, these amplitudes can always be simplified to The situation is different if N f is odd or if the numbers of positively-and negatively-charged HS quarks are 6 These three-body processes can also be used to test SO(Nc) and Sp(Nc) models, where the low-energy spectrum does not contain neutral pions, aside from the flavor singlet meson associated with the chiral anomaly [12]. We leave a detailed study of these models to future work. not equal. This is relevant for N f = 3 flavors as considered in this work. In this case, even if Q 2 ∝ 1, there is no analogous G-parity, and there exist chirally-invariant operators that facilitate the decays of singlet pions, e.g., where U = exp(2iπ a T a /f π ), T a are the generators of SU (N f ), and the overall coefficient is estimated using naive dimensional analysis [21]. At the parton level, traces of Q 2n−1 (with integer n) arise from quark loops with an odd number of A insertions. These diagrams have been claimed to vanish due to Furry's theorem [12]. However, this is not the case when an arbitrary number of HS gluons is attached to the loop. An example of such a parton-level diagram is shown in Fig. 3. The operator in Eq. (9) induces two-and four-body decays of the unstable singlet pions into SM fermions. Both partial widths are suppressed by α 2 D 4 and by loop or phase space factors, leading to extremely long lifetimes. We estimate that the four-body channel dominates with a partial width of while the two-body channel is additionally suppressed by m /m π . For α D = 10 −2 , 10 −4 , and GeV-scale masses, the lifetime of the singlet pions is comparable to the timescale of recombination, ∼ 10 13 seconds. Measurements of the CMB anisotropy spectrum are a powerful probe of such an energy injection and, for certain lifetimes, limit the fraction of DM that can undergo visible decays to 10 −11 [41,42]. Decays before recombination

FIG. 4:
The hidden sector pion abundance in the mπ − plane. The contours correspond to regions where the π abundance matches the observed dark matter energy density for mπ/fπ = 3, αD = 10 −2 , m A /mπ = 3, and various values of the mass ratio, mV /mπ. Multiple processes determine the dark matter abundance. For 10 −3 , the abundance is governed by WIMP-like annihilations into Standard Model particles, ππ → SM. For ∼ 10 −4 − 10 −3 , efficient semi-annihilations, ππ → πV , dictate freeze-out. If 10 −5 − 10 −4 , the ability of semi-annihilations to reduce the π density becomes limited by the vector meson decay rate, Γ(V → SM) ∼ 2 . If 10 −6 , then the relic abundance is determined by SIMP-like 3 → 2 reactions. For simplicity, we have assumed that the hidden and visible sector baths remain in kinetic equilibrium throughout freeze-out. This assumption will be relaxed in Sec. IV. Although this is valid for 10 −6 , we expect that kinetic decoupling will occur before dark matter freeze-out for 10 −6 . This case is analogous to models in which the abundance is governed by kinetic decoupling [39,40]. may be constrained by CMB spectral distortions [42]. Thus, unless the abundance of the unstable component is negligible, the CMB places an important constraint on a wide swath of otherwise viable parameter space. We emphasize that this is not a generic feature of strongly interacting dark sectors since the HS pions can be made absolutely stable, as discussed above. Models where some of the pions are unstable may be cosmologically viable, provided all unstable species are heavier than the lightest stable pion and decay at temperatures below the corresponding mass splitting. In this case, 2 → 2 scattering in the HS depletes the abundance of the unstable pions before they decay. CMB constraints on this possibility are discussed in Sec. IV I. The relevant mass splittings can arise from chiral symmetry breaking corrections from HS quark masses, through higher-order operators in the chiral Lagrangian such as [31] where B 0 is a dimensionful constant related to the HS quark condensate and defined in Eq. (A3). These operators, with coefficients α 6,7 ∼ O(10 −3 ), can lift the unstable pions to the degree needed for cosmology if the values of α 6,7 are chosen judiciously [31,43]. However, obtaining the cosmologically viable spectrum requires α 7 > 0 (of opposite sign from the expected η -mixing contribution) and is about 2σ away from SM measurements [31]. It is therefore unlikely (but not impossible) that this mechanism operates in SM-like theories. It is nonetheless plausible that this spectrum may be realized in less SM-like hidden sectors without G-parity.
We now turn to the problem of DM freeze-out in the early universe. The first examples of SIMPs were designed to realize the 3 → 2 annihilation mechanism [2]. This process arises from a dimension-9 operator in the Wess-Zumino-Witten term. The corresponding rate in the early universe is given by where x ≡ m π /T SM , T SM is the temperature of the SM bath, and n eq π is the equilibrium number density of π. Due to the strong exponential suppression in Eq. (12), the correct relic abundance is obtained for large values of m π /f π close to the perturbativity bound, i.e., m π /f π ∼ 4π [3]. This observation, combined with the effective field theory expectation for vector meson masses [19][20][21], suggests that the vector mesons play an important role in the DM cosmology. Indeed, the Wess-Zumino-Witten term gives rise to the pion semi-annihilation (2 → 1) process, ππ → πV , with V decaying to the SM. This of-fers an alternative mechanism to deplete the pion number density. 7 We calculate the relic abundance of DM pions by numerically solving the full system of Boltzmann equations, incorporating annihilations, 3 → 2 and 2 → 1 processes, and vector meson decays. A more detailed discussion is presented in Appendix C and will also appear in forthcoming work. In Fig. 4, we show regions in the m π − plane where the π abundance matches the observed DM energy density for α D = 10 −2 , m π /f π = 3, m A /m π = 3, and various values of m V /m π , assuming throughout that the HS and SM are in kinetic equilibrium. These parameter choices are motivated below.
The behavior of the relic abundance contours in Fig. 4 has a simple interpretation in terms of the dominant processes controlling freeze-out. The favored DM mass range is strongly dependent on . At large , direct annihilations into SM particles, ππ → A * → SM, determine the π relic abundance, favoring WIMP-like regions of parameter space. For smaller values of , semiannihilation (ππ → πV , V → SM) dictates the cosmological history. If the vector mesons decay to SM particles more rapidly than they downscatter into HS pions, i.e., Γ(V → SM) n eq π σv(πV → ππ), then semi-annihilation is independent of and is largely driven by the strength of ππ → πV alone. DM freezes out in this case when where H is the Hubble parameter. Alternatively, when Γ(V → SM) n eq π σv(πV → ππ), semi-annihilation is limited by the rate of vector meson decays, which depends on the strength of kinetic mixing, as in Eq. (6). Hence, freeze-out occurs when and the ability of semi-annihilation to deplete the pion number density is quenched by 1. Note that Eqs. (14) and (15) imply that semi-annihilation is dependent on n eq V /n eq π and is exponentially sensitive to the mass ratio, m V /m π . For even smaller values of , vector meson decay rates are extremely suppressed and semiannihilation cannot efficiently deplete the pion abundance. At this point, the only remaining viable mechanism is the 3 → 2 process, πππ → ππ. This region of parameter space corresponds to standard SIMP models [2].
In generating the contours of Fig. 4, we have assumed that DM-SM elastic scattering and vector meson decays maintain kinetic equilibrium during freeze-out. This assumption breaks down in Fig. 4 for 10 −6 , in which case the HS and SM kinetically decouple before 3 → 2 reactions freeze out. This scenario is analogous to elastically decoupling relics (ELDERs) of Refs. [39,40] but 7 In the SM, the same interaction is responsible for ω → 3π [37].
with a different process (vector meson decays to the SM) driving the kinetic equilibration of the two sectors. We emphasize that for most of the viable parameter space, semi-annihilations govern the cosmology of DM freezeout, while regions in which 3 → 2 processes or direct annihilations are important are mostly ruled out by existing constraints. This will be explored in more detail in Sec. IV.
In the semi-annihilation region of Fig. 4, the correct relic abundance is obtained for much larger values of m π compared to the pure 3 → 2 regime. Equivalently, smaller values of m π /f π are feasible for fixed DM mass. This significantly alleviates tensions with perturbativity and bounds on DM self-interactions. 8 We emphasize that the cosmological importance of processes involving vector mesons is a natural consequence of large values of If m V is fixed according to Eq. (13), then the correct relic abundance can be obtained for a wide range of DM masses without conflicting with perturbativity or limits on DM self-interactions. For instance, if is large enough for semi-annihilations to be efficient (vector mesons rapidly decay to the SM) then Eq. (13) implies that GeV-scale DM pions acquire an adequate relic density for 9 The mild dependence of m π /f π on the DM mass motivates m π /f π 3 as a good representative choice for models of thermal DM.
Previous literature has focused on maintaining kinetic equilibrium through DM-SM elastic scattering via A exchange [12]. Since such processes are controlled by the strength of kinetic mixing, demanding that the DM and SM sectors remain in kinetic equilibrium during freezeout places a lower limit on the size of . The presence of neutral vector mesons provides another means of maintaining kinetic equilibrium between the two sectors. The decays and inverse-decays, V ↔ SM, efficiently equilibrate the SM and vector mesons, while ππ ↔ πV and πV ↔ πV enforce equilibrium within the dark sector. For m V 2 m π and m π /f π 4π, equilibration between the HS and SM bath is dominantly governed by vector meson decays, which allows for significantly lower values of . As mentioned above, for values of near this lower bound, kinetic decoupling of the HS from the SM occurs before 3 → 2 freeze-out, leading to ELDER-like cosmology, albeit at a much smaller [39,40].

IV. EXPERIMENTAL SIGNALS
A rich experimental program is currently underway to search for light mediator particles, such as the minimal dark photon, with beam dump, fixed-target, and collider experiments [54,55]. These searches can also explore the strongly interacting models discussed in this paper, which give rise to new signals from HS vector meson production and decay, as shown in Fig. 1. In fact, existing data significantly constrains SIMP-like models. However, much of the cosmologically motivated parameter space remains out of reach. In this section, we show that searches for displaced visible decays of HS vector mesons can be used to explore the best motivated regions of masses and couplings. In Secs. IV A-IV C, we focus on existing data from past beam dump experiments and the currently running HPS and SeaQuest experiments. In Sec. IV D, we discuss signatures at the proposed LDMX detector. We also briefly comment on possible searches at colliders such as BaBar, Belle-II, and the LHC in Sec. IV E, the prospects for electron recoil signals at upcoming low-threshold detectors in Sec. IV F, and bounds from galaxy clusters, supernova cooling, and the CMB in Secs. IV G-IV I. Finally, in Sec. IV J, we give a detailed discussion of our main results. These are summarized in Fig. 5, where we show existing constraints and future projections as gray shaded regions and colored contours, respectively.

A. Past beam dumps
The SLAC E137 beam dump experiment searched for long-lived particles produced in e − -Al collisions. It utilized a 20 GeV electron beam and ∼ 30 C of current, equivalent to ∼ 10 20 electrons on target (EOT) and ∼ 100 ab −1 of integrated luminosity [46]. An electromagnetic calorimeter (ECAL) was placed ∼ 400 meters downstream of the target after ∼ 200 meters of natural shielding provided by a dirt hill and a ∼ 200 meter open air decay region. No signal events were observed. At this experiment, HS vector mesons may have been produced as shown in Fig. 1. Following Refs. [56,57], we use MadGraph5 [58] to calculate the detector acceptance and event yield.
If a HS vector meson is sufficiently long-lived, it can traverse the hill before decaying into electrons in the open region. As shown in the lower panel of Fig. 2 for m π /f π = 3, the lifetime of HS vector mesons that decay to the SM through three-body (two-body) processes is optimal for E137 for ∼ 10 −4 − 10 −3 (10 −5 − 10 −6 ) and GeV-scale vectors. Data from E137 is sensitive to such events as long as these electrons are energetic and directly pass through the ECAL. We define the signal region by requiring that an electron from such decays points back to the target within a resolution of ∼ O(10) mrad and deposits more than 1 − 3 GeV of energy in the calorimeter. This range of energies corresponds to the different analysis thresholds utilized in Ref. [46]. We also include the effect of energy loss as the electron traverses the open air, incorporating a radiation length of ∼ 304 meters [59,60]. This effect is typically negligible in searches for minimal dark photons or axions. However, in the models studied here, final state leptons are produced from a cascade of decays and inherit a smaller fraction of the initial beam energy. For comparable π and V masses, e ± from the decay V → e + e − each carry only about E beam /4 5 GeV of energy. These electrons then lose ∼ 50% of their energy on average while traversing the air, making the location of the decay and the precise value of the energy threshold important. These effects are more significant for three-body decays of HS vector mesons, as the final state electrons are comparatively softer. For m A GeV, roughly 90% and 5% − 90% of events pass the pointing and energy threshold cuts, respectively, where the lower part of the range corresponds to demanding a minimum energy deposition of 3 GeV from V ± → π ± e + e − . For concreteness, we choose a 1 GeV threshold. 10 The E137 regions in Fig. 5 correspond to 95% C.L. exclusions where at least 3 signal events are predicted.
Long-lived particles can also be produced at proton beam dumps [61]. Experiments like CHARM [62] and ν-Cal [63] are sensitive to a similar region of parameter space as E137 [56]. Shorter baseline beam dumps, such as the Orsay experiment [47], offer complementary sensitivity to smaller decay lengths. For m A O(100) MeV, decays of the A into DM pions are constrained by searches for A production and decay followed by DM scattering in a downstream detector such as at LSND [48,64], E137 [46,49], and MiniBooNE [50]. These regions are labeled as "π-scatt." in Fig. 5.

B. HPS
The Heavy Photon Search (HPS) experiment is well-suited to detect visible decays of the HS vector mesons [51]. The first science run was completed in 2016 with a 2.3 GeV electron beam, a 4 µm tungsten target, and 5.6 × 10 17 EOT. A future run is anticipated with a more energetic beam (6.6 GeV), a thicker target (8.8 µm), and a larger luminosity (6.8 × 10 18 EOT) [65]. We simulate vector meson production (shown in Fig. 1) using MadGraph5 [58]. The HPS spectrometer is a smallacceptance detector. We demand that the vector mesons decay to e + e − pairs within 1 − 8 cm of the target and assume that the geometric and kinematic efficiency is 5 − 10%. 11 The displaced vertex, along with the e + e − invariant mass peak (from V → e + e − ) and missing energy (due to the HS pions) can be used to suppress back- FIG. 5: Existing constraints (gray regions) and sensitivity of future searches (colored lines) to signals of strongly interacting hidden sectors. In the left column, we assume that all hidden sector vector mesons decay to Standard Model particles via two-body (V → + − ) or three-body (V → π + − ) processes. In the right column, we assume that the vector mesons that do not mix with the dark photon are heavier than 2 mπ. In this case, only the hidden sector ρ and φ decay into Standard Model particles via two-body processes, while the remaining vector mesons decay invisibly into dark matter pions. The shaded regions are excluded by BaBar [45], E137 [46], Orsay [47], and searches for dark matter scattering at LSND [48], E137 [49], and MiniBooNE [50], as described in the text. The colored contours correspond to the projected reach of HPS [51] (orange), an upgraded version of SeaQuest [52] (magenta), and the proposed LDMX experiment [53] (purple and green). In evaluating experimental exclusions and projected sensitivity, we have fixed αD = 10 −2 , m A /mπ = 3, mV /mπ = 1.8, and mπ/fπ = 3 (4π) in the top (bottom) row. The experimental results are insensitive to small variations in mV /mπ (except for values near thresholds and resonances). In contrast, the dark matter abundance strongly depends on the V -π mass splitting (see Fig. 4).
In each panel, hidden sector pions make up all of the dark matter along the solid (dashed) black contours for mV /mπ = 1. grounds to negligible levels. The future sensitivity of HPS, corresponding to a yield of ∼ 5 − 10 events, is shown as the orange contours in Fig. 5. The kinematic distribution of this signal is shown in the left panel of Fig. 6.

C. SeaQuest
The SeaQuest experiment is currently running at the Fermi National Accelerator Laboratory (FNAL) [52]. It is a nuclear physics spectrometer primarily designed to measure SM Drell-Yan production of muons. The detector trigger has recently been upgraded to retain both low-mass and long-lived dimuon pairs [66][67][68]. Future plans also include the installation of a recycled ECAL, which would allow for additional sensitivity to displaced electrons. SeaQuest utilizes the FNAL 120 GeV main injector proton beam and consists of a 5 meter magnetized iron target/shield (FMAG), a 3 meter focusing magnet (KMAG), and a series of triggering/tracking stations designed for accurate vertex and momentum reconstruction. A parasitic implementation of the displaced vertex trigger is expected to acquire 1.44×10 18 protons on target (POT) within the next two years, and a more dedicated run could obtain an integrated luminosity of up to 10 20 POT [69,70].
We consider three mechanisms for dark photon production at SeaQuest: Bremsstrahlung in the proton-nucleus collisions (see Fig. 1), Drell-Yan, and SM meson decays, e.g., π 0 SM → γA [71]. The kinematic distribution of the signal from proton Bremsstrahlung is shown in the right panel of Fig. 6. As in the previous sections, the dark photons decay promptly into HS vector mesons (A → πV ). For sufficiently small , V can easily traverse the iron shield before decaying into a lepton pair in the region downstream of FMAG. If V decays to e + e − , the level of SM backgrounds is expected to be small, as the magnetic field of FMAG effectively sweeps away soft SM radiation. The optimal ECAL location and form of the displaced electron trigger are uncertain. We will therefore investigate the conservative scenario in which decays occur immediately after FMAG and before the first existing tracking station, corresponding to a decay length of ∼ 5 − 6 meters (as measured from the front of FMAG). We simulate production and decays at SeaQuest using MadGraph5 [58] and PYTHIA 8.2 [72]. We estimate the reach of a future displaced electron search by demanding at least 10 events in which V decays to e + e − in the fiducial region of 5 − 6 meters. We further require that the electron pair remains within the geometric acceptance of the spectrometer, incorporating the transverse momentum kick of KMAG, |∆p T | 0.4 GeV. The magenta contours in Fig. 5 illustrate the potential capability of a future version of SeaQuest, with an upgraded sensitivity to electrons and 10 20 POT. More details on the modeling and capability of this instrument will be presented in the forthcoming work, Ref. [73]. We also note that similar sensitivity and an enhanced mass reach could be achieved with the proposed SHiP experiment (which utilizes the 400 GeV proton beam at CERN SPS) for dark photons heavier than 10 GeV [74].

D. LDMX
The proposed LDMX experiment will be able to test SIMP-like models in two distinct channels: visible displaced V decays and invisible dark photon decays. LDMX is designed to look for missing momentum by accurately measuring a beam electron's momentum before and after it scatters in a ∼ 10% radiation-length tungsten target [53][54][55]. Two experimental phases are planned. Phase 1 (2) will acquire 4 × 10 14 (10 16 ) EOT with a beam energy of E beam = 4 (8) GeV.
The missing momentum program at LDMX requires a large hermetic calorimeter in order to reliably measure the scattered electron energy. This opens the possibility of performing a search for the same visible signals discussed in the previous sections. If a HS vector meson produced in the collision decays to e + e − sufficiently far inside the calorimeter, it can give rise to a displaced electromagnetic shower. The size of the detector can then be used to range out normal electromagnetic backgrounds, such as the shower from the scattered electron. The majority of signal events originate from A production in the first radiation length of the ECAL and not from the thin target. We model this signal as described in Sec. IV B and estimate the sensitivity of such a search using the Phase 2 luminosity by requiring that the visible decays occur between 7 cm (20 radiation lengths in tungsten) and 300 cm (approximate length of the detector). We select candidate events such that the recoil electron energy is less than 0.3 E beam and apply an overall 50% "quality" cut to the yield to crudely model detector efficiencies. The estimated sensitivity, corresponding to a yield of 3 events, is shown as the purple contours in Fig. 5.
The other class of signals arises from invisible A decays. For the SIMP models discussed here, this signal occurs when the A decays to ππ or to vector mesons that decay outside of the detector. As shown in Fig. 2, the branching fraction for A → ππ is O(10)% for cosmologically motivated parameters, m π /f π 5, and falls rapidly for larger m π /f π . We estimate the signal yield for the missing momentum search by rescaling the results in Refs. [54,55], including decays to HS pions and vector mesons. The resulting reach is shown as the green contours in Fig. 5. The currently running NA64 experiment is sensitive to similar missing mass phenomenology [75,76]. However, the present dataset does not probe parameter space beyond existing beam dump constraints. MeV, αD = 10 −2 , and mπ/fπ = 3. We show the distribution in the energy of the visibly decaying hidden sector vector meson, EV = E e + e − , for both direct (V → e + e − ) and indirect production through the decay of an on-shell dark photon (A → πV → πe + e − ), the latter of which is the focus of this work. In the left (right) panel, the beam energy and decay length are fixed to those of the HPS (SeaQuest) experiment. In the right panel, we have additionally selected events that pass the geometric criteria as detailed in Sec. IV C. Direct V -Bremsstrahlung is suppressed by V − A mixing relative to production through A decay. As a result, for a fixed V lifetime, the indirect V -production mechanism dominates, leading to enhanced displaced vertex rates relative to the minimal dark photon-like scenario.

E. Colliders
Processes analogous to those presented in Fig. 1 can be looked for at low-energy e + e − colliders and at the LHC. Searches at BaBar for invisible decays of dark photons in association with a SM photon constitute the strongest probe of strongly interacting hidden sectors at present for 1 GeV m A 10 GeV [45,77]. In recasting this search, we have included decays of the A into HS pions and long-lived vector mesons, demanding that the latter decay outside of the detector (modeled as a cylinder with a 1.5 meter radius). A similar search at Belle-II is expected to improve these limits in the near future, probing values of up to a factor of ∼ 3 smaller compared to BaBar with only 20 fb −1 of data [55]. In Fig. 5, the Belle-II reach is shown as a solid (dotted) blue contour for m A GeV (m A GeV). We note that this sensitivity projection for m A GeV is relatively uncertain as it relies on the rejection and measurement of a peaking background [78].
Searches for new GeV-scale particles at the LHC are complicated by enormous SM backgrounds, which can be suppressed by looking for, e.g., resonances and displaced vertices. Refs. [79,80] used these techniques to demonstrate that LHCb will have impressive sensitivity to the minimal visible dark photon decays (A → + − ) in Run 3 of the LHC. In fact, an analysis has recently been performed with limited data in Ref. [81]. These studies rely on kinematic features that are unique to A → + − in key ways. For example, the D * → D 0 A search requires the lepton pair and D 0 momenta to reconstruct the D * ; this is not the case for the cascade decays considered in this work, A → πV (→ + − ), due to the presence of the final state HS pions [79]. The inclusive analysis of Ref. [80] relies on data-driven estimates of signal and therefore it is not simple to determine the signal and reach to the more complicated A decays considered here (for recent work on recasting LHCb dilepton searches, see Ref. [82]). It would be interesting to investigate if the above search strategies can be generalized to signals present in SIMP-like theories.
Other, more inclusive, searches at B-factories and the LHC can also be used to study the hidden sectors discussed here. For instance, BaBar has performed a search for pairs of tracks that reconstruct a displaced vertex, which may arise from the leptonic or hadronic decays of long-lived particles [83]. For momentum transfers below the HS confinement scale, these types of signals arise from the same visible decays of HS vector mesons, e.g., e + e − → πV (→ µ + µ − ) + · · · . Processes involving larger momentum transfers must be treated as production of HS quarks [84] that subsequently shower and hadronize, yielding a large multiplicity of pions and vector mesons [85][86][87][88]. At the LHC, similar searches at (or near) the ATLAS, CMS, and LHCb detectors could also offer greater sensitivity to heavier HS states, leveraging the additional presence of missing energy arising from the production of stable DM pions or the long lifetimes of HS vector mesons [89][90][91][92]. We leave a detailed exploration of these signals to future work.
Throughout this work, we have assumed that the SM-HS mediator can be produced on-shell. If the dark photon is above the kinematic threshold of certain accelerators, HS states can be produced directly through an offshell A . Standard constraints on visibly decaying dark photons can then be applied to the HS vector mesons, whose coupling to SM states is described by an effective kinetic mixing, eff ∼ (e D / g) (m 2 V / m 2 A ), generated by a heavy A . For instance, taking m A /m V 100 and α D ∼ 10 −1 , constraints from long-baseline beam dumps, which are typically sensitive to ∼ 10 −8 − 10 −6 for light visibly decaying dark photons [54,55], instead probe ∼ 10 −4 − 10 −2 when V 's are produced through an offshell A . However, these limits are always subdominant to the ones discussed throughout this work when on-shell dark photons can be produced and decay to long-lived vector mesons.

F. Direct Detection
Light DM can also be detected with electron or nuclear recoils in low-threshold detectors [55,93]. If the global U (1) D symmetry is preserved by the HS quark masses, i.e., [Q, M q ] = 0, then A exchange can lead to elastic scattering in semi-conductor direct detection experiments such as SENSEI and SuperCDMS [55]. Using the optimistic thresholds and exposures from Ref. [93], we find that SENSEI will be able to test regions of parameter space in the upper-left corner of each panel in Fig.  5 (see also Ref. [12]). Direct detection signals are absent if the A couples off-diagonally to mass eigenstates that are split by more than m π v 2 ∼ 0.05 keV (m π /50 MeV). This can occur if there is a single Higgs field that is responsible for both the A and HS quark masses, resulting in [Q, M q ] = 0. If such mass splittings are smaller than the temperature at DM freeze-out, T ∼ m π /20, the cosmological evolution, as well as signals at fixed-target and collider experiments, are not significantly affected.

G. Self-Interactions
The self-scattering cross-section per DM mass is bounded from observations of galaxy clusters, σ scatter /m π few × cm 2 /g [16][17][18]94]. As in Ref. [3], we calculate this bound assuming that a single species of pion is the lightest state in the theory. As a representative benchmark, we show contours of 1 and 5 cm 2 /g as vertical gray dotted lines in Fig. 5. This bound does not apply if the HS pions make up a subdominant fraction ( 10%) of the total DM energy density.

H. Supernovae
Light particles that are feebly coupled to the SM can lead to anomalous cooling in stellar systems. For instance, DM or long-lived states that are produced in supernovae from an off-shell intermediate A can diffuse out of the core, leading to energy loss. For m π, V 10 MeV, such processes could have lead to qualitative changes in the cooling rate and neutrino burst duration of supernova SN1987A. From the discussion in Ref. [53], we estimate that these bounds may apply for couplings of size α D 2 ∼ 10 −15 . Although not shown, this potential exclusion corresponds to the region of parameter space in the lower-left corner of each panel in Fig. 5. A careful determination of the supernova cooling bound will appear in Refs. [95,96].

I. CMB
For N f = odd flavors of light HS quarks, the singlet DM pions are cosmologically unstable. As discussed in Sec. III, reasonable coefficients for higher-order operators in the chiral Lagrangian lift the masses of the unstable pions above those of the stable species. In this case, 2 → 2 processes (such as π 0 π 0 → π ± π ± ) significantly deplete the fractional abundance of the unstable pions before they decay to SM particles.
Regardless, measurements of the CMB are extremely sensitive to energy injections at the time of recombination from a decaying DM subcomponent [41,42]. In estimating this bound, we take the limits from Ref. [42] and assume that the unstable and stable pion species are split in mass at the level of ∆ ∼ 5%. This choice does not significantly affect the freeze-out calculations discussed throughout this work since DM chemically decouples at T ∼ m π /20. The frozen-out fractional abundance of the unstable species is then proportional to exp(−∆ x 2→2 ), where x 2→2 ∼ O(100) is estimated by demanding that n π σv(2 → 2) ∼ H when T SM ∼ m π /x 2→2 .
Regions of parameter space that are in conflict with measurements of the CMB anisotropies are shown as the light gray dot-dashed contours in Fig. 5. This constraint is severely weakened for larger values of m π /f π since 2 → 2 scattering more efficiently depopulates the unstable species. We emphasize that this exclusion is not intrinsic to models of strongly interacting DM. For instance, as discussed in Sec. III, the DM pions can be made absolutely stable for N f = 4 light flavors, while considerations of DM freeze-out and the phenomenology in Fig. 5 are not qualitatively affected.

J. Discussion of Results
We show existing constraints and the projected reach of experiments in Fig. 5 as gray regions and colored contours, respectively, in the m A − plane. In each panel, we have fixed α D = 10 −2 , m V /m π = 1.8, and m A /m π = 3. In the top (bottom) row, we take m π /f π = 3 (4π). Since we do not expect small variations in the vector meson-pion mass ratio, m V /m π , to significantly impact the phenomenology, thermal freeze-out targets are shown for both m V /m π = 1.8 and m V /m π = 1.6 as solid and dashed black lines, respectively. The final abundance of the HS pions matches the observed DM energy density along these contours.
In the left panel of each row of Fig. 5, we assume that both the neutral and charged vector mesons are lighter than 2 m π . In this case, both two-and three-body decays to SM particles are present in the theory, corresponding to the processes shown in the bottom-left and bottomright diagrams of Fig. 1. The large A production rate makes E137 a powerful test of this mass spectrum. As shown in Fig. 5, E137 excludes a wide swath of models in which sub-GeV HS pions make up all of the DM. In this model, beam dump and fixed-target experiments are sensitive to two distinct ranges of , due to the differing lifetimes of vector mesons that undergo two-or threebody decays (see Fig. 2). For m π /f π = 3, this results in the two lobe-like features in the top-left panel of Fig. 5, corresponding to two-body (lower lobe) and three-body (upper lobe) vector meson decays.
A qualitatively different picture emerges if only twobody decays of vector mesons into the SM are allowed, as in the right panel of each row of Fig. 5. This scenario arises if only the vector mesons that directly mix with the A are lighter than 2 m π . The heavier vector mesons decay invisibly to pairs of HS pions, and the process shown in the bottom-right diagram of Fig. 1 does not occur with an appreciable rate. For this mass spectrum and 10 −4 , the decay lengths of vector mesons are much shorter than the E137 baseline, leaving open a large set of viable models that can be explored in the near future. The discovery prospects of HPS, SeaQuest, and LDMX are shown as the colored contours in Fig. 5. HPS can probe cosmologically interesting models in a future run given the mass spectrum shown in the right panels. On longer timescales, LDMX and especially SeaQuest will explore new territory.
For larger values of m π /f π , as shown in the bottom row of Fig. 5, the branching ratio for A → V π is enhanced (see Fig. 2). More importantly, the decay lengths corresponding to V → + − (V → π + − ) are significantly increased (decreased). As a result, the two baselines present in the theory (see Fig. 2) become comparable in length, and viable parameter space is reopened for ∼ 10 −3 and m A ∼ few × 100 MeV, as shown in the bottom-left panel of Fig. 5. However, larger values of m π /f π also favor heavier thermal DM, shifting the semi-annihilation and WIMP-like freeze-out regions to m π GeV. For m π /f π = 4π and m π GeV, freeze-out is instead controlled by kinetic decoupling, as in Refs. [39,40]. A detailed discussion of such cosmology will appear in forthcoming work. Large values of m π /f π lead to a decrease in the branching fraction of A → ππ, limiting the reach of experiments looking for missing energy signals, like BaBar, Belle-II, and LDMX. For m π /f π = 4π, the signals at these experiments arise primarily from long-lived vector mesons decaying outside of the detectors. In the bottom-left panel of Fig. 5, vector mesons often decay visibly within the detectors at Belle-II and LDMX for ∼ 10 −3 . As a result, the effective invisible A branching fraction and the corresponding sensitivities of searches at BaBar, Belle-II, and LDMX are reduced.
In Fig. 5, several parameters were held fixed. For the mass ratios considered here, variations in α D do not qualitatively change our results. For instance, larger values of α D lead to an overall decrease in the decay length of the vector mesons, which can be compensated by a decrease in the kinetic mixing, . As a result, the sensitivity contours of the fixed-target and beam dump experiments are correspondingly shifted down in the m A − plane, along with an overall suppression in the mass reach due to the suppressed production cross-section. Different choices for the mass ratios can have more significant effects. For example, if m A > 2 m V , then the new decay channel, A → V + V − → π + π − + 4 , becomes kinematically accessible. Furthermore, the visible signals we have considered are suppressed if m A m V + m π , although HS particle production is still possible through an offshell dark photon. In this case, direct detection, LDMX, and other missing momentum/energy experiments provide the most powerful probes of these models.

V. SUMMARY AND CONCLUSIONS
Strongly interacting hidden sectors offer a compelling possibility for accounting for the cosmological abundance of DM. This class of models requires the existence of a coupling between the hidden sector (HS) and the SM to be cosmologically viable. GeV-scale dark photons are ideal in this regard since kinetic mixing with SM hypercharge facilitates equilibration between the two sectors. Cosmology motivates scenarios in which the lightest spin-1 resonances of the HS are parametrically close in mass to the DM pions. In this case, visible decays of these vector mesons play an important role in the cosmological history of DM freeze-out. In this work, we showed that 2 → 1 processes (ππ → πV followed by V → SM) typically dominate over SIMP-like 3 → 2 annihilations and greatly expand the allowed ranges of masses and couplings for thermal DM.
The presence of visibly-decaying vector mesons also opens a new avenue for experimental searches. The experimentally viable range of HS-SM couplings spans several orders of magnitude, ∼ 10 −7 − 10 −3 , corresponding to a wide range of vector meson lifetimes. Thus, large regions of well-motivated parameter space can be probed by searches for production and displaced decays of HS vector mesons. We showed that existing data from past beam dump and collider experiments al-ready exclude some of these models. However, most of the cosmologically-motivated parameter space remains untested. Existing experiments like HPS will be able to cover new ground in this regard.
The wide range of viable vector meson lifetimes requires a multi-experiment approach to comprehensively cover the remaining low-mass parameter space. The next generation of fixed-target experiments in combination with long-lived particle searches at colliders will have unprecedented discovery reach for the vector mesons as well as the DM itself. For m A few × O(100) MeV, the proposed LDMX experiment will be able to search for invisible dark photon decays of the form A → ππ, as well as the displaced vertex signals discussed throughout this work. For m A GeV, Belle-II and an upgraded version of the existing SeaQuest spectrometer at FNAL will provide the most stringent tests. Above the muon threshold, similar searches at B-factories and the LHC can probe dark photons significantly heavier than a GeV.
In this paper, we have focused on a particular realization of a strongly interacting HS that mirrors SM QCD with three colors and three light flavors. However, the cosmology and the corresponding experimental signals are applicable to a much broader range of models with different gauge groups and flavor structures. This motivates ongoing and future experiments to test this interesting class of DM scenarios.
In this Appendix, we briefly review the building blocks that are necessary for specifying interactions in the HS. Throughout this work, we have investigated the dynamics of a confining SU (N c ) gauge theory, as discussed in Sec. II. We focus on N c = 3 colors and N f = 3 flavors of light HS quarks in the fundamental representation of SU (N c ), in analogy to SM QCD. The spontaneous breaking of the global chiral symmetry down to the diagonal subgroup proceeds as in the SM, i.e., SU (N f ) L × SU (N f ) R → SU (N f ) V . Below the confinement scale, the N 2 f − 1 broken generators correspond to the presence of pseudo-Nambu-Goldstone pions, π a , which are the relevant low-energy degrees of freedom in the HS.
The interactions of the HS pions are described by the chiral Lagrangian, analogous to that of the SM. We write the exponentiated pion fields as where T a are the broken generators of SU (N f ) and f π is the associated pion decay constant of the HS. We work in the convention where these generators are normalized to Tr T a T b = 1 2 δ ab . Since the number of colors and light flavors in the HS is identical to that of the light quark sector of the SM, we adopt the usual SM names for the HS mesons. We write the pion matrix as The superscripts indicate charge under an additional U (1) D gauge symmetry as discussed below. The kinetic and mass terms are where M q is the HS quark mass matrix and B 0 is a dimensionful scale. In the limit that M q ∝ 1, the mass of the degenerate pions is m 2 π = 2 B 0 M q to leading order in chiral perturbation theory. We follow the phenomenological approach of Ref. [32] in order to incorporate the spin-1 vector mesons into the low-energy theory. Adopting the SM naming convention, we parametrize the vector meson fields as where in the sum we have included the diagonal generator, T 0 = 1/ √ 6, corresponding to the U (1) subgroup of the U (N f ) flavor symmetry. The kinetic terms for the vector mesons are given by where the field strength is We assume that the vector mesons share a common mass, m V , unless stated otherwise. As mentioned in Sec. II, the KSRF relation implies that the vector meson coupling, g, can be estimated as g m V /( √ 2 f π ). We have utilized this parametrization throughout. This coupling also directly enters into the V ππ interaction by promoting the derivative of the pion field to the covariant form Interactions of the dark photon, A , with the HS mesons are incorporated by gauging a U (1) D subgroup of the unbroken diagonal flavor symmetry. We add to the covariant pion derivative, where e D is the U (1) D gauge coupling and Q is the U (1) D charge matrix of the HS quarks, as specified in Eq. (2). The vector mesons acquire a coupling to the A through kinetic mixing [12,32] Alternatively, interactions between the A and HS mesons can be implemented in the Vector Meson Dominance limit of the Hidden Local Symmetry framework [33]. In this case, the A only couples to π ± and V ± pairs through intermediate neutral vector mesons, as in Eq. (A8).

Appendix B: Decays
In this Appendix, we provide expressions for the decays of the A and HS vector mesons. This is especially relevant for the experimental signals outlined in Sec. IV. For the sake of brevity, we introduce the notation where i = , π, V denotes a SM lepton, HS pion, or HS vector meson, respectively. The A can decay to SM fermions through kinetic mixing with SM hypercharge. In the limit that m A m Z , the width for the decay of an A into a pair of SM leptons is approximated by We incorporate the decays of the dark photon to SM hadrons through the data-driven parameter R( √ s) ≡ σ(e + e − → hadrons)/σ(e + e − → µ + µ − ), as in Ref. [59]. In particular, If m A 2 m π , A can decay to pairs of HS pions that are directly charged under U (1) D . From Eq. (A7), the corresponding width is where α D ≡ e 2 D /4π. The interactions between the A and pairs of vector mesons follow in a similar manner, leading to For m A m π + m V , A also decays to πV final states. The calculation for such processes is analogous to determining π 0 SM → γγ in the SM and follows from the standard anomalous chiral current algebra. We find If m V 2 m π , then vector mesons can only decay to the SM. If Tr [Q T a ] = 0 as well, this decay proceeds through V a kinetically mixing with A , as in Eq. (A8), and V a dominantly decays to pairs of SM fermions. In the limit that m V , m A m Z , the corresponding decays to SM leptons are given by Similar to Eq. (B3), we can also express the decay of vector mesons into SM hadrons through On the other hand, if m V 2 m π but Tr [Q T a ] = 0, then V a decays through a three-body process involving an off-shell A , i.e., V a → π b A * → π b + − , as discussed in Sec. II. In the limit that m m V , m A m Z , the non-zero widths are approximated as − 39 r 2 π + r 2 V + 15 3r 4 π + 3r 4 V + 2r 2 π r 2 V + 12 × log 1 + −r 2 π +r 2 V +1 (r 2 π −1) 2 +r 4 V −2(r 2 π +1)r 2 V 1/2 1 − r 2 π −r 2 V +1 (r 2 π −1) 2 +r 4 V −2(r 2 π +1)r 2 V 1/2 1 − −r 2 π +r 2 V +1 (r 2 π −1) 2 +r 4 V −2(r 2 π +1)r 2 V 1/2 1 + r 2 π −r 2 V +1 (r 2 π −1) 2 +r 4 V −2(r 2 π +1)r 2 V 1/2 + r 4 π − 2r 2 π r 2 V + 1 + r 2 V − 1 2 1/2 17 r 6 π − r 6 V + 42 r 4 V − r 4 π + 39 r 2 π r 4 V − r 4 π r 2 V + 24 r 2 π − r 2 V + 36 r 4 π + r 4 V − 54 r 2 π + r 2 V − 6 r 4 π + r 4 V − 4r 2 π r 2 V r 2 π + r 2 V + 24 log and
Above, brackets denote thermal averaging, and we have been explicit at which temperature various quantities are to be evaluated. The equilibrium distributions are approximated assuming Maxwell-Boltzmann statistics, n eq g m 2 T 2 π 2 K 2 (m/T ) , ρ eq g m 2 T 2 π 2 (m K 1 (m/T ) + 3 T K 2 (m/T )) , P eq T n eq , where K 1,2 are modified Bessel functions of the second kind, and n eq π(V ) is evaluated with g = 8 (9) flavor degrees of freedom. We numerically solve Eq. (C1) for n π (T SM ), n V (T SM ), and T HS (T SM ) after substituting , ρ n n eq ρ eq , P n n eq P eq T n .
We now provide expressions for the effective decay and annihilation rates that enter into Eq. (C1). In the limit that m V m , the thermally-averaged decay rate into SM leptons is given by We also include decays to hadrons in Eq. (C1), as described in Appendix B. In the non-relativistic limit, the thermallyaveraged cross-section for DM annihilation into SM leptons is approximately, where we have again taken m 0. In our analysis, we numerically evaluate the general form of σv(ππ → + − ) , following the procedure outlined in Refs. [97,98]. We also incorporate direct annihilations to SM hadrons, as discussed in Ref. [89]. In evaluating σv 2 (πππ → ππ) and σv(π SM → π SM) ∆E , we utilize the results in Ref. [3] and Ref. [39], respectively. We follow Refs. [37,99] for the calculation of the semi-annihilation process. In the non-relativistic limit, we find σv(πV → ππ) 5 (m π /f π ) 8 m 2 V (m V − m π )(3 m π + m V ) 3/2 T 1179648 π 5 m 10 π (m π + 2 m V ) 2 (m 2 V − m 2 π + m π m V ) If the vector mesons that do not undergo two-body decays to SM particles are decoupled from the low-energy theory, as discussed in Sec. IV, then various quantities in Eq. (C1) are rescaled, n eq V → 2 9 n eq V , σv(πV → ππ) →